Bayesian linear regression

Bayesian statistics
Part of a series on Statistics

Theory
Admissible decision rule Bayesian efficiency Bayesian probability Probability interpretations Bayes' theorem Bayes' rule Bayes factor Bayesian inference Bayesian network Prior Posterior Likelihood Conjugate prior Posterior predictive Hyperparameter Hyperprior Principle of indifference Principle of maximum entropy Empirical Bayes method Cromwell's rule Bernstein–von Mises theorem Bayesian information criterion Credible interval Maximum a posteriori estimation
Techniques
Bayesian linear regression Bayesian estimator Approximate Bayesian computation
Statistics portal

Regression analysis
Part of a series on Statistics

Models
Linear regression Simple regression Ordinary least squares Polynomial regression General linear model
Generalized linear model Discrete choice Logistic regression Multinomial logit Mixed logit Probit Multinomial probit Ordered logit Ordered probit Poisson
Multilevel model Fixed effects Random effects Mixed model
Nonlinear regression Nonparametric Semiparametric Robust Quantile Isotonic Principal components Least angle Local Segmented
Errors-in-variables
Estimation
Least squares Ordinary least squares Linear (math) Partial Total Generalized Weighted Non-linear Non-negative Iteratively reweighted Ridge regression
Least absolute deviations Bayesian Bayesian multivariate
Background
Regression model validation Mean and predicted response Errors and residuals Goodness of fit Studentized residual Gauss–Markov theorem
Statistics portal

In statistics, Bayesian linear regression is an approach to linear regression in which the statistical analysis is undertaken within the context of Bayesian inference. When the regression model has errors that have a normal distribution, and if a particular form of prior distribution is assumed, explicit results are available for the posterior probability distributions of the model's parameters.

Model setup

Consider a standard linear regression problem, in which for $i=1,...,n$ we specify the conditional distribution of $y_{i}$ given a $k \times 1$ predictor vector $\mathbf {x} _{i}$ :

y_{i} = \mathbf{x}_{i}^{\rm T} \boldsymbol\beta + \epsilon_{i},

where ${\boldsymbol {\beta }}$ is a $k \times 1$ vector, and the $\epsilon _{i}$ are independent and identical normally distributed random variables:

\epsilon_{i} \sim N(0, \sigma^2).

This corresponds to the following likelihood function:

\rho(\mathbf{y}|\mathbf{X},\boldsymbol\beta,\sigma^{2}) \propto (\sigma^{2})^{-n/2} \exp\left(-\frac{1}{2{\sigma}^{2}}(\mathbf{y}- \mathbf{X} \boldsymbol\beta)^{\rm T}(\mathbf{y}- \mathbf{X} \boldsymbol\beta)\right).

The ordinary least squares solution is to estimate the coefficient vector using the Moore-Penrose pseudoinverse:

\hat{\boldsymbol\beta} = (\mathbf{X}^{\rm T}\mathbf{X})^{-1}\mathbf{X}^{\rm T}\mathbf{y}

where $\mathbf {X}$ is the $n \times k$ design matrix, each row of which is a predictor vector $\mathbf{x}_{i}^{\rm T}$ ; and $\mathbf {y}$ is the column $n$ -vector $[y_1 \; \cdots \; y_n]^{\rm T}$ .

This is a frequentist approach, and it assumes that there are enough measurements to say something meaningful about ${\boldsymbol {\beta }}$ . In the Bayesian approach, 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 belief about the parameters ${\boldsymbol {\beta }}$ and $\sigma$ . The prior can take different functional forms depending on the domain and the information that is available a priori.

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

(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})=(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})+({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} )({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}}).

The likelihood is now re-written as

\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 }}})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} )({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})\right),

where

vs^{2}=(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})\quad {\text{ and }}\quad v=n-k,

where $k$ is the number of regression coefficients.

This suggests a form for the prior:

\rho(\boldsymbol\beta,\sigma^{2}) = \rho(\sigma^{2})\rho(\boldsymbol\beta|\sigma^{2}),

where $\rho(\sigma^{2})$ is an inverse-gamma distribution

\rho (\sigma ^{2})\propto (\sigma ^{2})^{-{\frac {v_{0}}{2}}-1}\exp \left(-{\frac {v_{0}s_{0}^{2}}{2{\sigma }^{2}}}\right).

In the notation introduced in the inverse-gamma distribution article, this is the density of an $\text{Inv-Gamma}( a_0,b_0)$ distribution with $a_{0}={\tfrac {v_{0}}{2}}$ and $b_{0}={\tfrac {1}{2}}v_{0}s_{0}^{2}$ with $v_{0}$ and $s_{0}^{2}$ as the prior values of $v$ and $s^{2}$ , respectively. Equivalently, it can also be described as a scaled inverse chi-squared distribution, $\mbox{Scale-inv-}\chi^2(v_0, s_0^2).$

Further the conditional prior density $\rho(\boldsymbol\beta|\sigma^{2})$ is a normal distribution,

\rho ({\boldsymbol {\beta }}|\sigma ^{2})\propto (\sigma ^{2})^{-{\frac {k}{2}}}\exp \left(-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})^{\rm {T}}\mathbf {\Lambda } _{0}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})\right).

In the notation of the normal distribution, the conditional prior distribution is $\mathcal{N}\left(\boldsymbol\mu_0, \sigma^2\mathbf{\Lambda}_0^{-1}\right).$

Posterior distribution

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

{\begin{aligned}\rho ({\boldsymbol {\beta }},\sigma ^{2}|\mathbf {y} ,\mathbf {X} )&\propto \rho (\mathbf {y} |\mathbf {X} ,{\boldsymbol {\beta }},\sigma ^{2})\rho ({\boldsymbol {\beta }}|\sigma ^{2})\rho (\sigma ^{2})\\&\propto (\sigma ^{2})^{-{\frac {n}{2}}}\exp \left(-{\frac {1}{2{\sigma }^{2}}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})\right)(\sigma ^{2})^{-{\frac {k}{2}}}\exp \left(-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})^{\rm {T}}{\boldsymbol {\Lambda }}_{0}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})\right)(\sigma ^{2})^{-(a_{0}+1)}\exp \left(-{\frac {b_{0}}{{\sigma }^{2}}}\right)\end{aligned}}

With some re-arrangement,^[1] the posterior can be re-written so that the posterior mean $\boldsymbol\mu_n$ of the parameter vector ${\boldsymbol {\beta }}$ can be expressed in terms of the least squares estimator $\hat{\boldsymbol\beta}$ and the prior mean ${\boldsymbol {\mu }}_{0}$ , with the strength of the prior indicated by the prior precision matrix $\boldsymbol\Lambda_0$

\boldsymbol\mu_n = (\mathbf{X}^{\rm T}\mathbf{X}+\boldsymbol\Lambda_0)^{-1}(\mathbf{X}^{\rm T}\mathbf{X}\hat{\boldsymbol\beta}+\boldsymbol\Lambda_0\boldsymbol\mu_0) .

To justify that $\boldsymbol\mu_n$ is indeed the posterior mean, the quadratic terms in the exponential can be re-arranged as a quadratic form in $\boldsymbol\beta-\boldsymbol\mu_n$ .^[2]

(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})+({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})^{\rm {T}}{\boldsymbol {\Lambda }}_{0}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})=({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})+\mathbf {y} ^{\rm {T}}\mathbf {y} -{\boldsymbol {\mu }}_{n}^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}){\boldsymbol {\mu }}_{n}+{\boldsymbol {\mu }}_{0}^{\rm {T}}{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}.

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

\rho ({\boldsymbol {\beta }},\sigma ^{2}|\mathbf {y} ,\mathbf {X} )\propto (\sigma ^{2})^{-{\frac {k}{2}}}\exp \left(-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +\mathbf {\Lambda } _{0})({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})\right)(\sigma ^{2})^{-{\frac {n+2a_{0}}{2}}-1}\exp \left(-{\frac {2b_{0}+\mathbf {y} ^{\rm {T}}\mathbf {y} -{\boldsymbol {\mu }}_{n}^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}){\boldsymbol {\mu }}_{n}+{\boldsymbol {\mu }}_{0}^{\rm {T}}{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}}{2{\sigma }^{2}}}\right).

Therefore, the posterior distribution can be parametrized as follows.

\rho(\boldsymbol\beta,\sigma^{2}|\mathbf{y},\mathbf{X}) \propto \rho(\boldsymbol\beta |\sigma^{2},\mathbf{y},\mathbf{X}) \rho(\sigma^{2}|\mathbf{y},\mathbf{X}),

where the two factors correspond to the densities of $\mathcal{N}\left( \boldsymbol\mu_n, \sigma^2\boldsymbol\Lambda_n^{-1} \right)\,$ and $\text{Inv-Gamma}\left(a_n,b_n \right)$ distributions, with the parameters of these given by

\boldsymbol\Lambda_n=(\mathbf{X}^{\rm T}\mathbf{X}+\mathbf{\Lambda}_0), \quad \boldsymbol\mu_n = (\boldsymbol\Lambda_n)^{-1}(\mathbf{X}^{\rm T}\mathbf{X}\hat{\boldsymbol\beta}+\boldsymbol\Lambda_0\boldsymbol\mu_0) ,

a_{n}=a_{0}+{\frac {n}{2}},\qquad b_{n}=b_{0}+{\frac {1}{2}}({\mathbf {y}}^{{{\rm {T}}}}{\mathbf {y}}+{\boldsymbol \mu }_{0}^{{{\rm {T}}}}{\boldsymbol \Lambda }_{0}{\boldsymbol \mu }_{0}-{\boldsymbol \mu }_{n}^{{{\rm {T}}}}{\boldsymbol \Lambda }_{n}{\boldsymbol \mu }_{n}).

This can be interpreted as Bayesian learning where the parameters are updated according to the following equations.

\boldsymbol\mu_n=(\mathbf{X}^{\rm T}\mathbf{X}+\boldsymbol\Lambda_0)^{-1} (\boldsymbol\Lambda_0\boldsymbol\mu_0+\mathbf{X}^{\rm T}\mathbf{X}\hat{\boldsymbol\beta})=(\mathbf{X}^{\rm T}\mathbf{X}+\boldsymbol\Lambda_0)^{-1} (\boldsymbol\Lambda_0\boldsymbol\mu_0+\mathbf{X}^{\rm T}\mathbf{y}) ,

\boldsymbol\Lambda_n=(\mathbf{X}^{\rm T}\mathbf{X}+\boldsymbol\Lambda_0) ,

a_{n}=a_{0}+{\frac {n}{2}},

b_{n}=b_{0}+{\frac {1}{2}}({\mathbf {y}}^{{{\rm {T}}}}{\mathbf {y}}+{\boldsymbol \mu }_{0}^{{{\rm {T}}}}{\boldsymbol \Lambda }_{0}{\boldsymbol \mu }_{0}-{\boldsymbol \mu }_{n}^{{{\rm {T}}}}{\boldsymbol \Lambda }_{n}{\boldsymbol \mu }_{n}).

Model evidence

The model evidence $p(\mathbf{y}|m)$ is the probability of the data given the model $m$ . It is also known as the marginal likelihood, and as the prior predictive density. Here, the model is defined by the likelihood function $p(\mathbf{y}|\mathbf{X},\boldsymbol\beta,\sigma)$ and the prior distribution on the parameters, i.e. $p(\boldsymbol\beta,\sigma)$ . 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 Bayesian model comparison. 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 $p(\mathbf{y},\boldsymbol\beta,\sigma|\mathbf{X})$ over all possible values of ${\boldsymbol {\beta }}$ and $\sigma$ .

p(\mathbf{y}|m)=\int p(\mathbf{y}|\mathbf{X},\boldsymbol\beta,\sigma)\, p(\boldsymbol\beta,\sigma)\, d\boldsymbol\beta\, d\sigma

This integral can be computed analytically and the solution is given in the following equation.^[3]

p(\mathbf{y}|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)}

Here $\Gamma$ 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 ${\boldsymbol {\beta }}$ and $\sigma$ .

p(\mathbf{y}|m)=\frac{p(\boldsymbol\beta,\sigma|m)\, p(\mathbf{y}|\mathbf{X},\boldsymbol\beta,\sigma,m)}{p(\boldsymbol\beta,\sigma|\mathbf{y},\mathbf{X},m)}

Note that this equation is nothing but 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 inference method such as Monte Carlo sampling^[4] or variational Bayes.

The special case $\boldsymbol\mu_0=0, \mathbf{\Lambda}_0 = c\mathbf{I}$ 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.

Notes

↑ The intermediate steps of this computation can be found in O'Hagan (1994) at the beginning of the chapter on Linear models.
↑ The intermediate steps are in Fahrmeir et al. (2009) on page 188.
↑ The intermediate steps of this computation can be found in O'Hagan (1994) on page 257.
↑ Carlin and Louis(2008) and Gelman, et al. (2003) explain how to use sampling methods for Bayesian linear regression.

References

Box, G. E. P.; Tiao, G. C. (1973). Bayesian Inference in Statistical Analysis. Wiley. ISBN 0-471-57428-7.
Carlin, Bradley P. and Louis, Thomas A. (2008). Bayesian Methods for Data Analysis, Third Edition. Boca Raton, FL: Chapman and Hall/CRC. ISBN 1-58488-697-8.
O'Hagan, Anthony (1994). Bayesian Inference. Kendall's Advanced Theory of Statistics. 2B (First ed.). Halsted. ISBN 0-340-52922-9.
Gelman, Andrew, Carlin, John B., Stern, Hal S. and Rubin, Donald B. (2003). Bayesian Data Analysis, Second Edition. Boca Raton, FL: Chapman and Hall/CRC. ISBN 1-58488-388-X.
Walter, Gero; Augustin, Thomas (2009). "Bayesian Linear Regression—Different Conjugate Models and Their (In)Sensitivity to Prior-Data Conflict" (PDF). Technical Report Number 069, Department of Statistics, University of Munich.
Goldstein, Michael; Wooff, David (2007). Bayes Linear Statistics, Theory & Methods. Wiley. ISBN 978-0-470-01562-9.
Fahrmeir, L., Kneib, T., and Lang, S. (2009). Regression. Modelle, Methoden und Anwendungen (Second ed.). Heidelberg: Springer. doi:10.1007/978-3-642-01837-4. ISBN 978-3-642-01836-7.
Rossi, Peter E.; Allenby, Greg M.; McCulloch, Robert (2006). Bayesian Statistics and Marketing. John Wiley & Sons. ISBN 0470863676.
Thomas P. Minka (2001) Bayesian Linear Regression, Microsoft research web page

External links

Bayesian estimation of linear models (R programming wikibook). Bayesian linear regression as implemented in R.

Software

python
- Bayesian Type-II Linear Regression code, tutorial
- ARD Linear Regression code
- ARD Linear Regression with kernelized features code, tutorial

Least squares and regression analysis

Computational statistics

Correlation and dependence

Regression analysis

Regression as a
statistical model

Linear regression	Simple linear regression Ordinary least squares Generalized least squares Weighted least squares General linear model

Predictor structure	Polynomial regression Growth curve (statistics) Segmented regression Local regression

Non-standard	Nonlinear regression Nonparametric Semiparametric Robust Quantile Isotonic

Non-normal errors	Generalized linear model Binomial Poisson Logistic

Decomposition of variance

Model exploration

Background

Design of experiments

Numerical approximation

Applications

Regression analysis category
Statistics category
Statistics portal
Statistics outline
Statistics topics

Statistics

Descriptive statistics

Continuous data

Center	Mean arithmetic geometric harmonic Median Mode

Dispersion	Variance Standard deviation Coefficient of variation Percentile Range Interquartile range

Shape	Moments Skewness Kurtosis L-moments

Count data

Index of dispersion

Summary tables

Dependence

Graphics

Data collection

Study design	Population Statistic Effect size Statistical power Sample size determination Missing data

Survey methodology	Sampling Standard error stratified cluster Opinion poll Questionnaire

Controlled experiments	Design control optimal Controlled trial Randomized Random assignment Replication Blocking Interaction Factorial experiment

Uncontrolled studies	Observational study Natural experiment Quasi-experiment

Statistical inference

Statistical theory

Frequentist inference

Point estimation	Estimating equations Maximum likelihood Method of moments M-estimator Minimum distance Unbiased estimators Mean-unbiased minimum-variance Rao–Blackwellization Lehmann–Scheffé theorem Median unbiased Plug-in

Interval estimation	Confidence interval Pivot Likelihood interval Prediction interval Tolerance interval Resampling Bootstrap Jackknife

Testing hypotheses	1- & 2-tails Power Uniformly most powerful test Permutation test Randomization test Multiple comparisons

Parametric tests	Likelihood-ratio Wald Score

Specific tests

Z (normal) Student's t-test F

Goodness of fit	Chi-squared Kolmogorov–Smirnov Anderson–Darling Normality (Shapiro–Wilk) Likelihood-ratio test Model selection Cross validation AIC BIC

Rank statistics	Sign Sample median Signed rank (Wilcoxon) Hodges–Lehmann estimator Rank sum (Mann–Whitney) Nonparametric anova 1-way (Kruskal–Wallis) 2-way (Friedman) Ordered alternative (Jonckheere–Terpstra)

Bayesian inference

Correlation	Pearson product–moment Partial correlation Confounding variable Coefficient of determination

Regression analysis	Errors and residuals Regression model validation Mixed effects models Simultaneous equations models Multivariate adaptive regression splines (MARS)

Linear regression	Simple linear regression Ordinary least squares General linear model Bayesian regression

Non-standard predictors	Nonlinear regression Nonparametric Semiparametric Isotonic Robust Heteroscedasticity Homoscedasticity

Generalized linear model	Exponential families Logistic (Bernoulli) / Binomial / Poisson regressions

Partition of variance	Analysis of variance (ANOVA, anova) Analysis of covariance Multivariate ANOVA Degrees of freedom

Categorical / Multivariate / Time-series / Survival analysis

Categorical

Multivariate

Time-series

General	Decomposition Trend Stationarity Seasonal adjustment Exponential smoothing Cointegration Structural break Granger causality

Specific tests	Dickey–Fuller Johansen Q-statistic (Ljung–Box) Durbin–Watson Breusch–Godfrey

Time domain	Autocorrelation (ACF) partial (PACF) Cross-correlation (XCF) ARMA model ARIMA model (Box–Jenkins) Autoregressive conditional heteroskedasticity (ARCH) Vector autoregression (VAR)

Frequency domain	Spectral density estimation Fourier analysis Wavelet

Survival

Survival function	Kaplan–Meier estimator (product limit) Proportional hazards models Accelerated failure time (AFT) model First hitting time

Hazard function	Nelson–Aalen estimator

Test	Log-rank test