Evaluate posterior predictive distribution in Bayesian linear regression

Question

I'm confused on how to evaluate the posterior predictive distribution for Bayesian linear regression, past the basic case described here on page 3, and copied below.

$$ p(\tilde y \mid y) = \int p(\tilde y \mid \beta, \sigma^2) p(\beta, \sigma^2 \mid y) $$

The basic case is this linear regression model:

$$ y = X \beta + \epsilon, \hspace{10mm} y \sim N(X \beta, \sigma^2)$$

If we use either a uniform prior on $\beta$, with a scale-Inv$\chi^2$ prior on $\sigma^2$, OR the normal-inverse-gamma prior (see here) the posterior predictive distribution is analytic and is student t.

What about for this model? $$ y = X \beta + \epsilon, \hspace{10mm} y \sim N(X \beta, \Sigma)$$

When $y \sim N(X\beta, \Sigma)$, but $\Sigma$ is known, the posterior predictive distribution is multivariate Gaussian. Usually, you don't know $\Sigma$, but have to estimate it. Maybe you say its diagonal and make the diagonal a function of the covariates in some way. This is discussed in the linear regression chapter of Gelman's Bayesian Data Analysis.

Is there an analytic form for the posterior predictive distribution in this case? Can I just plug my estimate of it into a multivariate student t? If you estimate more than one variance, is the distribution still multivariate student t?

I am asking because say I have some $\tilde y$ already on hand. I want to know whether it is more likely to have been predicted by e.g. linear regression A, linear regression B

Answer 1

If you assume a uniform prior on $\beta$, then the posterior for $\beta$ is $$ \beta|y \sim N(\hat{\beta},V_\beta). $$ with $$\hat{\beta} = [X'\Sigma^{-1}X]X'y \quad \mbox{and} \quad V_\beta =[X'\Sigma^{-1}X]^{-1}.$$ To find the predictive distribution, we need more information. If $\tilde{y} \sim N(\tilde{X}\beta,\tilde{\Sigma})$ and is conditionally independent of $y$ given $\beta$, then $$ \tilde{y}|y \sim N(\tilde{X}\hat{\beta}, \tilde{\Sigma} + V_\beta).$$ But typically for these types of models, $y$ and $\tilde{y}$ are not conditionally independent, instead, we typically have $$ \left( \begin{array}{c} y \\ \tilde{y} \end{array} \right) \sim N\left( \left[\begin{array}{c} X\beta \\ \tilde{X}\beta \end{array}\right], \left[ \begin{array}{cc} \Sigma & \Sigma_{12} \\ \Sigma_{21} & \tilde{\Sigma} \end{array} \right]\right). $$ If this is the case, then $$ \tilde{y}|y \sim N(\tilde{X}\hat{\beta} + \Sigma_{21}\Sigma^{-1}(y-X\hat{\beta}), \tilde{\Sigma} - \Sigma_{21}\Sigma^{-1}\Sigma_{12}). $$ This assumes $\Sigma, \Sigma_{12},$ and $\tilde{\Sigma}$ are all known. As you point out, typically they are unknown and need to be estimated. For the common models that have this structure, e.g. time series and spatial models, there generally won't be a closed form for the predictive distribution.

Answer 2

Under non-informative or multivariate Normal-Wishart priors, you have the analytical form as a multivariate Student's distribution, for a classical mutlivariate, multiple regression. I guess the developments in this document are related to your question (you may like the Appendix A :-) ). I typically compared the outcome with a posterior predictive distribution obtained using WinBUGS and the analytical form: they are exactly equivalent. The problem only becomes difficult when you have additional random effects in mixed-effect models, especially in unbalanced design.

In general, with classical regressions, y and ỹ are conditionally independent (residuals are i.i.d) ! Of course if it is not the case, then the proposed solution here is not correct.

In R,(here, solution for uniform priors), assuming you made a lm model (named "model") of one of the response in your model, and called it "model", here is how to obtain the multivariate predictive distribution

library(mvtnorm)
Y = as.matrix(datas[,c("resp1","resp2","resp3")])
X =  model.matrix(delete.response(terms(model)), 
           data, model$contrasts)
XprimeX  = t(X) %*% X
XprimeXinv = solve(xprimex)
hatB =  xprimexinv %*% t(X) %*% Y
A = t(Y - X%*%hatB)%*% (Y-X%*%hatB)
F = ncol(X)
M = ncol(Y)
N = nrow(Y)
nu= N-(M+F)+1 #nu must be positive
C_1 =  c(1  + x0 %*% xprimexinv %*% t(x0)) #for a prediction of the factor setting x0 (a vector of size F=ncol(X))
varY = A/(nu) 
postmean = x0 %*% hatB
nsim = 2000
ysim = rmvt(n=nsim,delta=postmux0,C_1*varY,df=nu) 

Now, quantiles of ysim are beta-expectation tolerance intervals from the predictive distribution, you can of course directly use the sampled distribution to do whatever you want.