I am working with GBM models using the caret package and looking to find a method to solve the prediction intervals for my predicted data. I have searched extensively but only come up with a few ideas to find prediction intervals for Random Forest. Any help/R code would be greatly appreciated!


1 Answer 1


EDIT: As pointed out in the comments below this gives the confidence intervals for predictions and not strictly the prediction intervals. Was a bit trigger happy with my reply and should have given this some extra thought.

Feel free to ignore this answer or try to build on the code to get the prediction intervals.

I have used the simple bootstrap for creating prediction intervals a few times but there may be other (better) ways.

Consider the oil data in the caret package and suppose we want to generate partial dependencies and 95% intervals for the effect of Stearic on Palmitic. Below is just a simple example but you can play around with it to suit your needs. Make sure the gbm package is update to allow the grid.points argument in plot.gbm

#train the gbm using just the defaults.
tr <- train(Palmitic ~ ., method = "gbm" ,data = fattyAcids, verbose = FALSE)

#Points to be used for prediction. Use the quartiles here just for illustration
x.pt <- quantile(fattyAcids$Stearic, c(0.25, 0.5, 0.75))

#Generate the predictions, or in this case, the partial dependencies at the selected points. Substitute plot() for predict() to get predictions
p <- plot(tr$finalModel, "Stearic", grid.levels = x.pt, return.grid = TRUE)

#Bootstrap the process to get prediction intervals

bootfun <- function(data, indices) {
  data <- data[indices,]

  #As before, just the defaults in this example. Palmitic is the first variable, hence data[,1]
  tr <- train(data[,-1], data[,1], method = "gbm", verbose=FALSE)

  # ... other steps, e.g. using the oneSE rule etc ...
  #Return partial dependencies (or predictions)

  plot(tr$finalModel, "Stearic", grid.levels = x.pt, return.grid = TRUE)$y
  #or predict(tr$finalModel, data = ...)

#Perform the bootstrap, this can be very time consuming. Just 99 replicates here but we usually want to do more, e.g. 500. Consider using the parallel option
b <- boot(data = fattyAcids, statistic = bootfun, R = 99)

#Get the 95% intervals from the boot object as the 2.5th and 97.5th percentiles
lims <- t(apply(b$t, 2, FUN = function(x) quantile(x, c(0.025, 0.975))))

This is one way to do it which at least try to account for the uncertainties arising from tuning the gbm. A similar approach has been used in http://onlinelibrary.wiley.com/doi/10.2193/2006-503/abstract

Sometimes the point estimate is outside the interval, but modifying the tuning grid (i.e., increasing the number of trees and/or the depth) usually solves that.

Hope this helps!

  • 2
    $\begingroup$ If I understand your code correctly, what you have there is a 95% confidence interval for predictions. That's not the same as a 95% prediction interval, which adds the residual (random) error. $\endgroup$
    – Hong Ooi
    Commented Sep 30, 2014 at 8:32
  • 1
    $\begingroup$ D'oh! You are correct. Was a bit too quick in replying. Thanks, I'll edit my answer. $\endgroup$
    – ErikL
    Commented Sep 30, 2014 at 8:44
  • $\begingroup$ thanks for help! I am having an issue though with the bootstrap function. I posted that issue at stats.stackexchange.com/questions/117329/… . I am not exactly sure how to set up the bootstrap function properly with my dataset. $\endgroup$ Commented Sep 30, 2014 at 13:21
  • 1
    $\begingroup$ I guess at this point this is not what I am looking for, so I am still looking for an answer! $\endgroup$ Commented Oct 1, 2014 at 19:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.