In Bayesian linear regression Advantages of predictive posterior compared to posterior of model coefficients

Question

In Bayesian linear regression, if we want to get confidence intervals for predictions of a new observation. I was thinking of the following two options.

1. Use the quantiles from samples sampled from the coefficents in the posterior $\beta | \mathbf{y}; \mathbf{X}$

2. Use the predictive posterior distribution, $y^* \, | \beta, \mathbf{y}; \mathbf{x}^*, \mathbf{X}$

Where $\mathbf{x}^*$ is a new observation. What are the advantages of 2. over 1.

---

If we have the setup:

\begin{gather} \mathbf{y} \, | \beta ; \mathbf{X} \sim \mathcal{N}( \mathbf{X} \beta, \mathbf{I} \sigma^2) \\ \beta \sim \mathcal{N}(0, \mathbf{I}\sigma_{\beta}^2) \nonumber \end{gather}

Then we can get a closed form solution for the posterior of $\beta$

\begin{equation} \beta \, | \mathbf{y}; \mathbf{X} \sim \mathcal{N}\big((\mathbf{X}^T \mathbf{X} + \lambda \mathbf{I})^{-1} \mathbf{X} \mathbf{y}^T , \; (\mathbf{X}^T \mathbf{X} + \lambda \mathbf{I})^{-1} \sigma^2\big) \end{equation}

where $\lambda := \frac{\sigma^2}{\sigma_{\beta}^2}$.

Therefore can see that for a new point $\mathbf{x}^*$ the variance of $(\mathbf{x}^*)^T\beta$ will be

$$ (\mathbf{x}^*)^T((\mathbf{X}^T \mathbf{X} + \lambda \mathbf{I})^{-1} \sigma^2)\mathbf{x}^* $$

But, if we compute the predictive posterior:

\begin{align*} p(y^* \, | \beta, \mathbf{y}; \mathbf{x}^*, \mathbf{X}) = \int p(y^* \, | \beta; \mathbf{x}^*, \sigma^2) p(\beta \, | \mathbf{y}; \mathbf{X}, \sigma^2, \sigma_{\beta}^2) \; d\beta \end{align*}

we can compute this integral and get that the variance is:

$$ \sigma^2 + (\mathbf{x}^*)^T((\mathbf{X}^T \mathbf{X} + \lambda \mathbf{I})^{-1} \sigma^2)\mathbf{x}^* $$

Answer 1

Your question is more involved than it seems. Let me explain.

First, you have to decide on a tentative model for regression. Well, every Bayesian model involves a likelihood function and a prior probability. The likelihood function usually is rather well understood. Broadly speaking, the likelihood function describes how the data agree with the tentative model.

But what about the prior? The prior requires you to specify everything you know about the solution. For example: images cannot have negative pixel values, and physical processes cannot change abruptly. All of this is part of your professional knowledge. Or should be. And YOU have to decide on your prior!

Second, after multiplying the prior probability distribution and the likelihood function, you have obtained an unnormalized posterior probability distribution. But your model is still tentative. There is no alternative than to compute the normalization constant. This constant is sometimes called evidence or partition function.

Third, the next step is "model selection." The evidence is a numerical indication of the quality of the tentative model. But with the evidence values of several alternative models at hand, you can make a well-motivated selection of the best model. Only after you have selected the best model can you start with predictions and the corresponding error bars.

I have written a tutorial about this subject What is your model?, which describes the many pitfalls. You don't need to buy it as you can read 80% for free on Amazon. Good luck.

Answer 2

When using MCMC to sample the predictions you would first take a sample from the posterior for the parameters $\beta$, and then plug in the parameters and the data $\tilde{\mathbf{x}}$.

Let's use the Stan code example to illustrate it.

// [...]
model {
  y ~ poisson_log(alpha + beta * x);
  { alpha, beta } ~ normal(0, 1);
}
generated quantities {
  array[N_tilde] int<lower=0> y_tilde
    = poisson_log_rng(alpha + beta * x_tilde);
}

What the code does, is it first samples the alpha and beta parameters, then it predicts the mean of the Poisson distribution in alpha + beta * x_tilde, and for each x_tilde it samples one value from the Poisson distribution with the predicted mean.

So it does use the samples from the posterior distribution but the samples alone do not produce the posterior predictive samples. I'm not sure how you would like to use them for this purpose other than by plugging them into the model to sample the predictions.