Overestimation of the noise precision in Bayesian linear regression when $n\gtrsim p$

Question

I would like to infer the regression coefficients and the noise precision of a standard linear regression problem defined by

$$ y=X\theta + \epsilon, $$

where $X$ is a $n\times p$ design matrix, $\theta$ is a length $p$ vector of regression coefficients, $y$ is a length $n$ vector of observations of the dependent variable, and $\epsilon$ is a length $n$ vector of independent, zero-mean Gaussian noise with precision $\tau$.

Choosing a joint conjugate prior on the regression coefficients and the noise precision, the posterior distribution is multivariate-normal-gamma with posterior hyperparameters (see wikipedia for details)

\begin{align} \Lambda_n&=X'X + \Lambda_0\\ \mu_n&=\Lambda_n^{-1}\left(\Lambda_0\mu_0 + X^Ty\right)\\ a_n&=a_0 + \frac n2\\ b_n&=b_0 + \frac 12\left(y^Ty + \mu_0\Lambda_0\mu_0 - \mu_n\Lambda_n\mu_n\right), \end{align}

where $\Lambda$, $\mu$ denote the precision matrix and mean for $\theta$, and $a$, $b$ denote the shape and scale for $\tau$. Zero subscripts denote prior hyperparameters and $n$ subscripts denote posterior hyperparameters.

Whenever $n\gg p$m the inference works well. But when $n\gtrsim p$, the inference for the precision parameter fails and the support of the posterior is at larger values than the true value. E.g. for $n=500$ and $p=100$, I get the following results. Do you have any idea why?

The code is here.

Answer 1

This problem turns out to be well-known in the frequentist literature. In particular, if we use an impropr prior $\Lambda_0=b_0=0$, the posterior scale hyperparameter for the distribution on $\tau$ is

$$\begin{align} b_n&=\frac{1}{2}\left(y^Ty - \mu_n^T\Lambda_n\mu_n\right)\\ &=\frac 12\left(y^Ty -\mu_n^TX^Ty-y^TX\mu_n+ \mu_n^T\Lambda_n\mu_n\right)\\ &=\frac 12\left(y-X\mu_n\right)^T\left(y-X\mu_n\right), \end{align}$$

where we have used the fact that $\Lambda_n=X^TX$ and $\mu_n=\left(X^TX\right)^{-1}X^Ty$. Thus, $b_n$ is $n/2$ times the sample variance of the residuals $y-X\mu_n$. But we already used the data to estimate the regression coefficients and the sample variance of the residuals is a biased estimator of the population variance. In particular, we have $\nu=n-p$ degrees of freedom and an unbiased estimate of the population variance is $$ \frac{1}{n-p}\left(y-X\mu_n\right)^T\left(y-X\mu_n\right). $$

Whenever $n\gg p$ just using the sample variance is fine because $\frac{n}{n-p}\approx 1$. However, as soon as $p$ becomes comparable with $n$ the population variance is underestimated by the sample variance. The inference fails. The problem is discussed in the context of maximum likelihood inference on page 388 in "Data Analysis Using Regression and Multilevel/Hierarchical Models".

I like conjugate priors. But in this case they caused me some trouble. Moral of the story: don't just pick the functional form of your priors because they are conjugate.

As an aside: using a variational mean-field approximation such that the posterior factorises with respect to the regression coefficients and the noise precision works much better than the closed form solution provided by the conjugate priors (better being defined as the posterior being consistent with the true value that was used to generate the data).