Considering Uncertainty

Let’s consider a pretend offer of $180,000 for our example house. You might first think you ought to reject the offer due to it being ~$4,500 less than the expected price. However a less naive approach would be to compare the $180,000 offer in the context of the model’s observed accuracy when making predictions on a holdout set¹⁰.

We need to be careful regarding how we consider performance and the errors of our model. For this particular problem, the variability in our model’s errors increase with the magnitude of the predictions for Sale_Price¹¹. We want an error metric that will not inflate with Sale_Price. Hence rather than review the difference between Sale_Price and Expected[Sale_Price], it is more appropriate to review either:

• \(\log_{10}(Sale\_Price)\) against \(Expected[\log_{10}(Sale\_Price)]\) OR
• an error metric that is in terms of percent of Sale_Price

I will use the latter method and focus on the percentage error¹². The offer of $180,000 corresponds with a 2.5% difference from expectations. This falls well within the distribution of errors seen by our model on a holdout dataset¹³.

offer <- tibble(PE = (184503 - 180000) / 180000, offer = "offer")

data_preds %>% 
  mutate(PE = (Sale_Price - .pred) / Sale_Price) %>% 
  relocate(PE) %>% 
  ggplot(aes(x = PE))+
  geom_histogram(bins = 50)+
  geom_vline(aes(xintercept = PE, color = offer), data = offer)+
  scale_x_continuous(limits = c(-1, 1), labels = scales::percent)+
  labs(x = "Percent Error")+
  theme_bw()

On average our model is off by around ~12% of the actual Sale_Price observed.

data_preds %>% 
  yardstick::mape(Sale_Price, .pred) %>% 
  gt::gt() %>% 
  gt::fmt(gt::everything(), fns = function(x) fmt_if_number(x, digits = 1))

90% of errors on the holdout dataset were between -27.1% and 21.3% of the actual Sale_Price.

data_preds %>% 
  mutate(PE = (Sale_Price - .pred) / Sale_Price,
         APE = abs(PE)) %>% 
  summarise(quant_05 = quantile(PE, 0.05) * 100,
            quant_95 = quantile(PE, 0.95) * 100) %>% 
  gt::gt() %>% 
  gt::fmt(gt::everything(), fns = function(x) fmt_if_number(x, digits = 1))

All of which suggests the $180,000 offer does not represent a substantial outlier from the typical variability of observed from predicted Sale_Price¹⁴.

Is your model performant enough to be useful?

Before using the model for predictive inference, one should have reviewed overall performance on a holdout dataset to ensure the model is sufficiently accurate for the business context. For example, for our problem is an average error of ~12% and 90% prediction intervals of +/- ~25% of Sale_Price useful? If the answer is “no,” that suggests the need for more effort in improving the accuracy of the model (e.g. trying other transformations, features, model types)¹⁵. For our examples we are assuming the answer is ‘yes,’ our model is accurate enough (so it is appropriate to move-on and focus on prediction intervals).

Observation Specific Intervals

While a good starting point, a limitation with using aggregate error metrics to estimate intervals indiscriminately across observations¹⁶ is that different attributes may associate with different levels of uncertainty in the prediction. Specifically, predictions further from the centroid of the data generally have more uncertainty in the expected price.

A toy seesaw is a good analogy for this phenomenon:

image source

As the angle of the bench changes… the further you are from the center, the more distance you will move up/down →

The angle of the seesaw represents variability in the model function being estimated. The distance from the seesaw’s pivot point corresponds with the distance from the centroid of the data.

This variability in estimating the model’s expected value is represented by dashed blue lines in the chart below of a generic linear model:

Image source

Because of this varying uncertainty of the expected value of the model (represented by the differing interval widths encapsulated by the blue lines at any value of x), it is preferable to determine a plausible price range that is specific to the attributes of the observation¹⁷.

The uncertainty described here corresponds with the confidence intervals and not quite as directly with the prediction intervals (the focus of this post). The distinction and relationship between these will be discussed in A Few Things to Know About Prediction Intervals.

Prediction Intervals and Confidence Intervals

Prediction intervals and confidence intervals¹⁸ are often confused. Confidence intervals generally refer to making inferences on averages – this is most useful for evaluating parameter estimates, performance metrics, relationships with covariates, etc. However if you are interested in price ranges on individual observations (as in our case), prediction intervals are what you want¹⁹. Rob Hyndman has a helpful post where he describes the differences in more detail: The difference between prediction intervals and confidence intervals²⁰.

In linear regression, “prediction intervals” refer to a type of confidence interval²¹, namely the confidence interval for a single observation (a “predictive confidence interval”). Confidence intervals have a specific statistical interpretation. In later posts on this topic, the intervals I create do not quite mirror the interpretations that go with a predictive confidence interval. I will use the term “prediction interval” somewhat loosely to refer to a plausible range of values for an observation²².

Analytic Method of Calculating Prediction Intervals

Linear regression models can produce prediction intervals analytically. The equation below is for simple linear regression (meaning just one ‘x’ input) but is helpful for gaining an intuition on the key parts that contribute to the width of a prediction interval:

\(\hat{y}_h \pm t_{(1-\alpha/2, n-2)} \times \sqrt{MSE \times \left( 1+\dfrac{1}{n} + \dfrac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2}\right)}\)

Pseudo-translation of key parts of equation:

• {\(\hat{y}_h\) : prediction}
  
• {\(t_{(1-\alpha/2, n-2)}\) : multiplier for desired level of confidence, e.g. 95%}
  
• {MSE: multiplier for average error}
  
• {\(\dfrac{MSE}{n}\) : multiplier based on number of observations – smaller n contributes to greater variability.}
• {\(MSE\dfrac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2}\) : multiplier for distance from centroid of data}

You will notice the equation for the prediction interval is very similar to that of the equation for a confidence interval²³:

\(\hat{y}_h \pm t_{(1-\alpha/2, n-2)} \times \sqrt{MSE \left(\dfrac{1}{n} + \dfrac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2}\right)}\)

The variance for the prediction interval just has an extra Mean Squared Error (MSE) term²⁴. The prediction interval is essentially the variance in estimating the model²⁵ combined with the variability of individual observations in the sample.

Paraphrasing, the uncertainty in predictions can be thought of as coming from two sources:

1. uncertainty in fitting the model, i.e. variability in determining the expected value of the target (‘target’ just means the variable of interest, the thing we are predicting)
2. uncertainty in the sample, i.e. inherent variability of the target in observations

The first source may vary in part depending on position relative to the centroid (think back to the seesaw analogy in Observation Specific Intervals and the influence of position on amount of movement i.e. uncertainty). It will also vary depending on the number of observations – the more observations you have in your training data, the smaller the general variability due to model estimation will be. The latter part is assumed to be constant across observations²⁶ and does not change as the number of observations increases.

Do not think of these sources of uncertainty as being compounded on one another to produce the prediction interval. They are added together²⁷ and then the square root is taken. In most cases the variability due to the sample will be greater than the variability in estimating the expected value. Hence, the interval’s width due to uncertainty in estimating the model is somewhat ‘watered down’ by the variability in the sample. Therefore (for prediction intervals produced analytically at least) we will generally not see large differences in the widths of prediction intervals across observations (even at points that are relatively far from the centroid of the data). Therefore our prediction intervals in later sections of this post will not actually be that much different from those that might have been produced using aggregated error metrics similar to those discussed in Considering Uncertainty.

However understanding the intuition on these sources of uncertainty in prediction intervals is still helpful. In follow-up posts I use more flexible methods for building prediction intervals. These other methods can produce substantive differences in the width of prediction intervals between observations (in ways that depend on more than distance of an observation from the centroid of the data).

Visual Comparison of Prediction Intervals and Confidence Intervals

The chart below shows the key outputs concerning prediction that can be outputted by a generic linear model:

• expected value, i.e. the predictions, the average of the value of the target given a set of attributes (black line)
• confidence intervals, i.e. range for the expected values (red dashed line)
• prediction intervals (our primary interest in this post), i.e. range for the predictions on an individual observation (green dotted line)

Image source

Confidence intervals are more narrow (as they concern only uncertainty in the model’s estimation of the expected value). Prediction intervals are wider as they concern both model uncertainty and the sampling uncertainty of any observation (the latter of which contributes far greater variance). Both confidence and prediction intervals are wider the further the prediction is from the centroid of the data – however the effect is far greater for the confidence interval (as the uncertainty in the prediction interval is dominated by the random variance of the sample, which is assumed to be constant across observations). This explains why the curve of the confidence intervals is far more pronounced.

Inference or Prediction?

Predictive modeling is typically used either for making predictions or inferences. The assumptions associated with your model usually matter more when doing inference compared to prediction²⁸. Therefore, to get a sense of how careful you should be regarding model assumptions for our example you might ask:

“Are prediction intervals an example of prediction or inference?”

The answer is, while it (kind of) depends…

really, prediction intervals should be thought of as a kind of inference.

For this post the prediction interval is explicitly an inference on the predictions. In a later post I will actually just be making predictions for quantiles at upper and lower bounds of interest. However, generally, intervals should be thought of as a kind of inference²⁹. Therefore you should be thoughtful regarding the assumptions you are making about the nature of the data and your model specification³⁰.

For example, an important assumption when doing predictive inference is heteroskedasticity of errors – which means that the variability of our errors should be constant across our predictions[^puzzle]. As mentioned in Considering Uncertainty, if building a model for Sale Price we would find that errors tend to increase with our predictions. Hence instead of building a model to predict Sale Price we build a model to predict \(\log_{10}(Sale\:Price)\) – this helps our model’s outputted prediction intervals to be more appropriate across observations³¹.

In follow up posts the methods I describe do not depend as strongly on assumptions to produce reliable prediction intervals³².

[puzzle]: You may be puzzled why differing uncertainty depending on the distance of a point is from the centroid of the data is not thought of as a kind of heteroskedasticity. My understanding is this is becuase that difference has to do with uncertainty in estimating the model, not uncertainty in the sample (which is assumed to be constant).

Cautions With Overfitting

You want the inferences and predictions of your model to be generalizable, i.e. usable outside of the context of your training dataset. This is why we holdout data to review performance: we want to review the model’s performance on data it has never seen. If your model overfits³³ the training data, the prediction interval will often underestimate the expected variability in price and your prediction intervals may be too narrow. Two strategies you might use to correct for this:

1. Consider a more simple model³⁴
2. Tune your prediction intervals based on coverage measured on holdout data (coverage represents the proportion of observations that actually fall within their prediction interval).

A quick heuristic to check if your model is overfitting is to compare the model’s performance on a holdout dataset against performance on the data that was used to train the model. If performance is substantially better on the training dataset, it may suggest the model has overfit.

Simpler model

If your model seems to be overfitting, you may want to try fitting a more simple model³⁵. If your model has a small difference in performance metrics between training and holdout datasets, you likely will have more faith in the ranges given by the prediction intervals.

Problem with comparing difference in performance between training & holdout data

However, for the most part, your performance is going to always be better on the training data than on the holdout data³⁶. With regard to overfitting, you really care about whether performance is worse on the holdout dataset compared to an alternative simpler model’s performance on the holdout set. You don’t really care if a model’s performance on training and holdout data is similar, just that performance on a holdout dataset is as good as possible.

Let’s take an example when building random forest models and tuning on the min_n attribute (smaller values for min_n³⁷ represent more complicated models.):

See my gist for code

The gap in train-holdout (AKA, analysis-assessment³⁸) performance is greater for more complicated models³⁹ (which, as mentioned above, may signal overfitting). However more complicated models also perform better on the holdout data. Hence, we would likely select the complicated model⁴⁰ even though the performance estimates on the training data will be the most unrealistic.

While a difference in performance on training and holdout datasets may provide a helpful indicator for overfitting, it is not really what you care about during model selection⁴¹.

Tune your intervals based on a holdout set

Pretend that your selected model has a substantial difference in performance when evaluated on training versus holdout data. You are worried that the width of your model’s outputted prediction intervals will be optimistically narrow. However you still want to keep this model because it has the best performance when evaluated on holdout data⁴². An alternative to fitting a more simple model is to adjust the confidence level of your prediction intervals, tuning them based on the coverage level on a holdout dataset. To describe this solution, I will define three ways of thinking about coverage:

• target coverage: The level of coverage you want to attain on a holdout dataset (i.e. the proportion of observations you want to fall within your prediction intervals).
• expected coverage: The level of confidence in the model for the prediction intervals, e.g. asking the model to predict 90% prediction intervals.
• empirical coverage: The level of coverage actually observed when evaluated on a dataset, typically a holdout dataset not used in training the model.

Ideally all three align. However let’s say in our hypothetical example you have a 90% target coverage. If our model is slightly overfit, you might see that a 90% expected coverage leads to an 85% empirical coverage on a holdout dataset. To align your target and empirical coverage at 90%, may require setting expected coverage at something like 93%⁴³.

This is called “adaptive coverage” and is discussed briefly by Emmanuel Candes in the video below:

In a follow-up post, Quantile Regression for Prediction Intervals I will walk through a similar example where I tune the expected coverage.

Coverage

Let’s review what proportion of our holdout observations are actually “covered” by their 90% prediction intervals⁴⁶.

preds_intervals %>%
  mutate(covered = ifelse(Sale_Price >= .pred_lower & Sale_Price <= .pred_upper, 1, 0)) %>% 
  summarise(n = n(),
            n_covered = sum(
              covered
            ),
            stderror = 100 * sd(covered) / sqrt(n),
            coverage_prop = 100 * n_covered / n) %>% 
  gt::gt() %>% 
  gt::fmt(gt::everything(), fns = function(x) fmt_if_number(x, digits = 2))

~93% of our holdout observations were actually covered⁴⁷. We can likely trust these intervals aren’t too narrow⁴⁸. It is usually better to be on the slightly conservative than optimistic side regarding coverage⁴⁹.

In addition to reviewing overall coverage, it may be helpful to review coverage across predictions. This may be helpful in providing evidence of whether we have “marginal” or “conditional” coverage⁵⁰.

The chart below splits the holdout data into five segments ordered by predicted price with equal number of observations⁵¹ and checks the proportion covered over each quintile.

preds_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()

This suggests there may be slightly different levels of coverage at different quintiles (e.g. the cheapest predicted houses have less coverage than those in the 2nd quintile.)⁵².

A test for checking whether variability in the categorical variable covered (covered / not-covered) does not vary across the quintile of predicted Sale_Price (1st, 2nd, 3rd, 4th, 5th) could be done using a Chi-square statistical test:

preds_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()

The results provide evidence that whether an observation is covered or not-covered may depend (to some extent) on which quintile of predicted Sale_Price the observation falls within⁵³.

I will show in a follow-up post on Simulating Prediction Intervals intervals that seem to have more consistent coverage rates across predicted Sale_Price.

Interval Width

Another helpful way to evaluate prediction intervals is by their width. More narrow bands indicate a more precise model.

Remember from the section on Considering Uncertainty that the errors on raw Sale_Price increase with the size of the prediction. To adjust for this, an appropriate metric is interval width as a percentage of predicted Sale_Price⁵⁴. I will review these across quintiles of the predictions.

lm_interval_widths <- preds_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 = stdev / 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_interval_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)

The 90% interval represents an average of ~54% of the predicted Sale_Price. I.e. the difference between the upper and lower bounds of the prediction interval for Sale_Price will be about half of the predicted Sale_Price (to eye-ball this, see Prediction Intervals on Raw Sale Price).

The relative interval width seems to be consistent across values for predicted Sale_Price (~54% for all quintiles). There is also little difference in uncertainty between observations – varying on average less than half a percent between observations as indicated by the standard deviation.

Recall from Analytic Method of Calculating Prediction Intervals that the only factor contributing to differences in interval width across predictions is the distance of an observations from the centroid of the data. This has a small impact on relative interval width in our problem (given the number of observations and sample variance – which are treated as constant across observations)⁵⁵. Methods used in future posts will show a greater diversity of interval widths.

Bayesian Inference

In this and the immediate follow-up posts I approach prediction intervals from a more frequentist perspective. However when you really care about thinking through the nature of uncertainty in your model, Bayesian methods are often more appropriate. I may cover Bayesian approaches in a future post, in the meantime, here are a few resources:

• The stan computing platform for Bayesian modeling and inference has various excellent R package interfaces. rstan, rstanarm, brms and tidybayes are the ones of greatest interest to me⁵⁶.
  
• Michael Clark provides a helpful applied introduction to Bayesian analysis in R.
• Note also that you can interact with stan through the parsnip package by setting parsnip::set_engine(engine = "stan") for linear models.
• There is also a python interface, PyStan. The most popular Bayesian library in python however is PyMC3⁵⁷. You might also check-out Tensorflow Probability⁵⁸.

Pros & Cons

Upsides of prediction intervals using parametric methods with linear regression:

• Most common (largest number of people are familiar)
• Straightforward calculation with low computation costs

Downsides:

• Less flexible than alternative approaches.
• Prediction intervals produced in this way have a variety of strong assumptions that they depend on (generally more heavily than other approaches).
• Assumes the correct model (If your model overfits, the associated prediction intervals will likely be overly optimistic⁵⁹).

Prediction Intervals on Raw Sale Price

Because the model is on \(\log_{10}{(Sale\:Price)}\) when we transform our scale to be in terms of Sale_Price, the errors would increase with the predicted price of the house – hence why the transformed scales in Review Prediction Intervals are appropriate. The figure below shows the prediction intervals on the scale of raw Sale_Price.

p + 
  scale_x_continuous(labels = scales::dollar)+
  scale_y_continuous(labels = scales::dollar)

parsnip Support for Prediction Intervals

Run the following code to check which parsnip supported models have modules for prediction intervals or quantiles (cross-posted here)".

library(tidyverse)

envir <- parsnip::get_model_env()

ls(envir) %>% 
  tibble(name = .) %>% 
  filter(str_detect(name, "_predict")) %>% 
  mutate(prediction_modules  = map(name, parsnip::get_from_env)) %>% 
  unnest(prediction_modules) %>% 
  filter(str_detect(type, "pred_int|quantile"))

In future posts, I will use some approaches that are not directly (currently) supported by predict methods in parsnip.