How to represent a complete picture of uncertainty in linear regression predictions [rather than intervals]?

Question

Say I have a set of data for three quantities $x$, $z$, and $y$ at various locations, $i$. I then fit a linear model,

$$ y_i = a x_i + b z_i + \epsilon_i $$

to this data, assuming $\epsilon_i \sim N(0, \sigma)$, with the value of $\sigma$ unknown.

I'm then interested in predicting $y$ at some new locations $j \in \{1,2,...n\}$, and most importantly I'd like to quantify uncertainty in that prediction. We have $x$ and $z$ values in these new locations. The question is how do I generate something like a histogram of $y$ at each of these locations that account for all uncertainty. Presumably, simply saying

$$y_j \sim N( \hat a x_j + \hat b z_j, \hat \sigma ),$$

where $\hat a$ and $\hat b$ are the estimates of the coefficients from the linear regression model and $\hat \sigma$ is the residual standard error, is not account for all of the uncertainty, because it assumes the estimated coefficients are true.

So for example how might I plot a histogram for each $y_j$, given the model that accounts for all uncertainty. I think this is related to prediction, vs. confidence intervals for which there are a lot of tutorials online, but I'm interested in the full distribution. The R code below does roughly what I want to do, but I don't think it accounts for enough uncertainty.

## data for fitting
n = 100;
x = rnorm(n, 1, .5);
z = rnorm(n, .5, .1);
y = rnorm(n, mean = 2*x - 6*z, sd = 1)

mod = lm(y~x+z)

## new data
xn = rnorm(n, 1, .5);
zn = rnorm(n, .5, .1);
ytrue = rnorm(n, mean = 2*xn - 6*zn, sd = 1)

## predict with uncertainty
j=2; #choose a site to plot uncertainty
y_pred_samp = rnorm(500, 
              mean = mod$coef[1]*xn[j] + mod$coef[2]*zn[j], 
              sd = sigma(mod))
  
hist(y_pred_samp, main = paste('true value is: ', ytrue[j]))

Answer 1

Your R code treats the parameters as if the estimates were the true values, and only covers the uncertainty in residual noise. This is typically not what we want in a prediction interval. (You are right that this is the term you want to learn about.)

In order to account for the uncertainty in the parameter estimates, we need two new ingredients:

1. The covariance matrix of the parameter estimates (viewed as random variables, since they are transformations of the observations, which are in themselves realizations of a random variable)
2. The $t$ distribution instead of the normal distribution that we started with.

Ordinary least squares is special in that we have a neat closed form solution of the predictive distribution (assuming, of course, that our initial assumptions about the model form and the normal distribution of the error are correct). Specifically, the predictive distribution is a scaled and mean-shifted $t$ distribution, with central $1-\alpha$ level prediction intervals given by

$$ \hat{y}\pm t_{n-p}^{(\alpha/2)}\hat{\sigma}\sqrt{1+x^t(X^tX)^{-1}x}, $$

where

• $x$ contains the predictor values for the prediction (assumed known and non-random)
• $X$ is the design matrix of the predictors used in training
• $\hat{y}=x^t\hat{\beta}$ contains the expectation prediction
• $\hat{\sigma}$ contains the estimate of the residual standard deviation
• $t_{n-p}^{(\alpha/2)}$ contains the $\alpha/2$ quantile of the $t$ distribution with $n-p$ degrees of freedom (where $n$ is the number of historical observations, and $p$ is the number of parameters estimated)

You can find this in any slightly more advanced book on regression, e.g., Faraday, 2002, Practical Regression and Anova using R, section 3.5. Proofs are usually to be found in standard textbooks on mathematical statistics.

So in your specific situation, you could create a histogram of predicted outcomes for a specific value of your predictors as follows:

x_new <- 0.2
z_new <- 0.4
newdata <- matrix(c(1,x_new,z_new),nrow=1,dimnames=list(NULL,c("Intercept","x","z")))
y_pred <- predict(mod,newdata=as.data.frame(newdata))

sims <- y_pred + 
    rt(500,n-3) * summary(mod)$sigma * 
    sqrt(1+newdata%*%solve(crossprod(model.matrix(mod)))%*%t(newdata))[1,1]
hist(sims)

Note that inverting a matrix using solve(), as we do above for clarity, may be unstable, depending on the matrix. It works fine here, though.

Alternatively, if all you wanted was a predictive interval, you can use predict.lm() (as above), with the parameter interval="prediction". This will give you prediction intervals as above, which is something different from confidence intervals (which in turn you could get using interval="confidence").