\documentstyle[11pt,fullpage,psfig,times]{article}

\parskip=0.5em
\parindent=0pt

\unitlength=0.75mm

\def\tight{\parsep=0pt\parskip=0pt\itemsep=0pt}

\frenchspacing

\long\def\comment#1{}

\let\bar\overline

\def\tp{^T}

\newcommand{\eq}[1]{(\ref{#1})}

\def\Var{\mbox{Var}\,}
\let\var\Var
\def\sigsq{{\hat\sigma^2_\epsilon}}

\def\Varsigma{{\mathbf{\not{\hspace*{-0.33em}\Sigma}}}}
\def\tr{\mbox{tr\,}}

\def\diag{\mbox{\it diag\,}}

\let\vc\Varsigma
\let\phi\varphi

\def\Corr{\mbox{Corr\,}}

\let\ul\mb
\let\by\times

\def\ul#1{\underline{#1}}
\def\ulem#1{\ul{\em #1}}

% \def\mb#1{{\mathbf #1}}
\let\mb\ul

\def\xtx{X\tp X}

\def\Re{{\cal R}}

\newenvironment{boxed}{\begin{center}\begin{tabular}{|p{5.5in}|}
\hline}{\\ \hline\end{tabular}\end{center}}

\newtheorem{definition}{Definition}
\newtheorem{theorem}{Theorem}
\newtheorem{corollary}{Corollary}
\newtheorem{lemma}{Lemma}

\begin{document}

\begin{center}\Large\bf 
        A Bayesian Approach to Regression Models
\end{center}

\section{References}

Once again it seems to make sense to refer you to the sources for this
material before presenting the material.  Sections
\ref{noninformative}, \ref{gls}, and \ref{priors} below follow in
spirit the presentation in portions of Chapter 8 of
\begin{itemize}
\item Gelman, A., Carlin, J. B., Stern, H. S. and Rubin,
D. B. (1995). {\em Bayesian data analysis.} New York: Chapman-Hall.
\end{itemize}
A much briefer account, that emphasizes the ``shrinkage'' aspect of
Bayesian estimation, may be found on pp. 38--50 of
\begin{itemize}
\item Carlin, B. P. and Louis, T. A. (1996). {\em Bayes and empirical
Bayes methods for data analysis}. New York: Chapman-Hall. 
\end{itemize}
The Carlin-Louis material follows closely that of the classic (but not
classical!)
reference
\begin{itemize}
\item Lindley, D. V. and Smith, A. F. M. (1972).  Bayes estimates for
the linear model (with discussion). {\em Journal of the Royal
Statistical Society, Series B, 34,} 1--41.
\end{itemize}
Finally, much more algebraic detail may be found in Chapters 8 and 9 of 
\begin{itemize}
\item Box, G. E. P. and Tiao, G. C. (1973). {\em Bayesian inference in
statistical analysis}.  Reading MA: Addison-Wesley.
\end{itemize}

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.
 
\section{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 
\[
        \mb y = X \mb \beta + \mb \epsilon, 
        ~~~\mb\epsilon \sim N(\mb 0, \sigma^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 $\mb \theta$:
\[
        \mb\theta\tp ~ = ~ (\mb \beta\tp, \sigma^2)\tp ~ = ~
        (\beta_0, \beta_1, \ldots, \beta_k, \sigma^2)\tp
\]
So far I have specified enough for you to write down the (or ``a'', if
you are fussy about it) {\em likelihood} of the data
\begin{eqnarray*}
        p(\mb y| X, \mb\theta) 
        & = & \prod _1^n \frac{1}{\sqrt{2\pi\sigma^2}}
                \exp \left[ - \frac12  (y_i-\mb x_i\tp\mb\beta)^2 / 
                \sigma^2 \right] \\
        & = & \frac{1}{(2\pi\sigma^2)^{n/2}} 
                \exp \left[ - \frac12 \sum_1^n (y_i-\mb x_i\tp\mb\beta)^2 /
                \sigma^2 \right] \\ 
        & = & \frac{1}{(2\pi)^{n/2}|\sigma^2 I|^{1/2}} 
                \exp \left[ - \frac12 (\mb y - X\mb \beta)\tp 
                (\sigma^2 I)^{-1} (\mb y - X\mb\beta) \right] ,
\end{eqnarray*}
where as usual $\mb x_i\tp$ is the $i^{th}$ {\em row} of the matrix
$X$. 

\ul{\em Classical} statistical analysis is driven by the
interpretation of $p(\mb y| X, \mb\theta)$ as a probability density
for the {\em random quantities} $\mb y$, given fixed and known $X$
\ul{\em and} fixed and {\em unknown} $\mb\theta$.  Inferences for
$\mb\theta$ are driven by ``what if'' style questions and answers:
\begin{itemize}
\item {\em Point estimation:} What value of $\mb\theta$ maximizes the
probability of observing the data $\mb y$ that we saw, rather than
some other data? (MLE's)
\item {\em Interval estimation:} Can we develop recipies $A(\mb y,
  X)$, $B(\mb y, X)$ such that $P[$ $A(\mb y, X)$ $\leq$ $\theta_j$
  $\leq$ $B(\mb y, X)$ $|$ $\mb\theta$ $]$ $>$ $0.95$, in repeated
  samples of $\mb y$?  Then apply the recipies to the $\mb y$ actually
  observed.  (Confidence intervals; the probability statement is about
  the recipies, not the particular numbers calculated)
\item {\em Hypothesis testing:}
\begin{itemize}
\item Can we develop test functions $T(\mb y, X)$ such that the
distribution $q_0[T(\mb y, X)|\mb\theta=\mb\theta_0]$ is
well-understood?  If we calculate $T(\mb y, X)$ for the data $\mb y$
actually observed and it is an outlier for the $q_0$ distribution this
is evidence against $\mb\theta=\mb\theta_0$ (general hypothesis testing)
\item 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)
\end{itemize}
\end{itemize}
Note that all of the reasoning about distributions is done
prospectively, as if predicting $\mb y$ on the basis of hypothetical
values of $\mb\theta$ before $\mb y$ is actually observed.  Then $\mb
y$ is observed and compared with the predictions made to ``learn''
which hypothetical values of $\mb\theta$ make sense for the observed
$\mb y$.

\ul{\em Bayesian} statistical analysis supplies an additional set of
modeling assumptions, that $\mb\theta\sim p(\mb\theta|X)$ is also
random, with marginal, or ``prior'', distribution\footnote{Since the
particular values of the entries in $X$ are not usually thought of as
constraining the ``true value'' of $\mb\theta$ in a regression
problem, it usually makes sense to make $\mb\theta$ independent of
$X$: $p(\mb\theta|X) = p(\mb\theta)$.}  $p(\mb\theta|X)$; and it is
driven by the interpretation of
\begin{equation}
        p(\mb\theta|\mb y, X) = \frac{p(\mb y| X, \mb\theta) p(\mb\theta|X)}
        { p(\mb y | X)} \propto p(\mb y| X, \mb\theta) p(\mb\theta|X)
        \label{BayesRule}
\end{equation}
as a probability density for the random quantities $\mb\theta$,
\ul{\em after} conditioning on observing $\mb y$, the data actually
observed. 
Inferences for $\mb\theta$ are driven by direct calculation on the
``posterior'' distribution $p(\mb\theta|\mb y, X)$:

\clearpage

\begin{itemize}
\item {\em Point estimation:} What value of $\mb\theta$ does
$p(\mb\theta|\mb y, X)$ assign the highest probability to? (posterior
mode) What is the average $\mb\theta$, given $\mb y$ and $X$?
(posterior mean)
\item {\em Interval estimation:} Can we develop recipies $A(\mb y,
X)$, $B(\mb y, X)$ such that in the posterior distribution $P[A(\mb y,
X) \leq \theta_j \leq B(\mb y, X)~|~\mb y, X]~>~0.95$?  (credible
intervals's; the probability statement is about the actual numbers calculated)
\item {\em Hypothesis testing:} 
\begin{itemize}
\item Can we calculate the distribution of additional data $\mb y^*$
that follows the same model we are fitting?  If so, does the test
function $T(\mb y, X)$ evaluated on the data actually observed look
like an outlier in this ``predictive'' distribution? (Posterior
predictive checks)
\item Is the probability of the data $P(\mb y|X, M_1)$ under model 1
greater than that under model 2, $P(\mb y |X, M_2)$?  (Bayes factors)
\end{itemize}
\end{itemize}

As indicated above, a common method of point estimation for Bayesian
inference is to guess that $\mb\theta$ is in the region of the
parameter space where the posterior $p(\mb\theta|\mb y, X) \propto
p(\mb y| X, \mb\theta) p(\mb\theta|X)$ is highest.  The likelihood
imposes one set of constraints on where $\mb\theta$ should be (near
the MLE, where the likelihood $p(\mb y| X, \mb\theta)$ is highest) and
the prior imposes another set of constraints (near the prior mode,
where the prior distribution $p(\mb\theta|X)$ is highest).  The
Bayesian point estimate is a compromise\footnote{This 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.}  between these two
sets of constraints.

For this reason, especially strong features of the likelihood $p(\mb
y|X,\mb\theta)$ or of the prior $p(\mb\theta|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(\mb y|X,\mb\theta)$ 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(\mb\theta|\mb y, X)$ rather than
$p(\mb y|X,\mb\theta)$).  As we move away from flat priors, the
special features of the prior will pull the inferences away from the
classical results.

\section{Analysis of the Ordinary Linear Regression Model with
Non-informative Priors}\label{noninformative}

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 {\em noninformative prior
distribution}\footnote{Noninformative 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(\mb y, \mb\theta|X)$, when the prior is improper, is not
possible.}  which is usually so featureless that every problem turns
into a likelihood-dominated problem.

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(\mb\beta,\log \sigma) \propto {\rm constant}
\]
or, after a calculus change-of-variables
\[
        p(\mb\beta,\sigma^2) \propto \sigma^{-2}
\]
Note that this is flat (i.e. constant, or uniform) in the regression
coefficients $\mb\beta$, so no constraints beyond those imposed by the
likelihood are placed on them.  It is not flat in $\sigma^2$, however:
values of $\sigma^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
$\mb\theta=(\mb\beta\tp,\sigma^2)\tp$.

\subsection{Setting up the joint posterior distribution}

Following Equation \eq{BayesRule}, and plugging in the likelihood and
default prior for linear regression, we then have
\begin{eqnarray}
        \log p(\mb\theta|\mb y, X) & = & {\rm (Constant)} + \log
        p(\mb y| X, \mb\theta) + \log p(\mb\theta|X) 
        \label{log-BR} \\
        & = & {\rm (Constant)} + 
         \left[ - n\log\sigma -  \frac12 (\mb y - X\mb \beta)\tp 
                (\mb y - X\mb\beta) / \sigma^2 \right] +
         \left[ \vphantom{\frac12} - 2\log\sigma \right] \nonumber\\
        & = & {\rm (Constant)} 
         - (n+2)\log\sigma - \frac12 (\mb y - X\mb \beta)\tp 
                (\mb y - X\mb\beta) / \sigma^2 , \nonumber
\end{eqnarray}
where the ``constant'' terms in each line are not necessarily the
same, but none depend on $\mb\theta$.  One can readily
verify\footnote{cf. Box and Tiao, p. 439.} (multiply it out, and use
the definition of $\hat{\mb\beta}\,$!) that
\[  
        (\mb y - X\mb \beta)\tp 
                (\mb y - X\mb\beta) / \sigma^2 = 
        (\mb y - X\hat{\mb \beta})\tp 
                (\mb y - X\hat{\mb\beta}) / \sigma^2 +
        (\mb \beta - \hat{\mb\beta})\tp X\tp X
                (\mb\beta - \hat{\mb\beta}) / \sigma^2
\]
where
\[
        \hat{\mb\beta} = (X\tp X)^{-1}X\tp \mb y
\]
is the usual OLS estimator of $\mb\beta$, and hence
\begin{eqnarray}
        \log p(\mb\theta|\mb y, X) & = & {\rm (Constant)} 
         - (n+2)\log\sigma - \frac12 (\mb y - X\hat{\mb \beta})\tp 
                (\mb y - X\hat{\mb\beta}) / \sigma^2 \nonumber \\
        & & ~~~~~~~~
         - \frac12 (\mb \beta - \hat{\mb\beta})\tp X\tp X
                (\mb\beta - \hat{\mb\beta}) / \sigma^2
        \label{posterior-kernel}
\end{eqnarray}

\subsection{Identifying the pieces}

The part of Equation \eq{posterior-kernel} that depends on $\mb\beta$
gives the shape\footnote{This is an {\em extremely} useful trick in
general for determining the so-called ``complete conditinals''
$p(\theta_1|\theta_2, \ldots, \theta_p, \mbox{\em data})$ in setting up
Markov Chain Monte Carlo (MCMC) calculations.} of $\log
p(\mb\beta|\sigma^2, \mb y, X)$; since this is a quadratic form in
$\mb\beta$ it follows that $p(\mb\theta|\sigma^2, \mb y, X)$ is
multivariate normal:
%\[
$$
        \mb\beta~|~\sigma^2, \mb y, X ~ \sim ~ N(\hat{\mb\beta},
        \sigma^2(X\tp X)^{-1})
        \eqno{(Bayes)}
$$
%\]
Note how this corresponds to the classical result
%\[
$$
        \hat{\mb\beta}~|~\mb\beta,\sigma^2, X ~\sim~ 
        N(\mb\beta,\sigma^2(X\tp X)^{-1})
        \eqno{(Classical)}
$$
%\]

Thinking now of 
\[
        p(\mb\beta,\sigma^2|\mb y, X) ~ = ~ p(\mb\beta|\sigma^2, \mb y, X)
        p(\sigma^2|\mb y, X) ,
\]
the marginal density for $p(\sigma^2|\mb y, X)$ could be obtained by
integrating $\mb\beta$ out of the joint density $p(\mb\beta|\sigma^2,
\mb y, X) p(\sigma^2|\mb y, X)$, but a cleverer thing to do is to
realize that this also says that
\[
        p(\sigma^2|\mb y, X) = 
        \frac{p(\mb\beta,\sigma^2|\mb y, X)}{p(\mb\beta|\sigma^2, \mb y, X)}
\]
for {\em any} value of $\mb\beta$.  After some algebra (plugging in
$\mb\beta\equiv\hat{\mb\beta}$ simplifies things I think) one obtains
that
\[
        \sigma^2 ~|~ \mb y, X ~ \sim ~ \mbox{Inv-$\chi^2$}\,(n-k-1,s^2)
\]
a scaled, inverse chi-squared\footnote{Definition: If $V$ $\sim$
$\chi^2_d$ then $ds^2/V$ $\sim$ Inv-$\chi^2\,(d,s^2)$.  See Appendix A
of Gelman et al. for formulas.} random variable,
where $s^2$ is the usual squared residual standard error,
\[
        s^2 = \frac{1}{n-k-1} (\mb y - X\hat{\mb\beta})\tp(\mb y -
        X\hat{\mb\beta}) \equiv \frac{RSS}{n-k-1}.
\]
This is equivalent to saying that
%\[
$$
        \left. \frac{(n-k-1)s^2}{\sigma^2} ~\right|~ 
        \mb y, X ~ \sim ~ \chi^2_{n-k-1} ,
        \eqno{(Bayes)}
$$
%\]
which corresponds to the classical result
%\[
$$
        \left. \frac{(n-k-1)s^2}{\sigma^2} ~\right|~ 
        \mb\beta,\sigma^2, X ~ \sim ~ \chi^2_{n-k-1} 
        \eqno{(Classical)}
$$
%\]

To summarize, if we use the default noninformative prior
$p(\mb\beta,\sigma^2) \propto \sigma^{-2}$, we obtain a posterior for
$\mb\beta$ and $\sigma^2$ that can be broken down as follows:
\begin{itemize} 
\item $p(\mb\beta|\sigma^2,\mb y, X)$ is $N(\hat{\mb\beta},
        \sigma^2(X\tp X)^{-1})$; 
\item $p(\sigma^2|\mb y, X)$ is $\mbox{Inv-$\chi^2$}\,(n-k-1,s^2)$.
\end{itemize}
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$.

\subsection{Further inferences}

It is also possible to show that the usual $t_{n-k-1}$ distributions
for each $\hat\beta_j$ in classical OLS regression carry over to the
same $t_{n-k-1}$ posterior distributions for each $\beta_j|\mb 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 $(\mb \beta, \sigma^2)$ and then ``throw away''
the $\sigma^2$'s simulated:
\begin{enumerate}
\item Simulate $\sigma^2$ from $p(\sigma^2|\mb y, X)$, i.e. from
        $\mbox{Inv-$\chi^2$}\,(n-k-1,s^2)$; 
\item Then simulate $\mb\beta$ from $p(\mb\beta|\sigma^2,\mb y, X)$,
        i.e. from $N(\hat{\mb\beta}, \sigma^2(X\tp X)^{-1})$.
\end{enumerate}
This simulation is based on Monte Carlo estimation of the integral 
\begin{equation}
        \int \int T(\mb\beta,\sigma^2) \, p(\mb\beta,\sigma^2|\mb y, X) 
        \, d(\sigma^2) \, d(\mb\beta) ~ = ~
        \int \int T(\mb\beta,\sigma^2) \, p(\mb\beta|\sigma^2,\mb y, X)
        \, p(\sigma^2|\mb y, X) \, d(\sigma^2) \, d(\mb\beta) ,
        \label{test-function}
\end{equation}
where $T(\mb\beta,\sigma^2)$ is any test function of interest.

Finally we turn to the problem of predicting new $y$ values.  Suppose
\[
        \tilde X = \left[\begin{array}{rrrr}
        1 & \tilde x_{11} & \cdots & \tilde x_{1k} \\
        1 & \tilde x_{11} & \cdots & \tilde x_{1k} \\
        \vdots & \vdots & \ddots & \vdots \\
        1 & \tilde x_{m1} & \cdots & \tilde x_{mk} \\
        \end{array}\right]
\]
are new values of the regressors, and
\[
        \tilde{\mb y} = \left[\begin{array}{cccc}
        \tilde y_1 & \tilde y_2 & \cdots & \tilde y_m
        \end{array}\right]\tp
\] 
are the corresponding new $y$ values.  The linear regression model
says as usual that $p(\tilde{\mb y} | \mb\beta, \sigma^2, \tilde X)$
will be $N(\tilde X\mb\beta,\sigma^2 I_{m\times m})$.  We can then use
Equation \eq{test-function} with $T(\mb\beta,\sigma^2) = p(\tilde{\mb y} |
\mb\beta, \sigma^2, \tilde X)$, 
\[
        p(\tilde{\mb y} | \mb y, X, \tilde X) = \int\int p(\tilde{\mb y} |
        \mb\beta, \sigma^2, \tilde X, \mb y, X) p(\mb\beta|\sigma^2, \mb y, X)
        p(\sigma^2|\mb y, X) \, d(\sigma^2) d(\mb\beta)
\]
in two ways:
\begin{itemize}
\item To derive the exact ``predictive'' distribution of $\tilde{\mb
y}$.  It turns out for example that (omitting $X$ and $\tilde X$ from
the conditioning to save space)
\begin{eqnarray*}
        E[\tilde{\mb y} ~|~ \sigma^2, \mb y] 
        & = & E[E(\tilde{\mb y} ~|~ 
                \mb\beta,\sigma^2,\mb y)~|~\sigma^2,\mb y] \\
        & = & E[X\mb\beta ~|~\sigma^2,\mb y] \\
        & = & \tilde X\hat{\mb\beta}
\end{eqnarray*}
and similarly
\begin{eqnarray*}
        \Var(\tilde{\mb y} | \sigma^2, \mb y) 
        & = & E[\Var(\tilde{\mb y}~|~\mb\beta,\sigma^2,\mb y) ~|~ 
        \sigma^2, \mb y] + \Var(E[\tilde{\mb y}~|~\mb\beta,\sigma^2,\mb y] ~|~ 
        \sigma^2, \mb y) \\
        & = & E[\sigma^2 I_{m\times m}~|~\sigma^2, \mb y] 
        + \Var(\tilde X\mb\beta~|~\sigma^2, \mb y) \\
        & = & \sigma^2(I_{m\times m} + \tilde X (X\tp X)^{-1} \tilde X\tp)
\end{eqnarray*}
Note again how these results correspond with the classical ones: if
$\sigma^2$ is known, $\tilde{\mb y}$ will follow a multivariate normal
distribution with this mean and variance; if $\sigma^2$ is not known,
then the corresponding $t$ distribution results hold.
\item To obtain simulated values of $\tilde{\mb y}$ by simulating from
the joint posterior density of $(\tilde{\mb y}, \mb\beta, \sigma^2)$
and throwing away the simulated $\mb\beta$ and $\sigma^2$ values:
\begin{enumerate}
\item Simulate $\sigma^2$ from $p(\sigma^2|\mb y, X)$, i.e. from
        $\mbox{Inv-$\chi^2$}\,(n-k-1,s^2)$; 
\item Then simulate $\mb\beta$ from $p(\mb\beta|\sigma^2,\mb y, X)$,
        i.e. from $N(\hat{\mb\beta}, \sigma^2(X\tp X)^{-1})$;
\item Finally simulate $\tilde{\mb y}$ from $p(\tilde{\mb
y}|\mb\beta,\sigma^2, \tilde X)$, i.e. from $N(\tilde
X\mb\beta,\sigma^2 I_{m\times m})$. 
\end{enumerate}
\end{itemize}

\subsection{Example: Radon measurements}

Gelman et al. (1995, pp. 189 and 263) consider the data in Table
\ref{radon-table} 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
\begin{verbatim}
> 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
\end{verbatim}

\begin{table}
\begin{center}
\begin{tabular}{@{}c|c@{}}
\begin{tabular}[t]{rlc}
Radon   & County        & Floor \\
\hline
\hline
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 \\
\hline
\end{tabular}
&
\begin{tabular}[t]{rlc}
Radon   & County        & Floor \\
\hline
\hline
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 \\
\hline
\end{tabular} 
\end{tabular}
\end{center}
\caption{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.}\label{radon-table}
\end{table}

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 $\sigma^2$ and $\mb\beta$ 
\begin{verbatim}
> 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.
> 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
\end{verbatim}

Now we could give analytic results here but it is more informative to
simulate:
\begin{verbatim}
> 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
\end{verbatim}
Plotting the results takes more work than the simulations:
\begin{verbatim}
> 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)))
     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
\end{verbatim}
The results are shown in Figure \ref{ols-vs-post}.

\begin{figure}
\centerline{\psfig{figure=ols-vs-post.ps,height=3.5in}}
\caption{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 $\beta$'s) or point estimate (for $\sigma^2$) from OLS
analysis.}\label{ols-vs-post}
\end{figure}

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
$\tilde y$ as follows
\begin{verbatim}
> 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
\end{verbatim}
If the measurement is taken on the first floor, the simulated
distribution can be obtained as follows
\begin{verbatim}
> 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
\end{verbatim} 
The results are shown in Figure \ref{pred-radon}.

\begin{figure}
\centerline{\psfig{figure=pred-radon.ps,height=3.5in}}
\caption{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.}\label{pred-radon}
\end{figure}

\section{Bayesian GLS}\label{gls}

A Bayesian analysis corresponding to GLS (generalized least squares)
proceeds as follows\footnote{This follows very closely Section 8.7 of
Gelman et al.}.  We will assume that
\begin{eqnarray*}
        \mb y~|~\mb \beta, \Sigma_y, X & \sim & N(X\mb\beta,\Sigma_y) \\
        \mbox{and } p(\mb\beta|\Sigma_y,X) & \propto & \mbox{constant} .
\end{eqnarray*}
Let us assume $\Sigma_y$ is known; this lets us calculate
$p(\mb\beta|\Sigma_y, \mb y, X)$.  Let $\Sigma_y^{1/2}$ be a Cholesky
square root of $\Sigma_y$ (i.e., a triangular matrix such that
$[\Sigma_y^{1/2}]\tp\Sigma_y^{1/2} = \Sigma_y$).  Then as with our
discussion of WLS, we can premultiply by $\Sigma_y^{-1/2}$ and obtain
\[
        \Sigma_y^{-1/2} \mb y ~|~  \beta, \Sigma_y
        ~ \sim ~ N(\Sigma_y^{-1/2} X\mb\beta,I)
\]
Applying the ordinary regression results above to $\mb z =
\Sigma_y^{-1/2} \mb y$ we immediately obtain
\[
        \mb\beta~|~\Sigma_y, \mb y, X ~ \sim ~ N(\hat{\mb\beta},
        V_\beta)
\]
where
\begin{eqnarray*}
        \hat{\mb\beta} 
        & = & (X\tp \Sigma_y^{-1} X)^{-1} X\tp \Sigma_y^{-1} \mb y \\
        V_\beta & = & (X\tp \Sigma_y^{-1} X)^{-1} 
\end{eqnarray*}
These calculations turn out to be very useful in Section \ref{priors}
below, where we are considering Bayesian linear regression with
informative, conjugate prior information.

For completness we also indicate the form of $p(\Sigma_y|\mb y, X)$,
but we don't intend to use it in the rest of these notes.  As before,
we use
\[
        p(\mb\beta, \Sigma_y~|~\mb y, X) ~=~ p(\mb\beta~|~\Sigma_y,\mb y, X)
        ~ p(\Sigma_y~|~\mb y, X)
\]
to see that, for {\em any} value of $\mb\beta$,
\[
        p(\Sigma_y~|~\mb y, X) ~ = ~ \frac{p(\mb\beta, \Sigma_y~|~\mb y, X)}
        {p(\mb\beta~|~\Sigma_y,\mb y, X)} 
        ~ \propto ~ \frac{p(\Sigma_y) N(\mb y~|~X\mb\beta, \Sigma_y)}
        {N(\mb\beta~|~X\hat{\mb\beta},V_\beta)}
\]
where $p(\Sigma_y)$ is now the prior on $\Sigma_y$ (the prior on
$\mb\beta$ is still flat).  Substituting
$\mb\beta\equiv\hat{\mb\beta}$ again and simplifying, we see that
\begin{equation}
        p(\Sigma_y|\mb y, X) ~ \propto ~
        p(\Sigma_y) |V_\beta|^{1/2}\exp\left[ - \frac12
                (\mb y - X\hat{\mb\beta})\tp \Sigma_y^{-1} 
                (\mb y - X\hat{\mb\beta}) \right] .
        \label{post-sig}
\end{equation}
This is relatively easy to write down but difficult to work
with---both because $p(\Sigma_y)$ is difficult to specify (or
``elicit'') and because $\Sigma_y$ is wrapped up in the computation of
both $\hat{\mb\beta}$ and $V_\beta$---except in very special
circumstances, such as WLS (weighted least squares) where at least
$\Sigma_y = \sigma^2 \mbox{diag}(1/w_1, \ldots, 1/w_n)$ only depends
on a single parameter.

\comment{

The predictive distribution of new data $\tilde{\mb y}$ is more
complicated but follows a similar pattern to what we have seen before;
for example if we specify
\[
        \left.\left[\begin{array}{c}
        \mb y \\ \tilde{\mb y}
        \end{array}\right] ~ \right| ~ \tilde X, X,\mb\beta,\sigma^2
        ~ \sim ~ N\left(
        \left[\begin{array}{c}
        X\mb\beta \\ \tilde X\mb\beta
        \end{array}\right],~
        \left[\begin{array}{cc}
        \Sigma_y & \Sigma_{y,\tilde y} \\
        \Sigma_{\tilde y,y} & \Sigma_{\tilde y}
        \end{array}\right]\right)
\]
then $(\tilde{\mb y} ~|~ \mbox{[everything else]})$ will be
multivariate normal with mean and variance
\begin{eqnarray*}
        E[\tilde{\mb y}~|~ \mb y, \mb\beta, \mbox{ the $\Sigma$'s}] & = &
        \tilde X\mb\beta + \Sigma_{y, \tilde y}\Sigma_y^{-1}(\mb y-X\mb\beta)\\
        \Var(\tilde{\mb y}~|~ \mb y, \mb\beta, \mbox{ the $\Sigma$'s}) & = &
        \Sigma_{\tilde y} - \Sigma_{y,\tilde y}\Sigma_y^{-1}\Sigma_{y,\tilde y}
\end{eqnarray*}

}

\section{Prior Information}\label{priors}

\subsection{Prior information for $\mathbf\mb\beta; ~\Sigma_y$ known}

We start by assuming that the variance-covariance matrix $\Sigma_y$
is known, and that we have multivariate normal prior information on
$\mb\beta$: 
\begin{eqnarray*}
        \mb y~|~\mb \beta, \Sigma_y & \sim & N(X\mb\beta,\Sigma_y) \\
        \mbox{and } \mb\beta & \sim & N(\mb\beta_0,\Sigma_\beta)
\end{eqnarray*}
where $\Sigma_y$, $\mb\beta_0$, and $\Sigma_\beta$ are known and fixed
in advance.  

\begin{quote}
{\em Aside:} Another way to say the same thing is that
\begin{equation}
\left.
\begin{array}{rcl}
        \mb y & = & X\mb\beta + \mb\epsilon, ~~~ \mb\epsilon \sim
        N(\mb 0, \Sigma_y) \\
        \mb \beta & = & \mb\beta_0 + \mb\delta, ~~~ \mb\delta \sim
        N(\mb 0, \Sigma_\beta)
\end{array}
\right\}
\label{hlm}
\end{equation}
from which it immediately follows that
\begin{eqnarray*}
        \mb y & = & X\mb\beta_0 + X\mb\delta + \mb\epsilon \\
        \mbox{i.e. } \mb y & \sim & N(X\mb\beta_0, X\Sigma_\beta X\tp
        + \Sigma_y)
\end{eqnarray*}
\begin{itemize}
\item
If $\mb\beta_0$ were {\em not} known, it could now be estimated by
(non-Bayesian) WLS.  If we then proceed with a Bayesian analysis for
$\mb\beta$ pretending that $\mb\beta_0$ was fixed in advance to be
$\hat{\mb\beta}_0$ from the WLS analysis, then the resulting analysis
would be an example of an {\em empirical Bayes} analysis.

\item
If $\mb\beta_0$ had the further linear structure
$\mb\beta_0=Z\mb\gamma$, the resulting model \eq{hlm} would be called
a (two-level) {\em hierarchical linear model} (HLM).
\end{itemize}
{\em End of vocabulary lesson.}
\end{quote}

Continuing under the assumption that $\Sigma_y$, $\mb\beta_0$, and
$\Sigma_\beta$ are known and fixed in advance, we know that 
\[
        p(\mb\beta|\Sigma_y,\mb y,X) ~ \propto ~ N(y|X\mb\beta,\Sigma_y)
        \cdot N(\mb\beta|\mb\beta_0,\Sigma_\beta)
\]
and so
\begin{eqnarray*}
        \log p(\mb\beta|\Sigma_y,\mb y, X) & = & ({\rm Constant}) 
        - \frac12 (y - X\mb\beta)\tp \Sigma_y^{-1} (y - X\mb\beta)
        - \frac12 (\mb\beta - \mb\beta_0)\tp \Sigma_\beta^{-1}
        (\mb\beta - \mb\beta_0)  \\
        & = & ({\rm Constant})  - \frac12
        \left(
        \left[\begin{array}{ll}
        \mb y \\ \mb\beta_0
        \end{array}\right]
        -
        \left[\begin{array}{rr}
        X\mb\beta \\ I\mb\beta
        \end{array}\right]
        \right)\tp
        \left(
        \begin{array}{cc}
        \Sigma_y^{-1} & 0 \\
        0 & \Sigma_\beta^{-1}
        \end{array}
        \right)
        \left(
        \left[\begin{array}{ll}
        \mb y \\ \mb\beta_0
        \end{array}\right]
        -
        \left[\begin{array}{rr}
        X\mb\beta \\ I\mb\beta
        \end{array}\right]
        \right) \\
        & = & ({\rm Constant})  - \frac12
        (\mb y_* - X_*\mb\beta)\tp \Sigma_*^{-1} (\mb y_* - X_*\mb\beta)
\end{eqnarray*}
where the ``constant's'' do not depend on $\mb\beta$, where $I$ is
$(k+1)\times(k+1)$ identity matrix, and where
\[
        \mb y_* =       \left[\begin{array}{ll}
        \mb y \\ \mb\beta_0
        \end{array}\right]
        ~~~~~~~~
        X_* = 
        \left[\begin{array}{cc}
        X \\ I
        \end{array}\right]
        ~~~~~~~~
        \Sigma_* = 
        \left(
        \begin{array}{cc}
        \Sigma_y & 0 \\
        0 & \Sigma_\beta
        \end{array}
        \right)
\]
are ``augmented'' data and design matrices.  In other words {\em
regression with prior information on $\mb\beta$ is like regression
with the ``augmented'' data $\mb y_*$ and design matrix $X_*$ and
noninformative prior on $\mb\beta$, and known variance-covariance
matrix $\Sigma_*$.}  In other words, the ``prior'' information can be
interpreted as additional ``observations'' in a noninformative-prior
Bayesian analysis of the linear regression model
\[
        \mb y_* = X_*\mb\beta + \mb\epsilon_*, ~~~ \mb\epsilon_* \sim
        N(\mb 0, \Sigma_*)
\]

We can immediately apply the results on ``Bayesian GLS'' analysis (Section
\ref{gls}) to conclude that 
\[
        \mb\beta|\Sigma_y, \mb y, X ~ \sim ~ N(\hat{\mb\beta},V_\beta)
\]
with 
\begin{eqnarray*}
        \hat{\mb\beta} & = & (X_*\tp \Sigma_*^{-1} X_*)^{-1} 
                X_*\tp \Sigma_*^{-1} y_* \\
        & = & \left(
        \left[\begin{array}{cc}
        X\tp & I\tp
        \end{array}\right]
        \left[\begin{array}{cc}
        \Sigma_y^{-1} & 0 \\
        0 & \Sigma_\beta^{-1}
        \end{array}\right]
        \left[\begin{array}{c}
        X \\ I
        \end{array}\right]
        \right)^{-1}
        \left[\begin{array}{cc}
        X\tp & I\tp
        \end{array}\right]
        \left[\begin{array}{cc}
        \Sigma_y^{-1} & 0 \\
        0 & \Sigma_\beta^{-1}
        \end{array}\right]
        \left[\begin{array}{l}
        \mb y \\ \mb\beta_0
        \end{array}\right] \\
        & = & \left(X\tp \Sigma_y^{-1} X + \Sigma_\beta^{-1}\right)^{-1}
                \left[X\tp\Sigma_y^{-1}\mb y 
                + \Sigma_\beta^{-1} \mb\beta_0\right]\\
        \\
        \mbox{and } V_\beta & = &  (X_*\tp \Sigma_*^{-1} X_*)^{-1} \\
        & = & \left(X\tp \Sigma_y^{-1} X + \Sigma_\beta^{-1}\right)^{-1}
\end{eqnarray*}
If some of the $\beta_j$'s should have noninformative prior
distributions, we set the corresponding rows and columns of
$\Sigma_*^{-1}$ (or equivalently, of $\Sigma_\beta^{-1}$) equal to all
zeros.

The case of ordinary regression with iid priors on $\mb\beta$ gives
some insight into the ``shrinkage'' interpretation of Bayesian
estimates.  This case occurs when
\[
        \Sigma_y = \sigma^2 I_{n\times n} ~~~~~~~ \Sigma_\beta =
        \tau^2 I_{(k+1)\times(k+1)}
\]
so that the posterior mean and variance reduce to
\begin{eqnarray*}
        \hat{\mb\beta} 
        & = & \left(X\tp X + \frac{\sigma^2}{\tau^2} I\right)^{-1} 
                \left[X\tp \mb y + \frac{\sigma^2}{\tau^2} \mb\beta_0\right]\\
        \mbox{and } V_\beta & = & \sigma^2
                \left(X\tp X + \frac{\sigma^2}{\tau^2} I\right)^{-1} 
\end{eqnarray*}
There are two interesting things to note about this formulation:
\begin{itemize}
\item As the prior information gets weaker
($\tau^2\rightarrow\infty$), $\hat{\mb\beta}\rightarrow (X\tp
X)^{-1}X\tp\mb y$, the classical estimate, as in Section
\ref{noninformative}; and as the prior information gets stronger
($\tau^2\rightarrow0$), $\hat{\mb\beta}\rightarrow\mb\beta_0$, the
prior mean.  This is the usual ``shrinkage'' behavior.
\item When $\mb\beta_0 = \mb0$, the estimate is of the form
\[
        \hat{\mb\beta} ~ = ~ \left(X\tp X + \lambda I\right)^{-1} 
                X\tp \mb y
\]
(with $\lambda=\sigma^2/\tau^2$), i.e. it is the {\em ridge regression}
estimate.
\end{itemize}

\subsection{Prior information for $\Sigma_y$}

Now that we have beaten $p(\mb\beta|\Sigma_y,\mb y, X) =
N(\mb\beta|\hat{\mb\beta},V_\beta)$ into the ground, let's complete
the formulation of a joint posterior by specifying $p(\Sigma_y|\mb y,
X)$.  We will only do this in the case of ordinary regression,
$\Sigma_y = \sigma^2 I$, with a general normal prior $\mb\beta \sim
N(\mb\beta_0,\Sigma_\beta)$.  For this case, we saw in Section
\ref{noninformative} that the natural posterior with a non-informative
prior was an Inv-$\chi^2$.  Therefore it will not be a surprise that
the conjugate prior family for $\sigma^2$ in Bayesian ordinary
regression is
\[
        \mbox{Inv-$\chi^2$}\,(\sigma^2|n_0,\sigma_0^2) .
\]
Following the same calculations as before [see especially Equation
\eq{post-sig}], using the joint prior
\[
        p(\mb\beta,\sigma^2|X) ~ = ~ N(\mb\beta|\mb\beta_0,\Sigma_\beta) \cdot
                \mbox{Inv-$\chi^2$}\,(\sigma^2|n_0,\sigma_0^2)
\]
it follows that 
\begin{itemize}
\item $p(\mb\beta ~|~ \sigma^2,\mb y, X)$ $=$
$N(\mb\beta|\hat{\mb\beta},V_\beta)$ where
\begin{eqnarray}
        \hat{\mb\beta} & = & (X_*\tp \Sigma_*^{-1} X_*)^{-1} 
                X_*\tp \Sigma_*^{-1} y_* 
        ~ = ~ \left(X\tp X/\sigma^2 + \Sigma_\beta^{-1}\right)^{-1}
                \left[X\tp \mb y/\sigma^2 + \Sigma_\beta \mb\beta_0\right]
        \label{mean} \\
        V_\beta & = & (X_*\tp \Sigma_*^{-1} X_*)^{-1} 
        ~ = ~ \left(X\tp X/\sigma^2 + \Sigma_\beta^{-1}\right)^{-1}
        \label{var}
\end{eqnarray}
and
\item $p(\sigma^2|\mb y, X)$ is given by
\begin{equation}
        p(\sigma^2|\mb y, X) 
        ~ = ~ \mbox{Inv-$\chi^2$}\,\left(\sigma^2\left|n_0+n',
        \frac{n_0\sigma_0^2 + n's^2}{n_0 + n'}\right.\right)
        \label{sig}
\end{equation}
where $s^2 = (\mb y_* - X_*\hat{\mb\beta})\tp \Sigma_*^{-1} (\mb y_* -
X_*\hat{\mb\beta})/n'$; $n' = n_* + k-1$; and $n_* = rank(\Sigma_*)
= n + rank(\Sigma_\beta^{-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.
\end{itemize}

\subsection{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
%
% oldfit obtained by taking a bootstrap sample of the radon data:
% > perm _ sample(41,replace=T)
% > perm
%  [1]  1  9 23 37 17 18 15 27 19 30 34 19 36 31 13  4 32 14 34 40 22 36 18 29 10  2
% [27] 41 13 22 40 28 14  1  9 29 10  6 31  9  2 19
% > sort(perm)
%  [1]  1  1  2  2  4  6  9  9  9 10 10 13 13 14 14 15 17 18 18 19 19 19 22 22 23 27
% [27] 28 29 29 30 31 31 32 34 34 36 36 37 40 40 41
% > oldrad _ radon[perm,]
% > formula(olsfit)
% LogRadon ~ Floor + Clay + Goodhue
% > oldfit _ lm(LogRadon ~ Floor + Clay + Goodhue,data=oldrad)
%
\begin{verbatim}
> 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
\end{verbatim}
For our current data we might like to specify prior distributions
based on the earlier study:
\begin{eqnarray*}
        \mb\beta & \sim & N(\hat{\mb\beta}_{old},s^2_{old}(X_{old}\tp
        X_{old})^{-1})\\ 
        \sigma^2 & \sim & \mbox{Inv-$\chi^2$}\,(41,s_{old}^2)
\end{eqnarray*}
These can be obtained as follows
\begin{verbatim}
> 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
\end{verbatim}
Following Equations \eq{mean}, \eq{var} and \eq{sig}, 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):
\begin{verbatim}
> 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)
+ }
\end{verbatim}
The equal-tail 95\% posterior credible intervals here can now be
compared with the corresponding intervals for the non-informative
prior analysis.  For $\sigma^2$ we obtain
\begin{verbatim}
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
\end{verbatim}
and for the regression coefficients $\mb\beta$ the equal-tail 95\%
credible intervals look like this:
\begin{verbatim}
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
\end{verbatim}

\end{document}









