View source for Bayesian linear regression

{{Short description|Method of statistical analysis}}
{{Bayesian statistics}}
{{Regression bar}}
{{Distinguish|Bayes linear statistics}}
'''Bayesian linear regression''' is a type of [[conditional model]]ing in which the mean of one variable is described by a [[linear combination]] of other variables, with the goal of obtaining the [[posterior probability]] of the regression coefficients (as well as other parameters describing the [[probability distribution|distribution]] of the regressand) and ultimately allowing the [[out-of-sample]] prediction of the [[:wikt:regressand|regressand]] (often labelled <math>y</math>) ''[[Conditional expectation|conditional on]]'' observed values of the regressors (usually <math>X</math>). The simplest and most widely used version of this model is the ''normal linear model'', in which <math>y</math> given <math>X</math> is distributed [[normal distribution|Gaussian]]. In this model, and under a particular choice of [[prior probability|prior probabilities]] for the parameters—so-called [[conjugate prior]]s—the posterior can be found analytically. With more arbitrarily chosen priors, the posteriors generally have to be approximated.

==Model setup==
Consider a standard [[linear regression]] problem, in which for <math>i = 1, \ldots, n</math> we specify the mean of the [[Conditional probability distribution|conditional distribution]] of <math>y_i</math> given a <math>k \times 1</math> predictor vector <math>\mathbf{x}_i</math>:
<math display="block">y_{i} = \mathbf{x}_i^\mathsf{T} \boldsymbol\beta  + \varepsilon_i,</math>

where <math>\boldsymbol\beta</math> is a <math>k \times 1</math> vector, and the <math>\varepsilon_i</math> are [[i.i.d.|independent and identically]] [[normally distributed]] random variables:
<math display="block">\varepsilon_{i} \sim N(0, \sigma^2).</math>

This corresponds to the following [[likelihood function]]:

<math display="block">\rho(\mathbf{y}\mid\mathbf{X},\boldsymbol\beta,\sigma^{2}) \propto (\sigma^2)^{-n/2} \exp\left(-\frac{1}{2\sigma^2} (\mathbf{y}- \mathbf{X} \boldsymbol\beta)^\mathsf{T}(\mathbf{y}- \mathbf{X} \boldsymbol\beta)\right).</math>

The [[ordinary least squares]] solution is used to estimate the coefficient vector using the [[Moore–Penrose pseudoinverse]]:
<math display="block"> \hat{\boldsymbol\beta} = (\mathbf{X}^\mathsf{T}\mathbf{X})^{-1}\mathbf{X}^\mathsf{T}\mathbf{y}</math>

where <math>\mathbf{X}</math> is the <math>n \times k</math> [[design matrix]], each row of which is a predictor vector <math>\mathbf{x}_i^\mathsf{T}</math>; and <math>\mathbf{y}</math> is the column <math>n</math>-vector <math>[y_1 \; \cdots \; y_n]^\mathsf{T}</math>.

This is a [[frequentist]] approach, and it assumes that there are enough measurements to say something meaningful about <math>\boldsymbol\beta</math>.  In the [[Bayesian inference|Bayesian]] approach,<ref>{{cite journal | last1 = Huang | first1 = Yunfei |  last2 = Gompper | first2 =  Gerhard |  last3 = Sabass | first3 = Benedikt |year = 2020 | title = A Bayesian traction force microscopy method with automated denoising in a user-friendly software package | journal = Computer Physics Communications | volume = 256 | page = 107313 | doi = 10.1016/j.cpc.2020.107313 | arxiv = 2005.01377 | bibcode = 2020CoPhC.25607313H }}</ref> the data are supplemented with additional information in the form of a [[prior probability distribution]]. The prior belief about the parameters is combined with the data's likelihood function according to [[Bayes' theorem]] to yield the [[posterior probability|posterior belief]] about the parameters <math>\boldsymbol\beta</math> and <math>\sigma</math>. The prior can take different functional forms depending on the domain and the information that is available ''a priori''.

Since the data comprise both <math>\mathbf{y}</math> and <math>\mathbf{X}</math>, the focus only on the distribution of <math>\mathbf{y}</math> conditional on <math>\mathbf{X}</math> needs justification. In fact, a "full" Bayesian analysis would require a joint likelihood <math>\rho(\mathbf{y},\mathbf{X}\mid\boldsymbol\beta,\sigma^{2},\gamma)</math> along with a prior <math>\rho(\beta,\sigma^{2},\gamma)</math>, where <math>\gamma</math> symbolizes the parameters of the distribution for <math>\mathbf{X}</math>. Only under the assumption of (weak) exogeneity can the joint likelihood be factored into <math>\rho(\mathbf{y}\mid\boldsymbol\mathbf{X},\beta,\sigma^{2})\rho(\mathbf{X}\mid\gamma)</math>.<ref>See Jackman (2009), p. 101.</ref> The latter part is usually ignored under the assumption of disjoint parameter sets. More so, under classic assumptions <math>\mathbf{X}</math> are considered chosen (for example, in a designed experiment) and therefore has a known probability without parameters.<ref>See Gelman et al. (2013), p. 354.</ref>

==With conjugate priors==

===Conjugate prior distribution===
For an arbitrary prior distribution, there may be no analytical solution for the [[posterior distribution]]. In this section, we will consider a so-called [[conjugate prior]] for which the posterior distribution can be derived analytically.

A prior <math>\rho(\boldsymbol\beta,\sigma^{2})</math> is [[conjugate prior|conjugate]] to this likelihood function if it has the same functional form with respect to <math>\boldsymbol\beta</math> and <math>\sigma</math>. Since the log-likelihood is quadratic in <math>\boldsymbol\beta</math>, the log-likelihood is re-written such that the likelihood becomes normal in <math>(\boldsymbol\beta-\hat{\boldsymbol\beta})</math>.  Write

<math display="block">\begin{align}
(\mathbf{y}- \mathbf{X} \boldsymbol\beta)^\mathsf{T}(\mathbf{y}- \mathbf{X} \boldsymbol\beta) 
&= [(\mathbf{y}- \mathbf{X} \hat{\boldsymbol\beta}) + (\mathbf{X} \hat{\boldsymbol\beta} - \mathbf{X} \boldsymbol\beta)]^\mathsf{T} [(\mathbf{y}- \mathbf{X} \hat{\boldsymbol\beta}) + (\mathbf{X} \hat{\boldsymbol\beta} - \mathbf{X} \boldsymbol\beta)] \\
&= (\mathbf{y}- \mathbf{X} \hat{\boldsymbol\beta})^\mathsf{T}(\mathbf{y}- \mathbf{X} \hat{\boldsymbol\beta}) + (\boldsymbol\beta - \hat{\boldsymbol\beta})^\mathsf{T}(\mathbf{X}^\mathsf{T}\mathbf{X})(\boldsymbol\beta - \hat{\boldsymbol\beta})
+ \underbrace{2(\mathbf{X} \hat{\boldsymbol\beta} - \mathbf{X} \boldsymbol\beta)^\mathsf{T} (\mathbf{y}- \mathbf{X} \hat{\boldsymbol\beta})}_{= \ 0}\\
&= (\mathbf{y}- \mathbf{X} \hat{\boldsymbol\beta})^\mathsf{T}(\mathbf{y}- \mathbf{X} \hat{\boldsymbol\beta}) + (\boldsymbol\beta - \hat{\boldsymbol\beta})^\mathsf{T}(\mathbf{X}^\mathsf{T}\mathbf{X})(\boldsymbol\beta - \hat{\boldsymbol\beta})\,.
\end{align}</math>

The likelihood is now re-written as
<math display="block">\rho(\mathbf{y}|\mathbf{X},\boldsymbol\beta,\sigma^{2}) \propto (\sigma^2)^{-\frac{v}{2}} \exp\left(-\frac{vs^{2}}{2{\sigma}^{2}}\right)(\sigma^2)^{-\frac{n-v}{2}} \exp\left(-\frac{1}{2{\sigma}^{2}}(\boldsymbol\beta - \hat{\boldsymbol\beta})^\mathsf{T}(\mathbf{X}^\mathsf{T}\mathbf{X})(\boldsymbol\beta - \hat{\boldsymbol\beta})\right),</math>
where
<math display="block">vs^2 =(\mathbf{y}- \mathbf{X} \hat{\boldsymbol\beta})^\mathsf{T}(\mathbf{y}- \mathbf{X} \hat{\boldsymbol\beta}) \quad  \text{ and } \quad v = n-k,</math>
where <math>k</math> is the number of regression coefficients.

This suggests a form for the prior:
<math display="block">\rho(\boldsymbol\beta,\sigma^2) = \rho(\sigma^2)\rho(\boldsymbol\beta\mid\sigma^2),</math>
where <math>\rho(\sigma^2)</math> is an [[inverse-gamma distribution]]
<math display="block"> \rho(\sigma^2) \propto (\sigma^2)^{-\frac{v_0}{2}-1} \exp\left(-\frac{v_0 s_0^2}{2\sigma^2}\right).</math>

In the notation introduced  in the [[inverse-gamma distribution]] article, this is the density of an <math> \text{Inv-Gamma}( a_0, b_0)</math> distribution with <math>a_0=\tfrac{v_0}{2}</math> and <math>b_0=\tfrac{1}{2} v_0s_0^2 </math> with <math>v_0</math> and <math>s_0^2</math> as the prior values of <math>v</math> and <math>s^{2}</math>, respectively.  Equivalently, it can also be described as a [[scaled inverse chi-squared distribution]], <math>\text{Scale-inv-}\chi^2(v_0, s_0^2).</math>

Further the conditional prior density <math>\rho(\boldsymbol\beta|\sigma^{2})</math> is a [[normal distribution]],

<math display="block"> \rho(\boldsymbol\beta\mid\sigma^2) \propto (\sigma^2)^{-k/2} \exp\left(-\frac{1}{2\sigma^2}(\boldsymbol\beta - \boldsymbol\mu_0)^\mathsf{T} \mathbf{\Lambda}_0 (\boldsymbol\beta - \boldsymbol\mu_0)\right).</math>

In the notation of the [[normal distribution]], the conditional prior distribution is <math> \mathcal{N}\left(\boldsymbol\mu_0, \sigma^2 \boldsymbol\Lambda_0^{-1}\right).</math>

===Posterior distribution===

With the prior now specified, the posterior distribution can be expressed as

<math display="block"> \begin{align}
\rho(\boldsymbol\beta,\sigma^2\mid\mathbf{y},\mathbf{X}) &\propto \rho(\mathbf{y}\mid\mathbf{X},\boldsymbol\beta,\sigma^2)\rho(\boldsymbol\beta\mid\sigma^2)\rho(\sigma^2)  \\
& \propto (\sigma^2)^{-n/2} \exp\left(-\frac{1}{2{\sigma}^2}(\mathbf{y}- \mathbf{X} \boldsymbol\beta)^\mathsf{T}(\mathbf{y}- \mathbf{X} \boldsymbol\beta)\right) (\sigma^2)^{-k/2} \exp\left(-\frac{1}{2\sigma^2}(\boldsymbol\beta -\boldsymbol\mu_0)^\mathsf{T} \boldsymbol\Lambda_0 (\boldsymbol\beta - \boldsymbol\mu_0)\right)  (\sigma^2)^{-(a_0+1)} \exp\left(-\frac{b_0}{\sigma^2}\right)
\end{align}</math>

With some re-arrangement,<ref>The intermediate steps of this computation can be found in O'Hagan (1994) at the beginning of the chapter on Linear models.</ref> the posterior can be re-written so that the posterior mean <math>\boldsymbol\mu_n</math> of the parameter vector <math>\boldsymbol\beta</math> can be expressed in terms of the least squares estimator <math>\hat{\boldsymbol\beta}</math> and the prior mean <math>\boldsymbol\mu_0</math>, with the strength of the prior indicated by the prior precision matrix <math>\boldsymbol\Lambda_0</math>

<math display="block">\boldsymbol\mu_n = (\mathbf{X}^\mathsf{T}\mathbf{X}+\boldsymbol\Lambda_0)^{-1}(\mathbf{X}^\mathsf{T} \mathbf{X}\hat{\boldsymbol\beta}+\boldsymbol\Lambda_0\boldsymbol\mu_0) .</math>

To justify that <math>\boldsymbol\mu_n</math> is indeed the posterior mean, the quadratic terms in the exponential can be re-arranged as a [[quadratic form (statistics)|quadratic form]] in <math>\boldsymbol\beta - \boldsymbol\mu_n</math>.<ref>The intermediate steps are in Fahrmeir et al. (2009) on page 188.</ref>

<math display="block"> (\mathbf{y}- \mathbf{X} \boldsymbol\beta)^\mathsf{T}(\mathbf{y}- \mathbf{X} \boldsymbol\beta) + (\boldsymbol\beta - \boldsymbol\mu_0)^\mathsf{T}\boldsymbol\Lambda_0(\boldsymbol\beta - \boldsymbol\mu_0) =(\boldsymbol\beta-\boldsymbol\mu_n)^\mathsf{T}(\mathbf{X}^\mathsf{T}\mathbf{X}+\boldsymbol\Lambda_0)(\boldsymbol\beta-\boldsymbol\mu_n)+\mathbf{y}^\mathsf{T}\mathbf{y}-\boldsymbol\mu_n^\mathsf{T}(\mathbf{X}^\mathsf{T}\mathbf{X}+\boldsymbol\Lambda_0)\boldsymbol\mu_n+\boldsymbol\mu_0^\mathsf{T} \boldsymbol\Lambda_0\boldsymbol\mu_0 .</math>

Now the posterior can be expressed as a [[normal distribution]] times an [[inverse-gamma distribution]]:

<math display="block">\rho(\boldsymbol\beta,\sigma^2\mid\mathbf{y},\mathbf{X}) \propto (\sigma^2)^{-k/2} \exp\left(-\frac{1}{2{\sigma}^{2}}(\boldsymbol\beta - \boldsymbol\mu_n)^\mathsf{T}(\mathbf{X}^\mathsf{T} \mathbf{X}+\mathbf{\Lambda}_0)(\boldsymbol\beta - \boldsymbol\mu_n)\right) (\sigma^2)^{-\frac{n+2a_0}{2}-1} \exp\left(-\frac{2 b_0+\mathbf{y}^\mathsf{T}\mathbf{y}-\boldsymbol\mu_n^\mathsf{T}(\mathbf{X}^\mathsf{T} \mathbf{X}+\boldsymbol\Lambda_0)\boldsymbol\mu_n+\boldsymbol\mu_0^\mathsf{T} \boldsymbol\Lambda_0 \boldsymbol\mu_0}{2\sigma^2}\right) .</math>

Therefore, the posterior distribution can be parametrized as follows.
<math display="block">\rho(\boldsymbol\beta,\sigma^2\mid\mathbf{y},\mathbf{X}) \propto  \rho(\boldsymbol\beta \mid \sigma^2,\mathbf{y},\mathbf{X}) \rho(\sigma^2\mid\mathbf{y},\mathbf{X}), </math>
where the two factors correspond to the densities of <math> \mathcal{N}\left( \boldsymbol\mu_n, \sigma^2\boldsymbol\Lambda_n^{-1} \right)\,</math> and <math> \text{Inv-Gamma}\left(a_n,b_n \right) </math> distributions, with the parameters of these given by

<math display="block">\boldsymbol\Lambda_n=(\mathbf{X}^\mathsf{T}\mathbf{X}+\mathbf{\Lambda}_0), \quad \boldsymbol\mu_n = (\boldsymbol\Lambda_n)^{-1}(\mathbf{X}^\mathsf{T} \mathbf{X} \hat{\boldsymbol\beta} + \boldsymbol\Lambda_0 \boldsymbol\mu_0) ,</math>
<math display="block">a_n= a_0 + \frac{n}{2}, \qquad b_n=b_0+\frac{1}{2}(\mathbf{y}^\mathsf{T} \mathbf{y} + \boldsymbol\mu_0^\mathsf{T} \boldsymbol\Lambda_0\boldsymbol\mu_0-\boldsymbol\mu_n^\mathsf{T} \boldsymbol\Lambda_n \boldsymbol\mu_n) .</math>

which illustrates Bayesian inference being a compromise between the information contained in the prior and the information contained in the sample.

===Model evidence===
The [[model evidence]] <math>p(\mathbf{y}\mid m)</math> is the probability of the data given the model <math>m</math>. It is also known as the [[marginal likelihood]], and as the ''prior predictive density''. Here, the model is defined by the likelihood function <math>p(\mathbf{y}\mid\mathbf{X},\boldsymbol\beta,\sigma)</math> and the prior distribution on the parameters, i.e. <math>p(\boldsymbol\beta,\sigma)</math>. The model evidence captures in a single number how well such a model explains the observations. The model evidence of the Bayesian linear regression model presented in this section can be used to compare competing linear models by [[Bayes factors]]. These models may differ in the number and values of the predictor variables as well as in their priors on the model parameters. Model complexity is already taken into account by the model evidence, because it marginalizes out the parameters by integrating <math>p(\mathbf{y},\boldsymbol\beta,\sigma\mid\mathbf{X})</math> over all possible values of <math>\boldsymbol\beta</math> and <math>\sigma</math>.
<math display="block">p(\mathbf{y}|m)=\int p(\mathbf{y}\mid\mathbf{X},\boldsymbol\beta,\sigma)\, p(\boldsymbol\beta,\sigma)\, d\boldsymbol\beta\, d\sigma</math>
This integral can be computed analytically and the solution is given in the following equation.<ref>The intermediate steps of this computation can be found in O'Hagan (1994) on page 257.</ref>
<math display="block">p(\mathbf{y}\mid m)=\frac{1}{(2\pi)^{n/2}}\sqrt{\frac{\det(\boldsymbol\Lambda_0)}{\det(\boldsymbol\Lambda_n)}} \cdot \frac{b_0^{a_0}}{b_n^{a_n}} \cdot \frac{\Gamma(a_n)}{\Gamma(a_0)}</math>

Here <math>\Gamma</math> denotes the [[gamma function]]. Because we have chosen a conjugate prior, the marginal likelihood can also be easily computed by evaluating the following equality for arbitrary values of <math>\boldsymbol\beta</math> and <math>\sigma</math>.<ref>{{cite journal |last=Chib |first=Siddhartha |year=1995 |title=Marginal Likelihood from the Gibbs Output |journal=Journal of the American Statistical Association |volume=90 |issue=432 |pages=1313–1321 |doi=10.2307/2291521 }}</ref>
<math display="block">p(\mathbf{y}\mid m)=\frac{p(\boldsymbol\beta,\sigma|m)\, p(\mathbf{y} \mid \mathbf{X}, \boldsymbol\beta,\sigma,m)}{p(\boldsymbol\beta, \sigma \mid \mathbf{y},\mathbf{X},m)}</math>
Note that this equation follows from a re-arrangement of [[Bayes' theorem]]. Inserting the formulas for the prior, the likelihood, and the posterior and simplifying the resulting expression leads to the analytic expression given above.

==Other cases==
In general, it may be impossible or impractical to derive the posterior distribution analytically. However, it is possible to approximate the posterior by an [[approximate Bayesian computation|approximate Bayesian inference]] method such as [[Monte Carlo sampling]],<ref>Carlin and Louis (2008) and Gelman, et al. (2003) explain how to use sampling methods for Bayesian linear regression.</ref> [[Integrated nested Laplace approximations |INLA]] or [[variational Bayes]].

The special case <math>\boldsymbol\mu_0=0, \mathbf{\Lambda}_0 = c\mathbf{I}</math> is called [[ridge regression]].

A similar analysis can be performed for the general case of the multivariate regression and part of this provides for Bayesian [[estimation of covariance matrices]]: see [[Bayesian multivariate linear regression]].

==See also==
* [[Bayes linear statistics]]
* [[Constrained least squares]]
* [[Regularized least squares]]
* [[Tikhonov regularization]]
* [[Spike and slab variable selection]]
* [[Bayesian interpretation of kernel regularization]]

{{More footnotes|date=August 2011}}

==Notes==
{{Reflist}}

==References==
* {{cite book |author-link1=George E. P. Box |last1=Box |first1=G. E. P. |author2-link=George Tiao|last2=Tiao |first2=G. C. |year=1973 |title=Bayesian Inference in Statistical Analysis |publisher=Wiley |isbn=0-471-57428-7 }}
* {{cite book|author=Carlin, Bradley P. |author2=Louis, Thomas A.| title = Bayesian Methods for Data Analysis |edition=Third | publisher = Boca Raton, FL: Chapman and Hall/CRC | year = 2008 | isbn = 978-1-58488-697-6}}
* {{cite book|author=Fahrmeir, L. |author2=Kneib, T. |author3=Lang, S.| title = Regression. Modelle, Methoden und Anwendungen |edition=Second |publisher = Springer |location=Heidelberg | year = 2009 | isbn = 978-3-642-01836-7 | doi = 10.1007/978-3-642-01837-4}}
* {{cite book |last1=Gelman |first1=Andrew |author-link1=Andrew Gelman |last2=Carlin |first2=John B. |last3=Stern |first3=Hal S. |last4=Dunson |first4=David B. |last5=Vehtari |first5=Aki |last6=Rubin |first6=Donald B. |display-authors=1 |chapter=Introduction to regression models |pages=353–380 |title = Bayesian Data Analysis |edition=Third | publisher = Boca Raton, FL: Chapman and Hall/CRC | year = 2013 | isbn = 978-1-4398-4095-5 }}
* {{cite book |first=Simon |last=Jackman |chapter=Regression models |pages=99–124 |title=Bayesian Analysis for the Social Sciences |location= |publisher=Wiley |year=2009 |isbn=978-0-470-01154-6 }}
* {{cite book |first1=Peter E. |last1=Rossi |first2=Greg M. |last2=Allenby |first3=Robert |last3=McCulloch |title=Bayesian Statistics and Marketing |publisher=John Wiley & Sons |year=2006 |isbn=0470863676 }}
* {{cite book|author=O'Hagan, Anthony| title = Bayesian Inference| volume= 2B |series = Kendall's Advanced Theory of Statistics| year = 1994 | edition= First | isbn = 0-340-52922-9| publisher = Halsted}}

==External links==
* [[b:en:R Programming/Linear Models#Bayesian estimation|Bayesian estimation of linear models (R programming wikibook)]]. Bayesian linear regression as implemented in [[R (programming language)|R]].

{{Least Squares and Regression Analysis}}
{{Statistics|correlation}}

[[Category:Bayesian inference|Linear regression]]
[[Category:Single-equation methods (econometrics)]]