A Bayesian Approach to Regression Models

1 References

Once again it seems to make sense to refer you to the sources for this material before presenting the material. Sections 3, 4, and 5 below follow in spirit the presentation in portions of Chapter 8 of

• Gelman, A., Carlin, J. B., Stern, H. S. and Rubin, D. B. (1995). Bayesian data analysis. New York: Chapman-Hall.

A much briefer account, that emphasizes the “shrinkage” aspect of Bayesian estimation, may be found on pp. 38–50 of

• Carlin, B. P. and Louis, T. A. (1996). Bayes and empirical Bayes methods for data analysis. New York: Chapman-Hall.

The Carlin-Louis material follows closely that of the classic (but not classical!) reference

• Lindley, D. V. and Smith, A. F. M. (1972). Bayes estimates for the linear model (with discussion). Journal of the Royal Statistical Society, Series B, 34, 1–41.

Finally, much more algebraic detail may be found in Chapters 8 and 9 of

• Box, G. E. P. and Tiao, G. C. (1973). Bayesian inference in statistical analysis. Reading MA: Addison-Wesley.

In this handout I will give a thumbnail sketch of the “essential” difference between classical and Bayesian inference, and show you the ingredients for the ordinary linear regression model. Then I will describe Bayesian inference with non-informative priors for the ordinary linear regression model, and show how this can be immediately extended to analyses with informative conjugate priors. I will close with an example.

2 Introduction

It is sometimes useful to lay out the computational ingredients of Bayesian and classical inference side-by-side to see what differences to expect in practical data analysis problems. Here is my thumbnail sketch of the difference.

Let us suppose that we are interested in the usual linear regression model

\underline{y} = X\underline{β} + \underline{ϵ}, \underline{ϵ} ∼ N(\underline{0}, σ^{2}I)

where X is fixed an known in advance. To make the notation more compact I will collect together all of the parameters of this model into one parameter vector \underline{θ}:

\underline{θ}^{T} = (\underline{β}^{T} , σ^{2})^{T} = (β_{0}, β_{1}, …, β_{k}, σ^{2})^{T}

So far I have specified enough for you to write down the (or “a”, if you are fussy about it) likelihood of the data

p(\underline{y}|X, \underline{θ}) = ∏_{1}^{n} \frac{1}{\sqrt{2πσ^{2}}} exp [−\frac{1}{2}(y_{i} − \underline{x}_{i}^{T} \underline{β})^{2}/σ^{2}]

= \frac{1}{(2πσ^{2})^{n/2}} exp [−\frac{1}{2} ∑_{1}^{n}(y_{i} − \underline{x}_{i}^{T} \underline{β})^{2}/σ^{2}]
= \frac{1}{(2π)^{n/2}|σ^{2}I|^{1/2}} exp [−\frac{1}{2}(\underline{y} − X\underline{β})^{T} (σ^{2}I)^{−1}(\underline{y} − X\underline{β})] ,

where as usual \underline{x}_{i}^{T} is the i^{th} row of the matrix X.

\underline{Classical} statistical analysis is driven by the interpretation of p(\underline{y}|X, \underline{θ}) as a probability density for the random quantities \underline{y}, given fixed and known X \underline{and} fixed and unknown \underline{θ}. Inferences for \underline{θ} are driven by “what if” style questions and answers:

• Point estimation: What value of \underline{θ} maximizes the probability of observing the data \underline{y} that we saw, rather than some other data? (MLE’s)

• Interval estimation: Can we develop recipies A(\underline{y}, X), B(\underline{y}, X) such that P [ A(\underline{y}, X) ≤ θ_{j} ≤ B(\underline{y}, X) | \underline{θ} ] > 0.95, in repeated samples of \underline{y}? Then apply the recipies to the \underline{y} actually observed. (Confidence intervals; the probability statement is about the recipies, not the particular numbers calculated)

• Hypothesis testing:

– Can we develop test functions T (\underline{y}, X) such that the distribution q_{0}[T (\underline{y}, X)|\underline{θ} = \underline{θ}_{0}] is wellunderstood? If we calculate T (\underline{y}, X) for the data \underline{y} actually observed and it is an outlier for the q_{0} distribution this is evidence against \underline{θ} = \underline{θ}_{0} (general hypothesis testing)

– Is the likelihood of the data under the “small model” much smaller than the likelihood of the data under a “big model” in which it is nested? (Likelihood ratio testing)

Note that all of the reasoning about distributions is done prospectively, as if predicting \underline{y} on the basis of hypothetical values of \underline{θ} before \underline{y} is actually observed. Then \underline{y} is observed and compared with the predictions made to “learn” which hypothetical values of \underline{θ} make sense for the observed \underline{y}.

\underline{Bayesian} statistical analysis supplies an additional set of modeling assumptions, that \underline{θ} ∼ p(\underline{θ}|X) is also random, with marginal, or “prior”, distribution^{1} p(\underline{θ}|X); and it is driven by the interpretation of

p(\underline{θ}|\underline{y}, X) = \frac{p(\underline{y}|X, \underline{θ})p(\underline{θ}|X)}{p(\underline{y}|X)} ∝ p(\underline{y}|X, \underline{θ})p(\underline{θ}|X)  (1)

as a probability density for the random quantities \underline{θ}, \underline{after} conditioning on observing \underline{y}, the data actually observed. Inferences for \underline{θ} are driven by direct calculation on the “posterior” distribution p(\underline{θ}|\underline{y}, X):

1Since the particular values of the entries in X are not usually thought of as constraining the “true value” of \underline{θ} in a regression problem, it usually makes sense to make \underline{θ} independent of X: p(\underline{θ}|X) = p(\underline{θ}).

• Point estimation: What value of \underline{θ} does p(\underline{θ}|\underline{y}, X) assign the highest probability to? (posterior mode) What is the average \underline{θ}, given \underline{y} and X? (posterior mean)

• Interval estimation: Can we develop recipies A(\underline{y}, X), B(\underline{y}, X) such that in the posterior distribution P [A(\underline{y}, X) ≤ θ_{j} ≤ B(\underline{y}, X) | \underline{y}, X] > 0.95? (credible intervals’s; the probability statement is about the actual numbers calculated)

• Hypothesis testing:

– Can we calculate the distribution of additional data \underline{y}^{∗} that follows the same model we are fitting? If so, does the test function T (\underline{y}, X) evaluated on the data actually observed look like an outlier in this “predictive” distribution? (Posterior predictive checks)

– Is the probability of the data P (\underline{y}|X, M_{1}) under model 1 greater than that under model 2, P (\underline{y}|X, M_{2})? (Bayes factors)

As indicated above, a common method of point estimation for Bayesian inference is to guess that \underline{θ} is in the region of the parameter space where the posterior p(\underline{θ}|\underline{y}, X) ∝ p(\underline{y}|X, \underline{θ})p(\underline{θ}|X) is highest. The likelihood imposes one set of constraints on where \underline{θ} should be (near the MLE, where the likelihood p(\underline{y}|X, \underline{θ}) is highest) and the prior imposes another set of constraints (near the prior mode, where the prior distribution p(\underline{θ}|X) is highest). The Bayesian point estimate is a compromise^{2} between these two sets of constraints.

For this reason, especially strong features of the likelihood p(\underline{y}|X, \underline{θ}) or of the prior p(\underline{θ}|X) will dominate the Bayesian analysis. For example a common situation for iid or exchangeable data is “likelihood dominated” or “data dominated” inference, in which there are so many observations that the likelihood p(\underline{y}|X, \underline{θ}) has a single very sharp, very narrow peak, and the prior looks flat or featureless by contrast. In this case, the prior will exert very little influence in a Bayesian analysis (since the shape of the posterior will be essentially the same as the shape of the likelihood), and the Bayesian results will be very similar to the classical results. One way to ensure likelihood dominance is to choose “flat” or “noninformative” priors. Then we expect the Bayesian answers to be numerically and computationally parallel the classical ones (though the probability interpretations of the answers will be different, in terms of p(\underline{θ}|\underline{y}, X) rather than p(\underline{y}|X, \underline{θ})). As we move away from flat priors, the special features of the prior will pull the inferences away from the classical results.

3 Analysis of the Ordinary Linear Regression Model with Non-informative Priors

A central difficulty of a purist approach to Bayesian statistics is the careful specification of prior distributions for the parameters. As indicated above, in likelihood-dominated problems, careful specification of the prior is not very important unless there is tension between the prior information and the data. An alternative is to use a noninformative prior distribution^{3} which is usually so featureless that every problem turns into a likelihood-dominated problem.

2This compromise is also called “shrinkage”: the Bayesian answer “shrinks” the classical answer toward the prior mean or mode. In typical location models, when the prior mean is 0, the effect is to reduce, or shrink, the unsigned magnitude of the classical answer.

3Noninformative priors are usually flat (i.e. constant) or nearly flat, and if the parameter space is unbounded this leads to “improper” prior distributions, that is, a prior distribution that assigns infinite probability (instead of probability 1) to the parameter space. Often the posterior distribution is still “proper” so that formally we can work with the posterior as indicating probable locations for the paramter, as in any other Bayesian problem. However, a careful interpretation of the posterior as a conditional density derived from the joint density p(\underline{y}, \underline{θ}|X), when the prior is improper, is not possible.

In regression analysis, starting with a non-informative prior has two additional advantages: (1) it shows clearly the parallels between classical and Bayesian analysis; (2) it is a good foundation on which to build Bayesian analyses when we have prior information that can be specified through a conjugate prior distribution.

The standard, or default, non-informative prior for regression analysis is

p(\underline{β}, log σ) ∝ constant

or, after a calculus change-of-variables

p(\underline{β}, σ^{2}) ∝ σ^{−2}

Note that this is flat (i.e. constant, or uniform) in the regression coefficients \underline{β}, so no constraints beyond those imposed by the likelihood are placed on them. It is not flat in σ^{2}, however: values of σ^{2} near zero are preferred over values near infinity. It is also easy to verify that this prior is improper in each coordinate of the parameter vector \underline{θ} = (\underline{β}^{T} , σ^{2})^{T}.

3.1 Setting up the joint posterior distribution

Following Equation (1), and plugging in the likelihood and default prior for linear regression, we then have

log p(\underline{θ}|\underline{y}, X) = (Constant) + log p(\underline{y}|X, \underline{θ}) + log p(\underline{θ}|X)  (2)
= (Constant) + [−n log σ − \frac{1}{2}(\underline{y} − X\underline{β})^{T} (\underline{y} − X\underline{β})/σ^{2}] + [ − 2 log σ]
= (Constant) − (n + 2) log σ − \frac{1}{2}(\underline{y} − X\underline{β})^{T} (\underline{y} − X\underline{β})/σ^{2},

where the “constant” terms in each line are not necessarily the same, but none depend on \underline{θ}. One can readily verify^{4} (multiply it out, and use the definition of ˆ \underline{β}!) that

(\underline{y} − X\underline{β})^{T} (\underline{y} − X\underline{β})/σ^{2} = (\underline{y} − X ˆ\underline{β})^{T} (\underline{y} − X ˆ\underline{β})/σ^{2} + (\underline{β} − ˆ \underline{β})^{T} X^{T} X(\underline{β} − ˆ \underline{β})/σ^{2}

where

ˆ\underline{β} = (X^{T} X)^{−1}X^{T} \underline{y}

is the usual OLS estimator of \underline{β}, and hence

log p(\underline{θ}|\underline{y}, X) = (Constant) − (n + 2) log σ − \frac{1}{2}(\underline{y} − X ˆ\underline{β})^{T} (\underline{y} − X ˆ\underline{β})/σ^{2}
− \frac{1}{2}(\underline{β} − ˆ \underline{β})^{T} X^{T} X(\underline{β} − ˆ \underline{β})/σ^{2}  (3)

3.2 Identifying the pieces

The part of Equation (3) that depends on \underline{β} gives the shape^{5} of log p(\underline{β}|σ^{2}, \underline{y}, X); since this is a quadratic form in \underline{β} it follows that p(\underline{θ}|σ^{2}, \underline{y}, X) is multivariate normal:

\underline{β} | σ^{2}, \underline{y}, X ∼ N(ˆ\underline{β}, σ^{2}(X^{T} X)^{−1})  (Bayes)

4cf. Box and Tiao, p. 439.

5This is an extremely useful trick in general for determining the so-called “complete conditinals” p(θ_{1}|θ_{2}, …, θ_{p}, data) in setting up Markov Chain Monte Carlo (MCMC) calculations.

Note how this corresponds to the classical result

ˆ\underline{β} | \underline{β}, σ^{2}, X ∼ N(\underline{β}, σ^{2}(X^{T} X)^{−1})  (Classical)

Thinking now of

p(\underline{β}, σ^{2}|\underline{y}, X) = p(\underline{β}|σ^{2}, \underline{y}, X)p(σ^{2}|\underline{y}, X),

the marginal density for p(σ^{2}|\underline{y}, X) could be obtained by integrating \underline{β} out of the joint density p(\underline{β}|σ^{2}, \underline{y}, X)p(σ^{2}|\underline{y}, X), but a cleverer thing to do is to realize that this also says that

p(σ^{2}|\underline{y}, X) = \frac{p(\underline{β}, σ^{2}|\underline{y}, X)}{p(\underline{β}|σ^{2}, \underline{y}, X)}

for any value of \underline{β}. After some algebra (plugging in \underline{β} ≡ ˆ \underline{β} simplifies things I think) one obtains that

σ^{2} | \underline{y}, X ∼ Inv-χ^{2} (n − k − 1, s^{2})

a scaled, inverse chi-squared^{6} random variable, where s^{2} is the usual squared residual standard error,

s^{2} = \frac{1}{n − k − 1}(\underline{y} − X ˆ\underline{β})^{T} (\underline{y} − X ˆ\underline{β}) ≡ \frac{RSS}{n − k − 1}.

This is equivalent to saying that

\frac{(n − k − 1)s^{2}}{σ^{2}} | \underline{y}, X ∼ χ_{n−k−1}^{2},  (Bayes)

which corresponds to the classical result

\frac{(n − k − 1)s^{2}}{σ^{2}} | \underline{β}, σ^{2}, X ∼ χ_{n−k−1}^{2}  (Classical)

To summarize, if we use the default noninformative prior p(\underline{β}, σ^{2}) ∝ σ^{−2}, we obtain a posterior for \underline{β} and σ^{2} that can be broken down as follows:

• p(\underline{β}|σ^{2}, \underline{y}, X) is N(ˆ\underline{β}, σ^{2}(X^{T} X)^{−1});

• p(σ^{2}|\underline{y}, X) is Inv-χ^{2} (n − k − 1, s^{2}).

As with the classical results, the Bayesian results only make sense (and the posterior will only be proper) if X is of full rank k + 1 and n > k + 1.

3.3 Further inferences

It is also possible to show that the usual t_{n−k−1} distributions for each β̂_{j} in classical OLS regression carry over to the same t_{n−k−1} posterior distributions for each β_{j}|\underline{y}, X in this Bayesian framework, etc., but it is often simpler to obtain approximations to this distribution by Monte Carlo simulation from the joint posterior for (\underline{β}, σ^{2}) and then “throw away” the σ^{2}’s simulated:

6Definition: If V ∼ χ_{d}^{2} then ds^{2}/V ∼ Inv-χ^{2} (d, s^{2}). See Appendix A of Gelman et al. for formulas.

1. Simulate σ^{2} from p(σ^{2}|\underline{y}, X), i.e. from Inv-χ^{2} (n − k − 1, s^{2});

2. Then simulate \underline{β} from p(\underline{β}|σ^{2}, \underline{y}, X), i.e. from N(ˆ\underline{β}, σ^{2}(X^{T} X)^{−1}).

This simulation is based on Monte Carlo estimation of the integral

∫ ∫ T (\underline{β}, σ^{2}) p(\underline{β}, σ^{2}|\underline{y}, X) d(σ^{2}) d(\underline{β}) = ∫ ∫ T (\underline{β}, σ^{2}) p(\underline{β}|σ^{2}, \underline{y}, X) p(σ^{2}|\underline{y}, X) d(σ^{2}) d(\underline{β}), (4)

where T (\underline{β}, σ^{2}) is any test function of interest.

Finally we turn to the problem of predicting new y values. Suppose

X̃ = [ \begin{matrix}1 & x̃_{11} & ⋯ & x̃_{1k} \\ 1 & x̃_{11} & ⋯ & x̃_{1k} \\ ⋮ & ⋮ & ... & ⋮ \\ 1 & x̃_{m1} & ⋯ & x̃_{mk} \end{matrix} ]

are new values of the regressors, and

˜\underline{y} = [ ỹ_{1} ỹ_{2} ⋯ ỹ_{m} ]^{T}

are the corresponding new y values. The linear regression model says as usual that p(˜\underline{y}|\underline{β}, σ^{2}, X̃) will be N(X̃\underline{β}, σ^{2}I_{m×m}). We can then use Equation (4) with T (\underline{β}, σ^{2}) = p(˜\underline{y}|\underline{β}, σ^{2}, X̃),

p(˜\underline{y}|\underline{y}, X, X̃) = ∫ ∫ p(˜\underline{y}|\underline{β}, σ^{2}, X̃, \underline{y}, X)p(\underline{β}|σ^{2}, \underline{y}, X)p(σ^{2}|\underline{y}, X) d(σ^{2})d(\underline{β})

in two ways:

• To derive the exact “predictive” distribution of ˜ \underline{y}. It turns out for example that (omitting X and X̃ from the conditioning to save space)

E[˜\underline{y} | σ^{2}, \underline{y}] = E[E(˜\underline{y} | \underline{β}, σ^{2}, \underline{y}) | σ^{2}, \underline{y}]
= E[X\underline{β} | σ^{2}, \underline{y}]
= X̃ ˆ\underline{β}

and similarly

Var(˜\underline{y}|σ^{2}, \underline{y}) = E[Var(˜\underline{y} | \underline{β}, σ^{2}, \underline{y}) | σ^{2}, \underline{y}] + Var(E[˜\underline{y} | \underline{β}, σ^{2}, \underline{y}] | σ^{2}, \underline{y})
= E[σ^{2}I_{m×m} | σ^{2}, \underline{y}] + Var(X̃\underline{β} | σ^{2}, \underline{y})
= σ^{2}(I_{m×m} + X̃(X^{T} X)^{−1}X̃^{T})

Note again how these results correspond with the classical ones: if σ^{2} is known, ˜ \underline{y} will follow a multivariate normal distribution with this mean and variance; if σ^{2} is not known, then the corresponding t distribution results hold.

• To obtain simulated values of ˜ \underline{y} by simulating from the joint posterior density of (˜\underline{y}, \underline{β}, σ^{2}) and throwing away the simulated \underline{β} and σ^{2} values:

1. Simulate σ^{2} from p(σ^{2}|\underline{y}, X), i.e. from Inv-χ^{2} (n − k − 1, s^{2});

2. Then simulate \underline{β} from p(\underline{β}|σ^{2}, \underline{y}, X), i.e. from N(ˆ\underline{β}, σ^{2}(X^{T} X)^{−1});

3. Finally simulate ˜ \underline{y} from p(˜\underline{y}|\underline{β}, σ^{2}, X̃), i.e. from N(X̃\underline{β}, σ^{2}I_{m×m}).

3.4 Example: Radon measurements

Gelman et al. (1995, pp. 189 and 263) consider the data in Table 1 to learn about the effects of geographic location in Minnesota and floor of house on short term radon measurements in residential dwellings. An OLS regression analysis produces the following results

> radon _read.table("radon.txt",header=T)
> radon _cbind(LogRadon = log(radon$Radon),radon)
> names(radon)
[1] "LogRadon" "Radon"     "County"   "Floor"
> radon$Floor _ifelse(radon$Floor=="B",0,1)
> radon _cbind(radon,BE=ifelse(radon$County=="BlueEarth",1,0))
> radon _cbind(radon,Clay=ifelse(radon$County=="Clay",1,0))
> radon _cbind(radon,Goodhue=ifelse(radon$County=="Goodhue",1,0))
> names(radon)
[1] "LogRadon"  "Radon"  "County"  "Floor"  "BE"  "Clay"  "Goodhue"
> summary(olsfit _lm(LogRadon ~ Floor + Clay + Goodhue,data=radon))
Call: lm(formula = LogRadon ~ Floor + Clay + Goodhue, data = radon)
Residuals:
     Min       1Q   Median      3Q   Max
 -1.653 -0.5572 -0.01608 0.5933 1.855
Coefficients:
                Value Std. Error t value Pr(>|t|)
(Intercept)  1.9561  0.2158      9.0663  0.0000
       Floor -0.3283  0.3155     -1.0409  0.3047
        Clay -0.0798  0.3018     -0.2644  0.7929
     Goodhue -0.0379  0.3041     -0.1248  0.9014
Residual standard error: 0.7895 on 37 degrees of freedom
Multiple R-Squared: 0.0332
F-statistic: 0.4235 on 3 and 37 degrees of freedom, the p-value is 0.7373
Correlation of Coefficients:
         (Intercept)   Floor     Clay
  Floor -0.2089
   Clay -0.6526      -0.1493
Goodhue -0.6762      -0.0114  0.4868

The predictors do not look like they matter much, but we will forge ahead. We first obtain the basic ingredients needed to specify the posterior distributions of σ^{2} and \underline{β}

> X _model.matrix(olsfit)
> print(r.df _41 - 3 - 1)
[1] 37
> R _qr.R(qr(X))       # look at help(qr), help(qr.R), etc.

Radon County  Floor
5.0 BlueEarth B
13.0 BlueEarth B
7.2 BlueEarth B
6.8 BlueEarth B
12.8 BlueEarth B
5.8 BlueEarth F
9.5 BlueEarth B
6.0 BlueEarth B
3.8 BlueEarth B
14.3 BlueEarth F
1.8 BlueEarth B
6.9 BlueEarth B
4.7 BlueEarth B
9.5 BlueEarth B
0.9 Clay  F
12.9 Clay  B
2.6 Clay  B
3.5 Clay  F
26.6 Clay  B
1.5 Clay  B
13.0 Clay  B

Radon County Floor
8.8 Clay  B
19.5 Clay  B
2.5 Clay  F
9.0 Clay  B
13.1 Clay  B
3.6 Clay  B
6.9 Clay  F
14.3 Goodhue B
6.9 Goodhue F
7.6 Goodhue B
9.8 Goodhue F
2.6 Goodhue B
43.5 Goodhue B
4.9 Goodhue B
3.5 Goodhue B
4.8 Goodhue B
5.6 Goodhue B
3.5 Goodhue B
3.9 Goodhue B
6.7 Goodhue B

Table 1: Short-term measurements of radon concentration (in picocuries/liter) in a sample of houses in three Minnesota counties. Floor=B for basement measurements, and Floor=F for first-floor measurements.

> Q _qr.Q(qr(X))
> RI _solve(R)         # easy since R is upper-triangular
> XTXI _RI %*% t(RI)  # cheaper than XTXI _solve(t(X)%*%X) !!!
> print(bhat _solve(R,t(Q) %*% radon$LogRadon))
                          # solves R %*% bhat = t(Q) %*% y for bhat ...
                          # cheaper than XTXI %*% t(X) %*% radon$LogRadon
                      [,1]
(Intercept)  1.95613199
       Floor -0.32833983
        Clay -0.07980519
     Goodhue -0.03794035
> print(sqrt(s2 _sum(residuals(olsfit)^2)/r.df))
[1] 0.7894929

Now we could give analytic results here but it is more informative to simulate:

> V _rchisq(1000,r.df)
> sig2 _r.df*s2/V         #  1000 draws from scaled inv-chi-squared
> beta _matrix(rnorm(4000,0,sqrt(sig2)),nrow=4,ncol=1000,byrow=T)
> beta _RI %*% beta + matrix(bhat,nrow=4,ncol=1000)
                              #  1000 corresponding draws from post for beta

Plotting the results takes more work than the simulations:

> smoothplot _function(x,xlab="",ylab="") {
+ plot(density(x,width=sqrt(var(x))),type="l",xlab=xlab,ylab=ylab)
+ }
> par(mfrow=c(2,3))
> smoothplot(sig2,"sig2")
> abline(v=quantile(sig2,c(0.025,0.5,0.975)))
> abline(v=s2,lty=2)
> olsdata _summary(olsfit)$coefficients
> olsintervals _olsdata[,1] + cbind(olsdata[,2]) %*% c(-2,0,2)
> smoothplot(beta[1,],"Intercept")
> abline(v=quantile(beta[1,],c(0.025,0.5,0.975)))
> abline(v=olsintervals[1,c(1,3)],lty=2)
> smoothplot(beta[2,],"Floor")
> abline(v=quantile(beta[2,],c(0.025,0.5,0.975)))
> abline(v=olsintervals[2,c(1,3)],lty=2)
> smoothplot(beta[3,],"Clay")
> abline(v=quantile(beta[3,],c(0.025,0.5,0.975)))
> abline(v=olsintervals[3,c(1,3)],lty=2)
> smoothplot(beta[4,],"Goodhue")
> abline(v=quantile(beta[4,],c(0.025,0.5,0.975)))
> abline(v=olsintervals[4,c(1,3)],lty=2)
> c(sig2=s2,quantile(sig2,c(0.025,0.5,0.975)))

Figure 1: Simulated posteriors for the parameters of the radon regression problem, with noninformative priors. Solid lines indicate equal tail area 95% posterior credible intervals and posterior medians; dotted lines indicate approximate 95% confidence intervals (for β’s) or point estimate (for σ^{2}) from OLS analysis.

      sig2       2.5%      50.0%     97.5%
 0.623299 0.4170907 0.6317195 1.053781
> dimnames(olsintervals)[[2]] _c("bhat-2SE","bhat","bhat+2SE")
> cbind(post025 = apply(beta,1,quantile,0.025),olsintervals,
+ post975=apply(beta,1,quantile,0.975))
                 post025   bhat-2SE         bhat  bhat+2SE   post975
(Intercept)  1.5070074  1.5246130  1.95613199 2.3876510 2.4275849
       Floor -0.9767665 -0.9592412 -0.32833983 0.3025615 0.3038866
        Clay -0.7473966 -0.6833730 -0.07980519 0.5237626 0.5542808
     Goodhue -0.6708612 -0.6461490 -0.03794035 0.5702683 0.5657630

The results are shown in Figure 1.

Now suppose we wish to predict the radon measurement of another house in Blue Earth county. If the radon concentration is measured in the basement, we can simulate the distribution of the new measurement ỹ as follows

> X[1,]
 (Intercept) Floor Clay Goodhue
             1      0     0        0
> errors _rnorm(1000,0,sqrt(sig2))
> ynew _c(c(1,0,0,0) %*% beta + errors)
> smoothplot(ynew,"Predicted radon, Blue Earth Basement Meas")
> abline(v=quantile(ynew,c(0.025,0.5,0.975)))

> abline(v=c(1,0,0,0) %*% bhat +
+ c(-2,2)*sqrt(s2*(1+c(1,0,0,0)%*%XTXI%*%c(1,0,0,0)))
+  , lty=2)
> rbind(Bayes=quantile(ynew,c(0.025,0.975)),OLS= c(1,0,0,0) %*% bhat +
+ c(-2,2)*sqrt(s2*(1+c(1,0,0,0)%*%XTXI%*%c(1,0,0,0)))
+ )
             2.5%     97.5%
Bayes 0.2933040 3.635836
  OLS 0.3192433 3.593021

If the measurement is taken on the first floor, the simulated distribution can be obtained as follows

> ynew _c(c(1,1,0,0) %*% beta + errors)
> smoothplot(ynew,"Predicted radon, Blue Earth First Fl. Meas")
> abline(v=quantile(ynew,c(0.025,0.5,0.975)))
> abline(v=c(1,1,0,0) %*% bhat +
+ c(-2,2)*sqrt(s2*(1+c(1,1,0,0)%*%XTXI%*%c(1,1,0,0)))
+  , lty=2)
> rbind(Bayes=quantile(ynew,c(0.025,0.975)),OLS= c(1,1,0,0) %*% bhat +
+ c(-2,2)*sqrt(s2*(1+c(1,1,0,0)%*%XTXI%*%c(1,1,0,0)))
+ )
               2.5%     97.5%
Bayes -0.13944779 3.313950
  OLS -0.09375259 3.349337

The results are shown in Figure 2.

4 Bayesian GLS

A Bayesian analysis corresponding to GLS (generalized least squares) proceeds as follows^{7}. We will assume that

\underline{y} | \underline{β}, Σ_{y}, X ∼ N(X\underline{β}, Σ_{y})
and p(\underline{β}|Σ_{y}, X) ∝ constant.

Let us assume Σ_{y} is known; this lets us calculate p(\underline{β}|Σ_{y}, \underline{y}, X). Let Σ_{y}^{1/2} be a Cholesky square root of Σ_{y} (i.e., a triangular matrix such that [Σ_{y}^{1/2}]^{T} Σ_{y}^{1/2} = Σ_{y}). Then as with our discussion of WLS, we can premultiply by Σ_{y}^{−1/2} and obtain

Σ_{y}^{−1/2}\underline{y} | β, Σ_{y} ∼ N(Σ_{y}^{−1/2}X\underline{β}, I)

Applying the ordinary regression results above to \underline{z} = Σ_{y}^{−1/2}\underline{y} we immediately obtain

\underline{β} | Σ_{y}, \underline{y}, X ∼ N(ˆ\underline{β}, V_{β})

7This follows very closely Section 8.7 of Gelman et al.

Figure 2: Simulated predictive distributions for a new radon measurement in Blue Earth County: basement measurement (left) and first floor measurement (right). Solid lines indicate equal tail area 95% posterior predictive credible intervals and medians; dotted lines indicate approximate 95% predictive confidence intervals from OLS analysis.

where

ˆ\underline{β} = (X^{T} Σ_{y}^{−1}X)^{−1}X^{T} Σ_{y}^{−1}\underline{y}
V_{β} = (X^{T} Σ_{y}^{−1}X)^{−1}

These calculations turn out to be very useful in Section 5 below, where we are considering Bayesian linear regression with informative, conjugate prior information.

For completness we also indicate the form of p(Σ_{y}|\underline{y}, X), but we don’t intend to use it in the rest of these notes. As before, we use

p(\underline{β}, Σ_{y} | \underline{y}, X) = p(\underline{β} | Σ_{y}, \underline{y}, X) p(Σ_{y} | \underline{y}, X)

to see that, for any value of \underline{β},

p(Σ_{y} | \underline{y}, X) = \frac{p(\underline{β}, Σ_{y} | \underline{y}, X)}{p(\underline{β} | Σ_{y}, \underline{y}, X)} ∝ \frac{p(Σ_{y})N(\underline{y} | X\underline{β}, Σ_{y})}{N(\underline{β} | X ˆ\underline{β}, V_{β})}

where p(Σ_{y}) is now the prior on Σ_{y} (the prior on \underline{β} is still flat). Substituting \underline{β} ≡ ˆ \underline{β} again and simplifying, we see that

p(Σ_{y}|\underline{y}, X) ∝ p(Σ_{y})|V_{β}|^{1/2} exp [−\frac{1}{2}(\underline{y} − X ˆ\underline{β})^{T} Σ_{y}^{−1}(\underline{y} − X ˆ\underline{β})] .  (5)

This is relatively easy to write down but difficult to work with—both because p(Σ_{y}) is difficult to specify (or “elicit”) and because Σ_{y} is wrapped up in the computation of both ˆ \underline{β} and V_{β}—except in very special circumstances, such as WLS (weighted least squares) where at least Σ_{y} = σ^{2}diag(1/w_{1}, …, 1/w_{n}) only depends on a single parameter.

5 Prior Information

5.1 Prior information for \underline{β}; Σ_{y} known

We start by assuming that the variance-covariance matrix Σ_{y} is known, and that we have multivariate normal prior information on \underline{β}:

\underline{y} | \underline{β}, Σ_{y} ∼ N(X\underline{β}, Σ_{y})
and \underline{β} ∼ N(\underline{β}_{0}, Σ_{β})

where Σ_{y}, \underline{β}_{0}, and Σ_{β} are known and fixed in advance.

Aside: Another way to say the same thing is that

\begin{matrix}\underline{y} & = & X\underline{β} + \underline{ϵ}, \underline{ϵ} ∼ N(\underline{0}, Σ_{y}) \\ \underline{β} & = & \underline{β}_{0} + \underline{δ}, \underline{δ} ∼ N(\underline{0}, Σ_{β}) \end{matrix} \}  (6)

from which it immediately follows that

\underline{y} = X\underline{β}_{0} + X\underline{δ} + \underline{ϵ}
i.e. \underline{y} ∼ N(X\underline{β}_{0}, XΣ_{β}X^{T} + Σ_{y})

• If \underline{β}_{0} were not known, it could now be estimated by (non-Bayesian) WLS. If we then proceed with a Bayesian analysis for \underline{β} pretending that \underline{β}_{0} was fixed in advance to be ˆ \underline{β}_{0} from the WLS analysis, then the resulting analysis would be an example of an empirical Bayes analysis.

• If \underline{β}_{0} had the further linear structure \underline{β}_{0} = Z\underline{γ}, the resulting model (6) would be called a (two-level) hierarchical linear model (HLM).

End of vocabulary lesson.

Continuing under the assumption that Σ_{y}, \underline{β}_{0}, and Σ_{β} are known and fixed in advance, we know that

p(\underline{β}|Σ_{y}, \underline{y}, X) ∝ N(y|X\underline{β}, Σ_{y}) · N(\underline{β}|\underline{β}_{0}, Σ_{β})

and so

log p(\underline{β}|Σ_{y}, \underline{y}, X) = (Constant) − \frac{1}{2}(y − X\underline{β})^{T} Σ_{y}^{−1}(y − X\underline{β}) − \frac{1}{2}(\underline{β} − \underline{β}_{0})^{T} Σ_{β}^{−1}(\underline{β} − \underline{β}_{0})
= (Constant) − \frac{1}{2} ([ \underline{y} ; \underline{β}_{0} ] − [ \begin{matrix}X\underline{β} \\ I\underline{β} \end{matrix} ])^{T} ( \begin{matrix}Σ_{y}^{−1} & 0 \\ 0 & Σ_{β}^{−1} \end{matrix} ) ([ \underline{y} ; \underline{β}_{0} ] − [ \begin{matrix}X\underline{β} \\ I\underline{β} \end{matrix} ])
= (Constant) − \frac{1}{2}(\underline{y}_{∗} − X_{∗}\underline{β})^{T} Σ_{∗}^{−1}(\underline{y}_{∗} − X_{∗}\underline{β})

where the “constant’s” do not depend on \underline{β}, where I is (k + 1) × (k + 1) identity matrix, and where

\underline{y}_{∗} = [ \underline{y} ; \underline{β}_{0} ]  X_{∗} = [ \begin{matrix}X \\ I \end{matrix} ]  Σ_{∗} = ( \begin{matrix}Σ_{y} & 0 \\ 0 & Σ_{β} \end{matrix} )

are “augmented” data and design matrices. In other words regression with prior information on \underline{β} is like regression with the “augmented” data \underline{y}_{∗} and design matrix X_{∗} and noninformative prior on \underline{β}, and known variance-covariance matrix Σ_{∗}. In other words, the “prior” information can be interpreted as additional “observations” in a noninformative-prior Bayesian analysis of the linear regression model

\underline{y}_{∗} = X_{∗}\underline{β} + \underline{ϵ}_{∗}, \underline{ϵ}_{∗} ∼ N(\underline{0}, Σ_{∗})

We can immediately apply the results on “Bayesian GLS” analysis (Section 4) to conclude that

\underline{β}|Σ_{y}, \underline{y}, X ∼ N(ˆ\underline{β}, V_{β})

with

ˆ\underline{β} = (X_{∗}^{T} Σ_{∗}^{−1}X_{∗})^{−1}X_{∗}^{T} Σ_{∗}^{−1}y_{∗}
= ([ X^{T} I^{T} ] [ \begin{matrix}Σ_{y}^{−1} & 0 \\ 0 & Σ_{β}^{−1} \end{matrix} ] [ \begin{matrix}X \\ I \end{matrix} ])^{−1} [ X^{T} I^{T} ] [ \begin{matrix}Σ_{y}^{−1} & 0 \\ 0 & Σ_{β}^{−1} \end{matrix} ] [ \underline{y} ; \underline{β}_{0} ]
= (X^{T} Σ_{y}^{−1}X + Σ_{β}^{−1})^{−1} [X^{T} Σ_{y}^{−1}\underline{y} + Σ_{β}^{−1}\underline{β}_{0}]
and V_{β} = (X_{∗}^{T} Σ_{∗}^{−1}X_{∗})^{−1}
= (X^{T} Σ_{y}^{−1}X + Σ_{β}^{−1})^{−1}

If some of the β_{j}’s should have noninformative prior distributions, we set the corresponding rows and columns of Σ_{∗}^{−1} (or equivalently, of Σ_{β}^{−1}) equal to all zeros.

The case of ordinary regression with iid priors on \underline{β} gives some insight into the “shrinkage” interpretation of Bayesian estimates. This case occurs when

Σ_{y} = σ^{2}I_{n×n} Σ_{β} = τ^{2}I_{(k+1)×(k+1)}

so that the posterior mean and variance reduce to

ˆ\underline{β} = (X^{T} X + \frac{σ^{2}}{τ^{2}} I)^{−1} [X^{T} \underline{y} + \frac{σ^{2}}{τ^{2}} \underline{β}_{0}]

and V_{β} = σ^{2} (X^{T} X + \frac{σ^{2}}{τ^{2}} I)^{−1}

There are two interesting things to note about this formulation:

• As the prior information gets weaker (τ^{2} → ∞), ˆ \underline{β} → (X^{T} X)^{−1}X^{T} \underline{y}, the classical estimate, as in Section 3; and as the prior information gets stronger (τ^{2} → 0), ˆ \underline{β} → \underline{β}_{0}, the prior mean. This is the usual “shrinkage” behavior.

• When \underline{β}_{0} = \underline{0}, the estimate is of the form

ˆ\underline{β} = (X^{T} X + λI)^{−1} X^{T} \underline{y}

(with λ = σ^{2}/τ^{2}), i.e. it is the ridge regression estimate.

5.2 Prior information for Σ_{y}

Now that we have beaten p(\underline{β}|Σ_{y}, \underline{y}, X) = N(\underline{β}|ˆ\underline{β}, V_{β}) into the ground, let’s complete the formulation of a joint posterior by specifying p(Σ_{y}|\underline{y}, X). We will only do this in the case of ordinary regression, Σ_{y} = σ^{2}I, with a general normal prior \underline{β} ∼ N(\underline{β}_{0}, Σ_{β}). For this case, we saw in Section 3 that the natural posterior with a non-informative prior was an Inv-χ^{2}. Therefore it will not be a surprise that the conjugate prior family for σ^{2} in Bayesian ordinary regression is

Inv-χ^{2} (σ^{2}|n_{0}, σ_{0}^{2}).

Following the same calculations as before [see especially Equation (5)], using the joint prior

p(\underline{β}, σ^{2}|X) = N(\underline{β}|\underline{β}_{0}, Σ_{β}) · Inv-χ^{2} (σ^{2}|n_{0}, σ_{0}^{2})

it follows that

• p(\underline{β} | σ^{2}, \underline{y}, X) = N(\underline{β}|ˆ\underline{β}, V_{β}) where

ˆ\underline{β} = (X_{∗}^{T} Σ_{∗}^{−1}X_{∗})^{−1}X_{∗}^{T} Σ_{∗}^{−1}y_{∗} = (X^{T} X/σ^{2} + Σ_{β}^{−1})^{−1} [X^{T} \underline{y}/σ^{2} + Σ_{β}\underline{β}_{0}]  (7)
V_{β} = (X_{∗}^{T} Σ_{∗}^{−1}X_{∗})^{−1} = (X^{T} X/σ^{2} + Σ_{β}^{−1})^{−1}  (8)

and

• p(σ^{2}|\underline{y}, X) is given by

p(σ^{2}|\underline{y}, X) = Inv-χ^{2} (σ^{2} |n_{0} + n′, \frac{n_{0}σ_{0}^{2} + n′s^{2}}{n_{0} + n′} )  (9)

where s^{2} = (\underline{y}_{∗}−X_{∗} ˆ\underline{β})^{T} Σ_{∗}^{−1}(\underline{y}_{∗}−X_{∗} ˆ\underline{β})/n′; n′ = n_{∗}+k−1; and n_{∗} = rank(Σ_{∗}) = n+rank(Σ_{β}^{−1}) is typically n + k + 1. Once again, the effect of the informative conjugate prior is similar to adding n_{0} additional observations’ worth of data to the non-informative analysis.

5.3 Continuing the radon example

Suppose our radon measurement data were a replication of an earlier study of 41 houses in the same three counties. Suppose moreover that an OLS regression had been performed on the earlier data set, with the results

> summary(oldfit)
Call: lm(formula = LogRadon ~ Floor + Clay + Goodhue, data = oldrad)
Residuals:
     Min       1Q   Median   3Q   Max
 -1.853 -0.6117 -0.08181 0.48 1.592
Coefficients:
                Value Std. Error t value Pr(>|t|)

(Intercept)  1.9987  0.2218      9.0128  0.0000
       Floor -0.3447  0.3161     -1.0905  0.2826
        Clay  0.0937  0.3216      0.2915  0.7723
     Goodhue  0.1815  0.3065      0.5924  0.5572
Residual standard error: 0.8232 on 37 degrees of freedom
Multiple R-Squared: 0.04275
F-statistic: 0.5508 on 3 and 37 degrees of freedom, the p-value is 0.6508
Correlation of Coefficients:
         (Intercept)   Floor     Clay
  Floor -0.2851
   Clay -0.5961      -0.1311
Goodhue -0.6816       0.0589  0.4507

For our current data we might like to specify prior distributions based on the earlier study:

\underline{β} ∼ N(ˆ\underline{β}_{old}, s_{old}^{2}(X_{old}^{T}X_{old})^{−1})
σ^{2} ∼ Inv-χ^{2} (41, s_{old}^{2})

These can be obtained as follows

> Xold _model.matrix(oldfit)
> XTXIO _solve(t(Xold)%*%Xold)
> s2old _sum(residuals(oldfit)^2)/37
> beta0 _summary(oldfit)$coefficients[,1]
> Sigbeta _XTXIO*s2old
> SigbetaI _solve(Sigbeta)
> n0 _41

Following Equations (7), (8) and (9), and noting that n′ = n^{∗} − (k + 1) = n + (k + 1) − (k + 1) = n, we may set up a simulation for the posterior distribution of the parameters again as follows (though this time, not thinking much about computational efficiency):

> V _rchisq(1000,41+41)
> sig2 _(41*s2old + 41*s2)/V # inv-chisq(41+41,(41*s2old + 41*s2)/(41+41))
> X _model.matrix(olsfit)
> XTX _t(X)%*%X
> XTy _t(X)%*%radon$LogRadon
> SbIbeta0 _SigbetaI%*%beta0
> beta _matrix(0,ncol=1000,nrow=4)
> for (i in 1:1000) {
+ Vbeta _solve(XTX/sig2[i] + SigbetaI)
+ Vbeta _(Vbeta + t(Vbeta))/2       # guarantee symmetry (numerical instabilities)
+ Ch _t(chol(Vbeta))                 # look at help(chol)
+ bmean _Vbeta %*% (XTy/sig2[i] + SbIbeta0)
+ beta[,i] _bmean + Ch%*%rnorm(4)
+ }

The equal-tail 95% posterior credible intervals here can now be compared with the corresponding intervals for the non-informative prior analysis. For σ^{2} we obtain

NONINFORMATIVE:
> quantile(sig2,c(0.025,0.5,0.975))
       2.5%      50.0%     97.5%
 0.4170907 0.6317195 1.053781
INFORMATIVE:
> quantile(sig2,c(0.025,0.5,0.975))
       2.5%      50.0%      97.5%
 0.4879442 0.6552671 0.9014055

and for the regression coefficients \underline{β} the equal-tail 95% credible intervals look like this:

NONINFORMATIVE:
> t(apply(beta,1,quantile,c(0.025,0.975)))
                     2.5%      97.5%
(Intercept)  1.5070074 2.4275849
       Floor -0.9767665 0.3038866
        Clay -0.7473966 0.5542808
     Goodhue -0.6708612 0.5657630
INFORMATIVE:
> t(apply(beta,1,quantile,c(0.025,0.975)))
                     2.5%       97.5%
(Intercept)  1.6968706 2.28816000
       Floor -0.7624844 0.09801322
        Clay -0.4491324 0.43696194
     Goodhue -0.3421766 0.53141934