Bayesian multivariate linear regression

In statistics, Bayesian multivariate linear regression is a Bayesian approach to multiple linear regression.

Details

Consider a collection of m linear regression problems for n observations, related through a set of common predictor variables $\{x_{c}\}$ , and a jointly normal errors $\{\epsilon _{c}\}$ :

y_{1}=\beta _{1}x_{1}+\epsilon _{1},\,

y_{c}=\beta _{c}x_{c}+\epsilon _{c},\,

y_{m}=\beta _{m}x_{m}+\epsilon _{m},\,

where the subscript c denotes a column vector of k observations for each measurement ( $n=k+m$ ).

The noise terms are jointly normal over each collection of k observations. That is, each row vector $\{r\}$ represents an m length vector of correlated observations on each of the dependent variables:

y_{r}=B^{T}x_{r}+\epsilon _{r},\,

where the noise $\epsilon _{r}$ is i.i.d. and normally distributed for all rows $\{r\}$ .

\epsilon _{r}\sim N(0,\Sigma _{\epsilon }^{2}).\,

where B is an $k\times m$ matrix

B=[\beta _{1},\cdots ,\beta _{c},\cdots ,\beta _{m}]\,

We can write the entire regression problem in matrix form as:

Y=XB+E,\,

where Y and E are $n\times m$ matrices.

The classical, frequentists linear least squares solution is to simply estimate the matrix of regression coefficients ${\hat {B}}$ using the Moore-Penrose pseudoinverse:

{\hat {B}}=(X^{T}X)^{-1}X^{T}Y

.

To obtain the Bayesian solution, we need to specify the conditional likelihood and then find the appropriate conjugate prior. As with the univariate case of linear Bayesian regression, we will find that we can specify a natural conditional conjugate prior (which is scale dependent).

Let us write our conditional likelihood as

\rho (E|\Sigma _{\epsilon })\propto (\Sigma _{\epsilon }^{2})^{-n/2}\exp(-{\frac {1}{2}}tr(E^{T}\Sigma _{\epsilon }^{-1}E)),\,

writing the error E in terms Y,X, and B yields

\rho (Y|X,B\Sigma _{\epsilon })\propto (\Sigma _{\epsilon }^{2})^{-n/2}\exp(-{\frac {1}{2}}tr((Y-XB)^{T}\Sigma _{\epsilon }^{-1}(Y-XB))),\,

We seek a natural conjugate prior—a joint density $\rho (B,\sigma _{\epsilon })$ which is of the same functional form as the likelihood. Since the likelihood is quadratic in $B$ , we re-write the likelihood so it is normal in $(B-{\hat {B}})$ (the deviation from classical sample estimate)

Using the same technique as with Bayesian linear regression, we decompose the exponential term using a matrix-form of the sum-of-squares technique. Here, however, we will also need to use the Matrix Differential Calculus (Kronecker product and vectorization transformations).

First, let us apply sum-of-squares to obtain new expression for the likelihood:

\rho (Y|X,B,\Sigma _{\epsilon })\propto \Sigma _{\epsilon }^{-(n-k)/2}\exp(-tr(-{\frac {1}{2}}S^{T}\Sigma _{\epsilon }^{-1}S))(\Sigma _{\epsilon }^{2})^{-k/2}\exp(-{\frac {1}{2}}tr((B-{\hat {B}})^{T}X^{T}\Sigma _{\epsilon }^{-1}X(B-{\hat {B}}))),\,

S=Y-{\hat {B}}X

We would like to develop a conditional form for the priors:

\rho (B,\Sigma _{\epsilon })=\rho (\Sigma _{\epsilon })\rho (B|\Sigma _{\epsilon }),\,

where $\rho (\Sigma _{\epsilon })$ is an inverse-Wishart distribution and $\rho (B|\Sigma _{\epsilon })$ is some form of normal distribution in the matrix $B$ . This is accomplished using the vectorization transformation, which converts the likelihood from a function of the matrices $B,{\hat {B}}$ to a function of the vectors $\beta =vec(B),{\hat {\beta }}=vec({\hat {B}})$ .