Simulate Prediction Interval

Steps 1 & 2 from Procedure:

prep_interval() takes in a workflow (model specification + pre-processing recipe¹¹) along with a training dataset, and outputs a named list containing bootstrapped model fits + prepped recipes (model_uncertainty) and the resulting residuals from cross-validation (sample_uncertainty).

# load custom functions `prep_interval()` and ``predict_interval()`
devtools::source_gist("https://gist.github.com/brshallo/4053df78265ab9d77f753d95f5faaf5b")

set.seed(1234)
prepped_for_interval <- prep_interval(lm_rec_mod, ames_train, n_boot = 200)

prepped_for_interval 
## $model_uncertainty
## # A tibble: 200 x 2
##    fit      recipe  
##    <list>   <list>  
##  1 <fit[+]> <recipe>
##  2 <fit[+]> <recipe>
##  3 <fit[+]> <recipe>
##  4 <fit[+]> <recipe>
##  5 <fit[+]> <recipe>
##  6 <fit[+]> <recipe>
##  7 <fit[+]> <recipe>
##  8 <fit[+]> <recipe>
##  9 <fit[+]> <recipe>
## 10 <fit[+]> <recipe>
## # ... with 190 more rows
## 
## $sample_uncertainty
## # A tibble: 2,197 x 1
##    .resid
##     <dbl>
##  1 -0.347
##  2 -0.321
##  3 -0.304
##  4 -0.291
##  5 -0.278
##  6 -0.266
##  7 -0.257
##  8 -0.247
##  9 -0.239
## 10 -0.232
## # ... with 2,187 more rows

Steps 3-6 in Procedure:

predict_interval() takes in the output from prep_interval(), along with the dataset for which we want to produce predictions (as well as the quantiles associated with our confidence level of interest), and returns a dataframe containing the specified prediction intervals.

pred_interval <- predict_interval(prepped_for_interval, ames_holdout, probs = c(0.05, 0.50, 0.95)) 

lm_sim_intervals <- pred_interval %>% 
  mutate(across(contains("probs"), ~10^.x)) %>% 
  bind_cols(ames_holdout) %>% 
  select(Sale_Price, contains("probs"), Lot_Area, Neighborhood, Years_Old, Gr_Liv_Area, Overall_Qual, Total_Bsmt_SF, Garage_Area)

lm_sim_intervals <- lm_sim_intervals %>% 
  rename(.pred = probs_0.50, .pred_lower = probs_0.05, .pred_upper = probs_0.95) %>% 
  relocate(c(.pred_lower, .pred_upper, .pred)) 

See gist for more documentation on prep_interval() and predict_interval().

Interval Width

The relative widths of prediction intervals from the simulation based approach vary between observations by more than 10x what we saw in the Interval Widths using analytic methods¹² (where interval width had a standard deviation of less than half a percentage point).

We also see greater differences in interval width between buckets of predictions (a range of about 5.8 percentage points vs 0.4).

This suggests the simulation based approach allows for greater differentiation in estimated levels of uncertainty by attributes¹³.

lm_sim_widths <- lm_sim_intervals %>% 
  mutate(interval_width = .pred_upper - .pred_lower,
         interval_pred_ratio = interval_width / .pred) %>% 
  mutate(price_grouped = ggplot2::cut_number(.pred, 5)) %>% 
  group_by(price_grouped) %>% 
  summarise(n = n(),
            mean_interval_width_percentage = mean(interval_pred_ratio),
            stdev = sd(interval_pred_ratio),
            stderror = sd(interval_pred_ratio) / sqrt(n)) %>% 
  mutate(x_tmp = str_sub(price_grouped, 2, -2)) %>% 
  separate(x_tmp, c("min", "max"), sep = ",") %>% 
  mutate(across(c(min, max), as.double)) %>% 
  select(-price_grouped) 

lm_sim_widths %>% 
  mutate(across(c(mean_interval_width_percentage, stdev, stderror), ~.x*100)) %>% 
  gt::gt() %>% 
  gt::fmt_number(c("stdev", "stderror"), decimals = 2) %>% 
  gt::fmt_number("mean_interval_width_percentage", decimals = 1) 

Simultaneously, the average relative interval width is slightly more narrow: about 52% for the simulation based approach against 54% for the analytic method¹⁴.

lm_sim_intervals %>% 
  mutate(interval_width = .pred_upper - .pred_lower,
         interval_pred_ratio = interval_width / .pred) %>% 
  summarise(n = n(),
            mean_interval_width_percentage = mean(interval_pred_ratio),
            stderror = sd(interval_pred_ratio) / sqrt(n)) %>% 
  mutate(across(c(mean_interval_width_percentage, stderror), ~.x * 100)) %>% 
  gt::gt() %>% 
  gt::fmt_number(c("mean_interval_width_percentage", "stderror"), decimals = 2) %>% 
  gt::fmt_number("mean_interval_width_percentage", decimals = 1) 

Coverage

Coverage is essentially the same: 92.3% (vs 92.7% for Coverage with the analytic method)¹⁵.

lm_sim_intervals %>%
  mutate(covered = ifelse(Sale_Price >= .pred_lower & Sale_Price <= .pred_upper, 1, 0)) %>% 
  summarise(n = n(),
            n_covered = sum(
              covered
            ),
            stderror = sd(covered) / sqrt(n),
            coverage_prop = n_covered / n) %>% 
  mutate(across(c(coverage_prop, stderror), ~.x * 100)) %>% 
  gt::gt() %>% 
  gt::fmt_number("stderror", decimals = 2) %>% 
  gt::fmt_number("coverage_prop", decimals = 1) 

Coverage rates across quintiles:

lm_sim_intervals %>% 
  mutate(price_grouped = ggplot2::cut_number(.pred, 5)) %>% 
  mutate(covered = ifelse(Sale_Price >= .pred_lower & Sale_Price <= .pred_upper, 1, 0)) %>% 
  group_by(price_grouped) %>% 
  summarise(n = n(),
            n_covered = sum(
              covered
            ),
            stderror = sd(covered) / sqrt(n),
            n_prop = n_covered / n) %>% 
  mutate(x_tmp = str_sub(price_grouped, 2, -2)) %>% 
  separate(x_tmp, c("min", "max"), sep = ",") %>% 
  mutate(across(c(min, max), as.double)) %>% 
  ggplot(aes(x = forcats::fct_reorder(scales::dollar(max), max), y = n_prop))+
  geom_line(aes(group = 1))+
  geom_errorbar(aes(ymin = n_prop - 2 * stderror, ymax = n_prop + 2 * stderror))+
  coord_cartesian(ylim = c(0.70, 1.01))+
  # scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  labs(x = "Max Predicted Price for Quintile",
       y = "Coverage at Quintile",
       title = "Coverage by Quintile of Predictions",
       subtitle = "On a holdout Set",
       caption = "Error bars represent {coverage} +/- 2 * {coverage standard error}")+
  theme_bw() 

The overall pattern appears broadly similar to that seen in the analytic method. However coverage rates seem to be slightly more consistent across quintiles in the simulation based approach – an improvement over the analytic method (where coverage levels had appeared to be slightly more dependent upon quintile of predicted Sale_Price)¹⁶.

For example, a chi-squared test in the previous post had shown significant variation in coverage rates across quintiles when applied to the analytic method. However the same test applied to the results from the simulation based prediction intervals does not show a significant difference¹⁷.

lm_sim_intervals %>% 
  mutate(price_grouped = ggplot2::cut_number(.pred, 5)) %>% 
  mutate(covered = ifelse(Sale_Price >= .pred_lower & Sale_Price <= .pred_upper, 1, 0)) %>% 
  with(chisq.test(price_grouped, covered)) %>% 
  pander::pander() 

Keep in mind we are looking at just one model specification on one set of data¹⁸.

See Adjusting Procedure in the Appendix for a couple examples of circumstances where you might want to make changes to Procedure.

Conformal Inference

Highly related to what I presented in this post, the academic discipline of using out-of-sample data to create prediction intervals is known as conformal inference (much of the research in this field comes from Carnegie Mellon University and Royal Holloway University, London). I may do a follow-up post where I walk through a more formal example of using conformal inference. In the meantime, here are a few resources I glanced through²²:

• ryantibs/conformal: github repo with conformalInference R package and links to relevant articles on distribution-free predictive inference. The way model types are specified by conformalInferene seems to be not too dissimilar from the approach I took in this post²³
• donlnz/nonconformist: python package for conformal inference
• Conformal Prediction: Link to Royal Holloway University website by creators of method – Vladimir Vovk and Alex Gammerman.
• Assumption-free prediction intervals for black-box regression algorithms - Aaditya Ramdas (YouTube): professor at CMU giving overview of problem, approaches, and current “state-of-the-art”.
• Tutorial on conformal inference, Dataiku article, Analytics Vidhya article

Other Examples Using Simulation

Here is a video walk-through of a simple example using a linear model:

A slightly modified version of this approach that adjusts for the leverage of the observations (leverage²⁴ has to do with distance of a point from the centroid of the data²⁵) can be found at this Cross Validated Thread.

Limitation with these examples (as encoded):

• re-estimate some component of their prediction intervals based on performance on the data used to train the model.
• are set-up for linear models

Confusion With Confidence Intervals

It seems to be common for people to seem to set-out to simulate prediction intervals but end-up simulating confidence intervals²⁶ – i.e. they account for uncertainty in estimating the model while forgetting to account for the uncertainty of the sample. I wrote more on this in Prediction Intervals and Confidence Intervals from Part 1²⁷.

Adjusting Procedure

Below are just a few examples. (Procedure is more just a toy set-up and by no means optimized.)

Coverage level depends on attribute(s)

If we see that coverage rates depend on predicted value of Sale_Price (or on some other attribute) we might:

• try altering the pre-processing recipe
• adjust Procedure to segment²⁸ the residuals by predicted Sale_Price and shuffle the errors within segments (rather than across the entire dataset). This would allow the uncertainty due to the sample to vary according to the predicted Sale_Price²⁹.

Random Splits Not Appropriate

The procedure I walk through assumes a random split in your data is fine and that your observations are iid. This is generally a requirement for using the bootstrap, however there are some approaches to using modified bootstrap that can adjust for this. In many cases it may make sense to have group or time based sampling schemes (see my presentation that discusses time-based resampling schemes). If a random sampling scheme is inappropriate, you will likely end-up creating too small of prediction intervals. I had considered making prep_interval() capable of taking in custom resampling specifications but did not set this up… but could be potentially edited.

Alternative Procedure With CV

Given the model training dataset with n number of observations…

1. Build b number of models using b bootstrap resamples on the training dataset (default of b is the square root of n)³⁰.
2. Build another set of k models³¹ using k-fold cross-validation³² (default of k is 10) and extract the residuals, composed of n elements – providing distribution for uncertainty due to the sample³³.

Given an observation (or set of observations) you would like to produce prediction intervals for…

3. For each new observation, produce a prediction using all of the b models created in step 1.
4. Take the difference between each model’s prediction (from step 3) and the mean of all the model’s predictions – the resulting b differences for each observation provides the distribution for the variability due to model estimation at each point.
5. For each observation, repeatedly sample from “residuals” (step 2) and “model error / differences” (step 4) and add together – do this h times (with replacement) OR create all possible combinations from step 2 and step 5. Rather than h, the latter approach will produce a distribution composed of n x b elements³⁴ for each new observation (default behavior³⁵).
6. Pull quantiles (from the distribution created in step 5) according to the desired level of coverage of your intervals – e.g. 0.05 and 0.95 for a 90% prediction interval.

See gist for documentation on implementation of this method.

Differences from Procedure

This approach is influenced by suggestions from Max Kuhn. For example, I …

• only use out-of-sample estimates to produce the interval
• estimate the uncertainty of the sample using the residuals from a separate set of models built with cross-validation³⁶