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

\parskip=0.5em
\parindent=0pt

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

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

\def\bzero{\beta_0}
\def\bone{\beta_1}
\def\btwo{\beta_2}
\def\bhat{\hat\beta}

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

\def\inner#1#2{#1^T #2}

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

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

%\slides

\begin{document}

\begin{center}\Large\bf
A More Careful Development of the Linear Regression Model
\end{center}

\section{Estimating the Coefficients in the Multiple Regression Model}

The {\em multiple regression\/} model that we have been working with
may be written as
\[
        y = \beta_0 + \beta_1 x_1 + \beta_2 x_2
                + \cdots + \beta_k x_k + \epsilon
\]
where $\epsilon\sim N(0,\sigma^2)$.  As usual, this stands for
\begin{eqnarray*}
        y_1 & = & \beta_0 + \beta_1 x_{11} + \beta_2 x_{12}
                + \cdots + \beta_k x_{1k}  + \epsilon_1 \\
        y_2 & = & \beta_0 + \beta_1 x_{21} + \beta_2 x_{22}
                + \cdots + \beta_k x_{2k}  + \epsilon_2 \\
         & \vdots & \\
        y_n & = & \beta_0 + \beta_1 x_{n1} + \beta_2 x_{n2}
                + \cdots + \beta_k x_{nk}  + \epsilon_n
\end{eqnarray*}
or, in matrix notation,
\[
        \vc{y} = X \vc{\beta} +
        \vc{\epsilon}, \mbox{ \ \  } \vc{\epsilon} \sim N
        (\vc{0}, \sigma^{2} I_{n \times n})
\]
i.e.
\[
        \left[\begin{array}{l}
        y_1 \\ y_2 \\ \vdots \\ y_n
        \end{array}\right]
        =
        \left[\begin{array}{ccccc}
        1 &  x_{11} &  x_{12} & \cdots &  x_{1k}   \\
        1 &  x_{21} &  x_{22} & \cdots &  x_{2k}   \\
          &  \vdots &  \vdots & \ddots &  \vdots   \\
        1 &  x_{n1} &  x_{n2} & \cdots &  x_{nk}  
        \end{array}\right]
        \left[\begin{array}{l}
        \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_k
        \end{array}\right]
        +
        \left[\begin{array}{l}
        \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n
        \end{array}\right]
\]
We know that
\begin{eqnarray*}
        \mbox{Mean: } E[y_i] & = & E[\beta_0 +
                \beta_1 x_{i1} + \ldots \beta_k x_{ik} + \epsilon_i] \\
                             & = & \beta_0 +
                \beta_1 x_{i1} + \ldots \beta_k x_{ik} + E[\epsilon_i] \\
                             & = & \beta_0 +
                \beta_1 x_{i1} + \ldots \beta_k x_{ik} + 0 \\
                             & = & \beta_0 +
                \beta_1 x_{i1} + \ldots \beta_k x_{ik} , \mbox{ and}\\
                             \\
 \mbox{Variance: } \Var(y_i) & = & \Var(\beta_0 +
                \beta_1 x_{i1} + \ldots \beta_k x_{ik} + \epsilon_i) \\
                             & = & \Var(\epsilon_i) \\
                             & = & \sigma^2 .
\end{eqnarray*}
which we write compactly in matrix notation as
\begin{eqnarray*}
        E[\vc{y}] & = & X\vc{\beta} \\
        \Var(\vc{y}) & = & \sigma^2 I_{n\times n} \mbox{, and so} \\
        \vc{Y} & \sim & N (X\vc{\beta},~\sigma^2 I_{n\times n})
\end{eqnarray*}
If we have some guesses
\[
        \hat\vc{\beta} = (\bhat_0, \bhat_1, \ldots, \bhat_k)^T
\]
for $\vc{\beta}$, then the fitted values are 
\[
        \hat\vc{y} = X \hat\vc{\beta} 
\]
and the residuals are
\[
        \vc{E} = \vc{y} - \hat\vc{y}
\]
To find good estimates of $\vc{\beta}$, we may minimize the sum of
squares of residuals,
\begin{eqnarray*}
      Q(\bhat_0,\bhat_1,\ldots,\bhat_k) & = & \sum_{i=1}^n E_i^2 
        ~ = ~ \sum_{i=1}^n (y_i - [\hat\beta_0 +
        \hat\beta_1 x_{i1} + \ldots \hat\beta_k x_{ik}])^2 \\
        & = & \sum_{i=1}^n (y_i - \hat y_i)^2 \\
        & = & (\vc{y} - \hat\vc{y})^T(\vc{y} - \hat\vc{y}) \\
        & = & (\vc{y} - X\hat\vc{\beta})^T(\vc{y} - X\hat\vc{\beta}) 
\end{eqnarray*}

One way to proceed is to take derivatives ($\frac{\partial}{\partial
\vc{x}} A\vc{x} = \vc{x}$, and $\frac{\partial}{\partial \vc{x}} \vc{x}^T
A\vc{x} = 2A\vc{x}$) and equate to
zero:
\begin{eqnarray*}
\lefteqn{  \frac{\partial}{\partial \hat\vc{\beta}} \left [
    {(\vc{y} - X \hat{\vc{\beta}})}^T {(\vc{y} - X \hat{\vc{\beta}})}
   \right ] } \\
 & = &  \frac{\partial}{\partial \hat\vc{\beta}}
        ( \vc{y} - X \hat{\vc{\beta}})^T 
        \cdot 2( \vc{y} - X \hat{\vc{\beta}})  \\
 & = & -2 X^T (\vc{y} -  X \hat{\vc{\beta}} )  
\end{eqnarray*}
Setting this last expression equal to zero produces the {\em normal
equations}, 
\[
          X^T (\vc{y} -  X \hat{\vc{\beta}} ) = X^TE = \vc{0}
\]
because this asserts that the columns of $X$ are orthogonal to, or in
older mathematical language ``normal to'', the residuals $E =
\vc{y}-X\hat\vc{\beta}$.  Distributing $X^T$ across the
$(\vc{y}-X\hat\vc{\beta})$ term we obtain
\[
          X^T \vc{y} =  X^TX \hat{\vc{\beta}} 
\]
and---assuming that $X^TX$ is an invertible matrix---we can premultiply
this equation by $(X^TX)^{-1}$ to obtain
\[
  \hat{\vc{\beta}} = { ( X^T X ) }^{-1} X^{T} \vc{y}
\]

Note that it is just as easy to obtain the fitted values $\hat\vc{y}$:
\[
        \hat\vc{y} = X\hat\vc{\beta} = X(X^TX)^{-1}X^{T}\vc{y}
\]
The matrix $H=X(X^TX)^{-1}X^{T}$ is often called the {\em hat matrix}
for linear regression, because the effect of multiplying $y$ by $H$ is
to ``put the hat on'':
\[
        \hat\vc{y} = H y
\]

{\bf Remarks:}
\begin{itemize}
\item The fact that the fitted values $\hat\vc{y}$ are a \ulem{linear
transformation of the original responses} $\vc{y}$ is \ulem{really}
what makes linear regression ``linear''.  Nonparametric regression
methods such as 
\begin{itemize}
\item Kernel smoothing
\item Spline smoothing
\item \ulem{but not loess}
\end{itemize}
can be shown to be linear smoothers in exactly this sense.  Linearity
makes obtaining distribution theory, large sample theory, etc.,
particularly easy.
\item Note that
\[
        \frac{\partial}{\partial\vc{y}} \hat\vc{y} = H
\]
and in particular
\[
        \frac{\partial}{\partial{y_i}} \hat y_i = h_{ii},
\]
the $i^{th}$ diagonal element of the hat matrix $H$.  Thus, the
larger $h_{ii}$ is, the more sensitive $\hat y_i$ is to changes in
$y_i$; this further explains why $h_{ii}$ is a good place to start for
influence diagnostics.
\item We can also obtain the residuals using $H$:
\[
        \vc{E} = \vc{y} - \hat\vc{y} = \vc{y} - H\vc{y} = (I-H)\vc{y}
\]
\item $H$ has many remarkable properties.  Let 
\begin{eqnarray*}
        P & = & H \\
        Q & = & 1-H 
\end{eqnarray*}
Then it is not hard to show, using the definition of $H$ and/or the
normal equations, that
\[
        \begin{array}{rrl}
        P^T = P, & Q^T = Q: & \mbox{symmetric} \\
        P^2 = P, & Q^2 = Q: & \mbox{idempotent} \\
        PQ  = 0: &         & \mbox{orthogonal} \\
        P + Q = I: &       & \mbox{complete}
        \end{array}
\]
and indeed $P$ and $Q$ are ``projection matrices'' in $n$-dimensional
real space: every vector $x$ in $\Re^n$ can be written uniquely as an
orthogonal decomposition
\[
        x = x_P + x_Q
\]
where $x_P = Px \in \mbox{ColSpace}\,(P)$, $x_Q = Qx \in
\mbox{ColSpace}\,(Q)$, and $x_P\perp x_Q$ (they are orthogonal).
\item {\em Calculating the projection and properties of the
projection matrix are important, but the algebraic formula
$H=X(X^TX)^{-1}X^{T}$ involving matrix inversion need not be the best
computational formula.  See Thisted's Statistical Computing text for
linear regression estimation formulae that are actually useful
computationally.}
\end{itemize}



\section{Some Geometry}
Consider a dependent variable $\vc{y} = {( y_1, \ldots , y_n ) }^T$
and $k=2$ independent variables,
\[
  \vc{X_1} =  \left ( 
  \begin{array}{c}
           x_{11} \\
           x_{21} \\
           \vdots \\
           x_{n1}
   \end{array}
   \right )
\;\;\;\; 
  \vc{X_2} =  \left ( 
  \begin{array}{c}
           x_{12} \\
           x_{22} \\
           \vdots \\
           x_{n2}
   \end{array}
   \right )
\]

\begin{figure}
\psfig{figure=geometry.eps}
\end{figure}

We wish to consider the model $y_i = \beta_0 + \beta_1 x_{i1} +
\beta_2 x_{i2} + e_i$.  Letting
\[
  X = \left (
  \begin{array}{lll}
        1 & x_{11} & x_{12} \\
        1 & x_{21} & x_{22} \\
       \vdots & \vdots & \vdots \\
        1 & x_{n1} & x_{n2} 
  \end{array}
  \right )
= ( \vc{1}, \vc{X_1}, \vc{X_2} ),
\]
$\vc{\epsilon} = {( \epsilon_1, \ldots , \epsilon_n)}^T $, and
$\vc{\beta} = {(\beta_1, \beta_2, \beta_3)}^T $, we may again write
this in matrix form as
\[
  \vc{y} = X \vc{\beta} + \vc{\epsilon}
\]
The least-squares fitted values $\hat\vc{y}$ minimize the squared
distance $\Vert \vc{y} - \hat\vc{y}\Vert^2$, and so, looking at the
picture and thinking about elementary geometry, it's clear that the
vector $\vc{E} = \vc{y} - \hat\vc{y}$ should be perpendicuar to ${\cal
M} = \mbox{Span}\,(\vc{1},\vc{X_1},\vc{X_2})=\mbox{ColSpace}\,(X)$,
the vector space spanned by the columns of $X$.

On the other hand, if $ X^T X $ is of full rank, then the algebra we
did above asserts that the least squares estimates satisfy
\[
  \hat{\vc{\beta}} = {(X^T X)}^{-1} X^T \vc{y}
\]
and the fitted or predicted values for the $y$'s satisfy
\[      
        \hat\vc{y} = H y = X(X^TX)^{-1}X^T \vc{y} = X \hat\vc{\beta}
\]
This equation shows immediately that $\hat\vc{y}\in
\mbox{ColSpace}\,(X)$, and since $\vc{E} = (I-H)\vc{y}$ we also know
that, for any $\vc{\beta}$ (least-squares estimates or not)
\[
        (X\vc{\beta})^T\vc{E} 
        = \vc{\beta}^TX^T(I-H)\vc{y} 
        = \vc{\beta}^T(X^T - X^TX(X^TX)^{-1}X^T)
        = \vc{\beta}^T(X^T - X^T) 
        = \vc{0}
\]
so that the residuals $\vc{E}=\vc{y} - \hat\vc{y}$ really are
orthogonal to every possible linear combination of the columns of $X$.





\section{Assumptions and Standard Error Calculations} 

All of the preceding was just algebra; lets do some statistics.
The model is
\[
  \vc{y} = X \vc{\beta} + \vc{\epsilon}
\]
and we assume that
\begin{enumerate}
\item $X$ is of full rank (this is easily fixed)
\item $ E(\vc{\epsilon}) = \vc{0} $
\item $ \mbox{Var} ( \vc{\epsilon} ) = \sigma^2 \vc{I}_{n \times n} $ which
incorporates the equal variance and zero covariance assumptions
\item $ \vc{\epsilon} \sim N(\vc{0}, \sigma^2 \vc{I} ) $ a multivariate normal
distribution.
\end{enumerate}

\subsection{Digression --- Multivariate Normal Distribution} 

\paragraph{Definition} \ \\

A random vector
\[
  \vc{x} = \left ( \begin{array}{c}
                    x_1 \\
                    x_2 \\
                    \vdots \\
                    x_n
                    \end{array}
            \right )
\]
is said to have a multivariate normal distribution $ N( \vc{\mu}, \Sigma ) $
if 
\begin{eqnarray*}
  \vc{\mu} = \left ( \begin{array}{c}
                    \mu_1 \\
                    \mu_2 \\
                    \vdots \\
                    \mu_n
                    \end{array}
            \right )
\\
  \Sigma = \left ( \begin{array}{cccc}
                   \sigma_{11} & \sigma_{12} & \ldots & \sigma_{1n} \\
                   \sigma_{21} & \sigma_{22} & \ldots & \sigma_{2n} \\
                   \vdots   &   \ddots & \ddots & \vdots \\
                   \sigma_{n1} & \sigma_{n2} & \ldots & \sigma_{nn} 
                    \end{array}
              \right )
\mbox{, and}
\\
  E( \vc{x} ) = \vc{\mu}
\mbox{, and}
\\
  \mbox{Cov}( x_i, x_j ) = \sigma_{i,j}
\end{eqnarray*}
{\bf AND }
the distribution of 
\[
  \vc{a}^T \vc{x} \;\; \mbox{is} \;\; N ( \vc{a}^T \vc{\mu} , \vc{a}^T
       \Sigma \vc{a} )
\]
for {\em all } vectors $ \vc{a} $.

In other words {\em it is multivariate normal iff every possible linear
combination is univariate normal.}



\paragraph{Transformation Properties} \ \\

If $\vc{x} \sim N_n ( \vc{\mu} , \Sigma ) $ and $ W_{p \times n} $ is
a matrix of constants then
\[
  W \vc{x} \sim N_p ( W \vc{\mu} , W \Sigma W^T )
\]

Now let $\Sigma_{n \times n} $ be decomposed as $ V_{n \times n} V_{n
\times n}^T $, where each $V$ is of full rank (called {\em a} square
root of $W$.  This can always be done (in many ways).  Then
\[
  V^{-1} \vc{x} \sim N_n (V^{-1} \vc{\mu},~
  \underbrace { V^{-1} ( V V^T ) V^{-T} }_{I_{n \times n}} )
\]
Thus, if $\vc{x}$ is a random vector with mean $\mu$ and variance
$\Sigma$ and $ \Sigma = V V^T $, then $ V^{-1} \vc{x} $ has mean $
V^{-1} \vc{\mu} $ and variance $I$.



\subsection{Sampling Distributions of $\hat\vc{\beta}$ and $\hat\vc{y}$}

We are now in the position to say something about $\hat\vc{\beta} = { ( X^T X )
}^{-1} X^{T} \vc{y}$ and $\hat\vc{y} = H \vc{y}$.

For starters, note that
\begin{eqnarray*}
        E( \vc{y} ) & = & X \vc{\beta}, \mbox{ and} \\
       \mbox{Var} ( \vc{y} ) & = & \sigma^2 \cdot I . 
\end{eqnarray*}
Thus
\[
  \begin{array}{lclcl}
  E ( \hat\vc{\beta} ) & = &   {( X^T X )}^{-1} X^{T} E (\vc{y}) & & \\
        & = &   {( X^T X )}^{-1} X^{T} X \vc{\beta} & =  & 
                  \vc{\beta} 
  \end{array}
\]
and
\begin{eqnarray*}
\mbox{Var}( \hat\vc{\beta} ) & = & { (X^{-T}X) }^{-1} X^{T} \mbox{Var}(\vc{y}) X {(X^{-T}X)}^{-1} \\
             & = & { (X^{-T}X) }^{-1} X^{T} ( \sigma^2 \cdot I ) X {(X^{-T}X)}^{-1} \\
             & = & \sigma^2 {( X^{T}X )}^{-1}
\end{eqnarray*}
The fitted values $\hat\vc{y}$ and residuals $\vc{E}$ are even easier
\begin{eqnarray*}
        E[\hat\vc{y}] & = & E[H\vc{y}] ~ = ~ HX\vc{\beta} ~ = ~ X\vc{\beta}
        ~ = ~ E[\vc{y}] \\
        E[\hat\vc{E}] & = & E[(I-H)\vc{y}] ~ = ~ E[\vc{y}] - E[H\vc{y}] 
        ~ = ~ E[\vc{y}] - E[\vc{y}] ~ = ~ \vc{0} 
\end{eqnarray*}
and
\begin{eqnarray*}
        \Var(\hat\vc{y}) & = & \Var(H\vc{y}) ~ = ~ H \Var(\vc{y}) H^T
        ~ = ~ H \sigma^2 I H^T ~ = ~ \sigma^2 HH^T  ~ = ~ \sigma^2 H \\
        \Var(\vc{E}) & = & \Var((I-H)\vc{y}) ~ = ~ \cdots ~ = ~ 
        \sigma^2 (I-H)
\end{eqnarray*}
Together with the assumption that $\vc{\epsilon}\sim
N(\vc{0},\sigma^2I)$, or $\vc{y}\sim
N(X\vc{\beta},\sigma^2I)$, these facts show that

\hfil\begin{tabular}{|p{6in}|}
\hline
\[
        \hat\vc{\beta} \sim N(\vc{\beta}, \sigma^2 (X^TX)^{-1}) ~~~~~~
        \hat\vc{y} \sim N(X\vc{\beta},\sigma^2H) ~~~~~~
        \vc{E} \sim N(\vc{0},\sigma^2(I-H))
\]
\\ \hline
\end{tabular}

How do we estimate $\sigma^2$?  From our earlier work we expect that
we can get an unbiased estimate as a multiple of $RSS$.  To calculate
$E[RSS]$ we need a few preliminaries.

First we need a fact from linear algebra, namely that if both
$AB$ and $BA$ are legal matrix multiplications, then
\[
        \mbox{trace}\,(AB) = \mbox{trace}\,(BA)
\]
where trace$\,(A) = \sum_i a_{ii}$ is the sum of the diagonal
elements of $A$.  We immediately put this to good use, to see that
\[
        \mbox{trace}\,(H) = \mbox{trace}\,(X(X^TX)^{-1}X^T)
        = \mbox{trace}\,(X^TX(X^TX)^{-1}) 
        = \mbox{trace}\,(I_{(k+1)\times(k+1)}) 
        = k + 1
\]
Next we need a result about the expected value of a quadratic form
\[
        E[\vc{y}^TA\vc{y}] 
        = \mbox{trace}\,(A\Var(\vc{y})) + E[\vc{y}]^TAE[\vc{y}]
\]
(this is analogous to the univariate fact that $E[cy^2] = c^2\sigma^2
+c^2E[y]^2$).  Finally we note that
\begin{eqnarray*}
RSS & = & {( \vc{y} - \hat{\vc{y}})}^T {( \vc{y} - \hat{\vc{y}})} \\
 & = & \vc{y}^T [ I - H ]^T [ I - H ] \vc{y} \\
 & = & \vc{y}^T [ I - H ] \vc{y}
\end{eqnarray*}
With these three bits of information we calculate
\begin{eqnarray*}
E( RSS ) & = & E [ \vc{y}^T [ I - H ] \vc{y} ] \\
 & = & \mbox{trace}\,([I-H]\sigma^2I) + (X\beta)^T(I-H)(X\beta)\\
 & = & \sigma^2 \mbox{trace}\,(I-H) + \vc{0}\\
 & = & \sigma^2 [ n - (k+1) ]
\end{eqnarray*}
Consequently,
\[
        E\left[\frac{RSS}{n-(k+1)}\right] = \sigma^2
\]
so that 
\[
        \hat\sigma^2 = \frac{RSS}{n-(k+1)}
\]
is an unbiased estimate of $\sigma^2$, as we expected.

\subsection{Variances of Fitted Values and Predicted Values}

We have done part of this already.  Note that
\begin{eqnarray*}
\hat{y}_i & = & \vc{X_i}^T \cdot \hat\vc{\beta} \\
          & = & \vc{X_i}^T {( X^T X )}^{-1} X^{T} \vc{y}
\end{eqnarray*}
so that
\[
  \mbox{Var} ( \hat{y}_i ) = \sigma^2   \vc{X_i}^T {( X^T X )}^{-1} \vc{X_i}
\]
Notice that the variance of the fitted value depends on the vector $X$
where $ y_i $ was observed.

At a new value $\vc{x}^*$ that was not in the original data set, the
fitted value is still
\[
        \hat y^* = \vc{x^*}^T \,\hat\vc{\beta}
\]
with variance
\[
        \Var(\hat y^*) = \sigma^2 \vc{x^*}^T {( X^T X )}^{-1} \vc{x^*}
\]
Note that $\vc{x^*}$ does not have to be part of the original $X$ matrix.
The variance calculated here is for estimating $E[y^*]$ in the equation
\[
        y^* = \vc{x^*}^T\beta + \epsilon^*
\]
Sometimes it is desirable to try to predict the actual value $y^*$
itself, and not just the expected value $E[y^*]$.  Assuming as usual
that $\epsilon^*$ is independent of the other $\epsilon$'s in the
problem, and comes from the same $N(0,\sigma^2)$ distribution, we have
simply
\begin{eqnarray*}
        \Var(\hat y^*_{pred}) & = & \Var(\hat y^*) + \Var(\epsilon^*) \\
        & = & \sigma^2 \vc{x^*}^T {( X^T X )}^{-1} \vc{x^*} + \sigma^2 \\
        & = & \sigma^2\left[1 + \vc{x^*}^T {( X^T X )}^{-1} \vc{x^*}\right]
\end{eqnarray*}

\section{Some final remarks}

Aside from being a clean way to derive and generalize all of the
standard formulas from simple regression, there are several reasons
for bothering with the linear algebra formulation:
\begin{itemize}\small\itemsep=0pt\parsep=0pt\parskip=0pt

\item The geometry suggested above is sometimes quite helpful in
figuring out what vague notions in regression mean.  For example,
``collinearity'' of $\vc{X_1}$ and $\vc{X_2}$ simply means that the
angle between $\vc{X_1}$ and $\vc{X_2}$ is so small that the position of
the hyperplane Span$\,(\vc{1},\vc{X_1},\vc{X_2})$ is very sensitive to
small changes in $\vc{X_1}$ or $\vc{X_2}$,

\item Standard error calculations for the $\beta$'s are clearly
straightforward (though messy calculations are hidden in the matrix
multiplications!).  Since we also have expressions for the
covariances, we can easily obtain standard errors for any linear
combination (contrast!) of the parameters.

\item The geometry is also useful for suggesting what is the ``same''
model---any parametrization that generates the same $H$ matrix is the
``same'' model as far as linear regression is concerned.  

\item Finally the geometry, mixed with a little distribution theory,
will be very helpful in creating and understanding standard errors
that are corrected for multiple inference, such as in Scheff\'e's method.

\end{itemize}

You may also be interested in knowing that, just as in the simple
regression case, the least-squares estimates are the same as the
maximum likelihood estimates of the $\beta$'s; and as before the MLE
of $\sigma^2$ is the biased estimate $RSS/n$.  

\end{document}

