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

\parskip=0.5em
\parindent=0pt

\frenchspacing

\long\def\comment#1{}

\def\tp{^T}

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

\def\vc#1{{\underline{#1}}}

\def\Var{\mbox{Var}\,}

\def\Varsigma{\mb{\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}}

\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 
        Odds and Ends from the Theory of Linear Regression
\end{center}

\section{Introduction}

In this handout I will say a few words about some topics that were
missed earlier in the ``theoretical'' part of the course:

\begin{itemize}
\item Reparametrizing, and Testing Contrasts 
\item Gauss-Markov and MLE stuff 
% \item Random Regressors and Specification Error 
\end{itemize}

\section{Reparametrizing, and Testing Contrasts}

We have seen that the ``dummy variables'' linear regression version of
an ANOVA model is always overparametrized.  To see how ANOVA
reparametrizations overcome this, we now examine the one-way ANOVA
model from a linear regression point of view.

For simplicity, consider a one-way ANOVA problem for a factor with
three levels, $a$, $b$ and $c$, and two observations per cell.  Let 
\[
X_1 = \left\{\begin{array}{rr}1 & \mbox{at level $a$} \\
                              0 & \mbox{else} \end{array}\right.
~~~~~~~~~~~~
X_2 = \left\{\begin{array}{rr}1 & \mbox{at level $b$} \\
                              0 & \mbox{else} \end{array}\right.
~~~~~~~~~~~~
X_3 = \left\{\begin{array}{rr}1 & \mbox{at level $c$} \\
                              0 & \mbox{else} \end{array}\right.
\]
so the linear regression model is
\begin{equation}
        y_{ij} = \beta_0 + X_{1j}\beta_1 + X_{2j}\beta_2 
         + X_{3j}\beta_3 + \epsilon_{ij}
        \label{basic-1}
\end{equation}
or in matrix form
\[
        y = 
        \left[\begin{array}{c}
        y_{11} \\ y_{12} \\ y_{21} \\ y_{22} \\ y_{31} \\ y_{32}
        \end{array}\right]
        =
        X \beta + \epsilon
        =
        \left[\begin{array}{cccc}
        1 & 1 & 0 & 0 \\
        1 & 1 & 0 & 0 \\
        1 & 0 & 1 & 0 \\
        1 & 0 & 1 & 0 \\
        1 & 0 & 0 & 1 \\
        1 & 0 & 0 & 1 
        \end{array}\right]
        \left[\begin{array}{c}
        \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 
        \end{array}\right]
        +
        \left[\begin{array}{c}
        \epsilon_{11} \\ \epsilon_{12} \\ \epsilon_{21} \\ 
        \epsilon_{22} \\ \epsilon_{31} \\ \epsilon_{32}
        \end{array}\right]
\]
Basically we want to know how to introduce a linear constraint, to get
a new set of parameters that can be estimated uniquely (technical term:
the new parameters are {\em identifiable}).  

Another common application of reparametrization is to estimate a {\em
contrast}, or linear combination of the parameters of the model:
whether or not the contrast equals zero corresponds to whether or not
a null hypothesis of some hypothesis test of interest is true.


\subsection{Matrix algebra details}

We want to consider reparametrizations of the model $y=X\beta +
\epsilon$, that is, we want to consider new sets of parameters
\begin{equation}
        \left[\begin{array}{c}
        \gamma_1 \\ \gamma_2 \\ \gamma_3 
        \end{array}\right]
        =
        L \cdot
        \left[\begin{array}{c}
        \beta_1 \\ \beta_2 \\ \beta_3 \\ \beta_4
        \end{array}\right]
        \label{reparam}
\end{equation}
where 
\[
        L = \left[\begin{array}{cccc}
        \ell_{11} & \ell_{12} & \ell_{13} & \ell_{14} \\
        \ell_{21} & \ell_{22} & \ell_{23} & \ell_{24} \\
        \ell_{31} & \ell_{32} & \ell_{33} & \ell_{34} 
        \end{array}\right] .
\]
(Note that contrast matrices are often presented in transposed form,
$C=L\tp$.)
Our goal is to choose a new design matrix $Z$ so that
\[
        y = X\beta + \epsilon = Z\gamma + \epsilon
\]
i.e. the two regression problems are exactly the same (more formally:
$X$ and $Z$ have the same column space).

How do we find $Z$? Note that, by \eq{reparam}, $X\beta = ZL\beta$, so
it must be that
\[
        X ~ = ~ ZL .
\]
Now, we would like to solve this by ``multiplying on the right by
$L^{-1}$'', but unfortunately $L$ is not square, so it doesn't have an
inverse.  Instead, notice what happens if we multiply on the right by
$L\tp(LL\tp)^{-1}$ (this is called a ``generalized inverse''):
\begin{eqnarray*}
        X & = & ZL \\
        X\cdot L\tp(LL\tp)^{-1} & = & ZL\cdot L\tp(LL\tp)^{-1} \\
        X\cdot L\tp(LL\tp)^{-1} & = & Z 
\end{eqnarray*}
so now we have a way of finding $Z$ from $X$ and $L$, namely $Z = X
L\tp(LL\tp)^{-1}$.  Finally, since
\[
        y = X\beta + \epsilon = ZL\beta + \epsilon = Z\gamma + \epsilon
\]
it follows that testing whether the coefficient $\gamma_j=0$ in the model
$y = Z\gamma + \epsilon$ is the same as testing whether the corresponding
linear combination $\ell_{j0}\beta_0 + \ell_{j1}\beta_1 + \cdots +
\ell_{jk'}\beta_{k'} = 0$ in the model $y = X\beta + \epsilon$.

To summarize, our goal in understanding reparametrization is:

\hfil\fbox{\parbox{5in}{
\begin{description}
\item[Step 1.] Figure out what ``constraint(s)'' are being imposed on the
$\beta$'s to get rid of the linear dependency (or code the contrast);
\item[Step 2.] Figure out how to write that as something like $\gamma
= L\beta$;
\item[Step 3.] Compute the design matrix for the $\gamma$'s as $Z =
XL\tp(LL\tp)^{-1}$.
\item[Step 4.] The new regression is then $y = Z\gamma + \epsilon$.
\end{description}
}}

Two notes before going on:
\begin{itemize}
\item In many simple problems, you can intuitively get from {\bf Step
1} to {\bf Step 4} without fooling around with the $L$ matrix in {\bf
Steps 2} and {\bf 3}.
\item After {\bf Step 4} it is often helpful to write the columns of
$Z$ in terms of the 0/1 dummy variables $X_1$, $X_2$, $X_3$, and
re-expressing the cell means $\mu_1$, $\mu_2$ and $\mu_3$ in terms of
the new parameters.
\end{itemize}

\subsection{Some illustrations}

\paragraph{Cell means model (delete the first column).} This is easy.
Setting $\beta_0=0$ is the same as deleting the first column in the
$X$ matrix; the matrix form of the regression becomes

this becomes
\[
        \left[\begin{array}{c}
        y_{11} \\ y_{12} \\ y_{21} \\ y_{22} \\ y_{31} \\ y_{32}
        \end{array}\right]
        =
        \left[\begin{array}{ccc}
        1 & 0 & 0 \\
        1 & 0 & 0 \\
        0 & 1 & 0 \\
        0 & 1 & 0 \\
        0 & 0 & 1 \\
        0 & 0 & 1 
        \end{array}\right]
        \left[\begin{array}{c}
        \beta_1 \\ \beta_2 \\ \beta_3 
        \end{array}\right]
        +
        \left[\begin{array}{c}
        \epsilon_{11} \\ \epsilon_{12} \\ \epsilon_{21} \\ 
        \epsilon_{22} \\ \epsilon_{31} \\ \epsilon_{32}
        \end{array}\right]
\]
and \eq{basic-1} becomes
\[
        y_{ij} = X_{1j}\beta_1 + X_{2j}\beta_2 
         + X_{3j}\beta_3 + \epsilon_{ij}
\]
and it is easy to see in the resulting model, $\mu_1 = \beta_1$,
$\mu_2 = \beta_2$ and $\mu_3 = \beta_3$, the cell means model.

Let's see how we would do this by setting up the $L$ matrix and
calculating $Z = XL\tp(LL\tp)^{-1}$:

We want $\beta_0=0$ and
\begin{eqnarray*}
        \mu_1 & = & \beta_1 \\
        \mu_2 & = & \beta_2 \\
        \mu_3 & = & \beta_3
\end{eqnarray*}
or 
\[
        \left[\begin{array}{c}
        \mu_1 \\ \mu_2 \\ \mu_3 
        \end{array}\right]
        = 
        \left[\begin{array}{cccc}
        0 & 1 & 0 & 0 \\
        0 & 0 & 1 & 0 \\
        0 & 0 & 0 & 1
        \end{array}\right]
        \left[\begin{array}{c}
        \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 
        \end{array}\right]
        = L\beta
\]
Hence
\[
        Z = X L\tp(LL\tp)^{-1} = 
        \left[\begin{array}{cccc}
        1 & 1 & 0 & 0 \\
        1 & 1 & 0 & 0 \\
        1 & 0 & 1 & 0 \\
        1 & 0 & 1 & 0 \\
        1 & 0 & 0 & 1 \\
        1 & 0 & 0 & 1 
        \end{array}\right]
        \left[\begin{array}{ccc}
        0  & 0  & 0\\
        1  & 0  & 0\\
        0  & 1  & 0\\
        0  & 0  & 1
        \end{array}\right]
        =
        \left[\begin{array}{ccc}
        1 & 0 & 0 \\
        1 & 0 & 0 \\
        0 & 1 & 0 \\
        0 & 1 & 0 \\
        0 & 0 & 1 \\
        0 & 0 & 1 
        \end{array}\right]
\]
which is the same design matrix that we got ``by eye'' above.

\paragraph{Delete the last column.} Some computer packages set
$\beta_3=0$ instead of $\beta_0=0$. If we set $\beta_3=0$, this is
like deleting the last column in the matrix formulation.
\[
        \left[\begin{array}{c}
        y_{11} \\ y_{12} \\ y_{21} \\ y_{22} \\ y_{31} \\ y_{32}
        \end{array}\right]
        =
        \left[\begin{array}{ccc}
        1 & 1 & 0  \\
        1 & 1 & 0  \\
        1 & 0 & 1  \\
        1 & 0 & 1  \\
        1 & 0 & 0  \\
        1 & 0 & 0  
        \end{array}\right]
        \left[\begin{array}{c}
        \beta_0 \\ \beta_1 \\ \beta_2 
        \end{array}\right]
        +
        \left[\begin{array}{c}
        \epsilon_{11} \\ \epsilon_{12} \\ \epsilon_{21} \\ 
        \epsilon_{22} \\ \epsilon_{31} \\ \epsilon_{32}
        \end{array}\right]
\]
and \eq{basic-1} becomes
\[
        y_{ij} = \beta_0 + X_{1j}\beta_1 + X_{2j}\beta_2 
          + \epsilon_{ij}
\]
We can also do this by setting up $L$ and computing
$Z=XL\tp(LL\tp)^{-1}$.  We want to set $\beta_3=0$ and to define
\[
        \left[\begin{array}{c}
        \gamma_0 \\ \gamma_1 \\ \gamma_2 
        \end{array}\right]
        = 
        \left[\begin{array}{cccc}
        1 & 0 & 0 & 0 \\
        0 & 1 & 0 & 0 \\
        0 & 0 & 1 & 0
        \end{array}\right]
        \left[\begin{array}{c}
        \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 
        \end{array}\right]
        = L\beta
\]
It then follows that
\[
        Z = X L\tp(LL\tp)^{-1} = 
        \left[\begin{array}{cccc}
        1 & 1 & 0 & 0 \\
        1 & 1 & 0 & 0 \\
        1 & 0 & 1 & 0 \\
        1 & 0 & 1 & 0 \\
        1 & 0 & 0 & 1 \\
        1 & 0 & 0 & 1 
        \end{array}\right]
        \left[\begin{array}{ccc}
        1  & 0  & 0\\
        0  & 1  & 0\\
        0  & 0  & 1\\
        0  & 0  & 0
        \end{array}\right]
        =
        \left[\begin{array}{ccc}
        1 & 1 & 0  \\
        1 & 1 & 0  \\
        1 & 0 & 1  \\
        1 & 0 & 1  \\
        1 & 0 & 0  \\
        1 & 0 & 0  
        \end{array}\right]
\]
as expected.  What do the parameters mean here?
\begin{eqnarray*}
        \mu_1 & = & \gamma_0 + \gamma_1 \\
        \mu_2 & = & \gamma_0 + \gamma_2 \\
        \mu_3 & = & \gamma_0.
\end{eqnarray*}
Hence, $\gamma_0 = \mu_3$ is the mean of the last cell, and $\gamma_1$ and
$\gamma_2$ are the differences between the third cell and the first
and second (respectively).

\paragraph{Grand mean plus treatment effects.} If we set the constraint
\[
        \beta_1 + \beta_2 + \beta_3 = 0
\]
then we can express $\beta_3 = -\beta_1 - \beta_2$ and \eq{basic-1}
becomes
\[
        y_{ij} = \beta_0 + X_{1j}\beta_1 + X_{2j}\beta_2 -
        X_{3j}(\beta_1 + \beta_2) + \epsilon_{ij}
\]
or equivalently
\begin{equation}
        y_{ij} = \beta_0 + (X_{1j}-X_{3j})\beta_1 + (X_{2j}-X_{3j})\beta_2 
        + \epsilon_{ij}
        \label{sums-1}
\end{equation}
and the matrix form becomes
\[
        \left[\begin{array}{c}
        y_{11} \\ y_{12} \\ y_{21} \\ y_{22} \\ y_{31} \\ y_{32}
        \end{array}\right]
        =
        \left[\begin{array}{rrr}
        1 &  1 &  0  \\
        1 &  1 &  0  \\
        1 &  0 &  1  \\
        1 &  0 &  1  \\
        1 & -1 & -1  \\
        1 & -1 & -1  
        \end{array}\right]
        \left[\begin{array}{c}
        \beta_0 \\ \beta_1 \\ \beta_2  
        \end{array}\right]
        +
        \left[\begin{array}{c}
        \epsilon_{11} \\ \epsilon_{12} \\ \epsilon_{21} \\ 
        \epsilon_{22} \\ \epsilon_{31} \\ \epsilon_{32}
        \end{array}\right]
\]

Now let's try setting up $L$ and $Z = XL\tp(LL\tp)^{-1}$:
\begin{description}

\item[Step 1.] This parametrization can be thought of as doing two
things: (a) we set $\beta_0=0$; and (b) we define:
\begin{eqnarray*}
        \mu & = & (\beta_1 + \beta_2 + \beta_3)/3 \\
        \alpha_1 & = & \beta_1 - \mu ~ = ~ \phantom{\,+\,} \frac23 \beta_1 -
                        \frac13\beta_2 - \frac13\beta_3 \\
        \alpha_2 & = & \beta_2 - \mu ~ = ~ -\frac13 \beta_1 +
                        \frac23\beta_2 - \frac13\beta_3 \\
        \alpha_3 & = & \beta_3 - \mu ~ = ~ -\frac13 \beta_1 -
                        \frac13\beta_2 + \frac23\beta_3 
\end{eqnarray*}
(as a check, note that $\alpha_1 + \alpha_2 + \alpha_3=0$ as desired).

\item[Step 2.] This means that $\alpha_3 = -\alpha_1 - \alpha_2$, so
we can omit it from the reparametrization, and so $L$ is the matrix in
the transformation
\[
        \left[\begin{array}{c}
        \mu \\ \alpha_1 \\ \alpha_2 
        \end{array}\right]
        = 
        \left[\begin{array}{rrrr}
        0 &  \frac13 &  \frac13 &  \frac13 \\
        0 &  \frac23 & -\frac13 & -\frac13 \\
        0 & -\frac13 &  \frac23 & -\frac13
        \end{array}\right]
        \left[\begin{array}{c}
        \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 
        \end{array}\right]
        = L\beta
\]
\item[Step 3.] Therefore
\[
        Z = X L\tp(LL\tp)^{-1} =
        \left[\begin{array}{cccc}
        1 & 1 & 0 & 0 \\
        1 & 1 & 0 & 0 \\
        1 & 0 & 1 & 0 \\
        1 & 0 & 1 & 0 \\
        1 & 0 & 0 & 1 \\
        1 & 0 & 0 & 1 
        \end{array}\right]
        \left[\begin{array}{rrr}
        0 &  0 &  0  \\
        1 &  1 &  0  \\
        1 &  0 &  1  \\
        1 & -1 & -1  
        \end{array}\right]
        =
        \left[\begin{array}{rrr}
        1 &  1 &  0  \\
        1 &  1 &  0  \\
        1 &  0 &  1  \\
        1 &  0 &  1  \\
        1 & -1 & -1  \\
        1 & -1 & -1  
        \end{array}\right]
\]
so the new model is
\[
        y = 
        \left[\begin{array}{rrr}
        1 &  1 &  0  \\
        1 &  1 &  0  \\
        1 &  0 &  1  \\
        1 &  0 &  1  \\
        1 & -1 & -1  \\
        1 & -1 & -1  
        \end{array}\right]
        \left[\begin{array}{c}
        \mu \\ \alpha_1 \\ \alpha_2  
        \end{array}\right]
        + \epsilon
\]

\item[Step 4.] Now if we define 
\[
Z_1 = X_1 - X_3 =
        \left\{\begin{array}{rr}1 & \mbox{at level $a$} \\
                              0 & \mbox{at level $b$} \\
                             -1 & \mbox{at level $c$}
                                \end{array}\right.
~~~~~~~~~~~~
Z_2 = X_2 - X_3 =
        \left\{\begin{array}{rr}0 & \mbox{at level $a$} \\
                              1 & \mbox{at level $b$} \\
                             -1 & \mbox{at level $c$} 
                                \end{array}\right.
\]
we may express the model as
\[
        y_{ij} = \mu + Z_{1j}\alpha_1 + Z_{2j}\alpha_2 
        + \epsilon_{ij}
        ~ = ~
        \mu + (X_1 - X_3)\alpha_1 + (X_2 - X_3)\alpha_2 
        + \epsilon_{ij}
\]
In this model, the cell means are
\begin{eqnarray*}
        \mu_1 & = & \mu + \alpha_1 \\
        \mu_2 & = & \mu + \alpha_2 \\
        \mu_3 & = & \mu + \alpha_3 
\end{eqnarray*}
so that the $\alpha$'s really are offsets from the grand mean $\mu$ to
each cell mean.

\end{description}

\comment{

NOTHING WRONG WITH THIS, JUST SAVING SPACE

\paragraph{SPLUS's default ``Helmert contrasts''.}  The default
parametrization for ANOVA models in SPLUS is in terms of the ``Helmert
contrasts.''  It is easiest to explain how to set this model up
directly in terms of $L$ and $Z = XL\tp(LL\tp)^{-1}$:

\begin{description}

\item[Step 1.] Like the ``grand mean plus treatment effects''
parametrization we do two things: (a) we set $\beta_0=0$; and (b) we
define:
\begin{eqnarray*}
        \mu & = & (\beta_1 + \beta_2 + \beta_3)/3 \\
        \gamma_1 & = & \{-\beta_1 + \beta_2\}/2 \\
        \gamma_2 & = & \{-(\beta_1 +\beta_2)/2 + \beta_3\}/3 \\
        \gamma_3 & = & \{-(\beta_1 + \beta_2 + \beta_3)/3 + \beta_4\}/4
\end{eqnarray*}
{\em (Here, $\beta_4\equiv0$ [since there is no fourth category] and
is just included to show you the pattern, should you ever have to
construct Helmert contrasts for more than three categories).}

\item[Step 2.] Again, $\gamma_3$ is expressible in terms of the other
parameters ($\gamma_3 = 3\mu$), so we omit it from the
reparametrization, and $L$ is the matrix in the transformation
\[
        \left[\begin{array}{c}
        \mu \\ \alpha_1 \\ \alpha_2 
        \end{array}\right]
        = 
        \left[\begin{array}{rrrr}
        0 &  \frac13 &  \frac13 &  \frac13 \\
        0 & -\frac12 &  \frac12 & 0 \\
        0 & -\frac16 & -\frac16 &  \frac13
        \end{array}\right]
        \left[\begin{array}{c}
        \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 
        \end{array}\right]
        = L\beta
\]
\item[Step 3.] Therefore
\[
        Z = X L\tp(LL\tp)^{-1} =
        \left[\begin{array}{cccc}
        1 & 1 & 0 & 0 \\
        1 & 1 & 0 & 0 \\
        1 & 0 & 1 & 0 \\
        1 & 0 & 1 & 0 \\
        1 & 0 & 0 & 1 \\
        1 & 0 & 0 & 1 
        \end{array}\right]
        \left[\begin{array}{rrr}
        0 &  0 &  0  \\
        1 & -1 & -1  \\
        1 &  1 & -1  \\
        1 &  0 &  2  
        \end{array}\right]
        =
        \left[\begin{array}{rrr}
        1 & -1 & -1  \\
        1 & -1 & -1  \\
        1 &  1 & -1  \\
        1 &  1 & -1  \\
        1 &  0 &  2  \\
        1 &  0 &  2  
        \end{array}\right]
\]
so the new model is 
\[
        y = 
        \left[\begin{array}{rrr}
        1 & -1 & -1  \\
        1 & -1 & -1  \\
        1 &  1 & -1  \\
        1 &  1 & -1  \\
        1 &  0 &  2  \\
        1 &  0 &  2  
        \end{array}\right]
        \left[\begin{array}{c}
        \mu \\ \gamma_1 \\ \gamma_2  
        \end{array}\right]
        + \epsilon
\]

\item[Step 4.] Now if we define 
\[
Z_1 = X_2 - X_1 =
        \left\{\begin{array}{rr}-1 & \mbox{at level $a$} \\
                              1 & \mbox{at level $b$} \\
                              0 & \mbox{at level $c$}
                                \end{array}\right.
~~~~~~~~~~~~
Z_2 = 2X_3 - X_2 - X_1 =
        \left\{\begin{array}{rr}-1 & \mbox{at level $a$} \\
                             -1 & \mbox{at level $b$} \\
                              2 & \mbox{at level $c$} 
                                \end{array}\right.
\]
we may express the model as
\[
        y_{ij} = \mu + Z_{1j}\gamma_1 + Z_{2j}\gamma_2 
        + \epsilon_{ij}
        ~ = ~
        \mu + (X_2 - X_1)\gamma_1 + (2X_3 - X_2 - X_1)\gamma_2 
        + \epsilon_{ij}
\]
In this model, the cell means are
\begin{eqnarray*}
        \mu_1 & = & \mu - \gamma_1 - \gamma_2 \\
        \mu_2 & = & \mu + \gamma_1 - \gamma_2 \\
        \mu_3 & = & \mu + 2\gamma_2 
\end{eqnarray*}
Note that it is easy to read off the ``treatment effects'' also:
\begin{eqnarray*}
        \alpha_1 & = & - \gamma_1 - \gamma_2 \\
        \alpha_2 & = & \phantom{\,+} \gamma_1 - \gamma_2 \\
        \alpha_3 & = & \,\,2\gamma_2 
\end{eqnarray*}


\end{description}


}


\paragraph{Coding Contrasts to Test Hypotheses.}  In this example we 
consider data from Friendly and Franklin (1980) on
the effects of presentation format on learning and memory, we examined
two hypotheses,
\begin{itemize}
\item Is mean recall in {\bf SFR} condition the same as the average
mean recall in the {\bf Meshed} and {\bf Before} conditions?
\begin{eqnarray*}
        H_0:&& \mu_{\em SFR}=(\mu_{\em Before}+\mu_{\em Meshed})/2  \\
        H_1:&& \mu_{\em SFR}\neq(\mu_{\em Before}+\mu_{\em Meshed})/2 
\end{eqnarray*}
\item Is mean recall the same under the {\bf Meshed} and {\bf Before}
conditions?
\begin{eqnarray*}
        H_0:&& \mu_{\em Before}=\mu_{\em Meshed} \\
        H_1:&& \mu_{\em Before}\neq\mu_{\em Meshed}
\end{eqnarray*}
\end{itemize}
which we coded as the contrast matrix $C$,
\begin{verbatim}
> Cmatrix
     [,1] [,2] 
[1,] -0.5    1
[2,] -0.5   -1
[3,]  1.0    0
\end{verbatim}
where the rows represent the cell means for (in alphabetical order)
{\bf Before}, {\bf Meshed} and {\bf SFR}.  The overparametrized 
cell means version of the model is
\[
        \ul y = X \left[\begin{array}{c}
                \beta_0 \\ \beta_1\\ \beta_2 \\ \beta_3 
                \end{array}\right] + \ul\epsilon
\]
where
\[
        X = \left[\begin{array}{cccc}
                1 & 1 & 0 & 0 \\
                \vdots & \vdots & \vdots & \vdots \\
                1 & 1 & 0 & 0 \\
                1 & 0 & 1 & 0 \\
                \vdots & \vdots & \vdots & \vdots \\
                1 & 0 & 1 & 0 \\
                1 & 0 & 0 & 1 \\
                \vdots & \vdots & \vdots & \vdots \\
                1 & 0 & 0 & 1 
            \end{array}\right]
            \begin{array}{@{}cc@{}}
              \left.\begin{array}{@{}c@{}} \ \\ \ \\ \ \end{array}\right\} &
              \begin{array}{@{}c@{}} \ \\
                \mbox{10 rows} \\
                \
              \end{array} \\
              \left.\begin{array}{@{}c@{}} \ \\ \ \\ \ \end{array}\right\} &
              \begin{array}{@{}c@{}} \ \\
                \mbox{10 rows} \\
                \
              \end{array} \\
              \left.\begin{array}{@{}c@{}} \ \\ \ \\ \ \end{array}\right\} &
              \begin{array}{@{}c@{}} \ \\
                \mbox{10 rows} \\
                \
              \end{array}
            \end{array}
\]
and the new parametrization is given by $L=C\tp$ (note that we have
added a column to $L$ for the intercept):
\[
        \left[\begin{array}{c}
        \gamma_0 \\
        \gamma_1 \\
        \gamma_2 
        \end{array}\right]
        =
        \left[\begin{array}{rrrr}
        1.0 & 0.0 &  0.0 &  0.0 \\
        0.0 & 0.5 &  0.5 & -1.0 \\
        0.0 & 1.0 & -1.0 &  0.0 
        \end{array}\right]
        \left[\begin{array}{c}
        \beta_0 \\ \beta_1\\ \beta_2 \\ \beta_3 
        \end{array}\right]
\]
The new model matrix $Z$ for the $\gamma_j$'s is then
\begin{eqnarray*}
        Z & = & XL\tp (LL\tp)^{-1} 
        ~ = ~ X
        \left[\begin{array}{rrr}
        1.0 &  0.0 &  0.0 \\
        0.0 &  0.5 &  1.0 \\
        0.0 &  0.5 & -1.0 \\
        0.0 & -1.0 &  0.0
        \end{array}\right]
        \left[\begin{array}{rrr}
        1.00 & 0.00 & 0.00 \\
        0.00 & 1.25 & 0.00 \\
        0.00 & 0.00 & 2.00
        \end{array}\right]^{-1} \\
        & = &  \left[\begin{array}{rrr}
                1 & 0.5 & 1 \\
                \vdots & \vdots & \vdots \\
                1 & 0.5 & 1 \\
                1 & 0.5 & -1 \\
                \vdots & \vdots & \vdots \\
                1 & 0.5 & -1 \\
                1 & -1 & 0 \\
                \vdots & \vdots & \vdots \\
                1 & -1 & 0  
            \end{array}\right]
        \left[\begin{array}{rrr}
        1.00 & 0.00 & 0.00 \\
        0.00 & 0.80 & 0.00 \\
        0.00 & 0.00 & 0.50
        \end{array}\right]
\end{eqnarray*}
Comparing this with our earlier notes where we fit the contrasts ``by
hand'', you may notice that in the earlier notes I only formed the
product $Z=XL\tp$ to get the new regression matrix for the contrasts,
not $Z=XL\tp(LL\tp)^{-1}$.  Since these contrasts were orthogonal,
$(LL\tp)^{-1}$ is a diagonal matrix (as can be seen above) and failing
to multiply by it affects the {\em scale} of the raw $\gamma_j$
parameters but not the magnitude or distributions of the
$t$-statistics
\[
        t = \frac{\hat\gamma_j}{SE_{\hat\gamma_j}}
\]
(or the $F$ statistics we typically use) since the missing scale
factor will cancel from both numerator and denominator of the
expressions defining them.


\section{Gauss-Markov and MLE}

This is about two distinct ``optimality'' results in the theory of
linear regression.  The first, the {\em Gauss-Markov theorem}, says
that the least-squares estimates are best linear unbiased estimates
(BLUE) of the regression coefficients.  The second says that the
least-squares estimates are also the maximum-likelihood estimates
(MLE's).

\subsection{The Gauss-Markov Theorem}

In the multiple regression model
\[
        \ul y_{n\by 1} = X_{n\by k'} \ul \beta_{k'\by 1} + \ul
        \epsilon_{n\by 1} 
\]
we will only assume that the error terms $\epsilon_i$ are independent,
mean 0, constant variance $\sigma^2_\epsilon$.  

Let us define $M = (X\tp X)^{-1}X\tp$, so that $\hat{\ul\beta} = M\ul y$
are the least-squares estimators of $\ul\beta$, and let
$\tilde{\ul\beta} = U\ul y$ be any unbiased linear estimate of
$\beta$.   Finally, let $A=U-M$.  

\begin{lemma}
        For any unbiased linear estimate $U\ul y$ of $\ul\beta$,
        $A=U-M$ is orthogonal to the rows of $X$; hence $AM\tp = MA\tp = 0$.
\end{lemma}

{\em Proof:}  Since both $\tilde{\ul\beta}$ and $\hat{\ul\beta}$ are
unbiased, 
\[
        \ul\beta = E[\tilde{\ul\beta}] = E[(M+A)\ul y] 
        = E[M\ul y] + E[M\ul y] = E[\hat{\ul\beta}] + AE[\ul y]
        = \ul\beta + AX\ul\beta ,
\]
and since this must be true for all $\ul\beta$ it follows that
$AX=0$. It immediately follows that $AM\tp = AX(X\tp X)^{-1} = 0$. $\Box$

\begin{theorem} {\em Gauss-Markov Theorem.}

        $\hat{\ul\beta}$ is the Best (minimum-variance) Linear Unbiased
        Estimate (BLUE) of $\ul\beta$.
\end{theorem}

{\em Proof:} We will do this by minimizing the variance of the
arbitrary unbiased linear estimate $\tilde{\ul\beta} = U\ul y =
(M+A)\ul y$.  Since $AM\tp = MA\tp = 0$, it follows that
\begin{eqnarray*}
        \Var(\tilde{\ul\beta}) & = & \Var[(M+A)\ul y] \\
        & = & (M+A) \Var(\ul y) (M+A)\tp ~ = ~ (M+A) \sigma^2_\epsilon
        I (M+A)\tp \\
        & = & \sigma^2_\epsilon (MM\tp + MA\tp + AM\tp + AA\tp) \\
        & = & \sigma^2_\epsilon (MM\tp + AA\tp)
\end{eqnarray*}
and hence, reading off the diagonal elements of this
variance-covariance matrix,
\[
        \Var(\tilde\beta_j) = \sigma^2_\epsilon \left(
        \sum_{i=1}^n m_{ij}^2 + \sum_{i=1}^n a_{ij}^2 
        \right) .
\]
Since the $m_{ij}$'s are fixed, we choose different unbiased linear
estimates $\tilde{\ul\beta}$ by choosing different $a_{ij}$'s; from
the above expression it is clear that the only way to minimize the
variance of all the $\tilde\beta_j$'s is to take $a_{ij} \equiv 0$;
so that $\tilde\beta_j=\hat\beta_j$ for every $j$. $\Box$

{\em Remark:} Essentially the same technique can be used to prove the
following result:
\begin{corollary}
        Let $\phi=L\ul\beta$ be a linear function of the $\beta$'s.
        The BLUE of $\phi$ is $L\hat{\ul\beta}$ (i.e. plug in the
        least-squares estimates).
\end{corollary}


\clearpage



\subsection{MLE Properties of Least-Squares}

\subsubsection{The Multivariate Normal Density}

The univariate normal density ($n=1$) we all remember is 
\[
        \phi(y) = \frac{1}{\sqrt{2\pi}\sigma} \exp \left[-\frac12 
        \left( \frac{y-\mu}{\sigma}\right)^2 \right]
\]
for $y \in \Re$.  We often write $Y \sim N(\mu,\sigma^2)$ to convey
that the random variable $Y$ has this density.

The multivariate normal density, for $\ul y \in \Re^n$, is
\begin{equation}
        \phi(\ul y) = \frac{1}{(2\pi)^{n/2}|\vc|^{1/2}}
        \exp [ -\frac12 (\ul y - \ul \mu)\tp \vc^{-1} (\ul y - \ul \mu) ]
\label{mvn}
\end{equation}
The mean vector $\ul \mu$ replaces the univariate mean $\mu$ and the
variance-covariance matrix $\vc$ replaces the univariate variance
$\sigma^2$. We often write $\ul Y \sim N_n(\ul\mu,\vc)$, or just $\ul
Y \sim N(\ul\mu,\vc)$ if the dimension $n$ is clear from context.
Sometimes we also say that the components of $\ul Y$ are ``jointly
normal'' or ``multivariate normal''.

One way to build up the multivariate normal is to start with a
$n$-dimensional vector of independent identically distributed $N(0,1)$
variables $\ul Z$ $=$ $(Z_1, \ldots, Z_n)$.  Since densities for
independent random variables multiply, the density for $\ul Z$ must be
\begin{eqnarray*}
        \phi(\ul z) & = & \prod _1^n \phi(z_i)  \\
        & = & \frac{1}{(2\pi)^{n/2}} 
        \exp [ - \frac12 \sum_1^n z_i^2 ] \\
        & = & \frac{1}{(2\pi)^{n/2}} 
        \exp [ - \frac12 \ul z\tp I^{-1} \ul z ]
\end{eqnarray*}
where $I$ is the $n\by n$ identity matrix.  

Then if $\ul Y = A\ul Z + \ul \mu$, it is clear that $E[\ul Y]$ $=$
$\ul \mu$ and $\Var(\ul Y)$ $=$ $E[(\ul Y - \ul \mu)(\ul Y - \ul
\mu)\tp]$ $=$ $AA\tp$ $=$ $\vc$ (so for example it is immediately
obvious that $\vc$ must by symmtric and positive definite).  

The multivariate change-of-variables formula from calculus will ensure
that the density for $\ul Y$ looks like \eq{mvn}.  In fact, it is
generally true that if $\ul Y = A\ul Y + \ul b$, where $A$ is $m\by n$
and $\ul b \in \Re^m$, and $\ul Y \sim N_n(\ul \mu, \vc)$, then $\ul Y
\sim N_m(A\ul \mu + \ul b, A\vc A\tp)$.

\comment{

NOTHING WRONG WITH THIS, JUST SAVING SPACE

As a more explicit example, in the case $k=2$ we often write
\[
        \ul \mu = \left[ \begin{array}{c}\mu_1 \\ \mu_2 \end{array} \right] ;
        \vc = \left[ \begin{array}{cc} \sigma_1^2 & \rho\sigma_1\sigma_2 \\
                                    \rho\sigma_1\sigma_2 & \sigma_2^2
                \end{array} \right]
\]
in which case you should be able to figure out from \eq{mvn} that
\small
\begin{eqnarray*}
        \lefteqn{\phi(y_1, y_2)} \\
         & = & \frac{1}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}
        \exp \left\{ -\frac1{2(1-\rho^2)} \left[ 
        \left(\frac{y_1 - \mu_1}{\sigma_1}\right)^2
        - 2\rho \left(\frac{y_1 - \mu_1}{\sigma_1}\right)\left(\frac{y_2 -
                \mu_2}{\sigma_2}\right)
        + \left(\frac{y_2 - \mu_2}{\sigma_2}\right) \right] \right\} \\
        & = & 
        \frac{1}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}
        \exp \left\{ -\frac1{2(1-\rho^2)} \left[ z_1^2 - 2z_1z_2 + z_2^2 
                \right] \right\}
\end{eqnarray*}
\normalsize
where $z_i = (y_i - \mu_i)/\sigma_i$; what is the matrix $A$ in the
transformation from $\ul Z$ to $\ul Y$ mentioned above?

}

\subsubsection{Maximum Likelihood Estimation}

Now let us derive maximum likelihood estimates for $\ul\beta$ and
$\sigma^2_\epsilon$ in the multiple regression model
\[
        \ul y_{n\by 1} = X_{n\by k'} \ul \beta_{k'\by 1} + \ul
        \epsilon_{n\by 1}    ~~~~ \ul\epsilon \sim
        N_n(\ul 0, \sigma^2_\epsilon I_{n\by n}) .
\]
Since this implies the $y_i$'s are independent
\[
        \ul y \sim N_n(X\beta, \sigma^2_\epsilon I_{n\by n})
\]
we can immediately write down the likelihood for $\ul\beta$ and
$\sigma^2_\epsilon$ by writing down the product of the individual
normals\footnote{{\em Warning:} $\ul x_i\tp$ is the $i^{th}$ \ul{\em row}
of the regression matrix $X$.},
\begin{eqnarray*}
        L(\ul y; \ul\beta,\sigma^2_\epsilon) 
	& = & \prod _1^n \frac{1}{\sqrt{2\pi\sigma^2_\epsilon}}
		\exp \left[ - \frac12  (y_i-\ul x_i\tp\ul\beta)^2 / 
		\sigma^2_\epsilon \right] \\
        & = & \frac{1}{(2\pi\sigma^2_\epsilon)^{n/2}} 
		\exp \left[ - \frac12 \sum_1^n (y_i-\ul x_i\tp\ul\beta)^2 /
		\sigma^2_\epsilon \right] \\ 
        & = & \frac{1}{(2\pi)^{n/2}|\sigma^2 I|^{1/2}} 
                \exp \left[ - \frac12 (\ul y - X\ul \beta)\tp 
		(\sigma^2_\epsilon I)^{-1} (\ul y - X\ul\beta) \right] ,
\end{eqnarray*}
and hence the log-likelihood is
\[
        \ell(\ul\beta,\sigma^2_\epsilon) 
        = \log L((\ul y; \ul\beta,\sigma^2_\epsilon)
        = - \frac n2 \log 2\pi - \frac n2 \log \sigma^2_\epsilon 
                - \frac{1}{2\sigma^2_\epsilon} (\ul y - X\ul \beta)\tp
                (\ul y - X\ul\beta)
\]
Now we differentiate with respect to $\ul\beta$ and
$\sigma^2_\epsilon$ to obtain
\begin{eqnarray*}
        \frac{\partial\ell}{\partial\ul\beta} & = &
        - \frac{1}{2\sigma^2_\epsilon}[-2X\tp(y - X\ul\beta)] \\
        \frac{\partial\ell}{\partial\sigma^2_\epsilon} & = &
        - \frac n2 \frac{1}{\sigma^2_\epsilon} 
        + \frac{1}{\sigma^4_\epsilon} (\ul y - X\ul \beta)\tp
                (\ul y - X\ul\beta)
\end{eqnarray*}
Setting the first set of partial derivatives equal to zero produces the
normal equations again, and hence
\[
        \hat{\ul\beta} = (X\tp X)^{-1} X\tp \ul y ,
\]
the usual least-squares regression estimates.  Setting the second set
of partial derivatives equal to zero produces
\[
        \hat\sigma^2_\epsilon = \frac{ (\ul y - X\hat{\ul\beta})\tp
                (\ul y - X\hat{\ul\beta})}{n} = \frac{RSS}{n}
\]
which we now know is a biased estimate of $\sigma^2_\epsilon$ (the
correct denominator should be $n-k'$).





\subsubsection{Some Geometry; Eigensystems; Square Roots; Principal Components}

A useful geometric property of the multivariate normal distribution is
that it is ``log-quadratic''; $\log \phi(\ul y)$ is a simple linear
function of the symmetric positive-definite quadratic form $(\ul y -
\ul \mu)\tp \vc^{-1} (\ul y - \ul \mu)$ in $(\ul y - \ul \mu)$.  This
means that the level surfaces $\{\ul y: \phi(\ul y) = K\}$, or
equivalently the level surfaces $\{\ul y: \log \phi(\ul y)= c\}$ will
be ellipsoids.

Finding the principal axes of the ellipsoidal level sets (contours of
$\phi(\ul y)$) is easy with Lagrange multipliers.  To find the first
principal axis for example, you need a point $\ul y$ that is a maximum
distance from $\ul \mu$ [maximize the squared distance $(\ul y - \ul
\mu)\tp(\ul y - \ul \mu)$], but is still on the level set [subject to
the constraint $(\ul y - \ul \mu)\tp\vc^{-1}(\ul y - \ul \mu) = c$].
Differentiating the Lagrange-multiplier objective function
\begin{equation}        
        f(\ul y) = (\ul y - \ul \mu)\tp(\ul y - \ul \mu) - 
        \lambda [(\ul y - \ul \mu)\tp\vc^{-1}(\ul y - \ul \mu) - c]
\label{pc-lagrange}
\end{equation}
with respect to $\ul y$ and setting these derivatives equal to zero
leads to the eigenvalue/eigenvector problem
\[
        (\vc - \lambda I) (\ul y - \ul \mu) = 0
\]
So there must be multiple critical points for \eq{pc-lagrange}, one
for one corresponding to each eigenvalue of $\sigma$.  It now a matter
of calculation to see that the eigenvalues will satisfy
\[
        \lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_p > 0
\]
with corresponding eigenvectors $\ul e_i$ $=$ $(\ul y_i - \ul \mu)$
lying along the principal axes of the ellipsoid $(\ul y - \ul
\mu)\tp\vc^{-1}(\ul y - \ul \mu) = c$; the half-lengths of the principal
axes will be $\sqrt{\lambda c}$.

Several interesting algebraic and geometric properties come out of
this analysis.  To every eigenvalue $\lambda_i$ there corresponds an
eigenvector $\ul e_i$, and it can be shown that all the $\ul e_i$ are
linearly independent.  If two eigenvectors $\ul e_i$ and $\ul e_j$
correspond to different eigenvalues $\lambda_i \neq \lambda_j$ then,
since $\ul e_j\tp \vc \ul e_i$ is a scalar (hence equal to its transpose),
\[
        \lambda_i \ul e_j\tp\ul e_i = \ul e_j\tp \vc \ul e_i = [\ul e_j\tp
        \vc \ul e_i]\tp = \ul e_i\tp \vc \ul e_j = \lambda_j \ul e_i\tp\ul e_j
\]
so it must be that $\ul e_i\tp\ul e_i = 0$: the eigenvectors from
different eigenvalues are orthogonal.  If $\lambda_i$ $=$ $\lambda_j$,
there will still be two eigenvectors, and using the Gram-Schmidt
procedure we can select two eigenvectors that are orthogonal to each
other and to all the other eigenvectors in the problem.  If the two
eigenvalues are equal, then the cross-section of the ellipse in the
plane containing the two linearly independent eigenvectors $\ul e_i$ and
$\ul y_j$ will be circular (it would always be an ellipse, and now the
major and minor principal axes will be of the same length
$2\sqrt{\lambda_i c}$); if three are equal the corresponding
cross-section will be spherical, etc.

We can assume all the $e_i$ are orthogonal; so let $\ul a_i$ = $\ul
e_i / \| \ul e_i \|$ and consider the matrix
\[
        A = \left[\begin{array}{cccc}
                \ul a_1 & \ul a_2 & \cdots & \ul a_p
        \end{array}\right]
\]
You should verify that $AA\tp = I$ (the matrix is orthogonal---it has
orthonormal columns); that $\vc A = A \, \diag(\lambda_1, \ldots,
\lambda_p)$ (the square $p\by p$ matrix with $\lambda$'s on the main
diagonal and zeros everywhere else); and therefore that $A\tp\vc A =
\diag(\lambda_1, \ldots, \lambda_p)$.  

Now consider the random vector $\ul W = A (\ul Y - \ul \mu)$.  This vector
$\ul W$ will be multivariate normal, and you should verify that it has
mean $\ul 0$ (the zero vector) and variance-covariance matrix
$\diag(\lambda_1, \ldots, \lambda_p)$.  Hence each $W_i$ is
uncorrelated with $W_j$ (in fact, independent); and $\Var(W_i) =
\lambda_i$.  The $W_i$'s are called the {\em principal components\/}
of $\ul Y$ and they play a large role in understanding the structure
of $\ul Y$.




\end{document}





