Bayes regression: how is it done in comparison to standard regression?

Question

I got some questions about the Bayesian regression:

1. Given a standard regression as $y = \beta_0 + \beta_1 x + \varepsilon$. If I want to change this into a Bayesian regression, do I need prior distributions both for $\beta_0$ and $\beta_1$ (or doesn't it work this way)?

2. In standard regression one would try to minimize the residuals to get single values for $\beta_0$ and $\beta_1$. How is this done in Bayes regression?

---

I really struggle a lot here:

$$ \text{posterior} = \text{prior} \times \text{likelihood} $$

Likelihood comes from the current dataset (so it's my regression parameter but not as a single value but as a likelihood distribution, right?). Prior comes from a previous research (let's say). So I got this equation:

$$ y = \beta_1 x + \varepsilon $$

with $\beta_1$ being my likelihood or posterior (or is this just totally wrong)?

I simply can't understand how the standard regression transforms into a Bayes one.

Answer 1

The simple linear regression model

$$ y_i = \alpha + \beta x_i + \varepsilon $$

can be written in terms of the probabilistic model behind it

$$ \mu_i = \alpha + \beta x_i \\ y_i \sim \mathcal{N}(\mu_i, \sigma) $$

i.e. dependent variable $Y$ follows normal distribution parametrized by mean $\mu_i$, that is a linear function of $X$ parametrized by $\alpha,\beta$, and by standard deviation $\sigma$. If you estimate such a model using ordinary least squares, you do not have to bother about the probabilistic formulation, because you are searching for optimal values of $\alpha,\beta$ parameters by minimizing the squared errors of fitted values to predicted values. On another hand, you could estimate such model using maximum likelihood estimation, where you would be looking for optimal values of parameters by maximizing the likelihood function

$$ \DeclareMathOperator*{\argmax}{arg\,max} \argmax_{\alpha,\,\beta,\,\sigma} \prod_{i=1}^n \mathcal{N}(y_i; \alpha + \beta x_i, \sigma) $$

where $\mathcal{N}$ is a density function of normal distribution evaluated at $y_i$ points, parametrized by means $\alpha + \beta x_i$ and standard deviation $\sigma$.

In the Bayesian approach instead of maximizing the likelihood function alone, we would assume prior distributions for the parameters and use the Bayes theorem

$$ \text{posterior} \propto \text{likelihood} \times \text{prior} $$

The likelihood function is the same as above, but what changes is that you assume some prior distributions for the estimated parameters $\alpha,\beta,\sigma$ and include them into the equation

$$ \underbrace{f(\alpha,\beta,\sigma\mid Y,X)}_{\text{posterior}} \propto \underbrace{\prod_{i=1}^n \mathcal{N}(y_i\mid \alpha + \beta x_i, \sigma)}_{\text{likelihood}} \; \underbrace{f_{\alpha}(\alpha) \, f_{\beta}(\beta) \, f_{\sigma}(\sigma)}_{\text{priors}} $$

"What distributions?" is a different question, since there is an unlimited number of choices. For $\alpha,\beta$ parameters you could, for example, assume normal distributions parametrized by some hyperparameters, or $t$-distribution if you want to assume heavier tails, or uniform distribution if you do not want to make many assumptions, but you want to assume that the parameters can be a priori "anything in the given range", etc. For $\sigma$ you need to assume some prior distribution that is bounded to be greater than zero since standard deviation needs to be positive. This may lead to the model formulation as illustrated below by John K. Kruschke.

(source: http://www.indiana.edu/~kruschke/BMLR/)

While in the maximum likelihood you were looking for a single optimal value for each of the parameters, in the Bayesian approach by applying the Bayes theorem you obtain the posterior distribution of the parameters. The final estimate will depend on the information that comes from your data and from your priors, but the more information is contained in your data, the less influential are priors.

Notice that when using uniform priors, they take form $f(\theta) \propto 1$ after dropping the normalizing constants. This makes Bayes theorem proportional to the likelihood function alone, so the posterior distribution will reach its maximum at exactly the same point as the maximum likelihood estimate. What follows, the estimate under uniform priors will be the same as by using ordinary least squares since minimizing the squared errors corresponds to maximizing the normal likelihood.

To estimate a model in the Bayesian approach in some cases you can use conjugate priors, so the posterior distribution is directly available (see example here). However, in the vast majority of cases, posterior distribution will not be directly available and you will have to use Markov Chain Monte Carlo methods for estimating the model (check this example of using Metropolis-Hastings algorithm to estimate parameters of linear regression). Finally, if you are only interested in point estimates of parameters, you could use maximum a posteriori estimation, i.e.

$$ \argmax_{\alpha,\,\beta,\,\sigma} f(\alpha,\beta,\sigma\mid Y,X) $$

For a more detailed description of logistic regression, you can check the Bayesian logit model - intuitive explanation? thread.

For learning more you could check the following books:

Kruschke, J. (2014). Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan. Academic Press.

Gelman, A., Carlin, J. B., Stern, H. S., and Rubin, D. B. (2004). Bayesian data analysis. Chapman & Hall/CRC.

Answer 2

Given a data set $D = (x_1,y_1), \ldots, (x_N,y_N)$ where $x \in \mathbb{R}^d, y \in \mathbb{R}$, a Bayesian Linear Regression models the problem in the following way:

Prior: $$w \sim \mathcal{N}(0, \sigma_w^2 I_d)$$

$w$ is vector $(w_1, \ldots, w_d)^T$, so the previous distribution is a multivariate Gaussian; and $I_d$ is the $d\times d$ identity matrix.

Likelihood: $$Y_i \sim \mathcal{N}(w^T x_i, \sigma^2)$$

We assume that $Y_i \perp Y_j | w, i \neq j$

For now we'll use the precision instead of the variance, $a = 1/\sigma^2$, and $b = 1/\sigma_w^2$. We'll also assume that $a,b$ are known.

The prior can be stated as $$p(w) \propto \exp \Big\{ -\frac{b}{2} w^t w \Big\}$$

And the likelihood $$p(D|w) \propto \exp \Big\{ -\frac{a}{2} (y-Aw)^T (y-Aw) \Big\}$$

where $y = (y_1,\ldots,y_N)^T$ and $A$ is a $n\times d$ matrix where the i-th row is $x_i^T$.

Then the posterior is $$p(w|D) \propto p(D|w) p(w)$$

After many calculations we discover that

$$p(w|D) \sim \mathcal{N}(w | \mu, \Lambda^{-1})$$

where ($\Lambda$ is the precision matrix)

$$\Lambda = a A^T A + b I_d $$ $$\mu = a \Lambda^{-1} A^T y$$

Notice that $\mu$ is equal to the $w_{MAP}$ of the regular linear regression, this is because for the Gaussian, the mean is equal to the mode.

Also, we can make some algebra over $\mu$ and get the following equality ($\Lambda = aA^TA+bI_d$):

$$\mu = (A^T A + \frac{b}{a} I_d)^{-1} A^T y$$

and compare with $w_{MLE}$:

$$w_{MLE} = (A^T A)^{-1} A^T y$$

The extra expression in $\mu$ corresponds to the prior. This is similar to the expression for the Ridge regression, for the special case when $\lambda = \frac{b}{a}$. Ridge regression is more general because the technique can choose improper priors (in the Bayesian perspective).

For the predictive posterior distribution:

$$p(y|x,D) = \int p(y|x,D,w) p(w|x,D) dw = \int p(y|x,w) p(w|D) dw$$

it is possible to calculate that

$$y|x,D \sim \mathcal{N}(\mu^Tx, \frac{1}{a} + x^T \Lambda^{-1}x)$$

Reference: Lunn et al. The BUGS Book

For using a MCMC tool like JAGS/Stan check Kruschke's Doing Bayesian Data Analysis