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

\parskip=0.5em
\parindent=0pt

\frenchspacing

\def\tp{^T}
% \def\tp{'}

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

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

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

\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}

\begin{document}

\begin{center}\Large\bf 
        Quadratic Forms and Geometry of Least Squares
\end{center}

\section{Basic Geometry}

\begin{figure}[h]
\unitlength=0.5in
\begin{picture}(12,7)(-0.2,2.5)
\input page1.tex
\end{picture}
\caption{Geometry of Linear Regression}
\label{fig0}
\end{figure}

We have seen that fitting the linear regression model
\[
        \ul{y} = X \ul{\beta} + \ul{\varepsilon},~~~
\ul{\varepsilon} \sim N (\ul{0}, \sigma^{2} I_{n \times n}),
\]
where
\[
X = \left[
\begin{array}{c c c c c}
| & | & | && |\\
\ul{1} &\ul{X}_{1} & \ul{X}_{2} & \cdots &
\ul{X}_{k}\\
| & | & | & & |
\end{array}
 \right], \mbox{ \ \ } \ul{\beta}\tp = (\beta_{0}, \beta_{1},
\ldots, \beta_{k}) ,
\]
by least squares was equivalent to decomposing the vector $\ul y\in \Re^n$
into orthogonal components
\[
        \ul y = P\ul y + Q\ul y = \hat{\ul y} + \ul E
\]
where
\begin{eqnarray*}
        \hat{\ul y} & = & P\ul y ~ \in ~ ColSpace(P) ~ = ~ ColSpace(X) \\
        \ul E & = & Q\ul y ~ \perp ColSpace(P)
\end{eqnarray*}
and
\[
\begin{array}{rclcl}
        P & = & X (X\tp X)^{-1} X\tp \\
        Q & = & I_{n \times n} - P \\
        \hat{\ul{\beta}} & = & (X\tp X)^{-1} X\tp \ul{y} \\
        \hat{\ul{y}} & = & X \hat{\ul{\beta}} & = & P \ul{y}\\
        \ul{E} & = & \ul{y} - \hat{\ul{y}} & = & Q \ul{y} 
\end{array}
\]
since
\[
\begin{array}{rclr}
        P^2 & = & P;~~ P\tp ~=~ P & {\em (idempotent, symmetric)} \\
        Q^2 & = & Q;~~ Q\tp ~=~ Q & {\em (idempotent, symmetric)} \\
        PQ  & = & QP ~ = ~ \mb 0 & {\em (orthogonal)} \\
        P+Q & = & I              & {\em (complete)} \\
\mbox{and } PX & = & X & {\em (projection~onto~}ColSpace(X){\em )}      
\end{array}     
\]
This gives rise to the picture in Figure~\ref{fig0}.  It is not hard
to show for example that the Pythagorean Theorem holds,
\[
        \Vert \ul y \Vert^2 = \ul y\tp \ul y = [(P+Q)\ul y]\tp[(P+Q)\ul y]
        = (\ul y \tp P + \ul y\tp Q )(Py + Qy) 
        = \ul y\tp P \ul y + \ul y\tp Q \ul y
        = \Vert \hat{\ul y}\Vert^2 + \Vert \ul E \Vert^2
\]
and it turns out, according to a theorem called ``Cochran's Theorem''
which we shall see below, that $\ul y\tp P \ul y$ and $\ul y\tp Q \ul
y$ are statistically independent of one another.  Much of the theory
of $F$-statistic testing in linear models is built on elaborations of
this Pythagorean Theorem, and on the corresponding statistical
indepedence facts about the component sums of squares.

\section{Some Results on Quadratic Forms}

Let \ul{$y$} $\in \Re^{n}$, and suppose $A$ is an $n
\times n$ matrix.


\begin{definition}
$\ul{y}\tp A \ul{y}$ is a {\em quadratic form} if
\[
\begin{array}{rll}
1. &  A\tp = A  & \mbox{(symmetric)}\\
2. & \ul{y}\tp A \ul{y} \geq 0, ~  \forall 
\ul{y} \in \Re^{n} & \mbox{(positive semidefinite)}
\end{array}
\] 
In this case, $A$ is the {\em defining matrix} of $\ul{y}\tp
A \ul{y}$.
\end{definition}

\begin{definition}
$A$ is a {\em projection matrix} if $A^2=A$ (we say $A$ is {\em idempotent}).
\end{definition}
Note that if $A$ is a projection matrix then $A(I_{n\times n}-A)=\ul
0$, and every $\ul y \in \Re^n$ has a unique decomposition 
\[
        \ul y = A\ul y + (I-A)\ul y = \ul y_A + \ul y_{A^{\perp}}
\]
into a vector $\ul y_A \in ColSpace(A)$ and a vector $\ul
y_{A^{\perp}} \perp ColSpace(A)$.  If $A$ is also symmetric, then
\[
        \mbox{trace$(A)$ $=$ rank$(A)$ $=$ dim[$ColSpace(A)$]} .
\]
Note that, therefore, trace$(A)$ is the ``number of dimensions'' that
$A\ul y$ can wander around in, and $\ul y\tp A \ul y = (A\ul y)\tp
(A\ul y) = \Vert A\ul y\Vert^2$ is the squared length of $A\ul y$.

\clearpage

{\bf Facts:}
\begin{enumerate}
\item
If $\ul{y}$ is a random vector with $E [\ul{y}] =
\ul{\mu}, \Var(\ul{y}) = V, \mbox{ and if }
\ul{y}\tp A \ul{y}$ is a quadratic form, then
\[
E [\ul{y}\tp A \ul{y}] = \mbox{ trace } A V +
\ul{\mu}\tp A \ul{\mu}
\]
\item 
If $y_{i} \sim N (\mu_{i}, \sigma^{2})$, independent, and $A$ is
symmetric and idempotent, then $\ul{y}\tp A \ul{y} / \sigma^{2}$ has a
{\em noncentral} $\chi^{2}$ {\em distribution on} $k$ {\em degrees of
freedom}, where $k = \mbox{ rank } A$. The quantity $\ul{\mu}\tp A
\ul{\mu} / 2 \sigma^{2}$ is called the {\em noncentrality parameter}
of the distribution.  Note: here, $V = \Var\ul{y} = \sigma^{2} I.$
\end{enumerate}


\begin{theorem}
{\em (Cochran's Theorem, Scheff\'e, 1959, pp 419--421)}

Let $y_{i} \sim N (\mu_{i}, 1), i = 1 \ldots n$, independent, and
suppose 
\[
\sum^{n}_{i = 1} y^{2}_{i} = \ul{y}\tp A_{1} \ul{y} +
\ul{y}\tp A_{2} \ul{y} + \ldots +
\ul{y}\tp A_{s} \ul{y}
\]
\noindent
where each $\ul{y}\tp A_{j} \ul{y}$ is a {\em quadratic
form}, with rank $A_{j} = n_{j}$.
\ul{Then:}
\begin{itemize}
\item[(a)] 
Each $A^{2}_{j} = A_{j}$ {\em (\ul{idempotent}, so rank $A_{j}
= \mbox{ trace } A_{j}$)}; and $A_{i} A_{j} = {\bf 0}$ for all $i\neq
j$ ($A_{i}$ and $A_{j}$ are {\em orthogonal}).
\item[(b)] 
The quadratic forms $\ul{y}\tp A_{j} \ul{y}$
are distributed as independent, noncentral $\chi^{2}$'s on $n_{j}$
degrees of freedom.
\end{itemize}
\noindent
\ul{if and only if} $n_{1} + n_{2} + \ldots + n_{s} = n$.
\end{theorem}

\begin{corollary}
If $y_{i} \sim N (\mu_{i}, \sigma^{2})$, the theorem still holds but
now $\ul{y}\tp A_{j} \ul{y} / \sigma^{2}$ will be
independent noncentral $\chi^{2}_{n_{j}}$'s.
\end{corollary}
 
Cochran's Theorem or its corollary tell us immediately, for example,
that since
\[
        {\sum}_i y_i^2 = \Vert \ul y \Vert^2 
        = \ul y\tp P \ul y + \ul y\tp Q \ul y
\]
and since trace$(P)$ $+$ trace$(Q)$ $=$ trace$(P+Q)$ $=$ trace$(I)$ $=$
$n$, then $\ul y\tp P \ul y$ and $\ul y\tp Q \ul y$ are (proportional
to) independent $\chi^2$'s on trace$(P)$ and trace$(Q)$ degrees of freedom.



\clearpage





\section{Overall ANOVA for regression with an intercept}

Let's continue now with the model
\begin{equation}
        \ul{y} = X \ul{\beta} + \ul{\varepsilon},~~~
        \ul{\varepsilon} \sim N (\ul{0}, \sigma^{2} I_{n \times n}),
        \label{whole}
\end{equation}
where
\[
X = \left[
\begin{array}{c c c c c}
| & | & | && |\\
\ul{1} &\ul{X}_{1} & \ul{X}_{2} & \cdots &
\ul{X}_{k}\\
| & | & | & & |
\end{array}
 \right], \mbox{ \ \ } \ul{\beta}\tp = (\beta_{0}, \beta_{1},
\ldots, \beta_{k}) , \mbox{ \ \  } k' = k + 1 ,
\]
with projection matrix $P$ and $Q=I-P$.  

We wish to compare this model with the simpler model 
\begin{equation}
        \ul y = \ul 1 \beta_0 + \ul\varepsilon
        \label{intercept}
\end{equation}
of regression onto just the intercept $\beta_0$.  You can verify that
the projection matrix for this regression is $J/n$ where $J$ is the
$n\times n$ matrix of all 1's.   It is easy to see that
\[
\begin{array}{rcccccc}
        \sum_{i=1}^n y_i^2 & = & \Vert \ul y \Vert^2 
        & = & \ul y\tp \,I\, \ul y  \\
%       & = & \ul y\tp (I - J/n) \ul y & + & \ul y\tp\, J/n\, \ul y \\
        & = & \ul y\tp (I - P) \ul y & + & \ul y\tp (P - J/n) \ul y
                & + & \ul y\tp\, J/n\, \ul y \\
        & = & \Vert Q\ul y \Vert^2 & + & \Vert (P-J/n)\ul y\Vert^2
                & + & \Vert (J/n) \ul y \Vert^2 \\
        & = & RSS & + & RegSS & + & n(\bar y)^2
\end{array}
\]
and we can readily verify\footnote{Use the facts that $P$ and $J/n$
are the projection matrices for the regressions \eq{whole} and
\eq{intercept} and the fact that $\ul 1$ is a column of $X$, so that
$(I-P)(J/n) = \mb 0$.} that the matrices $A_1 = (I-P)$, $A_2 =
  (P-J/n)$, and $A_3 = J/n$,
% \begin{eqnarray*}
%   A_1 & = & (I-P) \\
%   A_2 & = & (P-J/n) \\
%   A_3 & = & J/n
% \end{eqnarray*}
are symmetric ($A_i\tp = A$), idempotent ($A_i^2=A_i$) and orthogonal
($A_iA_j=\mb0$), and trace$(A_1)$ $+$ trace$(A_2)$ $+$ trace$(A_3)$ $=$ $n$.
So, by {\em Cochran's Theorem} and its corollary, the quadratic forms
\[
\ul{y}\tp (J/n) \ul{y} / \sigma^{2}, \mbox{ \ \ \ } \ul{y}\tp
(P - J/n) \ul{y} / \sigma^{2}, \mbox{ \ \ \ } \ul{y}\tp Q
\ul{y} / \sigma^{2}
\]
are independent noncentral $\chi^{2}$'s.
Moreover you can easily show that
\[
TSS= \sum^{n}_{i = 1} (y_{i} - \bar{y})^{2} = \sum y^{2}_{i} - n
(\bar{y})^{2} = \ul{y}\tp (I - J/n) \ul{y}
\]
so the RegSS and RSS terms are in fact independent $\chi^2$'s in the
overall ANOVA table
\begin{center}
\begin{tabular}{l|c|l}
Source & df & 
\multicolumn{1}{c}{SS}\\ \hline
Regression & $k' - 1$ & $RegSS = \ul{y}\tp (P - J/n)
\ul{y}$\\ 
Residual & $n - k'$ & $RSS = \ul{y}\tp Q \ul{y}$ \\
\hline
Total    & $n - 1$ & TSS = $\ul{y}\tp (I - J /n) \ul{y}$\\
\end{tabular}
\end{center}

The geometry for this ANOVA table is given in Figure~\ref{fig1}.

\begin{figure}[h]
\unitlength=0.5in
\begin{picture}(12,7)(-0.2,2.5)
\input page9.tex
\end{picture}
\caption{Geometry of the overall ANOVA table in linear regression}
\label{fig1}
\end{figure}

The overall $F$-statistic for the regression is
\[
        F = \frac{(TSS- RSS)([(n - 1) - (n - k')])}{RSS/ (n - k')} 
          = \frac{RegSS / (k' - 1)}{RSS/ (n - k')}
          = \frac{RegSS / k }{RSS/ (n - (k+1))} ,
\]
and tests $H_0:$ eq. \eq{intercept} holds, vs. $H_1:$ equation
\eq{whole} holds.  We note that
\begin{itemize}\parsep=0pt\parskip=0pt\itemsep=0pt
\item $ RegSS / \sigma^{2}  = \ul{y}\tp (P - J/n) \ul{y}
/ \sigma^{2} \ \ \sim \mbox{ noncentral } \chi^{2}_{k-1}:$
\begin{eqnarray*}
        df & = & \mbox{ rank } (P - J/n) = \mbox{ trace } \ (P - J/n) =
                \mbox{ trace } P - \mbox{ trace } (J/n) \\
        & = & \mbox{ trace } \left[ X (X\tp X)^{-1} X\tp \right] - 1 = 
                \mbox{ trace } (X\tp X (X\tp X)^{-1}){-1}\\
%       & & \mbox{(since trace$AB = $ trace$BA$ when $AB$ and $BA$ are
%               both defined).} \\
        & = & \mbox { trace } I_{k' \times k'} - 1 ~ = ~ k' - 1 ~ = ~ k
                \\[0.5em]
        E \left[ \ul{y}\tp (P - J/n) \ul{y} \right] & = & \mbox{ trace }
                \left[ (P - J/n) \sigma^{2} I \right] + (X \ul{\beta})\tp (P -
                J/n) (X \ul{\beta}) \\
        & = & \sigma^{2} (k' - 1) + \ul{\beta}\tp X\tp (P - J/n) X
                \ul{\beta}\\
        & = & \sigma^{2} (k' - 1) + \ul{\beta}\tp X\tp (I - J/n) X
                \ul{\beta}\\
         & & \mbox{since } PX = X.
\end{eqnarray*}
\item $RSS / \sigma^{2}  = \ul{y}\tp Q \ul{y} \ \
\sim \mbox{ noncentral } \chi^{2}_{n - k'}:$
\begin{eqnarray*}
        df & = & \mbox{ rank } Q = \mbox{ trace } Q = \mbox{ trace }
                (I - P) = \mbox{ trace } I - \mbox{ trace } P\\
        & = & n - k' ~ = ~ n - (k+1) 
                \\[0.5em]
        E \left[ \ul{y}\tp Q \ul{y} \right] & = & \mbox{ trace } \left[ Q
                (\sigma^{2} I) \right] + (X \ul{\beta})\tp Q (X \ul{\beta})\\
        & = & \sigma^{2} (n - k') + \ul{\beta}\tp X\tp Q X
                \ul{\beta} \\
        & = & \sigma^{2} (n - k') ~ = ~ \sigma^2 (n-(k+1))
\end{eqnarray*}
\end{itemize}

\noindent
(Note: $\ul{\beta}\tp X\tp Q X \ul{\beta} = 0 \mbox{
since } Q X = 0.$)

After some algebra we see that the ratio of the expected values in the
$F$ statistic is
\[
\frac{E \left[ RegSS / (k' - 1) \right]}{E \left[ RSS/ (n
- k') \right]} = 1 + \left[ \ul{\beta}\tp X\tp (I -
J/n) X \ul{\beta} \right] / (k' - 1) \sigma^{2} \geq 1
\]
The thing that prevents this from being exactly 1 is $\ul{\beta}\tp
X\tp (I - J/n) X \ul{\beta}$.  If you multiply out $ X\tp (I - J/n) X$
you will see\footnote{{\em Hint}: Apply the formula $\sum^{n}_{i = 1}
x^{2}_{i} - n (\bar{x})^{2} = \ul{x}\tp (I - J/n) \ul{x}$ to
situations where $\ul{x}\tp$ and $\ul{x}$ represent a column of $X$,
and recall that the first column of $X$ is all 1's.}  that
\[
\begin{array}{lll}
\ul{\beta}\tp X\tp (I - J/n) X \ul{\beta}
 & = &
(\beta_{0}, \beta_{1} \cdots \beta_{k})
\left[
\begin{array}{cccc}
0 & 0 & \ldots & 0\\
0 & \mbox{ stuff } & \ldots & \mbox{ stuff } \\
\vdots & \vdots & & \vdots \\
0 & \mbox{ stuff } & \ldots & \mbox{ stuff }
\end{array}
\right]
\left[
\begin{array}{c}
\beta_{0}\\
\beta_{1}\\
\vdots\\
\beta_{k}
\end{array}
\right]\\
& = &(\beta_{1} \cdots \beta_{k})
\left[
\begin{array}{ccc}
\mbox{ stuff } & \ldots & \mbox{ stuff }\\
\vdots & & \vdots\\
\mbox{ stuff } & \ldots & \mbox{ stuff }
\end{array}
\right]
\left[ 
\begin{array}{c}
\beta_{1}\\
\vdots\\
\beta_{k}
\end{array}
\right]
\end{array}
\]
\noindent
which will equal zero when $H_{0}: \beta_{1} = \beta_{2} = \cdots =
\beta_{k} = 0$ is true.


{\bf Distribution of $F$ = $\frac{RegSS / (k' - 1)}
{RSS/ (n - k')}$ under $H_{0}$:}

According to Cochran's theorem, $F$ is a ratio of independent {\em
noncentral} $\chi^{2}$'s divided by their respective degrees of freedom.
If they were {\em central} $\chi^{2}$'s (noncentrality parameter $= 0$), 
$F$ would follow the $F$ distribution we always use.

But they {\em are} central $\chi^{2}$'s:
\begin{itemize}
\item
Noncentrality param. of $RSS$ is $\frac{1}{2 \sigma^{2}} (X
\ul{\beta})\tp Q (X \ul{\beta})
= \ul{\beta}\tp X\tp Q X \ul{\beta} = 0$ since $Q X =
{\bf 0}$
\item
Noncentrality param. of $RegSS$ is $\frac{1}{2 \sigma^{2}}
\ul{\beta}\tp X\tp (I - J/n) X \ul{\beta} = 0$ under
$H_{0}$ (See above!).
\end{itemize}

\noindent
{\bf Distribution of $F$ under $H_{1}$}: 

\noindent
The numerator $\chi^{2}, RegSS / \sigma^{2}$, now has noncentrality
parameter
\[
\frac{1}{2 \sigma^{2}} \ \ul{\beta}\tp X\tp (I - J/n) X
\ul{\beta} \geq 0
\]
\noindent
The larger this parameter, the easier $H_{1}$ is to detect.






\section{Using sums of squares to test nested models}


\noindent
{\em Model:} $\ul{y} = X \ul{\beta} +
\ul{\varepsilon} \mbox{ }, \ul{\varepsilon} \sim N
(\ul{0}, \sigma^{2} I)$
\noindent
with some small changes of notation: Let $k' = k_{1} + k_{2}$, and
\[
\begin{array}{ll}
X = \left[
\begin{array}{cccc}
| & | & & | \\
\ul{X}_{1} & \ul{X}_{2} & \cdots & \ul{X}_{k'}\\
| & | & & |
\end{array}
\right] &
 = \left[
\underbrace{\begin{array}{cccc}
| & | & & | \\
\ul{X}_{1} & \ul{X}_{2} & \cdots &
\ul{X}_{k_{1}}\\
| & | & & |
\end{array}}
\underbrace{\begin{array}{cccc}
| & | & & |\\
\ul{X}_{k_{1} + 1} & \ul{X}_{k_{1} + 2} & \cdots &
\ul{X}_{k_{1} + k_{2}}\\
| & | & & |
\end{array}}
\right]\\
& = \left[
\hspace*{.75in} 
X_{1} \hspace*{1.5in} X_{2}
\hspace*{.9in} 
\right]\\
& \\
\noindent
\mbox{and}\\
& \\
\ul{\beta\tp} = (\beta_{1} \cdots \beta_{k'}) & = \left(
\underbrace{\begin{array}{cccc}
\beta_{1}, & \beta_{2}, & \cdots, & \beta_{k_{1}}, 
\end{array}}
\underbrace{\begin{array}{cccc}
\beta_{k_1 + 1},& \beta_{k_{1} + 2}, & \cdots, & \beta_{k_{1}
+k_{2}}
\end{array}} \right)\\
& = ( \hspace*{.7in} \ul{\beta}\tp_{1}, \hspace*{.6in} 
\hspace*{.9in}
\ul{\beta}\tp_{2} \hspace*{.9in} )
\end{array}
\]
\noindent
so that 
\[
\ul{y} = X \ul{\beta} + \ul{\varepsilon} = X_{1}
\ul{\beta}_{1} + X_{2} \ul{\beta}_{2} +
\ul{\varepsilon}
\]
\noindent
(We could treat ``regression with an intercept'' again by making the
column $\ul{X}_{1}$ just a column of 1's.)    
 
\noindent
We want to compare the two models

\begin{center}
\begin{tabular}{| l | l |} \hline
SMALL: $\ul{y} = X_{1} \ul{\beta}_{1} +
\ul{\varepsilon}$ & $\mbox{ Hat matrix } P_{1} = X_{1} (X\tp_{1}
X_{1})^{-1} X\tp_{1}$\\
LARGE: $\ul{y} = X_{1} \ul{\beta}_{1} + X_{2}
\ul{\beta}_{2} + \varepsilon$ & $\mbox{ Hat matrix } P_{12} = X (X\tp
X)^{-1} X\tp$\\ \hline
\end{tabular}
\end{center}

\[
\begin{array}{rlp{0.5in}l}
RSS(\ul \beta_1,\ul\beta_2) = RSS_12 & =  \sum^{n}_{i = 1} (y_{i} - \hat{y}_{i (12)})^{2}\\ 
& = ( \ul{y} - \ul{\hat{y}}_{(12)})\tp
(\ul{y} - \ul{\hat{y}}_{(12)}) \\
& = [ (I - P_{12}) \ul{y} ]\tp 
[(I - P_{12}) \ul{y}]\\
& =  \ul{y}\tp (I - P_{12}) \  \ul{y} = \ul{y}\tp
Q_{12} \ul{y} \  & & (Q_{12} = I - P_{12})\\
\mbox{ Similarly} &  \\ 
RSS(\ul \beta_1) = RSS_{1} & =  \sum^{n}_{i = 1} (y_{i} - y_{i (1)})^{2}
= \ul{y}\tp Q_{1} \ul{y} & & (Q_{1} = I - P_{1})
\end{array}
\]

\begin{figure}
\unitlength=0.5in
\begin{picture}(12,7)(-0.2,2.5)
\input page11.tex
\end{picture}
\caption{Geometry of the $F$-test comparing the nested models
$X=[X_1,\,X_2]$ and $X_1$}
\label{fig2}
\end{figure}

\noindent
We want a test for $H_{0}: \ul{\beta}_{2} = \ul{0}$ vs.
$H_{1}: \ul{\beta_{2}} \mbox{ free}$.  [See Figure~\ref{fig2}.]
This is a test for
\[
\begin{array}{rl}
H_{0}: \ul{y} & = X_{1} \ul{\beta}_{1} +
\ul{\varepsilon}\\
\mbox{ vs. } \ \ H_{1}: \ul{y} & = X_{1} \ul{\beta}_{1} + X_{2}
\ul{\beta}_{2} + \ul{\varepsilon}
\end{array}
\]
\noindent
Clearly we can write, thinking of the smaller model,
\[
\begin{array}{rl}
\sum y^{2}_{i} & = \sum \hat{y}^{2}_{i (1)} + \left( \sum y^{2}_{i} -
\sum \hat{y}^{2}_{i (1)} \right)\\
& = \ul{y}\tp P_{1} \ul{y} + \ul{y}\tp I
\ul{y} - \ul{y}\tp P_{1} \ul{y}\\
& = \ul{y}\tp P_{1} \ul{y} + \ul{y}\tp (I -
P_{1} ) \ul{y}\\
& = \ul{y}\tp P_{1} \ul{y} + \ul{y}\tp Q_{1}
\ul{y}\\
& = \ul{y}\tp P_{1} \ul{y} + \ul{y}\tp ( Q_{1} - Q_{12}) \ul{y}
        + \ul{y}\tp Q_{12} \ul{y} \\
& = \ul{y}\tp P_{1} \ul{y} + (RSS_1 - RSS_{12}) + RSS_{12} \\
& = \ul{y}\tp P_{1} \ul{y} + SS( \ul{\beta}_{2} | \ul{\beta}_{1}) + RSS_{12} 
\end{array}
\]
where the third-to-last inequality follows from $Q_1 =  (Q_1 - Q_{12}) +
Q_{12}$ and the fact that (check!) $Q_{12}(Q_1-Q_{12})=0$.

Now, if $SS(\ul{\beta}_{2} | \ul{\beta}_{1}) =
\ul{y}\tp (Q_{1} - Q_{12}) \ul{y}$:
\begin{enumerate}
\item[(1)]
is a quadratic form (positive-semidefinite)?
\item[(2)]
has noncentrality parameter $ = 0$ under $H_{0}$?
\end{enumerate}

\noindent
then according to Cochran's theorem, the $F$-statistic
\[
F = \frac{SS(\ul{\beta}_{2} | \beta_{1}) / \mbox{ trace } [Q_{1} -
Q_{12}]}
{RSS(\ul{\beta}_{1}, \ul{\beta}_{2}) / \mbox{ trace } Q_{12}}
\]
\noindent
can be used to test $H_{0}: \ul{\beta}_{2} = 0$ vs. $H_{1}:
(H_{0})^{c}$.
\noindent
In order to verify (1) and (2) for $SS(\ul{\beta}_{2}
| \ul{\beta}_{1}) = \ul{y}\tp (Q_{1} - Q_{12})
\ul{y}$ we will {\em reparametrize} the {\em larger} models to
make $Q_{12}$ easier to calculate.  The idea is similar to the
idea that $\hat{\ul\beta}_2$ can be obtained by regressing the residuals
$Q_1\ul y$ onto the residuals $Q_1X_2$ that we saw empirically earlier
in the course.
% idea behind {\em added variable plots}.

\noindent
Define
\[
\begin{array}{rl}
& ~~~~~~~~~~~~~~~~~~~~~~~ P_1 \\
X_{2 \cdot 1} ~ = ~ Q_{1} X_{2} & 
\rlap{$
\hphantom{~~= ~ X_{2} - }\overbrace{\vphantom{X_1^{-1}}
\hphantom{X_{1}\left( X\tp_{1} X_{1}\right)^{-1} X\tp_{1}}} 
$}
= ~ X_{2} -  X_{1}
\underbrace{ \left( X\tp_{1} X_{1} \right)^{-1} X\tp_{1} X_{2}}
\\
& ~~~~~~~~~~~~~~~~~~~~~~~~~~ A 
% \\
% & ~ = ~ X_{2} - P_{1} X_{2} ~ = ~ X_{2} - X_{1} A  \\
% & \\
% &  
\end{array}
\]
\renewcommand{\baselinestretch}{1.5} % changes spacing.
\renewcommand{\arraystretch}{1.5} \tiny\normalsize

\noindent
Clearly $X_{2 \cdot 1}$ is the {\em residuals} of (each column of) $X_{2}$
after regressing on $X_{1}$.
\begin{itemize}
\item 
$X_{2 \cdot 1}$ is the part of $X_{2}$ unexplained by regression on
$X_{1}$.
\item
$X_{2 \cdot 1}$ is orthogonal to $X_{1}:
X\tp_{1} X_{2 \cdot 1} = {\bf 0}$.
\end{itemize}

\noindent
With $X_{2 \cdot 1}$ and $A = (X\tp_{1} X_{1})^{-1} X\tp_{1} X_{2}$
we can write 
\[
\begin{array}{rl}
\ul{y} & = X_{1} \ul{\beta}_{1} + X_{2}
\ul{\beta}_{2} + \ul{\varepsilon} 
        ~ = ~ X_{1} \ul{\beta}_{1} + ( X_{1} A + X_{2} - X_{1} A)
\ul{\beta}_{2} + \ul{\varepsilon} \\
& = X_{1} ( \underbrace{ \ul{\beta}_{1} + A
\ul{\beta}_{2}})~~ + ~~( \underbrace{X_{2} - X_{1} A})
\ul{\beta}_{2} ~~ + ~~ \ul{\varepsilon}\\
& = ~~~ X_{1} ~~~ \ul{\alpha}_{1} ~~~~~~~ +
~~~~~~~ X_{2 \cdot 1} ~~~
\ul{\beta}_{2} ~~~~ + ~~ \ul{\varepsilon}
\end{array}
\]
\noindent
So the larger model is {\em identical} to 
\[
\ul{y} = X_{1} \ul{\alpha}_{1} + X_{2 \cdot 1}
\ul{\beta}_{2} + \ul{\varepsilon}
\]

\noindent
with model matrix $\tilde{X} = \left[ X_{1} \  X_{2 \cdot 1} \right]$
instead of $X = \left[ X_{1} \ X_{2} \right]$.
\noindent
Therefore
\[
\begin{array}{rl}
\ul{y}\tp Q_{12} \ul{y} & = RSS\left(
\ul{\beta}_{1}, \ul{\beta}_{2} \right) = RSS\left(
\ul{\alpha}_{1}, \ul{\beta}_{2} \right) =
\ul{y}\tp 
\left[ I - \tilde{X} \left( \tilde{X}\tp \tilde{X} \right)^{-1}
\tilde{X}\tp \right] \ul{y}\\
& = \ul{y}\tp \left\{
I - \left[ X_{1} \ X_{2 \cdot 1} \right]
\left[
\begin{array}{ll}
X\tp_{1} X_{1} & X\tp_{1} X_{2 \cdot 1}\\
X\tp_{2 \cdot 1} X_{1} & X\tp_{2 \cdot 1} X_{2 \cdot 1}
\end{array}
\right]^{-1}
\left[
\begin{array}{c}
X\tp_{1} \\
X\tp_{2 \cdot 1}
\end{array}
\right]
\right\}\ul{y}\\
& = \ul{y}\tp \left\{ I - \left[ X_{1} \  X_{2 \cdot 1} \right] 
\left[
\begin{array}{cc}
\left( X\tp_{1} X_{1} \right)^{-1} & {\bf 0} \\
{\bf 0 } & \left( X\tp_{2 \cdot 1} X_{2 \cdot 1} \right)^{-1}
\end{array}
\right]
\left[
\begin{array}{c}
X\tp_{1}\\
X\tp_{2 \cdot 1}
\end{array}
\right]
\right\} \ul{y}\\
& = \ul{y}\tp \left(I - P_{1} - P_{2 \cdot 1} \right) \ul{y}
\end{array}
\]
\noindent
Hence
\[
\begin{array}{rl}
\ul{y}\tp Q_{12} \ul{y} & = \ul{y}\tp \left(I -
P_{1} - P_{2 \cdot 1} \right) \ul{y}\\
& = \ul{y}\tp Q_{1} \ul{y} - \ul{y}\tp P_{2
\cdot 1} \ul{y}
\end{array}
\]
\noindent
Therefore
\[
SS\left( \ul{\beta}_{2} | \ul{\beta}_{1} \right) =
\ul{y}\tp \left(Q_{1} - Q_{12} \right) \ul{y} =
\ul{y}\tp P_{2 \cdot 1} \ul{y}
\]
\noindent
is a quadratic form with defining matrix $P_{2 \cdot 1} = X_{2 \cdot
1} \left(X\tp_{2 \cdot 1} X_{2 \cdot 1}\right)^{-1} X\tp_{2 \cdot
1}$. (It is positive definite because
\[
\ul{y}\tp_{1} P_{2 \cdot 1} \ul{y} = \sum \hat{y}^{2}_{i (2
\cdot 1)} \geq 0
\]
\noindent
from the model $\ul{y} = X_{2 \cdot 1}
\tilde{\ul{\beta}} + \ul{\varepsilon}$).

\noindent
Now, since $P_{2 \cdot 1}$ projects onto $X_{2 \cdot 1}$ which is {\em
orthogonal} to $X_{1}$, it follows that $P_{2 \cdot 1} X_{1} = {\bf
0}$. Therefore,
\[
\begin{array}{ll}
& E [SS( \ul{\beta}_{2} | \ul{\beta}_{1}
) ] 
= E \left[ \ul{y}\tp P_{2 \cdot 1}
\ul{y} \right] \\
= & \mbox{ trace } \left( P_{2 \cdot 1} \sigma^{2} I \right) + \left( \left[ X_{1}
\ X_{2} \right] \ul{\beta} \right)\tp P_{2 \cdot 1} \left(
\left[ X_{1} \  X_{2} \right] \ul{\beta} \right) \\
= & \sigma^{2} \mbox{ trace } \left[ X\tp_{2 \cdot 1} X_{2 \cdot 1} \left( X\tp_{2
\cdot 1} X_{2 \cdot 1} \right)^{-1} \right] + \left[
\ul{\beta}\tp_{1} \ \ul{\beta}\tp_{2} \right]
\left[
\begin{array}{c}
X\tp_{1}\\
X\tp_{2}
\end{array}
\right] P_{2 \cdot 1} \left[ X_{1} \  X_{2} \right]
\left[
\begin{array}{c}
\ul{\beta}_{1}\\
\ul{\beta}_{2}
\end{array}
\right]\\
= & \sigma^{2} k_{2} + \ul{\beta}\tp_{2} X\tp_{2} P_{2 \cdot
1} X_{2} \ul{\beta}_{2}
\end{array}
\]
This last term controls the noncentrality parameter, and equals 0 under $H_{0}:
\ul{\beta}_{2} = \ul{0}$.

Once again, $E \left[ RSS\left( \ul{\beta}_{1}, 
\ul{\beta}_{2} \right) \right] = E \left[ \ul{y}\tp
Q_{12} \ul{y} \right] = \sigma^{2} \left(n - k_{1} - k_{2}
\right)$, (verify this!) so
\[
\frac{E \left[ SS\left( \ul{\beta}_{2}| \ul{\beta}_{1}
\right) / k_{2} \right]}
{E \left[ RSS\left( \ul{\beta}_{1},  \ul{\beta}_{2} 
\right) /  \left( n - k_{1} - k_{2} \right) \right]} = 1 +
\ul{\beta}\tp_{2} X\tp_{2} P_{2 \cdot 1} X_{2}
\ul{\beta}_{2}
\]
equals 1 when $H_{0}: \ul{\beta}_{2} = \ul{0}$ is {\em
true}; and $\geq 1$ when $H_{0}$ is {\em false}.  This is suggestive
of the behavior of the $F$ statistic
\[
        F = \frac{ SS\left( \ul{\beta}_{2}| \ul{\beta}_{1}
        \right) / k_{2} }{ RSS\left( \ul{\beta}_{1},\ul{\beta}_{2} 
        \right) /  \left( n - k_{1} - k_{2} \right) }
\]


\section{Testing a General Linear Hypothesis}

Let us return to the general linear regression model
\[
        \ul y_{n\times1} = X_{n\times k'} \ul \beta_{k'\times 1} +
        \varepsilon_{n\times 1}
\]
and now consider inference for a set of general linear combinations of
the parameters
\[
        \ul\psi_{q\times1} = L_{q\times k'} \beta_{k'\times 1}
\]
Some examples of $L$'s that might be of interest include
\begin{itemize}
\item $L$ selects the last $k_2$ parameters (as in the previous section):
\[
        L_{k_2\times(k_1+k_2)} = 
        \underbrace{
        \left[\begin{array}{ccc}
        0 & \cdots & 0 \\
        0 & \cdots & 0 \\
        \vdots & \ddots & \vdots \\
        0 & \cdots & 0 
        \end{array}\right.
        }_{k_1~\rm columns}
        \underbrace{
        \left.\begin{array}{cccc}
        1 & 0 & \cdots & 0 \\
        0 & 1 & \cdots & 0 \\
        \vdots & \vdots & \ddots & \vdots \\
        0 & 0 & \cdots & 1 
        \end{array}\right]
        }_{k_2~\rm columns}
        ;~~~~~
        L\left[\begin{array}{c}\beta_1 \\ \vdots \\ \beta_{k_1 + k_2} 
        \end{array}\right] ~ = ~ \left[\begin{array}{c}\beta_{k_1+1} 
        \\ \vdots \\ \beta_{k_1 + k_2} \end{array}\right]
\]
\item $L$ codes the contrasts to test the hypotheses:
\[
        \begin{array}{@{H_0:~}l@{~vs.~ H_1:~}l}
        \mu_1=(\mu_2 + \mu_3)/2 & \mu_1\neq(\mu_2 + \mu_3)/2 \\
        \mu_2=\mu_3 & \mu_2 \neq \mu_3
        \end{array}
\]
in a one-way ANOVA model with three cells, i.e.
\[
        L_{2\times3} = \left[\begin{array}{rrr}
        1 & -1/2 & -1/2 \\
        0 & 1 & -1
        \end{array}\right]
        ; ~~~~~
        L\left[\begin{array}{c}\mu_1 \\ \mu_2 \\ \mu_3 \end{array}\right] 
        ~ = ~ \left[\begin{array}{c} \mu_1 - (\mu_2+\mu_3)/2 
        \\ \mu_2 - \mu_3 \end{array}\right]
\]
(note that $L$ is the transpose of the ``contrast matrix'' that we
have dealt with previously).
\end{itemize}

\renewcommand{\baselinestretch}{1} % changes spacing.
\renewcommand{\arraystretch}{1} \tiny\normalsize

Now, since $\ul\psi_{q\times1} = L_{q\times k'}\ul\beta_{k'\times1}$,
it follows that
\[
        \hat{\ul\psi} = L\hat{\ul\beta} = L(X\tp X)^{-1}X\tp y = A y
        \sim N(\psi,\not\Sigma_{\hat\psi})
\]
(either by least-squares, or more simply from a maximum likelihood
argument), where
\[
        \not\Sigma_{\hat\psi} = \Var Ay = A\Var y A\tp = AA\tp
        \sigma^2_\epsilon 
\]
Now let
\[
        B = \not\Sigma_{\hat\psi}/\sigma^2_\epsilon = AA\tp
        = L(X\tp X)^{-1}X\tp X (X\tp X)^{-1}L\tp = L(X\tp X)^{-1}L\tp
\]
It then follows (see for example Scheff\'e pp. 1--30) that 
\[
        (\hat{\ul\psi} - \ul\psi)\tp B^{-1} (\hat{\ul\psi} - \ul\psi) 
        \sim \sigma^2_\epsilon \chi^2_q, ~~~{\rm independent~of} ~RSS
\]
{\em Intuition:}
\begin{description}
\item[Chi-squared:] Let $R$ be a ``square root'' of $B^{-1}$; then
$R(\hat{\ul\psi} - \ul\psi)$ uncorrelates $(\hat{\ul\psi} - \ul\psi)$ and
the result is a vector of $q$ independent $N(0,1)'s$; the quadratic 
form $(\hat{\ul\psi} - \ul\psi)\tp B^{-1} (\hat{\ul\psi} - \ul\psi) =
\Vert R(\hat{\ul\psi} - \ul\psi)\Vert^2$ is just the sum of squares of 
these $q$ $N(0,1)$'s; 
\item[Indep. of $\mathbf RSS$:] Like our earlier arguments based on Cochran's
Theorem.
\end{description}
Thus, when the true value of $\ul\psi=\ul\psi_0$,
\[
        F = \frac{(\hat{\ul\psi} - \ul\psi_0)\tp B^{-1} 
        (\hat{\ul\psi} - \ul\psi_0)/q }{RSS/(n-k')}
        = \frac{(\hat{\ul\psi} - \ul\psi_0)\tp [L(X\tp X)^{-1}L\tp]^{-1}
        (\hat{\ul\psi} - \ul\psi_0)/q }{S^2_E}
        \sim F_{q,n-k'}
\]
% {\bf Warning: \em In Chapter 9 of the text, Fox always writes $[L(X\tp
% X)^{-1}L\tp]$ when he means $[L(X\tp X)^{-1}L\tp]^{-1}$.  There are
% errors on pp. 226 (three times) and 239 (once).}

This can immediately be ``inverted'' to provide direct tests of
hypotheses: 
\begin{itemize}
\item Under $H_0:~\ul\psi=\ul\psi_0$, the $F$ statistic above is
$F_{q,n-k'}$; 
\item Under $H_1:\ul\psi$ free, the $F$ statistic tends to be too
large.
\end{itemize}
One can show that when $L$ selects $k_2$ parameters and we test for
$\ul\psi=\ul 0$, we get exactly the same test that we derived
in connection with least-squares geometry in the last section.  In
this case $[L(X\tp X)^{-1}L\tp]^{-1}$ has the effect of ``selecting
out'' the rows and columns of $(X\tp X)^{-1}$ corresponding to the
parameters $\beta_{k_1+1}, \ldots, \beta_{k_1+k_2}$ and inverting the
resulting submatrix.  

\section{Confidence Regions}

The same basic fact, that 
\[
        F = \frac{(\hat{\ul\psi} - \ul\psi_0)\tp [L(X\tp X)^{-1}L\tp]^{-1}
        (\hat{\ul\psi} - \ul\psi_0)/q }{S^2_E}
        \sim F_{q,n-k'}
\]
when $\ul\psi=\ul\psi_0$, can be used to generate confidence regions
for parameters or combinations of parameters.  The idea of a
confidence region is the same as that of a confidence interval, except
that it is a region in $q$-space for all $q$ parameters $\ul\psi
=(\psi_1, \ldots, \psi_q)$.

\begin{definition}
A 100$\times(1-\alpha)$\% confidence ellipsoid for the parameters
$\ul\psi =(\psi_1, \ldots, \psi_q)$ is given by 
\[
        \{\ul\psi:~(\hat{\ul\psi} - \ul\psi)\tp [L(X\tp X)^{-1}L\tp]^{-1}
        (\hat{\ul\psi} - \ul\psi)\leq q S^2_E F_{\alpha,q,n-k'}\}
\]
\end{definition}
Once again, $L$ may be designed so that it selects a particular group
of parameters $\beta_j$ of interest, or a particular group of linear
combinations of $\beta_j$'s, etc.

Some examples will help to fix ideas here.  The {\tt prestige} data
% set was introduced in Fox Chapter 2, and 
contains various variables
that were thought to be predictive of occupational prestige in Canada,
taken largely from 1971 Canadian Census archives.  The variables are
fairly self-explanatory
{\small
\begin{verbatim}
> names(prestige)
[1] "education"    "income"       "pct.female"   "occ.prestige" "occ.code"    
[6] "occ.type"    
> dim(prestige)
[1] 102   6
\end{verbatim}
}
\noindent
Here we review a couple of the fits that we did when we first
encountered this data set.
{\small
\begin{verbatim}
> summary(fit.ed <- lm(occ.prestige ~ education, data = prestige))

Call: lm(formula = occ.prestige ~ education, data = prestige)
Residuals:
    Min     1Q Median    3Q   Max 
 -26.04 -6.523 0.6611 6.743 18.16

Coefficients:
               Value Std. Error  t value Pr(>|t|) 
(Intercept) -10.7320   3.6771    -2.9186   0.0043
  education   5.3609   0.3320    16.1478   0.0000

Residual standard error: 9.103 on 100 degrees of freedom
Multiple R-Squared: 0.7228 
F-statistic: 260.8 on 1 and 100 degrees of freedom, the p-value is 0 

Correlation of Coefficients:
          (Intercept) 
education -0.9695    
\end{verbatim}
}
\noindent
We see that neither the slope nor the intercept are likely to be zero,
according to their own $t$-tests; we could also see this with 95\%
$t$-intervals:
{\small
\begin{verbatim}
> qt((1-.05/2),100)
[1] 1.983972
> -10.732+c(-1,1)*1.983972*3.6771    # for beta_0
[1] -18.027263  -3.436737
> 5.3609  +c(-1,1)*1.983972*0.332    # for beta_1
[1] 4.702221 6.019579
\end{verbatim}
}
\noindent
on the other hand, the parameter estimates are very highly
correlated, and these intervals are not corrected for doing two
inferences at once.  

We can fix the multiple inference problem with the {\em Bonferroni}
correction, replacing $\alpha$ with $\alpha$/(number of tests),
{\small
\begin{verbatim}
> -10.732+c(-1,1)*qt((1-(.05/2)/2),100)*3.6771    # for beta_0
[1] -19.099801  -2.364199
> 5.3609  +c(-1,1)*qt((1-(.05/2)/2),100)*0.332    # for beta_1
[1] 4.605383 6.116417
\end{verbatim}
}
\noindent
but we still don't learn how the parameters covary.  Instead we
calculate the confidence ellipse
{\small
\begin{verbatim}
> X _ model.matrix(fit.ed)
> dim(X)
[1] 102   2
> X[1:5,]
                    (Intercept) education 
 GOV_ADMINISTRATORS           1     13.11
   GENERAL_MANAGERS           1     12.26
        ACCOUNTANTS           1     12.77
PURCHASING_OFFICERS           1     11.42
           CHEMISTS           1     14.62
> SE2 _ sum(residuals(fit.ed)^2)/100
> sqrt(SE2) # check against earlier output
[1] 9.103291 
> M _ t(X)%*%X/(2*SE2)
> b0 _ seq(-20,-2,length=50)
> b1 _ seq(3,8,length=50)
> bhat0 _ coefficients(fit.ed)[1]
> bhat1 _ coefficients(fit.ed)[2]
> qform _ function(psi, psi0, M)
+ {
+ #
+ # qform
+ #
+ # computes the quadratic form (psi - psi0)' %*% M %*% (Psi - Psi0)
+ # on every row of the matrix psi, using the single row in psi0
+ # 
+ #
+       qfn <- function(x, M)
+       {
+               x %*% M %*% x
+       }
+       n <- dim(psi)[1]
+       delta <- psi - matrix(rep(psi0, n), byrow = T, ncol = 2)
+       result <- apply(delta, 1, qfn, M)
+       return(result)
+ }
> qform(cbind(0,0),c(bhat0,bhat1),M)
[1] 1480.219
> qlev _ qform(expand.grid(b0,b1),c(bhat0,bhat1),M) # wait a while...
> contour(b0,b1,matrix(qlev,nrow=length(b0),ncol=length(b1)),
+ levels=qf(.95,2,100))
> abline(h=c(4.702221, 6.019579),v=c(-18.027263, -3.436737))
\end{verbatim}
}
\noindent
The resulting confidence ellipse is shown in Figure~\ref{contour1}.
Even though the individual confidence intervals are fairly wide, the
confidence ellipse shows that $\hat\beta_1$ can vary only a little,
for example, at each fixed value of $\hat\beta_0$.  Horizontal and
vertical lines indicating the original uncorrected $t$-based
confidence intervals for $\beta_0$ and $\beta_1$ are also overlaid on
the plot; one can see that the ``box'' formed by these two intervals
both misses some of the area covered by the ellipse---reflecting the
fact that the intervals are uncorrected for multiple inferences---and
covers much more area than the ellipse covers---reflecting the
inefficiency of using individual intervals.

\begin{figure}
\hfil\psfig{figure=contour1.ps,width=5.5in,height=3.5in,angle=0}
\caption{95\% Confidence Ellipse for $\beta_0$ and $\beta_1$ from {\tt
fit.ed}}\label{contour1}
\end{figure}

If one did want individual intervals for the $\beta$'s that were
guaranteed to cover the $\beta$'s with 95\% confidence, a natural
choice would be the ``shadow'' of the ellipse on each axis.  This
generates a larger box (not shown in Figure~\ref{contour1}); the two
intervals would be guaranteed to contain $\beta_0$ and $\beta_1$ in at
least (and in fact more than) 95\% of repeated samples from the
``model'' that generated the data.

{\bf Scheff\'e's method.}  An extension of the ``shadow'' idea is
Scheff\'e's method of generating simultaneous confidence intervals for
{\em all possible} linear combinations of the original parameters
$\ul\beta$.

Let $\psi = \sum_{j=1}^{k'} w_j\beta_j = \ul w\tp \ul\beta$ be any
linear combination of the original $k'$ model parameters; then
as usual
\[
        \hat\psi = \ul w\tp \hat{\ul\beta} = \ul w\tp (X\tp X)^{-1}X\tp \ul y
\]
and
\[
        \Var(\hat\psi) = \ul w\tp \Var(\hat{\ul\beta}) \ul w
        = \sigma^2_\epsilon \ul w\tp (X\tp X)^{-1} \ul w
\]
and hence
\[
        SE_{\hat\psi} = S_E \sqrt{\ul w\tp (X\tp X)^{-1} \ul w}
\]
A simple version\footnote{If you know that you will only be
considering $\psi$'s that are linear combinations of a subset of $q$
of the $\beta$'s then you may replace $(X\tp X)^{-1}$ with $L(X\tp
X)^{-1}L\tp$ where $L$ selects the subset, and use the cutoff
$qF(1-\alpha, q, n-k')$ instead.} and of Scheff\'e's method replaces
the $t$ cutoff in an uncorrected confidence interval
\[
        \hat\psi \pm t_{1-\alpha/2,n-k'} \cdot SE_{\hat\psi}
\]
with a scaled $F$ cutoff:
\[
        \hat\psi \pm \sqrt{k'F(1-\alpha, k', n-k')}
        \cdot SE_{\hat\psi}
\]
The $F$ cutoff $k'F(1-\alpha, k', n-k')$ comes from calculating a
confidence ellipsoid for the original set of $k'$ parameters
$\ul\beta$; the essential idea of Scheff\'e's method is that once you
have controlled these, you automatically control any linear
combination of them.  For a proof that these are really simultaneous
100$(1-\alpha)$\% intervals, see Scheff\'e, 1959, section 3.5.

Scheff\'e's intervals are very useful for 
\begin{itemize}
\item Data snooping, or any situation in which you are evaluating a
lot of linear contrasts or other linear combinations of the original
parameters, and you want to control the lengths of your intervals; and
\item Graphing simultanous confidence envelopes around regression
lines or surfaces.
\end{itemize}
We will illustrate the latter with the {\tt prestige} example.  We
wish to graph 
\begin{itemize}
\item The data
\item The fitted regression line $\hat y = \hat\beta_0 + x\hat\beta_1$
\item A simultaneous 95\% confidence envelope computed at every fitted value
\end{itemize}

\begin{figure}
\hfil\psfig{figure=envelope1.ps,width=5.5in,height=3.5in,angle=0}
\caption{95\% Envelope for Fitted Regression Line from from {\tt
fit.ed}}\label{envelope1}
\end{figure}

We recall that $\hat{\ul y} = X\hat{\ul\beta} = P\ul y = H\ul y$ and 
$\Var(y_i) = h_{ii} \sigma^2_\epsilon$, so that $SE_{\hat y_i} =
S_E\sqrt{h_{ii}}$.  Also $n=102$ and $k'=2$ for the simple regression
of {\tt occ.prestige} on {\tt education}.  
{\small
\begin{verbatim}
> attach(prestige)
> plot(education,occ.prestige)
> abline(fit.ed)
> SEy _ sqrt(lm.hat(fit.ed)*sum(residuals(fit.ed)^2)/100)
> l _ fitted(fit.ed) - sqrt(2*qf(.95,2,100))*SEy
> h _ fitted(fit.ed) + sqrt(2*qf(.95,2,100))*SEy
> ord _ order(education)
> lines(education[ord],l[ord],lty=3)
> lines(education[ord],h[ord],lty=3)
> detach()
\end{verbatim}
}
\noindent
The resulting graph is displayed in Figure~\ref{envelope1}.
Comparing the $t$ cutoff with the $F$ cutoff in Scheff\'e's method, 
{\small
\begin{verbatim}
> sqrt(2*qf(.95,2,100))
[1] 2.484873
> qt(1-.05/2,100)
[1] 1.983972
\end{verbatim}
}
\noindent
we see that there is a substantial difference between the corrected
and the uncorrected envelopes for the regression line here.

{\bf Remark:} There are other methods to generate confidence envelopes
that may do better.  For example, at each $x_i$ value, we could
calculate and plot the values
\begin{eqnarray*}
        \tilde u_i  & = & \max\{\beta_0 + x_i\beta_1:~(\beta_0,\beta_1)\in
        \mbox{ the 95\% confidence ellipsoid in Fig~\ref{contour1}}\} \\
        \tilde\ell_i & = &\min\{\beta_0 + x_i\beta_1:~(\beta_0,\beta_1)\in
        \mbox{ the 95\% confidence ellipsoid in Fig~\ref{contour1}}\} 
\end{eqnarray*}
along with the $\hat y_i$'s.

\end{document}



