Prediction intervals for cubist model of R caret package

Question

I am using model training function train from caret package to train cubist rule-based prediction model. How to generate the prediction intervals for a response variable of a new observation?

Here's a working example:

library(caret)
library(mlbench)
data(BostonHousing)

myControl <- trainControl(method = 'cv', 
                            repeats = 5, number = 10, 
                            returnResamp = 'none', 
                            returnData = FALSE, savePredictions = TRUE, 
                            verboseIter = TRUE, allowParallel = TRUE)

  inTrain <- createDataPartition(y = BostonHousing$medv,
                                 p = .66,
                                 list = FALSE)

  training <- BostonHousing[ inTrain,]
  testing <- BostonHousing[-inTrain,]                            

  cubistModel <- train(medv ~ ., data = training, method='cubist',  trControl = myControl, tuneLength = 10)

# New case:
predict(cubistModel, newdata = testing[1,])

Session info:

> sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: CentOS release 6.8 (Final)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
[7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] mlbench_2.1-1   Cubist_0.0.18   caret_6.0-71    ggplot2_2.1.0   lattice_0.20-33

loaded via a namespace (and not attached):
 [1] Rcpp_0.11.6        magrittr_1.5       splines_3.3.1      MASS_7.3-45        munsell_0.4.3      colorspace_1.2-6   foreach_1.4.3      minqa_1.2.4       
 [9] stringr_1.1.0      car_2.1-3          plyr_1.8.4         tools_3.3.1        nnet_7.3-12        pbkrtest_0.4-6     parallel_3.3.1     grid_3.3.1        
[17] gtable_0.2.0       nlme_3.1-128       mgcv_1.8-12        quantreg_5.26      MatrixModels_0.4-1 iterators_1.0.8    lme4_1.1-12        Matrix_1.2-6      
[25] nloptr_1.0.4       reshape2_1.4.1     codetools_0.2-14   stringi_1.1.1      compiler_3.3.1     scales_0.4.0       stats4_3.3.1       SparseM_1.7 

Answer 1

Basically, you can't (outside of resampling). While these are linear models in the terminal nodes, they are created by blending hierarchical linear models (i.e. the math is not tractable).