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

\parskip=0.5em
\parindent=0pt

\unitlength=0.75mm

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

\frenchspacing

\long\def\comment#1{}

\let\bar\overline

\def\tp{^T}

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

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

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

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

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

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

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

% \def\ul#1{\setbox1\hbox{$\underline{#1}$}\dp1=0pt\box1}
\def\ul#1{\underline{#1}}
\def\ulem#1{\ul{\em #1}}

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

\def\xtx{X\tp X}

\def\Re{{\cal R}}

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

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

\begin{document}

\begin{center}\Large\bf 
        Variable Selection and Collinearity in Linear Regression 
\end{center}

\section{The Trouble with Collinearity}

Consider the linear regression problem
\[
        \mb y_{n\times1} = X_{n\times k'} \mb \beta_{k'\times 1} +
        \mb\epsilon_{n\times 1}
\]
where $k'$ is the total number of regressors (for regression with an
intercept, $k'=k+1$), and of course $n>k'$.

One can show that the rank $r$ of the model matrix $X_{n\times k'}$
(number of linearly independent columns of $X$) is the same as the
rank of $X\tp X_{k'\times k'}$; if $r<k'$ then $X\tp X$ will not be
invertible and the usual equations for computing regression estimates
and their standard errors will not be usable.

This is the case of {\em perfect collinearity} and, as we have seen
before (especially in the context of ANOVA models), essentially the
only fix is to reparametrize the model and/or throw away columns of
$X$, so that the new model matrix $X_{new}$ has the same column space
as $X$, but $X_{new}$ has exactly $r$ columns; for then $(X_{new}\tp
X_{new})_{r\times r}$ is invertible and the usual regression
estimation recipies work.

{\em Partial collinearity}, or just {\em collinearity} for short,
occurs when some of the columns of the $X$ matrix are {\em nearly}
linearly dependent.  As the columns become more and more nearly
linearly dependent, the behavior of the estimated regression
coefficients, and especially their standard errors, degrades.

Consider a simple case where we have a $y$ variable and two
predictors, $x_1$ and $x_2$, with the model:
\[
  \mb{y} = \mb{x_1} \beta_1 + \mb{x_2} \beta_2 + \epsilon
\]
(i.e., no constant term).  Let
\[
  X = \left [ \begin{array}{cc}
                x_{11} & x_{21} \\
                \vdots & \vdots \\
                x_{1n} & x_{2n}
            \end{array} \right ]
\]
then we know that $\hat{\beta} = {( \xtx)}^{-1} X\tp\mb{y} $
and 
${\var}(\hat{\beta}) =  \sigma^2_\epsilon {(\xtx)}^{-1}$.
Now when $X$ is nearly singular (i.e. near collinearity), then
${(\xtx)}^{-1}$ get big (like $1/0$ ) and the variances explode.

Consider the following numerical example.  Let
\[
  \mb{x_1} = \left [ \begin{array}{c} 1 \\ 2 \\ 3 \\ 4 \end{array} \right ]
  \qquad \mbox{and} \qquad 
  \mb{x_2} = \left [ \begin{array}{c} 1 \\ 2 \\ 3 \\ 4+\epsilon 
                        \end{array} \right ]
\]
then
\[
  \xtx = \left [ \begin{array}{ll}
   30 & 30 + 4 \epsilon \\
  30 + 4 \epsilon & 30 + 8 \epsilon + \epsilon^2 
  \end{array} \right ]
\]
and
\[
  {( \xtx)}^{-1} =
  \frac{1}{900 + 240\epsilon + 30\epsilon^2 - 900 -240\epsilon -16\epsilon^2}
  \cdot \left [ \begin{array}{ll}
  30 + 8\epsilon + \epsilon^2 & -30 + 4\epsilon \\
  -30 + 4\epsilon & 30
   \end{array} \right ]
\]
so as $\epsilon \longrightarrow 0 $, 
\[
\var(\hat{\beta}_2) \propto \frac{1}{14\epsilon^2} \longrightarrow \infty .
\]
Of course the variance estimate for $\beta_1$ explodes in just the
same manner.

Similar behavior occurs in all cases where collinearity is a problem.
This behavior is mostly of concern for the parameter estimates---the
fitted values $\hat y_i$ are (or should be) quite stable.  Numerical
problems can (potentially) upset the stability of fitted values unless
the software is {\em very} careful.  

More broadly, one can show that, for regression {\em with} an intercept,
\[
        \Var(\hat\beta_j) = \mbox{[j$^{th}$ diag element of
        $\sigma^2_\epsilon (\xtx)^{-1}$]} = 
        \frac{1}{1-R^2_j} \cdot \frac{\sigma^2_\epsilon}{(n-1)S^2_j}
\]
(where $R^2_j$ is the squared multiple correlation between $X_j$ and
all the other columns of $X$, and $S^2_j$ is the sample variance of
$X_j$, the $j^{th}$ column of $X$), though the algebra can be quite
tedious if you are not careful (see Weisberg, 1985, Exercise 8.10,
pp. 223--224; or Rawlings et al, 1998, pp. 372--373).

Since $R^2_j$ measures the degree of collinearity between $X_j$ and
the other columns of $X$, the ``variance inflation factor''
\[
        VIF_j = \frac{1}{1-R^2_j} \geq 1
\]
measures the inflation of the variance of $\hat\beta_j$ due to
collinearity in the columns of $X$; note that $VIF_j\equiv1$ iff $X_j$
is orthogonal to the other columns of $X$ (after centering $X_j$ and
all the other columns of $X$ by their column means).  Sometimes VIF's
are a useful diagnostic tool to determine whether collinearity is a
problem.  More often, though, you don't worry too much about
collinearity until you realize that it's getting in the way of some
inference that you want to make.

Fixing partial collinearity is like fixing perfect collinearity: you
either have to reparameterize the model or throw out some columns of
the $X$ matrix.  As we have seen with reparametrizing ANOVA models,
the effect of this can force a non-trivial change in the
interpretation of the model parameters.

{\em Example}

The data are from an unpublished Master's Thesis in Civil Engineering
by Carl Hoffstedt, as reported by Weisberg (1985, p.  206).  They
relate the accident report rate, in accidents per million vehicle
miles ($y$) to 13 potential independent variables.  The data include
39 sections of large highways in the state of Minnesota in 1973.  The
variables are:
{\small\begin{verbatim}
Rate - 1973 accident rate per million vehicle miles
Len  - length of segment in miles
Adt  - average daily traffic count in thousands (estimated)
Trk  - truck volume as a percent of total volume
Sp   - speed limit (in 1973, before the 55 mph limits)
lw   - lane width in feet
shl  - width in feet of outer shoulder on the roadway
itg  - number of freeway-type interchanges per mile in the segment
sigs - number of signalized interchanges per mile in the segment
acpt - number of access points per mile in the segment
lan  - total number of lanes of traffic in both directions
type -  1 if Federal aid interstate highway,
        2 if principal arterial interstate highway,
        3 if minor arterial interstate highway,
        4 otherwise
\end{verbatim}}\noindent
(The full dataset is {\tt highway.dat} on the class website).

We start with a naive regression analysis:
\begin{verbatim}
> high _ read.table("highway.dat",header=T)
> high[1:5,]
  Rate   Len Adt Trk Sp lw shl  itg sigs acpt lan type 
1 4.58  4.99  69   8 55 12  10 1.20    0  4.6   8    1
2 2.86 16.11  73   8 60 12  10 1.43    0  4.4   4    1
3 3.02  9.75  49  10 60 12  10 1.54    0  4.7   4    1
4 2.29 10.65  61  13 65 12  10 0.94    0  3.8   6    1
5 1.61 20.01  28  12 70 12  10 0.65    0  2.2   4    1

> X _ high[,-1]
> y _ high[,1]

> rawlm _ lm(Rate ~ . , data=high)
> summary(rawlm)

Call: lm(formula = Rate ~ ., data = high)
Residuals:
    Min      1Q  Median     3Q   Max 
 -1.867 -0.8856 0.09851 0.6906 2.887

Coefficients:
               Value Std. Error  t value Pr(>|t|) 
(Intercept)  12.5190   6.5290     1.9175   0.0658
        Len  -0.0564   0.0328    -1.7222   0.0965
        Adt   0.0093   0.0326     0.2841   0.7785
        Trk  -0.1271   0.1057    -1.2026   0.2396
         Sp  -0.0685   0.0637    -1.0756   0.2916
         lw  -0.3433   0.4967    -0.6912   0.4954
        shl  -0.0947   0.1245    -0.7606   0.4535
        itg   0.2868   1.1426     0.2510   0.8037
       sigs   0.4908   0.3941     1.2453   0.2237
       acpt   0.0684   0.0241     2.8376   0.0085
        lan  -0.0131   0.2743    -0.0477   0.9623
       type   0.2338   0.4893     0.4778   0.6366

Residual standard error: 1.172 on 27 degrees of freedom
Multiple R-Squared: 0.7527 
F-statistic: 7.472 on 11 and 27 degrees of freedom, the p-value is
    1.051e-05 

[ ... ]



> source("collin.q") # This function is available in the data area of
                     # the course website.  It calculates R^2 for each
                     # X regressed on all the other X's, providing 
                     # essentially the same information as VIF=1/(1-R^2).
                     # (Thanks to Howard Seltman.)

> collin(X,cutoff=0.7)  # R^2 for regressing each X on the other X's;
                        # report only those R^2's above 0.7.
Adt 0.902 
Sp 0.739 
shl 0.747 
itg 0.836 
lan 0.741 
type 0.743 
\end{verbatim}
Clearly there is a lot of collinearity!

% But before we go too far, we should do at least a cursory look at some
% diagnostics; see Figure~\ref{plot-lm1}.
% 
% \begin{figure}
% \centerline{\psfig{figure=plot-lm1.ps,height=3.5in}}
% \caption{Standard diagnostic plots for the initial fit of the highway
% accident data.}\label{plot-lm1}
% \end{figure}


\section{Variable Selection and Specification Error}

A more pleasant-sounding name for ``throwing out some of the columns
of $X$'' is {\em variable selection} or {\em model selection}, so
automatic methods of variable selection almost always get discussed in
the same breath as collinearity problems.  However it is important to
note that variable selection is really a separate issue from fixing
collinearity.  As illustrated above, {\em collinearity is primarily a
technical problem}.  On the other hand, {\em variable selection is
often at least partially a scientific problem}, related to the process
of sifting through independent variables to find ones that are
causally related to the dependent variable in a particular scientific
context.

Variable selection and model selection must deal with two fundamental
obstacles, whether in a classical context like least-squares
regression, or in a Bayesian or other context:
\begin{itemize}
\item {\em There are a lot of models to consider.}  For a regression
  with $k'$ predictors, there are $2^{k'}-1$ additive models to
  consider, and many many more if we allow interactions between two or
  more predictors at a time.  For example, if you have $k'=8$
  predictors, there are 255 models to consider, and if you have
  $k'=50$ regressor variables, as with some of your end-of-semester
  projects, there are $2^{50} - 1 = 1.1259 \times 10^{15}$ models to
  consider.  Even with fast computers, we aren't going to find the
  best model by exhaustive search.
\item {\em There are not enough models to consider.}  Even if you were
  able to consider all $2^{50} - 1 = 1.1259 \times 10^{15}$ models in a
  model with 50 known regressors, you know or suspect now that the
  ``right'' variables to use to predict differences in yield and so
  forth are probably not even in the data set.
\end{itemize}

The first problem is primarily a computational one, and automatic
variable selection methods are primarily heuristic methods to reduce
the number of models considered, but still generate some reasonably
good models.  Alternatively, and often more effectively, you can
select one or more good models by hand if you have some subject area
knowledge to guide you.

The second problem is not a theoretical or computational problem; it
is fundamentally a scientific problem.  Statistical theory can shed
light on the consequences of this problem but not on how to solve it.

From our earlier work this semester (especially Exercise 7.3,
pp. 233-234 in Rawlings et al., 1998), we know:
\begin{itemize}
\item \ul{\em Overspecifying the model,} i.e. including $X$'s that aren't
  necessary, can {\em increase our uncertainty (higher SE's)} in
  parameter estimates and predictions, but will not affect the
  unbiasedness of estimates.

\item \ul{\em Underspecifying the model,} i.e. failing to include $X$'s
  that are needed and are collinear with $X$'s already in the model,
  can {\em increase bias} (including having completely the wrong
  sign!) and {\em reduce standard errors} for estimates of
  coefficients of $X$'s already in the model.
\end{itemize}
Variable selection tries to achieve a specification that is ``not too
hot'' and ``not too cold'' in the above senses.

We will consider two or three statistically-motivated approaches to
variable selection:
\begin{itemize}
\item {\em Automatic, ``stepwise'' procedures} that ``score'' each model and
  select the model with the highest (or lowest) score.
\item {\em Principal component transformations} of the $X$ matrix that
  produces a new regressor matrix $W$.  Variable selection is
  ``easier'' with $W$ because the columns of $W$ will be uncorrelated;
  and in many cases 
\item {\em Other methods} such as minimizing prediction error and
  so-called ``biased regression'' methods.
\end{itemize}
Note that {\em all} of these methods finess the
scientific question: we assume that the largest model we consider has
enough variables in it so that there is no bias in estimating
parameters due to omitting relevant regressors from the model.
If $X$ is missing some
important regressors, then the final model you get from these methods
will be underspecified, and your $\hat\beta_j$ estimates will be
biased, possibly severely (albeit of lower variance).


\section{Automatic Variable Selection}

In this section we will look at two or three automatic
methods for variable selection.  These are essentially all heuristic
methods for searching through the space of possible models in a
regression problem.

Each heuristic can be combined with several different ways of
``scoring'' the models, so we begin by considering some different ways
of scoring models.  For each criterion that follows, $p$ is the number
of regression parameters (regressors {\em plus intercept}) in the
(current) model.  The maximum value for $p$ is $k'=k+1$, the total
number of regressors (plus intercept) available in the $X$ matrix.
\begin{itemize}

\item $\mathbf R^2 = 1 - RSS/TSS$\\
  $R^2$ can be (and is) used to compare models with the {\em same
    number} of parameters $p$.  Thus for the highway data we
  could compare all models which include {\tt Len} and {\tt acpt} and
  (say) one extra variable, using the appropriate $R^2$ values.  The
  model with the largest value of $R^2$ would be considered as (a
  candidate for) the ``best''.  However, since $R^2$ must grow as
  regressors are added to the model, $R^2$ is not a good way to
  compare models with different numbers of regressors.

\item {\bf Extra-sum-of-squares F-test} \\
The F-test can be used to compare a model and a subset of that model.  For
example the model $(x_1, x_3, x_5 )$ and the model 
$(x_1,x_3,x_5,x_7,x_8)$ can be compared using an F-test.  
\begin{itemize}
\item Because the log-likelihood for normal-errors linear regression
  is proportional to the RSS for the model, the extra-sum-of-squares
  F-test 
\[
        F =  \frac{(RSS_{small}-RSS_{big})/(p_{big}-p_{small})}
        {RSS_{big}/(n-p_{big})}
\]
(in our current notation $p_{big}$ and $p_{small}$ count the intercept too)
makes the same decisions (and asymptotically with the same level) as
the likelihood ratio test for comparing nested models
\begin{eqnarray*}
  \lambda & = & 
  -2\left\{
     \log(L_{big}(\ul y; \hat{\ul\beta}_{big},\hat\sigma^2_{big})) 
  -  \log(L_{small}(\ul y; \hat{\ul\beta}_{small},\hat\sigma^2_{small}))
    \right\} 
   \\
  & = & n\log\frac{\hat\sigma^2_{small}}{\hat\sigma^2_{big}} + 
  \frac{ 
  (\ul y - X_{small}\hat{\ul\beta}_{small})\tp
  (\ul y - X_{small}\hat{\ul\beta}_{small})}
  {\hat\sigma^2_{small}} 
  - 
  \frac{ 
  (\ul y - X_{big}\hat{\ul\beta}_{big})\tp
  (\ul y - X_{big}\hat{\ul\beta}_{big})}
  {\hat\sigma^2_{big}} 
  \\
  & = & n\log\frac{RSS_{small}}{RSS_{big}} + 
  \frac{ 
  (\ul y - X_{small}\hat{\ul\beta}_{small})\tp
  (\ul y - X_{small}\hat{\ul\beta}_{small})}
  {RSS_{small}/n} 
  - 
  \frac{ 
  (\ul y - X_{big}\hat{\ul\beta}_{big})\tp
  (\ul y - X_{big}\hat{\ul\beta}_{big})}
  {RSS_{big}/n} 
  \\
  & = & n\log\frac{RSS_{small}}{RSS_{big}} + n - n
  ~ = ~ n\log\left(F\cdot\frac{p_{big}-p_{small}}{n-p_{big}} + 1\right) 
\end{eqnarray*}
\item Since the classical distribution theory for
  the likelihood ratio test and F-test of non-nested models is
  unwieldy, these tests are not directly useful for comparing
  non-nested models.
\item Comparing two models with $R^2=1-RSS/TSS$ is equivalent to
  looking at the difference $RSS_{small} - RSS_{big}$ (or the ratio of
  RSS's), i.e. it is like doing the $F$ or likelihood ratio test
  without the benefit of distribution theory to help decide when the
  difference in RSS's is big enough to care about.  All of the other
  methods below replace distribution theory with a numerical
  penalization for more parameters, so that model comparisons for
  which classical distribution theory isn't available (especially
  non-nested models) may be performed.
\end{itemize}

\item {\bf Adjusted $\mathbf R^2$} \\
\[
  R^2_{adj,p} = 1 - \left( \frac{n-1}{n-p} \right ) 
                \left ( 1 - R^2 \right )
\]
where $p$ is the number of regression parameters (including intercept)
included and $R^2$ is for the model with $p$ parameters.  This is an
adjustment which is designed to compensate for the fact that otherwise
as $p \nearrow $ then $R^2 \nearrow$.  Note that there are some
problems with this in that $R^2_{adj,p}$ can be negative, but this is
generally not a problem when the model fits reasonably well.  The idea
is to choose the model with the largest value of $R^2_{adj,p}$.
\item {\bf Mallows' $\mathbf C_p$} \\
\begin{equation}
  C_p = \frac{RSS_p}{\sigsq} + 2p - n
  \label{mallows}
\end{equation}
where $RSS_p$ is computed from the ``current'' model with $p$
parameters (including intercept), but $\sigsq$ is estimated
from the full model or the largest model which will be considered (and
which is consequently a superset of all the other models).  Clearly
this is related to the F-statistic, but it is adjusted.  Since
$E[RSS_p]=(n-p)\sigma^2$ under the ``correct'' model, it follows
that $E[C_p] \approx p $; and in fact $E[C_p]\approx p$ for every
model containing the ``correct'' model as well.  The idea then is to
select models with $C_p$ near $p$ as candidates for further
investigation.  A useful method for selecting such models is to plot
$C_p$ against $p$.

The full model always has
$C_{k+1} = k+1$, but notice some of the simpler models may also have
reasonable $C_p$'s.  We usually try to either minimize $C_p$, or
choose the smallest $p$ for which $C_p\approx p$.

\item {\bf Akaike's Information Criterion (AIC)} \\
This is one of the most popular ``penalized likelihood'' methods for
scoring goodness-of-fit.  
\begin{eqnarray*}
  AIC & = & -2\log(L_p(\ul y; \hat{\ul\beta}_p,\hat\sigma^2_p)) + 2p \\
  & = & n\log\hat\sigma^2_p + 
  \frac{ 
  (\ul y - X_p\hat{\ul\beta}_p)\tp
  (\ul y - X_p\hat{\ul\beta}_p)
  }
  {\hat\sigma^2_p} - 2p \\
  & = & n\log RSS_p - n\log(n) + n + 2p 
\end{eqnarray*}
We want to choose the model with the smallest AIC (sometimes -AIC is
reported, so that big values are for good models).

Note that for comparing different models on the same data set (so
additive terms involving $n$ only cancel out).  Also, for comparing
models with the same $n$, AIC and $Cp=RSS_p/\hat\sigma^2_\epsilon +
2p-n$ are clearly monotone functions of each other: smallest AIC
and smallest $C_p$ choose the same model.

For nested models
\[
        AIC_{small} - AIC_{big} = n\log \frac{RSS_{small}}{RSS_{big}}
        - 2(p_{big}-p_{small})
\]
so that this is effectively a likelihood ratio comparison, penalized
for extra parameters in the model.  It is known that AIC tends to
choose models with more parameters (regressors) than the true model;
but it is very good at minimizing prediction error.

\item {\bf Schwarz's Bayesian Information Criterion (BIC, a.k.a. SBC)} \\
This is the other popular ``penalized likelihood'' methods for
scoring goodness-of-fit.  
\begin{eqnarray*}
  BIC & = & -2\log(L_p(\ul y; \hat{\ul\beta}_p,\hat\sigma^2_p)) +
  p\log(n) \\
  & = & n\log RSS_p - n\log(n) + n + p\log(n) 
\end{eqnarray*}
Again the goal is to choose the model with the smallest BIC (sometimes
-BIC is reported, so that big values are for good models), and again
the terms involving only $n$ cancel out in model comparisons.  Note
that the penalty for extra parameters is now multiplied by $\log(n)$
rather than by 2 as in AIC, therefore BIC will tend to pick smaller
models than AIC.

The penalized likelihood ratio statistic is now
\[
        BIC_{small} - BIC_{big} = n\log \frac{RSS_{small}}{RSS_{big}}
        - (p_{big}-p_{small})\log(n)
\]
Perhaps more importantly, $BIC_{small} - BIC_{big}$ can be shown to be
a reasonable asymptotic approximation to the {\em Bayes factor}
comparing the two models, as $n$ grows, and therefore it can be shown
to select the ``true'' model with probability tending to one as $n$
grows (Schwarz, 1978, {\em Ann.  Stat.}; Kass and Raftery, 1995, {\em
  JASA}).

\end{itemize}


We now have some tools for comparing models, but they seem to require
that we fit all (or most) of the models.  This would not be much use,
unless there was a simple way of fitting models, or deciding which
models to think carefully about.  In the next several subsections we
describe and illustrate several such methods using the highway
accident data.

These automatic methods do not absolve you (the analyst) from looking
at the data and thinking about the substantive problem that the data
came from.  You should still
\begin{enumerate}
\item Examine the data before selecting model, to assure there are no
glaring errors
\item Carefully examine the best modesl for influential or outlying
observations
\item Make sure you can interpret and justify the best model (and
explain the results to others) in terms of the substantive/scientific
inquiry.
\end{enumerate}

{\em Continuation of Example}

Figure~\ref{fig:cp} gives an example of a $C_p$ plot, for the highway
data.  The solid line is $C_p = p $.  

The code for generating the plot follows.  It exploits the ``leaps()''
function in Splus (and R), which exhaustively\footnote{{\tt leaps()}
  does not actually have to fit all models, since it can sometimes
  guess that if a particular model fits badly, all submodels of that
  model will fit even worse.} fits all models and
reports the $n$ best (lowest $C_p$) models with $p=1$ predictor, $p=2$
predictors, $p=3$ predictors, etc.
\begin{verbatim}
> result _ leaps(X,y,nbest=2) # generates Cp for the 2 best (lowest Cp) 
                              # models with 1 predictor, 2 predictors, etc.
> names(result)
[1] "Cp"    "size"  "label" "which" "int"  
> plot(result$size,result$Cp)
> abline(0,1)
> title("Cp Plot for the Highway Data")
\end{verbatim} 
In this case, models of size $p$ $=$ 4 or 5 are suggested; a quick
look at {\tt help(leaps)} tells us that this means 3 or 4 regressors,
plus the intercept.  The suggested models are
\begin{verbatim}
> data.frame(terms=result$label,Cp=result$Cp,size=result$size)[5:8,]
             terms       Cp size 
5      Len,Sp,acpt 2.457994    4
6      Trk,Sp,acpt 3.106290    4
7 Len,Sp,sigs,acpt 1.563455    5
8  Len,Trk,Sp,acpt 1.930997    5
\end{verbatim} 
Note that $C_p/n \approx \hat\sigma^2_p +2p/n$, so one way to
think about $C_p$ is that it is a way of minimizing residual standard
error, with a penalty for including ``too many'' parameters in the
model.  It is also related to a penalized likelihood method, AIC, as
we shall see next.
\begin{figure}
\hfil\psfig{figure=highway-cp.ps,angle=-90,width=5in}
\caption{Plot of $C_p$ vs $p$.  Best model has least $C_p$ close to
 $p$.}\label{fig:cp}
\end{figure}

Another function (available only in Splus) called {\tt stepwise()},
can be used to generate the $n$ best (in terms of minimizing RSS,
which is equivalent to maximizing $R^2$) models with $p=1$ predictor,
$p=2$ predictors, $p=3$ predictors, etc.\footnote{Like {\tt leaps()},
  {\tt stepwise()} does not actually have to fit all models, since it
  can sometimes guess that if a particular model fits badly, all
  submodels of that model will fit even worse.}  The function {\tt
  sum.step.q} available on the course website (thanks to Howard
Seltman!), can summarize the output of {\tt stepwise} in terms of AIC,
BIC, etc., to choose the best model.
{\small\begin{verbatim}
> sw _ stepwise(X,y,meth="exhaustive")
> sum.step(sw,high,"Rate",includeAIC=T,includeR2=T)
         BIC   AIC    R2   R2a Len Adt Trk Sp lw shl itg sigs acpt lan type 
 3(#1) -66.0 -62.6 0.694 0.667   1   0   0  1  0   0   0    0    1   0    0
 4(#1) -66.1 -61.9 0.720 0.687   1   0   0  1  0   0   0    1    1   0    0
 4(#2) -66.3 -62.2 0.717 0.683   1   0   1  1  0   0   0    0    1   0    0
 3(#2) -66.3 -63.0 0.688 0.661   0   0   1  1  0   0   0    0    1   0    0
 4(#3) -67.0 -62.8 0.707 0.673   0   0   1  1  0   0   0    1    1   0    0
 5(#1) -67.1 -62.1 0.732 0.692   1   0   1  1  0   0   0    1    1   0    0
 3(#3) -67.3 -64.0 0.671 0.643   0   0   0  1  0   0   0    1    1   0    0
 5(#3) -67.7 -62.8 0.724 0.682   1   0   0  1  0   0   0    1    1   0    1
 5(#2) -67.7 -62.7 0.724 0.682   1   0   0  1  0   1   0    1    1   0    0
 2(#1) -68.2 -65.7 0.622 0.601   0   0   0  1  0   0   0    0    1   0    0
 6(#1) -68.6 -62.8 0.739 0.689   1   0   1  1  0   1   0    1    1   0    0
 6(#2) -68.8 -63.0 0.736 0.686   1   0   0  1  0   0   1    1    1   0    1
 6(#3) -68.9 -63.0 0.735 0.685   1   0   1  1  1   0   0    1    1   0    0
 2(#2) -68.9 -66.4 0.607 0.586   0   0   1  0  0   0   0    0    1   0    0
 2(#3) -69.3 -66.8 0.600 0.578   0   0   0  0  0   0   0    1    1   0    0
 7(#1) -70.0 -63.4 0.746 0.688   1   0   1  1  0   0   1    1    1   0    1
 7(#2) -70.1 -63.4 0.745 0.687   1   0   1  1  1   1   0    1    1   0    0
 7(#3) -70.2 -63.5 0.744 0.686   1   1   1  1  0   0   0    1    1   0    1
 1(#1) -71.1 -69.4 0.517 0.504   0   0   0  0  0   0   0    0    1   0    0
 8(#1) -71.6 -64.2 0.750 0.684   1   0   1  1  1   1   1    1    1   0    0
 8(#2) -71.7 -64.2 0.750 0.683   1   1   1  1  1   1   0    1    1   0    0
 8(#3) -71.8 -64.4 0.748 0.680   1   0   1  1  0   1   1    1    1   0    1
 1(#2) -73.1 -71.5 0.464 0.449   0   0   0  1  0   0   0    0    0   0    0
 9(#2) -73.5 -65.2 0.752 0.675   1   0   1  1  1   1   1    1    1   0    1
 9(#1) -73.5 -65.2 0.752 0.675   1   1   1  1  1   1   0    1    1   0    1
 9(#3) -73.6 -65.3 0.751 0.673   1   1   1  1  1   1   1    1    1   0    0
10(#3) -75.5 -66.4 0.752 0.663   1   0   1  1  1   1   1    1    1   1    1
10(#2) -75.5 -66.4 0.752 0.664   1   1   1  1  1   1   0    1    1   1    1
10(#1) -75.5 -66.3 0.753 0.664   1   1   1  1  1   1   1    1    1   0    1
11(#1) -77.5 -67.5 0.753 0.652   1   1   1  1  1   1   1    1    1   1    1
 1(#3) -77.8 -76.1 0.319 0.300   0   0   0  0  0   0   0    1    0   0    0
 0(#3) -83.4 -82.6 0.000 0.000   0   0   0  0  0   0   0    0    0   0    0
 1(#2) -83.7 -82.0 0.077 0.052   0   0   0  0  0   0   0    0    0   0    1
\end{verbatim}}\noindent
    Note that {\tt sum.step()} reports $-AIC/2$ and $-BIC/2$ (and
    omits terms like $n\log(n)$ that depend only on $n$) so large
    values (values near zero) are good; the output is sorted in order
    from best BIC (nearest zero) to worst BIC (farthest from zero).

The variables in the best several models here are:
\begin{verbatim}
Len  - length of segment in miles
Trk  - truck volume as a percent of total volume
Sp   - speed limit (in 1973, before the 55 mph limits)
sigs - number of signalized interchanges per mile in the segment
acpt - number of access points per mile in the segment
\end{verbatim}
Generally speaking a change of 2 to 6 is considered moderate evidence
of improved fit in BIC (this translates to a 1 to 3 difference in what
{\tt sum.step} reports), so there isn't much difference among the
first five or so models listed here.
% 
% but what are these variables telling us about the type of road that is
% accident-prone?  {\em More on this later\ldots}




\subsection{Stepwise Selection Procedures}

The exhaustive procedures used to illustrate $C_p$, AIC, BIC, $R^2$
and $R^2_{adj}$ above are not useful when there are many (more than 30
or so) regressors to consider in the problem.  In this case we need
one or another heuristic algorithm to search the model space for a few
good candidate models to explore further.

The standard search algorithms are variations on ``greedy''
strategies, that is, they search among models closely related to the
``current'' model that maximize the ``score'' (whether $F$-test,
$C_p$, or whatever), and jump to that model.  In statistics, these
methods are all called ``stepwise'' methods.

There are (at least) three basic strategies for stepwise selection of
models:
\begin{itemize}\tight
\item Start with a model an keep adding until you can't add
  anymore---called Forward Selection 
\item Start with a (big) model and take away variables while you 
can---called Backward Elimination
\item Do both at the same time---a true stepwise technique
\end{itemize}
Even within these there are several possible ways of choosing a
model.  Note that all of these methods are very heuristic; there is no
guarantee that the best model---even in terms of the specific
criterion used---will be found!

\paragraph{Forward Selection.} \

Start with the intercept-only straw-man model; and add one variable at
a time.
\begin{enumerate}\tight
\item Stop after a given number of steps or (equivalently) when a subset of
a certain size has been reached;  \ul{\bf or}
\item Use F-tests to decide which variable to add next, and stop when no
variable has a big F-test; \ul{\bf or}
\item Stop when adding the next variable would result in collinearity 
  $ \Longleftrightarrow $ Stop when $R^2_{new}$ (squared multiple
  correlation for regressing $X_{new}$ on the other $X$'s already in
  the model) gets too big.
\end{enumerate}

\paragraph{Backward Elimination.} \

Start with all the variables in the model; and remove one variable at
a time.
\begin{enumerate}\tight
\item Stop after a given number of steps;  \ul{\bf or}
\item Use F-tests to decide which variable to remove next, and stop
  when removing each variable produces a big F-test.
\end{enumerate}

\paragraph{True Stepwise.} \ 

One needs to provide criteria for including and deleting terms from
the model.  This can be done by specifying the values of the
F-statistic, F-IN and F-OUT, required for entry into or removal from
the the model.  Then the stepwise procedure repeatedly applies the
following rules, until no further changes in the model occur.
\begin{enumerate}\tight
\item Remove a variable is F-test is less than F-OUT
\item Remove if $<$ F-OUT and $R^2$ is bigger than the previous model with
this number of predictors
\item Swap two variables (in and out) if this improves $R^2$
\item Add a variable if F-test is greater than F-IN
\end{enumerate}
% Generally some combination of these criteria is used.  
%
A common choice is F-IN = F-OUT = 4.  Making F-IN large and F-OUT
small means that the current model is hard to get into and easy to get
out of.  By making F-OUT = 0 one obtains a forward selection method.
F-OUT should generally be (the same or) smaller than F-IN to stop the
method from getting stuck in a loop where the same variables are
repeatedly added in and taken out.  Since F-IN and F-OUT are fixed
numbers, the ``add'' and ``remove'' decisions are not $F$-tests with
fixed levels.

The stepwise procedure is probably preferable to just backwards
elimination or forward selection, as it allows one to enter and delete
variables. Due to correlations among the $X$ variables, it is quite
possible that the variable entered at the first step (i.e., the BEST
single variable) may be removed at some later step.

\comment{

Various programs offer different ways of proceeding.  The crucial
thing to decide is how to specify the enter and remove values.  One
approach is to specify the p-value (% point) of F for entering and
removing.  Thus at different stages the numerical value of F required
to enter or remove varies, depending on the degrees of freedom.  The
other approach (taken by BMDP, S, and Minitab) is to specify a fixed
value for entering and removing and forget about the significance of
the test.  As significance tests are probably meaningless here, the
second approach seems reasonable.

But what values should we actually use?  It {\em DEFINITELY} does make
a difference and is a matter of personal choice.  However one should
remember that these techniques are mainly used to find a handful of
plausible models which will be more carefully explored.

Some common recommendations are to take\\
F-IN $=$ 4.00 (which is about the 5\% level of F) \\
F-OUT $=$ 3.00 (which is about the 15\% level of F)


}

\paragraph{Stepwise regression with {\tt step()}} \

Splus has a function {\tt step} that will perform stepwise regression,
though it is based on minimizing AIC rather than using entry and
removal criteria with $F$ statistics.  When adding variables, {\tt
  step} chooses the variable that reduces AIC the most.  Similarly,
when removing variables, {\tt step} removes the variable that reduces
AIC the most.  It continues in adding and removing cycles until any
further move increases AIC.

{\tt step()} generates a lot of output, so I will abbreviate it here.
\begin{verbatim}
> step(rawlm,scope=list(upper= ~ ., lower ~ 1))  

Start: AIC= 70.0038 
 Rate ~ Len + Adt + Trk + Sp + lw + shl + itg + sigs + acpt + lan + type 

Step:  AIC= 67.2617 
 Rate ~ Len + Adt + Trk + Sp + lw + shl + itg + sigs + acpt + type 

Step:  AIC= 64.6104 
 Rate ~ Len + Adt + Trk + Sp + lw + shl + sigs + acpt + type 

Step:  AIC= 62.1651 
 Rate ~ Len + Adt + Trk + Sp + lw + shl + sigs + acpt 

Step:  AIC= 60.2145 
 Rate ~ Len + Trk + Sp + lw + shl + sigs + acpt 

Step:  AIC= 58.4081 
 Rate ~ Len + Trk + Sp + shl + sigs + acpt 

Step:  AIC= 56.5758 
 Rate ~ Len + Trk + Sp + sigs + acpt 

Step:  AIC= 55.6784 
 Rate ~ Len + Sp + sigs + acpt 
\end{verbatim}
In this case, it appears that we could never reduce AIC by adding
variables after some were already taken out.  [Note that AIC as
reported here is scaled yet another way.]

\comment{

OMITTED, SINCE I CAN'T EXPLAIN WHY STEP() DIDN'T ADD SIGS AT THE LAST
STEP BELOW


We can start with a small model in {\tt step()} and build up, e.g., by
saying
\begin{verbatim}
> small _ lm(Rate ~ 1,data=high)   
> as.formula(paste("~",paste(names(high),collapse="+"))
 ~ Rate + Len + Adt + Trk + Sp + lw + shl + itg + sigs + acpt + lan + type
> step(small, scope=
+ list(upper=as.formula(paste("~",paste(names(high),collapse="+"))),
+ lower=~1))        


Start:  AIC= 157.7748 
 Rate ~ 1 

Step:  AIC= 88.1151 
 Rate ~ acpt 

Step:  AIC= 80.2819 
 Rate ~ acpt + Sp 

Step:  AIC= 77.4802 
 Rate ~ acpt + Sp + Len 
\end{verbatim}
which shows that these heuristic algorithms don't converge to a unique
model.

END OF OMITTED TEXT.

}

The {\tt drop1()} and {\tt add1()} functions can be used to perform
individual steps of the {\tt step()} procedure, to see what happens
``by hand''.  For example,
\begin{verbatim}
> final  _ lm(formula = Rate ~ Len + Sp + sigs + acpt, data = high)
> as.formula(paste("~",paste(names(high),collapse="+"))
 ~ Rate + Len + Adt + Trk + Sp + lw + shl + itg + sigs + acpt + lan + type
> add1(final,scope=as.formula(paste("~",paste(names(high),collapse="+"))))
Single term additions

Model:
Rate ~ Len + Sp + sigs + acpt
       Df Sum of Sq      RSS       Cp 
<none>              41.95214 54.29100
   Adt  1  0.135480 41.81666 56.62330
   Trk  1  1.847865 40.10427 54.91091
    lw  1  0.319123 41.63302 56.43965
   shl  1  0.622958 41.32918 56.13582
   itg  1  0.351065 41.60107 56.40771
   lan  1  0.025793 41.92635 56.73298
  type  1  0.564620 41.38752 56.19416
\end{verbatim}
shows that there are no variables that can be added to the final model
to improve (lower) AIC (which is listed in the ``Cp'' column (!?) in
the output).

One nice feature of {\tt step()} is that it can also be used with
logistic regression and other {\tt glm()} objects.  The Venables and
Ripley MASS library ({\tt library(MASS)}) contains a function {\tt
  stepAIC()} which is supposed to be an improvement on {\tt step()},
especially for glm's where computation of AIC is not as
straightforward as in least-squares regression.  See {\tt
  help(stepAIC)} for details.


\paragraph{Stepwise regression with {\tt stepwise()}} \

The {\tt stepwise()} function that we used previously also has options
for forward, backwards, and stepwise regression.  The stepwise
method is Efroymson's, which adds variables when the F-statistic is
greater than F-IN (default: 2), and the criterion for
removing a variable is just heavy collinearity.  See {\tt
  help(stepwise)} for details.
{\small
\begin{verbatim}
> swe _ stepwise(X,y,method="efroymson")  
> sum.step(swe,high,"Rate",includeAIC=T,includeR2=T,sortcol=F)
        BIC   AIC    R2   R2a Len Adt Trk Sp lw shl itg sigs acpt lan type 
1(+9) -71.1 -69.4 0.517 0.504   0   0   0  0  0   0   0    0    1   0    0
2(+4) -68.2 -65.7 0.622 0.601   0   0   0  1  0   0   0    0    1   0    0
3(+1) -66.0 -62.6 0.694 0.667   1   0   0  1  0   0   0    0    1   0    0
4(+8) -66.1 -61.9 0.720 0.687   1   0   0  1  0   0   0    1    1   0    0
\end{verbatim}}
  Here we add first {\tt acpt}, then {\tt Sp}, then {\tt Len}, and
  finally {\tt sigs}; after this no more variables could be added
  according to the F-IN criterion, and at no stage could any variables
  be dropped based on collinearity.  Note that the final model (by
  AIC) is the same one achieved by {\tt step} above; the best of these
  by BIC omits {\tt sigs}.

\clearpage

\paragraph{Model summaries} \

Let's take a look at the parameter estimates and diagnostic plots of
the original model including all variables, as well as the two models
that were identified by our stepwise procedures.

{\small\begin{verbatim}
> summary(rawlm _ lm(Rate ~ ., data =high))$coefficients
                    Value Std. Error     t value    Pr(>|t|) 
(Intercept)  12.519017130 6.52897506  1.91745519 0.065816656
        Len  -0.056444936 0.03277522 -1.72218346 0.096474882
        Adt   0.009261974 0.03260618  0.28405580 0.778534668
        Trk  -0.127078208 0.10567014 -1.20259335 0.239574239
         Sp  -0.068470234 0.06365773 -1.07559977 0.291624232
         lw  -0.343297057 0.49669188 -0.69116704 0.495360764
        shl  -0.094711453 0.12451762 -0.76062688 0.453471734
        itg   0.286810535 1.14261457  0.25101250 0.803704939
       sigs   0.490779227 0.39409530  1.24533133 0.223709754
       acpt   0.068408278 0.02410784  2.83759426 0.008521754
        lan  -0.013096558 0.27430020 -0.04774535 0.962270618
       type   0.233826847 0.48934626  0.47783516 0.636613497
> plot(rawlm)
\end{verbatim}}\noindent

{\small\begin{verbatim}
> summary(mod1 _ lm(Rate ~ Len + Sp +sigs + acpt, data =high))$coefficients
                  Value Std. Error   t value     Pr(>|t|) 
(Intercept)  9.97563578 2.35441161  4.236997 0.0001633434
        Len -0.06188818 0.02535342 -2.441019 0.0200061400
         Sp -0.11573228 0.03881110 -2.981937 0.0052653982
       sigs  0.59434719 0.33121662  1.794436 0.0816397324
       acpt  0.06824634 0.02109316  3.235472 0.0027058274
> plot(mod1)
\end{verbatim}}\noindent

{\small\begin{verbatim}
> summary(mod2 _ lm(Rate ~ Len + Sp + acpt, data =high))$coefficients
                   Value Std. Error   t value      Pr(>|t|) 
(Intercept)  11.10415391 2.33971260  4.745948 0.00003447674
        Len  -0.07254915 0.02541716 -2.854338 0.00719988687
         Sp  -0.13122745 0.03901986 -3.363094 0.00187820573
       acpt   0.07591309 0.02130094  3.563837 0.00107925360
> plot(mod2)
\end{verbatim}}\noindent

Our task next time is to come up with a reasonable explanation for
what these variables mean.

\clearpage

\centerline{\psfig{figure=plot-lm1.ps,height=2.75in,width=5in}}

\hrule
\hrule

\centerline{\psfig{figure=plot-lm2.ps,height=2.75in,width=5in}}

\hrule
\hrule

\centerline{\psfig{figure=plot-lm3.ps,height=2.75in,width=5in}}

% {\small\begin{verbatim}
% > summary(mod3 _ lm(Rate ~ Trk + Sp + acpt, data =high))$coefficients
%                    Value Std. Error   t value      Pr(>|t|) 
% (Intercept)  11.63413664 2.39933275  4.848905 0.00002530127
%         Trk  -0.22984715 0.08491277 -2.706862 0.01043001186
%          Sp  -0.11930762 0.03979752 -2.997865 0.00497639882
%        acpt   0.07784186 0.02140489  3.636639 0.00088029306
% > plot(mod3)
% \end{verbatim}}\noindent
% 
% \centerline{\psfig{figure=plot-lm4.ps,height=3.5in}}



\section{Cross-validation and Prediction Error Sum of Squares (PRESS)}

We have seen earlier in the course (e.g. in the first course project)
that one way to evaluate a model is to hold out part of the data, do
model selection and fitting with what's left, and then calculate how
well the best model we can find does in predicting the hold-out data
(e.g. calculate the ``prediction'' RSS $\sum_{i\in\rm hold~out~data}
(y_i - \hat y_i)^2$, where the $\hat y_i$'s are the predictions of the
hold out data from the model fitted on the remainder.  We can even do
this several times with several different hold-out sets, to get a feel
for (a) how much the regression coefficients change; and (b) how much
the prediction RSS changes (both are kind of {\em bootstrapping}
ideas).

The most extreme version of this would be to hold one one observation
at a time, and predict it from the model fitted to the rest of the
data (this is an idea related to {\em jackknifing}).

By using the magic formula
\[
        (X_{(-i)}\tp X_{(-i)})^{-1} = (X\tp X)^{-1} +
        \frac{(X\tp X)^{-1} \mb x_i \mb x_i\tp  (X\tp X)^{-1}}
                {1 - h_{ii}}
\]
from our study of unusual or influential observations, we can
calculate 
\[
        E_{(-i)} = y_i - \mb x_i^T\hat{\mb\beta}_{(-i)} = 
        \frac{E_i}{1-h_{ii}}
\]
where $E_i$ and $h_{ii}$ are the residuals and diagonal entries in the
hat matrix of the full regression, and $\hat{\mb\beta}_{(-i)}$ os
calculated from the model without observation $i$, $(y_i,\mb x_i^T)$.
This suggests a criterion like
\[
        PRESS = \sum_{i=1}^n E_{(-i)}^2 
        = \sum_{i=1}^n \left[\frac{E_i}{1-h_{ii}}\right]^2
\]
PRESS has the aesthetic advantage of measuring how well $y_i$
is predicted from all the data {\em except} the $i^{th}$ observation,
for all $i$.

Clearly, different models can be compared based on the size of their
PRESS, so that PRESS becomes another way to ``score'' models in model
selection.  Often $C_p$ and AIC choose models that minimze PRESS
(since they are basically minimizing penalized RSS, which is like
PRESS, to begin with).

Outside of ordinary linear regression, where the ``magic formula''
above no longer holds, PRESS can be expensive to calculate.
Therefore, one often finds a few candidate models using a cheaper
criterion like $C_p$, AIC, or BIC, and then PRESS is used for final
model selection (if that is the goal).

\clearpage

\section{Other Ways of Coping with Collinearity: Biased Regression Methods}

Recall that the principal problem with collinearity in the regression
problem 
\begin{equation}
        \mb y = X \mb\beta + \mb \epsilon
        \label{ols-problem}
\end{equation}
is the technical difficulty that $(X\tp X)^{-1}$ becomes difficult to
compute.  If we still wish to work with this regression problem and
not some other one, then we have to modify our estimation procedure
somehow to obtain, for example, more stable estimates
$\hat{\mb\beta}$.  

\subsection{Principal Components}

The least-squares estimate of $\mb\beta$ is
\[
        \hat{\mb\beta} = (X\tp X)^{-1} X\tp \mb y
\]
but if $X$ has nearly collinear columns this is difficult to compute.
To make the problem easier to think about let us imagine that $X$ has
a column of 1's, to fit an intercept, and let us write
\[
        X = [~\mb 1 ~~ \mb x_1 ~~ \cdots ~~ \mb x_k~]
\]
where $\mb 1$ is the intercept column and the columns $\mb
x_1,\ldots,\mb x_k$ are the independent variables of interest.  We define
\[
        Z = [~\mb z_1 ~~ \cdots ~~ \mb z_k~] 
        = [~(\mb x_1 - \bar x_1)/s_1 ~~ \cdots ~~ (\mb x_k - \bar x_k)/s_k~]
\]
where $\bar x_j$ is the $j^{th}$ column mean and $s_j^2$ is the sample
variance of the $j^{th}$ column.  It is fairly easy to show that 
\begin{equation}
        \ul y = \mb 1 \beta_0 + Z\mb\alpha + \mb\epsilon
        \label{normalized}
\end{equation}
is the same regression problem as Equation \eq{ols-problem}, and 
\begin{eqnarray*}
        \hat{\mb\alpha} & = & (Z\tp Z) Z\tp \ul y \\
        \Var(\hat{\mb\alpha}) & = & \sigma^2_{\epsilon} (Z\tp Z)^{-1} 
\end{eqnarray*}
There will be collinearity in $X$ if and only if there is collinearity
in $Z$, and $Z$ has the pleasant property that
\[
        R_{XX} = \frac{1}{n-1} \cdot Z\tp Z
\]
is the matrix of sample correlations of the $\mb x_j$'s (with 1's down
the main diagonal).  

For any $\mb u$, $\mb u\tp Z\tp Z \mb u = \Vert Z\mb u\Vert^2\geq0$ so
$Z\tp Z$ is symmetric and positive-definite; we can then write the
eigenvalue-eigenvector decomposition
\[
        Z\tp Z = V L V\tp
\]
where $dim(V)=dim(Z) = n \times k$ and the columns of $V$ are {\em
orthonormal}, i.e. $V\tp V=I$, and
\[
        L = \left[\begin{array}{cccc}
        \lambda_1 & 0 & \cdots & 0 \\
        0 & \lambda_2 & \cdots & 0 \\
        \vdots & \vdots &\ddots & \vdots \\
        0 & 0 & \cdots & \lambda_k
        \end{array}\right]
\]
Because the columns of $V$ are orthonormal, it follows that
\[
        Z\tp Z V = VLV\tp V = VL
\]
so that the $j^{th}$ column of $V$ is the $j^{th}$ {\em eigenvector}
of $Z$, with associated {\em eigenvalue} $\lambda_j$; we can arrange
that $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_p$.

The {\em principal components} of $Z$ are the columns of the matrix
\[
        W = Z V
\]
Note that Equation \eq{normalized} can be written
\begin{eqnarray}
        \ul y & = & \mb 1 \beta_0 + Z\mb\alpha + \mb\epsilon 
        \nonumber \\
        & = & \mb 1 \beta_0 + (Z V) (V\tp \mb\alpha) + \mb\epsilon 
        \nonumber \\
        & = & \mb 1 \beta_0 + W\mb\gamma + \mb\epsilon
        \label{pc-problem}
\end{eqnarray}
and this is equivalent to our original Equation \eq{ols-problem}.

Finally, note that 
\[
        \Var(\hat{\mb\gamma}) = \sigma^2_{\epsilon} (W\tp W)^{-1}
        = \sigma^2_{\epsilon} (V\tp Z\tp Z V)^{-1} 
        = \sigma^2_{\epsilon} L^{-1} / (n-1)
\]
is a diagonal matrix, and $\Var(\hat\gamma_k) = \sigma^2\lambda_k^{-1}/(n-1)$,
so that (1) the $\hat\gamma$'s are uncorrelated; and (2) the data
tells us the least about the $\gamma_k$ for the smallest eigenvalue
$\lambda_k$.

This suggests that we can ``improve'' the model by dropping columns of
$W$ corresponding to high $\Var(\hat\gamma_k)$ i.e. low $\lambda_j$.
If the original model in Equation \eq{ols-problem} was the correct
model for the data, then the estimates $\hat{\mb\beta}$,
$\hat{\mb\alpha}$ and $\hat{\mb\gamma}$ in these problems should all
be unbiased.

To see what the effect of dropping columns of $W$ is, consider
minimizing the sum of mean square errors (SMSE) in the coefficients,
instead of in the predicted values,
\begin{eqnarray}
        SMSE & = & \sum_{j=1}^k E[(\hat\gamma_j - \gamma_j)^2] 
                \nonumber \\
        & = & \sum_{j=1}^k \{ \Var(\hat\gamma_j) + [bias(\hat\gamma_j)]^2\} 
                \nonumber \\
        & = & \frac{\sigma^2_{\epsilon}}{n-1} \sum_{j=1}^k \lambda_j^{-1}
                + \sum_{j=1}^k \{E[\hat\gamma_i - \gamma_j]\}^2
                \label{bias-variance}
\end{eqnarray}
where $\hat\gamma_j$ are the least-squares estimates.  When the full
set of least-squares estimates are used, $\hat\gamma_j$ will be
unbiased and the ``bias'' term on the right will be zero, so 
\[
        SMSE_{OLS} = \frac{\sigma^2_{\epsilon}}{n-1} \sum_{j=1}^k
        \lambda_j^{-1} 
\]
Because the columns of $W$ are orthogonal, when we drop columns we do
not change the least-squares estimates of the remaining coefficients.
Dropping the last $m$ columns for example we obtain
\[
        SMSE_{PC(k-m)} = \frac{\sigma^2_{\epsilon}}{n-1}
                \sum_{j=1}^{k-m}\lambda_j^{-1} 
                + \sum_{j=k-m+1}^k \gamma_j
\]
This is a classic case of the {\em bias/variance tradeoff} in
statistics: by dropping only a few terms we can reduce SMSE by
paying for a large drop in variance with a small increase in bias;
clearly we can minimize $SMSE_{PC(k-m)}$ over $m$ if desired.


\subsection{Ridge Regression}

Ridge regression is another way to play the bias-variance tradeoff
game.  The idea is that $Z\tp Z$ may have a numerically unstable
inverse, but 
\[
        (Z\tp Z + kI)^{-1}
\]
will be quite stable, for $k>0$.  The ridge regression estimator
for $\mb\beta$ is then
$$
        \hat{\mb\alpha}^{(k)} = (Z\tp Z + kI)^{-1}Z\tp \ul y \eqno{(*)}
$$
Again,
\[
        SMSE_{RR(k)} = \sum_{j=1}^k \{ \Var(\hat\alpha^{(k)}_j) 
                + [bias(\hat\alpha^{(k)}_j)]^2\} 
\]
and when $k=0$ this reduces to the least-squares expression, with no
bias term.  As $k$ grows (not necc. integer values!), the variance
term decreases since it is determined largely by $(Z\tp Z + kI)^{-1}$;
and of course the bias term is no longer zero.  We can if we wish
estimate the optimal $k$ by directly minimizing $SMSE_{RR(k)}$, or by
employing a criterion like PRESS, or Mallows' $C_p$ (where $p={\rm
  trace}\,(H_{RR(k)})$, the trace of the ridge regression hat matrix);
see the bullets below for a bit more information.

If we work with the coefficients $\mb\gamma$ in the principal components
regression \eq{pc-problem}, we can show that 
\[
        \hat\gamma_j^{(k)} = \frac{\lambda_j}{\lambda_j + k} \hat\gamma_j
\]
so we see that as $k$ grows the estimates of the $\gamma_j$'s are
being ``shrunk to zero''.  (This corresponds to a certain Bayesian
formulation of linear regression which we may be able to look at
briefly later in the semester).

\noindent
{\bf Notes:}

\begin{itemize}

\item In the formulation of ridge regression above we worked with the
  normalized matrix of predictors $Z$.  In fact there is no particular
  reason to use this normalized version, and we might prefer to use
\begin{itemize}
\item The matrix $X^*$ of regressors variables after subtracting the
  column means; or
\item The original, un-normalized matrix $X$ of regressors.
\end{itemize}
For example if the $X$ matrix has 0/1 dummy variables to code for
levels of some discrete factor, the normalized version is more
difficult to interpret and we'd better stick to the un-normalized $X$.
In that case the ridge regression estimates are just as in $(*)$, but
using $X$ in place of $Z$.

\item How do we choose the value of the ridge constant $k$?  
\begin{itemize}
\item
  One possibility is to plot the parameter estimates $\hat\beta_j$ as
  a function of $k$: in general $|\hat\beta_j|$ decreases with $k$.
  Then choose the largest value of $k$ that doesn't affect the values
  of the $\hat\beta_j$'s very much.  Such plots are called {\em ridge
    trace plots}.
\item Another possibility is to plot a measure of overall ``goodness''
  of the model, such as PRESS, as a function of $k$.  Then choose the
  value of $k$ that minimizes PRESS in the plot.  This is an example
  of choosing $k$ {\em by cross-validation}.
\end{itemize}

\end{itemize}

\end{document}









============================




\section{From Cross-Validation (Fox 16.2) to PRESS to Jackknifing}




\section{Detecting and Fixing Collinearity (Fox 13.1, 13.2)}


A seminal reference for this material is the book {\em Regression
Diagnostics} by Belsley, Kuh, and Welch (1980).


------------------------------------------------------------------

% \section{Collinearity (again) and Model Selection (at last)}



\section{Half Sample Analysis}
\label{sec:half}

There are 39 cases in the highway data.  Splitting that in half (to
only 20 cases) is a little extreme, but worth doing for an exercise.
It might be a better idea to look randomly select, say, 30
observations to analyze and the remaining 9 to verify against.

Anyway(s), here goes with a half sample.
{\small\begin{verbatim}
> ss <-sample(39,size=19)
> high1 <- high[ss,]     
> high2 <- high[-ss,]
\end{verbatim}}\noindent

The data {\tt high1} contains a sample of size 19 from the original
data, while the data {\tt high2} contains the other 20 observations.
Lets see what happens if we let the step function loose on these two
data sets.

{\small\begin{verbatim}
> high1.lm <- lm(formula(high.lm), data=high1)
> high1.lmb <- lm(Rate ~ 1, data=high1)
\end{verbatim}}\noindent

The simple forward stepwise gives a lot of output, but ends up with
{\small\begin{verbatim}
step(high1.lmb, formula(high1.lm), direction="both")
 .......
Model:
Rate ~ acpt

scale:  2.894905 

       Df Sum of Sq      RSS       Cp 
<none>              23.18474 34.76436
   Len  1  4.995871 18.18887 35.55830
   Adt  1  3.555030 19.62971 36.99914
   Trk  1  2.469646 20.71509 38.08452
    Sp  1  4.272275 18.91246 36.28189
    lw  1  1.114846 22.06989 39.43932
   shl  1  0.260603 22.92414 40.29357
   itg  1  1.292816 21.89192 39.26135
  sigs  1  4.960888 18.22385 35.59328
   lan  1  3.731358 19.45338 36.82281
  type  3  3.280835 19.90390 48.85296
Call:
lm(formula = Rate ~ acpt, data = high1)

Coefficients:
 (Intercept)      acpt 
    1.491071 0.1829118

Degrees of freedom: 19 total; 17 residual
Residual standard error (on weighted scale): 1.167822 
\end{verbatim}}\noindent

Going for backward elimination yields,
{\small\begin{verbatim}
step(high1.lm)

Single term deletions

Model:
Rate ~ Sp + acpt + lan

scale:  1.227786 

       Df Sum of Sq      RSS       Cp 
<none>              13.98973 23.81202
    Sp  1  5.463647 19.45338 26.82009
  acpt  1  9.404156 23.39389 30.76060
   lan  1  4.922730 18.91246 26.27918
Call:
lm(formula = Rate ~ Sp + acpt + lan, data = high1)

Coefficients:
 (Intercept)         Sp     acpt       lan 
    7.844092 -0.1326009 0.140156 0.5445826

Degrees of freedom: 19 total; 15 residual
Residual standard error (on weighted scale): 0.9657375 
\end{verbatim}}\noindent

This looks a little familiar, but not exactly anything we have seen
before.  What about the other half sample?
We again set up the {\tt high2} models and go with the fits.
Forward stepping leads to
{\small\begin{verbatim}
lm(formula = Rate ~ acpt + Trk, data = high2)
\end{verbatim}}\noindent
while backward elimination leads to
{\small\begin{verbatim}
lm(formula = Rate ~ Adt + itg + sigs + acpt, data = high2)
\end{verbatim}}\noindent

We seem to be getting deeper and deeper into trouble.

We might try predicting the ``other'' half from the best model for one
half.  Here is a part of this.
{\small\begin{verbatim}
> del1 <- predict(lm(formula = Rate ~ Sp + acpt + lan, data = high1), high2)
> del2 <- predict(lm(formula = Rate ~ Adt + itg + sigs + acpt, data = high2)
+ )
> plot(del2,del1)
> var( high2$Rate - del1)     
[1] 2.913679
> var( high2$Rate - del2)
[1] 1.491909
\end{verbatim}}\noindent

(The plot does not look very linear!).

Perhaps the half-sample was a little too extreme in this case.
Lets try a subsample of size 30 out of 39.
{\small\begin{verbatim}
> ss <- sample(39,30)
> high1 <- high[ss,] 
> high2 <- high[-ss,]
\end{verbatim}}\noindent

Fit the models to high1, but assess the fit based on high2.
To do this really properly would require some substantial programming.

Stepping down from the full model yields, 
{\small\begin{verbatim}
> step(high1.lm)

Call:
lm(formula = Rate ~ Len + lw + shl + acpt, data = high1)

Coefficients:
 (Intercept)       Len        lw        shl      acpt 
    27.43066 -0.131288 -1.840498 -0.1691469 0.1146749
\end{verbatim}}\noindent

So now I'll write a short PRESS function, and try dropping terms from
this model, just to see what happens.
{\small\begin{verbatim}
> mypress
function(pred, model)
{
        var(pred$Rate - predict(model, pred))
}
> 
> ####                                                                 
> mypress(high2, lm(formula = Rate ~ Len + lw + shl + acpt, data = high1))
[1] 4.436016
> mypress(high2, lm(formula = Rate ~ Len + lw + shl, data = high1))       
[1] 6.619988
> mypress(high2, lm(formula = Rate ~ Len + lw + acpt, data = high1))      
[1] 4.97908
> mypress(high2, lm(formula = Rate ~ Len + shl + acpt, data = high1))     
[1] 2.525569
> mypress(high2, lm(formula = Rate ~ lw + shl + acpt, data = high1))      
[1] 4.217933

\end{verbatim}}\noindent

Notice that the model where I dropped {\tt lw} is markedly better than
the base model, and the one where I dropped Len is slightly better.
This isn't exactly what one would have expected from the summary
{\small\begin{verbatim}
> summary(lm(formula = Rate ~ Len + lw + shl + acpt, data = high1))

Call: lm(formula = Rate ~ Len + lw + shl + acpt, data = high1)
Residuals:
    Min     1Q  Median    3Q   Max 
 -1.621 -0.562 -0.1297 0.686 2.103

Coefficients:
               Value Std. Error  t value Pr(>|t|) 
(Intercept)  27.4307   8.1390     3.3703   0.0024
        Len  -0.1313   0.0372    -3.5285   0.0016
         lw  -1.8405   0.6428    -2.8634   0.0084
        shl  -0.1691   0.0792    -2.1362   0.0426
       acpt   0.1147   0.0240     4.7691   0.0001
\end{verbatim}}\noindent

END DISCUSSION.  This must be very confusing.

{\small\begin{verbatim}
leaps(cbind(Len,Adt,Trk,Sp,lw,shl,itg,sigs,acpt,lan,type),Rate,nbest=2)
\end{verbatim}}\noindent

\end{document}



