Do the predictions of a Random Forest model have a prediction interval?

Question

If I run a randomForest model, I can then make predictions based on the model. Is there a way to get a prediction interval of each of the predictions such that I know how "sure" the model is of its answer. If this is possible is it simply based on the variability of the dependent variable for the whole model or will it have wider and narrower intervals depending on the particular decision tree that was followed for a particular prediction?

Answer 1

This is partly a response to @Sashikanth Dareddy (since it will not fit in a comment) and partly a response to the original post.

Remember what a prediction interval is, it is an interval or set of values where we predict that future observations will lie. Generally the prediction interval has 2 main pieces that determine its width, a piece representing the uncertainty about the predicted mean (or other parameter) this is the confidence interval part, and a piece representing the variability of the individual observations around that mean. The confidence interval is fairy robust due to the Central Limit Theorem and in the case of a random forest, the bootstrapping helps as well. But the prediction interval is completely dependent on the assumptions about how the data is distributed given the predictor variables, CLT and bootstrapping have no effect on that part.

The prediction interval should be wider where the corresponding confidence interval would also be wider. Other things that would affect the width of the prediction interval are assumptions about equal variance or not, this has to come from the knowledge of the researcher, not the random forest model.

A prediction interval does not make sense for a categorical outcome (you could do a prediction set rather than an interval, but most of the time it would probably not be very informative).

We can see some of the issues around prediction intervals by simulating data where we know the exact truth. Consider the following data:

set.seed(1)

x1 <- rep(0:1, each=500)
x2 <- rep(0:1, each=250, length=1000)

y <- 10 + 5*x1 + 10*x2 - 3*x1*x2 + rnorm(1000)

This particular data follows the assumptions for a linear regression and is fairly straight forward for a random forest fit. We know from the "true" model that when both predictors are 0 that the mean is 10, we also know that the individual points follow a normal distribution with standard deviation of 1. This means that the 95% prediction interval based on perfect knowledge for these points would be from 8 to 12 (well actually 8.04 to 11.96, but rounding keeps it simpler). Any estimated prediction interval should be wider than this (not having perfect information adds width to compensate) and include this range.

Let's look at the intervals from regression:

fit1 <- lm(y ~ x1 * x2)

newdat <- expand.grid(x1=0:1, x2=0:1)

(pred.lm.ci <- predict(fit1, newdat, interval='confidence'))
#        fit       lwr      upr
# 1 10.02217  9.893664 10.15067
# 2 14.90927 14.780765 15.03778
# 3 20.02312 19.894613 20.15162
# 4 21.99885 21.870343 22.12735

(pred.lm.pi <- predict(fit1, newdat, interval='prediction'))
#        fit      lwr      upr
# 1 10.02217  7.98626 12.05808
# 2 14.90927 12.87336 16.94518
# 3 20.02312 17.98721 22.05903
# 4 21.99885 19.96294 24.03476

We can see there is some uncertainty in the estimated means (confidence interval) and that gives us a prediction interval that is wider (but includes) the 8 to 12 range.

Now let's look at the interval based on the individual predictions of individual trees (we should expect these to be wider since the random forest does not benefit from the assumptions (which we know to be true for this data) that the linear regression does):

library(randomForest)
fit2 <- randomForest(y ~ x1 + x2, ntree=1001)

pred.rf <- predict(fit2, newdat, predict.all=TRUE)

pred.rf.int <- apply(pred.rf$individual, 1, function(x) {
  c(mean(x) + c(-1, 1) * sd(x), 
  quantile(x, c(0.025, 0.975)))
})

t(pred.rf.int)
#                           2.5%    97.5%
# 1  9.785533 13.88629  9.920507 15.28662
# 2 13.017484 17.22297 12.330821 18.65796
# 3 16.764298 21.40525 14.749296 21.09071
# 4 19.494116 22.33632 18.245580 22.09904

The intervals are wider than the regression prediction intervals, but they don't cover the entire range. They do include the true values and therefore may be legitimate as confidence intervals, but they are only predicting where the mean (predicted value) is, no the added piece for the distribution around that mean. For the first case where x1 and x2 are both 0 the intervals don't go below 9.7, this is very different from the true prediction interval that goes down to 8. If we generate new data points then there will be several points (much more than 5%) that are in the true and regression intervals, but don't fall in the random forest intervals.

To generate a prediction interval you will need to make some strong assumptions about the distribution of the individual points around the predicted means, then you could take the predictions from the individual trees (the bootstrapped confidence interval piece) then generate a random value from the assumed distribution with that center. The quantiles for those generated pieces may form the prediction interval (but I would still test it, you may need to repeat the process several more times and combine).

Here is an example of doing this by adding normal (since we know the original data used a normal) deviations to the predictions with the standard deviation based on the estimated MSE from that tree:

pred.rf.int2 <- sapply(1:4, function(i) {
  tmp <- pred.rf$individual[i, ] + rnorm(1001, 0, sqrt(fit2$mse))
  quantile(tmp, c(0.025, 0.975))
})
t(pred.rf.int2)
#           2.5%    97.5%
# [1,]  7.351609 17.31065
# [2,] 10.386273 20.23700
# [3,] 13.004428 23.55154
# [4,] 16.344504 24.35970

These intervals contain those based on perfect knowledge, so look reasonable. But, they will depend greatly on the assumptions made (the assumptions are valid here because we used the knowledge of how the data was simulated, they may not be as valid in real data cases). I would still repeat the simulations several times for data that looks more like your real data (but simulated so you know the truth) several times before fully trusting this method.

Answer 2

If you use R you can easily produce prediction intervals for the predictions of a random forests regression: Just use the package quantregForest (available at CRAN) and read the paper by N. Meinshausen on how conditional quantiles can be inferred with quantile regression forests and how they can be used to build prediction intervals. Very informative even if you don't work with R!

Answer 3

I realize this is an old post but I have been running some simulations on this and thought I will share my findings.

@GregSnow made a very detailed post about this but I believe when calculating the interval using predictions from individual trees he was looking at $[ \mu + \sigma, \mu - \sigma]$ which is only a 70% prediction interval. We need to look at $[\mu+1.96*\sigma,\mu-1.96*\sigma]$ to get the 95% prediction interval.

Making this change to @GregSnow code, we get the following results

set.seed(1)
x1 <- rep( 0:1, each=500 )
x2 <- rep( 0:1, each=250, length=1000 )
y <- 10 + 5*x1 + 10*x2 - 3*x1*x2 + rnorm(1000)

library(randomForest)
fit2 <- randomForest(y~x1+x2)
pred.rf <- predict(fit2, newdat, predict.all=TRUE)
pred.rf.int <- t(apply( pred.rf$individual, 1, function(x){ 
  c( mean(x) + c(-1.96,1.96)*sd(x), quantile(x, c(0.025,0.975)) )}))

pred.rf.int
                          2.5%    97.5%
1  7.826896 16.05521  9.915482 15.31431
2 11.010662 19.35793 12.298995 18.64296
3 14.296697 23.61657 14.749248 21.11239
4 18.000229 23.73539 18.237448 22.10331

Now, comparing these with the intervals generated by adding normal deviation to predictions with standard deviation as MSE like @GregSnow suggested we get,

pred.rf.int2 <- sapply(1:4, function(i) {
   tmp <- pred.rf$individual[i,] + rnorm(1000, 0, sqrt(fit2$mse))
   quantile(tmp, c(0.025, 0.975))
   })
t(pred.rf.int2)
          2.5%    97.5%
[1,]  7.486895 17.21144
[2,] 10.551811 20.50633
[3,] 12.959318 23.46027
[4,] 16.444967 24.57601

The interval by both these approaches are now looking very close. Plotting the prediction interval for the three approaches against the error distribution in this case looks as below

• Black lines = prediction intervals from linear regression,
• Red lines = Random forest intervals calculated on Individual Predictions,
• Blue lines = Random forest intervals calculated by adding normal deviation to predictions

Now, let us re-run the simulation but this time increasing the variance of the error term. If our prediction interval calculations are good, we should end up with wider intervals than what we got above.

set.seed(1)
x1 <- rep( 0:1, each=500 )
x2 <- rep( 0:1, each=250, length=1000 )
y <- 10 + 5*x1 + 10*x2 - 3*x1*x2 + rnorm(1000,mean=0,sd=5)

fit1 <- lm(y~x1+x2)
newdat <- expand.grid(x1=0:1,x2=0:1)
predict(fit1,newdata=newdat,interval = "prediction")
      fit       lwr      upr
1 10.75006  0.503170 20.99695
2 13.90714  3.660248 24.15403
3 19.47638  9.229490 29.72327
4 22.63346 12.386568 32.88035

set.seed(1)
fit2 <- randomForest(y~x1+x2,localImp=T)
pred.rf.int <- t(apply( pred.rf$individual, 1, function(x){ 
  c( mean(x) + c(-1.96,1.96)*sd(x), quantile(x, c(0.025,0.975)) )}))
pred.rf.int
                          2.5%    97.5%
1  7.889934 15.53642  9.564565 15.47893
2 10.616744 18.78837 11.965325 18.51922
3 15.024598 23.67563 14.724964 21.43195
4 17.967246 23.88760 17.858866 22.54337

pred.rf.int2 <- sapply(1:4, function(i) {
   tmp <- pred.rf$individual[i,] + rnorm(1000, 0, sqrt(fit2$mse))
   quantile(tmp, c(0.025, 0.975))
   })
t(pred.rf.int2)
         2.5%    97.5%
[1,] 1.291450 22.89231
[2,] 4.193414 25.93963
[3,] 7.428309 30.07291
[4,] 9.938158 31.63777

Now, this makes it clear that calculating the prediction intervals by the second approach is far more accurate and is yielding results quite close to the prediction interval from linear regression.

Taking the assumption of normality, there is another easier way to compute the prediction intervals from random forest. From each of the individual trees we have the predicted value ($\mu_i$) as well as the mean squared error ($MSE_i$). So prediction from each individual tree can be thought of as $ N(\mu_i,RMSE_i)$. Using the normal distribution properties our prediction from the random forest would have the distribution $N(\sum \mu_i/n, \sum RMSE_i/n)$. Applying this to the example we discussed above, we get the below results

mean.rf <- pred.rf$aggregate
sd.rf <- mean(sqrt(fit2$mse))
pred.rf.int3 <- cbind(mean.rf - 1.96*sd.rf, mean.rf + 1.96*sd.rf)
pred.rf.int3
1  1.332711 22.09364
2  4.322090 25.08302
3  8.969650 29.73058
4 10.546957 31.30789

These tally very well with the linear model intervals and also the approach @GregSnow suggested. But note that the underlying assumption in all the methods we discussed is that the errors follow a Normal distribution.

Answer 4

This is easy to solve with randomForest.

First let me deal with the regression task (assuming your forest has 1000 trees). In the predict function, you have the option to return results from individual trees. This means that you will receive 1000 column output. We can take the average of the 1000 columns for each row - this is the regular output RF would have produced any way. Now to get the prediction interval lets say +/- 2 std. deviations all you need to do is, for each row, from the 1000 values calculate +/-2 std. deviations and make these your upper and lower bounds on your prediction.

Second, in the case of classification, remember that each tree outputs either 1 or 0 (by default)and the sum over all 1000 trees divied by 1000 gives the class probablity (in the case of binary classification). In order to put a prediction interval on the probability you need to modify the min. nodesize option (see randomForest docuementation for the exact name of that option) once you set it a value >>1 then the individual trees will output numbers between 1 and 0. Now, from here on you can repeat the same process as described above for the regression task.

I hope that makes sense.

Answer 5

The problem of constructing prediction intervals for random forest predictions has been addressed in the following paper:

Zhang, Haozhe, Joshua Zimmerman, Dan Nettleton, and Daniel J. Nordman. "Random Forest Prediction Intervals." The American Statistician,2019.

The R package "rfinterval" is its implementation available at CRAN.

Installation

To install the R package rfinterval:

#install.packages("devtools")
#devtools::install_github(repo="haozhestat/rfinterval")
install.packages("rfinterval")
library(rfinterval)
?rfinterval

Usage

Quickstart:

train_data <- sim_data(n = 1000, p = 10)
test_data <- sim_data(n = 1000, p = 10)

output <- rfinterval(y~., train_data = train_data, test_data = test_data,
                     method = c("oob", "split-conformal", "quantreg"),
                     symmetry = TRUE,alpha = 0.1)

### print the marginal coverage of OOB prediction interval
mean(output$oob_interval$lo < test_data$y & output$oob_interval$up > test_data$y)

### print the marginal coverage of Split-conformal prediction interval
mean(output$sc_interval$lo < test_data$y & output$sc_interval$up > test_data$y)

### print the marginal coverage of Quantile regression forest prediction interval
mean(output$quantreg_interval$lo < test_data$y & output$quantreg_interval$up > test_data$y)

Data example:

oob_interval <- rfinterval(pm2.5 ~ .,
                            train_data = BeijingPM25[1:1000, ],
                            test_data = BeijingPM25[1001:2000, ],
                            method = "oob",
                            symmetry = TRUE,
                            alpha = 0.1)
str(oob_interval)

Answer 6

I've tried some options (this all WIP):

1. I actually made the dependent variable a classification problem with the results as ranges, instead of a single value. The results I got were poor, compared to using a plain value. I gave up this approach.

2. I then converted it to multiple classification problems, each of which was a lower-bound for the range (the result of the model being whether it would cross the lower bound or not) and then ran all the models (~20), and then combined the result to get a final answer as a range. This works better than 1 above but not as good as I need it to. I'm still working to improve this approach.

I used OOB and leave-one-out estimates to decide how good/bad my models are.