Instatistics,Bayesian multivariate linear regression is aBayesian approach tomultivariate linear regression, i.e.linear regression where the predicted outcome is a vector of correlatedrandom variables rather than a single scalar random variable. A more general treatment of this approach can be found in the articleMMSE estimator.

Consider a regression problem where thedependent variable to be predicted is not a singlereal-valued scalar but anm-length vector of correlated real numbers. As in the standard regression setup, there aren observations, where each observationi consists ofk−1explanatory variables, grouped into a vector${\displaystyle {\mathbf{x}}_i}$ of lengthk (where adummy variable with a value of 1 has been added to allow for an intercept coefficient). This can be viewed as a set ofm related regression problems for each observationi:where the set of errors${\displaystyle \{{ϵ}_{i,1},\dots ,{ϵ}_{i,m}\}}$ are all correlated. Equivalently, it can be viewed as a single regression problem where the outcome is arow vector${\displaystyle {\mathbf{y}}_i^{\mathsf{T}}}$ and the regression coefficient vectors are stacked next to each other, as follows:$${\displaystyle {\mathbf{y}}_i^{\mathsf{T}}={\mathbf{x}}_i^{\mathsf{T}}\mathbf{B}+{\mathit{ϵ}}_i^{\mathsf{T}}.}$$

Thecoefficient matrixB is a${\displaystyle k\times m}$ matrix where the coefficient vectors${\displaystyle {\mathit{\beta }}_1,\dots ,{\mathit{\beta }}_m}$ for each regression problem are stacked horizontally:$${\displaystyle \mathbf{B}=\left[\begin{array}{c}\left(\begin{array}{c}\\ {\mathit{\beta }}_1\\ \end{array}\right)\cdots \left(\begin{array}{c}\\ {\mathit{\beta }}_m\\ \end{array}\right)\end{array}\right]=\left[\begin{array}{c}\left(\begin{array}{c}{\beta }_{1,1}\\ ⋮\\ {\beta }_{k,1}\end{array}\right)\cdots \left(\begin{array}{c}{\beta }_{1,m}\\ ⋮\\ {\beta }_{k,m}\end{array}\right)\end{array}\right].}$$

The noise vector${\displaystyle {\mathit{ϵ}}_i}$ for each observationi isjointly normal, so that the outcomes for a given observation are correlated:$${\displaystyle {\mathit{ϵ}}_i\sim N(0,{\mathbf{\Sigma }}_{ϵ}).}$$

We can write the entire regression problem in matrix form as:$${\displaystyle \mathbf{Y}=\mathbf{X}\mathbf{B}+\mathbf{E},}$$whereY andE are${\displaystyle n\times m}$ matrices. Thedesign matrixX is an${\displaystyle n\times k}$ matrix with the observations stacked vertically, as in the standardlinear regression setup:

The classical, frequentistslinear least squares solution is to simply estimate the matrix of regression coefficients${\displaystyle \hat{\mathbf{B}}}$ using theMoore-Penrosepseudoinverse:$${\displaystyle \hat{\mathbf{B}}=({\mathbf{X}}^{\mathsf{T}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathsf{T}}\mathbf{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 oflinear 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^[1]$${\displaystyle \rho (\mathbf{E}|{\mathbf{\Sigma }}_{ϵ})\propto |{\mathbf{\Sigma }}_{ϵ}{|}^{-n/2}\exp \left(-\frac{1}{2}\mathrm{tr}\left({\mathbf{E}}^{\mathsf{T}}{\mathbf{\Sigma }}_{ϵ}^{-1}\mathbf{E}\right)\right),}$$writing the error${\displaystyle \mathbf{E}}$ in terms of${\displaystyle \mathbf{Y},\mathbf{X},}$ and${\displaystyle \mathbf{B}}$ yields$${\displaystyle \rho (\mathbf{Y}|\mathbf{X},\mathbf{B},{\mathbf{\Sigma }}_{ϵ})\propto |{\mathbf{\Sigma }}_{ϵ}{|}^{-n/2}\exp (-\frac{1}{2}\mathrm{tr}((\mathbf{Y}-\mathbf{X}\mathbf{B}{)}^{\mathsf{T}}{\mathbf{\Sigma }}_{ϵ}^{-1}(\mathbf{Y}-\mathbf{X}\mathbf{B}))),}$$

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

Using the same technique as withBayesian 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 andvectorization transformations).

First, let us apply sum-of-squares to obtain new expression for the likelihood:$${\displaystyle \rho (\mathbf{Y}|\mathbf{X},\mathbf{B},{\mathbf{\Sigma }}_{ϵ})\propto |{\mathbf{\Sigma }}_{ϵ}{|}^{-(n-k)/2}\exp (-\mathrm{tr}(\frac{1}{2}{\mathbf{S}}^{\mathsf{T}}\mathbf{S}{\mathbf{\Sigma }}_{ϵ}^{-1}))|{\mathbf{\Sigma }}_{ϵ}{|}^{-k/2}\exp (-\frac{1}{2}\mathrm{tr}((\mathbf{B}-\hat{\mathbf{B}}{)}^{\mathsf{T}}{\mathbf{X}}^{\mathsf{T}}\mathbf{X}(\mathbf{B}-\hat{\mathbf{B}}){\mathbf{\Sigma }}_{ϵ}^{-1})),}$$$${\displaystyle \mathbf{S}=\mathbf{Y}-\mathbf{X}\hat{\mathbf{B}}}$$

We would like to develop a conditional form for the priors:$${\displaystyle \rho (\mathbf{B},{\mathbf{\Sigma }}_{ϵ})=\rho ({\mathbf{\Sigma }}_{ϵ})\rho (\mathbf{B}|{\mathbf{\Sigma }}_{ϵ}),}$$where${\displaystyle \rho ({\mathbf{\Sigma }}_{ϵ})}$ is aninverse-Wishart distributionand${\displaystyle \rho (\mathbf{B}|{\mathbf{\Sigma }}_{ϵ})}$ is some form ofnormal distribution in the matrix${\displaystyle \mathbf{B}}$. This is accomplished using thevectorization transformation, which converts the likelihood from a function of the matrices${\displaystyle \mathbf{B},\hat{\mathbf{B}}}$ to a function of the vectors${\displaystyle \mathit{\beta }=\mathrm{vec}(\mathbf{B}),\hat{\mathit{\beta }}=\mathrm{vec}(\hat{\mathbf{B}})}$.

Write$${\displaystyle \mathrm{tr}((\mathbf{B}-\hat{\mathbf{B}}{)}^{\mathsf{T}}{\mathbf{X}}^{\mathsf{T}}\mathbf{X}(\mathbf{B}-\hat{\mathbf{B}}){\mathbf{\Sigma }}_{ϵ}^{-1})=\mathrm{vec}(\mathbf{B}-\hat{\mathbf{B}}{)}^{\mathsf{T}}\mathrm{vec}({\mathbf{X}}^{\mathsf{T}}\mathbf{X}(\mathbf{B}-\hat{\mathbf{B}}){\mathbf{\Sigma }}_{ϵ}^{-1})}$$

Let$${\displaystyle \mathrm{vec}({\mathbf{X}}^{\mathsf{T}}\mathbf{X}(\mathbf{B}-\hat{\mathbf{B}}){\mathbf{\Sigma }}_{ϵ}^{-1})=({\mathbf{\Sigma }}_{ϵ}^{-1}\otimes {\mathbf{X}}^{\mathsf{T}}\mathbf{X})\mathrm{vec}(\mathbf{B}-\hat{\mathbf{B}}),}$$where${\displaystyle \mathbf{A}\otimes \mathbf{B}}$ denotes theKronecker product of matricesA andB, a generalization of theouter product which multiplies an${\displaystyle m\times n}$ matrix by a${\displaystyle p\times q}$ matrix to generate an${\displaystyle mp\times nq}$ matrix, consisting of every combination of products of elements from the two matrices.

Thenwhich will lead to a likelihood which is normal in${\displaystyle (\mathit{\beta }-\hat{\mathit{\beta }})}$.

With the likelihood in a more tractable form, we can now find a natural (conditional) conjugate prior.

The natural conjugate prior using the vectorized variable${\displaystyle \mathit{\beta }}$ is of the form:^[1]$${\displaystyle \rho (\mathit{\beta },{\mathbf{\Sigma }}_{ϵ})=\rho ({\mathbf{\Sigma }}_{ϵ})\rho (\mathit{\beta }|{\mathbf{\Sigma }}_{ϵ}),}$$where$${\displaystyle \rho ({\mathbf{\Sigma }}_{ϵ})\sim {\mathcal{W}}^{-1}({\mathbf{V}}_0,{\mathit{\nu }}_0)}$$and$${\displaystyle \rho (\mathit{\beta }|{\mathbf{\Sigma }}_{ϵ})\sim N({\mathit{\beta }}_0,{\mathbf{\Sigma }}_{ϵ}\otimes {\mathbf{\Lambda }}_0^{-1}).}$$

Using the above prior and likelihood, the posterior distribution can be expressed as:^[1]where${\displaystyle \mathrm{vec}({\mathbf{B}}_0)={\mathit{\beta }}_0}$.The terms involving${\displaystyle \mathbf{B}}$ can be grouped (with${\displaystyle {\mathbf{\Lambda }}_0={\mathbf{U}}^{\mathsf{T}}\mathbf{U}}$) using:with$${\displaystyle {\mathbf{B}}_n={\left({\mathbf{X}}^{\mathsf{T}}\mathbf{X}+{\mathbf{\Lambda }}_0\right)}^{-1}\left({\mathbf{X}}^{\mathsf{T}}\mathbf{X}\hat{\mathbf{B}}+{\mathbf{\Lambda }}_0{\mathbf{B}}_0\right)={\left({\mathbf{X}}^{\mathsf{T}}\mathbf{X}+{\mathbf{\Lambda }}_0\right)}^{-1}\left({\mathbf{X}}^{\mathsf{T}}\mathbf{Y}+{\mathbf{\Lambda }}_0{\mathbf{B}}_0\right).}$$

This now allows us to write the posterior in a more useful form:

This takes the form of aninverse-Wishart distribution times aMatrix normal distribution:$${\displaystyle \rho ({\mathbf{\Sigma }}_{ϵ}|\mathbf{Y},\mathbf{X})\sim {\mathcal{W}}^{-1}({\mathbf{V}}_n,{\mathit{\nu }}_n)}$$and$${\displaystyle \rho (\mathbf{B}|\mathbf{Y},\mathbf{X},{\mathbf{\Sigma }}_{ϵ})\sim {\mathcal{M}\mathcal{N}}_{k,m}({\mathbf{B}}_n,{\mathbf{\Lambda }}_n^{-1},{\mathbf{\Sigma }}_{ϵ}).}$$

The parameters of this posterior are given by:$${\displaystyle {\mathbf{V}}_n={\mathbf{V}}_0+(\mathbf{Y}-\mathbf{X}{\mathbf{B}}_{\mathbf{n}}{)}^{\mathsf{T}}(\mathbf{Y}-\mathbf{X}{\mathbf{B}}_{\mathbf{n}})+({\mathbf{B}}_n-{\mathbf{B}}_0{)}^{\mathsf{T}}{\mathbf{\Lambda }}_0({\mathbf{B}}_n-{\mathbf{B}}_0)}$$$${\displaystyle {\mathit{\nu }}_n={\mathit{\nu }}_0+n}$$$${\displaystyle {\mathbf{B}}_n=({\mathbf{X}}^{\mathsf{T}}\mathbf{X}+{\mathbf{\Lambda }}_0{)}^{-1}({\mathbf{X}}^{\mathsf{T}}\mathbf{Y}+{\mathbf{\Lambda }}_0{\mathbf{B}}_0)}$$$${\displaystyle {\mathbf{\Lambda }}_n={\mathbf{X}}^{\mathsf{T}}\mathbf{X}+{\mathbf{\Lambda }}_0}$$

• Bayesian linear regression
• Matrix normal distribution

This article includes a list ofgeneral references, butit lacks sufficient correspondinginline citations. Please help toimprove this article byintroducing more precise citations.(November 2010) (Learn how and when to remove this message)

• Box, G. E. P.;Tiao, G. C. (1973). "8".Bayesian Inference in Statistical Analysis. Wiley.ISBN 0-471-57428-7.
• Geisser, S. (1965)."Bayesian Estimation in Multivariate Analysis".The Annals of Mathematical Statistics.36 (1):150–159.doi:10.1214/aoms/1177700279.JSTOR 2238083.
• Tiao, G. C.; Zellner, A. (1964). "On the Bayesian Estimation of Multivariate Regression".Journal of the Royal Statistical Society. Series B (Methodological).26 (2):277–285.doi:10.1111/j.2517-6161.1964.tb00560.x.JSTOR 2984424.

• Bayesian inference
• Single-equation methods (econometrics)

• Articles with short description
• Short description matches Wikidata
• Articles lacking in-text citations from November 2010
• All articles lacking in-text citations