\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{\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.

\section{Variable Selection (Fox 13.2) and Specification Error (Fox 9.6, 6.3)}

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$ variables, as
with the college admissions data, 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 the
college admissions data, 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.

To see what can go wrong (or right), consider the following argument
of Fox (Sect 9.3): We are interested in what happens to the
correctly-specified regression
\[
        \mb y = \mb 1 \beta_0 + X_1 \mb \beta_1 + X_2 \mb \beta_2 +
        \mb \epsilon 
\]
when the columns in $X_2$ are omitted.  To simplify the notation 
we consider the equivalent problem 
\[
        \mb y^* = \mb X_1^* \mb\beta_1 + X_2^* \mb\beta_2 + \mb \epsilon
\]
where $\mb y^* = \mb y - \mb 1 \beta_0$ and $X_1^*$ and $X_2^*$ are
obtained from $X_1$ and $X_2$ by subtracting column means.

Now imagine that we forget to include the columns $X_2^*$ in the
regression; so we are instead working with 
\[
        \mb y^* = X_1^*\mb\beta_1 + \widetilde{\mb\epsilon}
        \mbox{~~~ where ~} \widetilde{\mb\epsilon} = X_2^* \mb\beta_2
        + \mb \epsilon 
\]
The least squares estimator for $\mb\beta_1$ in this mis-specified
model is 
\begin{eqnarray*}
        \hat{\mb\beta}_1 & = & ({X_1^*}\tp X_1^*)^{-1} {X_1^*}\tp\mb y^*\\
        & = & (\frac1n {X_1^*}\tp X_1^*)^{-1} \frac1n {X_1^*}\tp\mb y^*\\
        & = & (\frac1n {X_1^*}\tp X_1^*)^{-1} \frac1n {X_1^*}\tp
                \mb X_1^* \mb\beta_1 + X_2^* \mb\beta_2 + \mb\epsilon\\
        & = & \mb\beta_1 + (\frac1n {X_1^*}\tp X_1^*)^{-1} \frac1n {X_1^*}\tp
                X_2^* \mb\beta_2
                + (\frac1n {X_1^*}\tp X_1^*)^{-1} \frac1n {X_1^*}\tp
                \mb\epsilon
\end{eqnarray*}
Under suitable regularity conditions (basically, joint random sampling
of $(y_i, x_{i1}, \ldots, x_{ik})$, it follows that, as $n\longrightarrow\infty$,
\begin{eqnarray*}
        \hat{\mb\beta}_1 & \stackrel{p}{\longrightarrow} &
                \mb\beta_1 + \vc_{11}^{-1} \vc_{12} \mb\beta_2
                + \vc^{-1}_{11}\mb\sigma_{1\mb\epsilon} \\
        & = & \mb\beta_1 + \vc_{11}^{-1} \vc_{12} \mb\beta_2
\end{eqnarray*}
where $\vc{ij}$ is the population variance-covariance matrix for $X_i$
and $X_j$, and $\mb\sigma_{1\mb\epsilon}$ is the vector of covariances
between $X_1$ and $\mb\epsilon$ (which must be zero, since
$\mb\epsilon$ is assumed to be independent of the $X$'s).  Moreover it
is also easy to see that 
\[
        \frac1n {X_1^*}\tp \widetilde{\mb\epsilon} \stackrel{p}{\longrightarrow}
        \vc_{12}\mb\beta_2
\]
We conclude that
\begin{itemize}
\item If $\vc_{12} = 0$, i.e. if $X_2$ is orthogonal to $X_1$,
 then the estimates $\hat{\mb\beta}_1$ will be
consistent for $\mb\beta_1$ as $n$ grows. 
\item If $\mb\beta_2=\mb 0$, i.e. if $X_2$ is irrelevant to the
problem, then the estimates $\hat{\mb\beta}_1$ will also be consistent
for $\mb\beta_1$ as $n$ grows.
\item Otherwise the estimates of $\mb\beta_1$ are biased, even if the
sample size grows to $\infty$.
\end{itemize}

How important is this bias?  If we attach causal meaning to the
regression, it depends on how $X_1$ and $X_2$ are related
``causally''.  Suppose we are interested in the effect of $X_1$ on
$Y$, as measured by $\mb\beta$.  
\begin{itemize}
\item In Figure~\ref{diag1}, the ``bias'' in estimating $\mb\beta_1$
without $X_2$ reflects how much regressing on $X_2$ first reduces the
effect of $X_1$ on $Y$ measured by $\mb\beta_1$.  In some applications
of regression to psychology and social science, finding such an $X_2$
that 
\begin{itemize}
\item is correlated with $Y$;
\item reduces the effect of $X_1$ on $Y$ when included in the
regression
\end{itemize}
is considered to be suggestive that $X_2$ is an intermediate step in
(at least one) causal chain leading from $X_1$ to $Y$.  Such an $X_2$
is called a {\em mediating variable} in the quant. psych. literature,
for example.

\begin{figure}
\begin{center}
\input diag1.tex
\end{center}
\caption{Assessing the causal effect of $X_1$ on $Y$: the bias in
estimating $\mb\beta_1$ without $X_2$ in the regression {\tt lm(Y \~{}
X$_1$)} does not affect a causal interpretation of the arrow from $X_1$
to $Y$, though it may suggest that $X_2$ is an intermediate variable
in (one) causal chain from $X_1$ to $Y$.}\label{diag1}
\end{figure}

\item In Figure~\ref{diag2}, $X_2$ is a common cause of $X_1$ and $Y$,
and the presence or absence of the arrow leading from $X_1$ to $Y$ may
be entirely due to $X_2$ (analogous to ``Simpson's paradox''); if the
``true picture'' does not have the arrow from $X_1$ to $Y$, but we
only regress on $X_1$, the common cause $X_2$ will induce a spurious
correlation between $X_1$ and $Y$ that the regression will pick up;
indeed this correlation is exactly the bias in estimating $\mb\beta_1$
without including $Y$ in the model.

\begin{figure}
\begin{center}
\input diag2.tex
\end{center}
\caption{Assessing the causal effect of $X_1$ on $Y$: The association
between $X_1$ and $Y$ picked up by {\tt lm(Y \~{} X$_1$)}  may be
``real'' or it may be ``spurious'', i.e. just caused by the bias in
estimating $\mb\beta_1$ without including $X_2$ in the
regression.}\label{diag2} 
\end{figure}

\item The real problem though is what happens if we have only
collected data on $X_1$, $X_2$ and $Y$ and we are looking at
Figure~\ref{diag2} but the real situtation is that we completely
forgot to measure $Z$ in Figure~\ref{diag3}.  Even if we condition on
$X_2$ (i.e. include $X_2$ in a regression model that is intended to
measure the effects of $X_1$) we may or may not be able to deduce that
$X_1$ and $Y$ are unrelated given $X_2$.

\begin{figure}
\begin{center}
\input diag4.tex
\end{center}
\caption{Assessing the causal effect of $X_1$ on $Y$: If $Z$ is
omitted, the marginal distribution of $(X_1,X_2,Y)$ will have
associations like that of Figure~\ref{diag2}.  Inclusing $X_2$ in the
regression {\tt lm(Y \~{} X$_1$ + X$_2$)} may or may not reveal the
true absence of an arrow from $X_1$ to $Y$, because of bias in
estimating $\mb\beta_1$ induced by not including $Z$ in the
regression.}\label{diag3}
\end{figure}

\end{itemize}

\section{Automatic Variable Selection (Fox 13.2.2)}

In this section we will look at two or three automatic, heuristic
methods for variable selection and compare them with what we would do
``by hand'', even if we didn't know much about the scientific context
of the data.  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
parameter due to omitting relevant regressors from the model.

\subsection{An Example}

The data in the Table \ref{tab:highway} 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.
\begin{table}\small
\hfil\begin{tabular}{rrrrrrrrrrrrrrr} \hline
 & $y$ & $x_1$ & $x_2$ & $x_3$ & $x_4$ & $x_5$ & $x_6$ & $x_7$ & $x_8$ & $x_9$ & $x_{10}$ & $x_{11}$
 & $x_{12}$ & $x_{13}$\\ \hline
1 &  4.58 & 4.99  & 69 &  8 & 55 & 12 & 10 & 1.20 & 0 & 4.60 & 8 & 1 & 0 & 0\\
2 &  2.86 & 16.11 & 73 &  8 & 60 & 12 & 10 & 1.43 & 0 & 4.40 & 4 & 1 & 0 & 0\\
3 &  3.02 & 9.75  & 49 & 10 & 60 & 12 & 10 & 1.54 & 0 & 4.70 & 4 & 1 & 0 & 0\\
4 &  2.29 & 10.65 & 61 & 13 & 65 & 12 & 10 & 0.94 & 0 & 3.80 & 6 & 1 & 0 & 0\\
5 &  1.61 & 20.01 & 28 & 12 & 70 & 12 & 10 & 0.65 & 0 & 2.20 & 4 & 1 & 0 & 0\\
6 &  6.87 & 5.97  & 30 &  6 & 55 & 12 & 10 & 0.34 & 1.84 & 24.80 & 4 & 0 & 1 & 0\\
7 &  3.85 & 8.57  & 46 &  8 & 55 & 12 &  8 & 0.47 & 0.70 & 11.00 & 4 & 0 & 1 & 0\\
6 &  6.12 & 5.24  & 25 &  9 & 55 & 12 & 10 & 0.38 & 0.38 & 48.50 & 4 & 0 & 1 & 0\\
9 &  3.29 & 15.79 & 43 & 12 & 50 & 12 &  4 & 0.95 & 1.39 & 7.50 & 4 & 0 & 1 & 0\\
10 & 5.88 & 8.26  & 23 &  7 & 50 & 12 &  5 & 0.12 & 1.21 & 8.20 & 4 & 0 & 1 & 0\\
11 & 4.20 & 7.03  & 23 &  6 & 60 & 12 & 10 & 0.29 & 1.85 & 5.40 & 4 & 0 & 1 & 0\\
12 & 4.61 & 13.28 & 20 &  9 & 50 & 12 &  2 & 0.15 & 1.21 & 11.20 & 4 & 0 & 1 & 0\\
13 & 4.80 & 5.40  & 18 & 14 & 50 & 12 &  8 & 0    & 0.56 & 15.20 & 2 & 0 & 1 & 0\\
14 & 3.85 & 2.96  & 21 &  8 & 60 & 12 & 10 & 0.34 & 0 & 5.40 & 4 & 0 & 1 & 0\\
15 & 2.69 & 11.75 & 27 &  7 & 55 & 12 & 10 & 0.26 & 0.60 & 7.90 & 4 & 0 & 1 & 0\\
16 & 1.99 & 8.86  & 22 &  9 & 60 & 12 & 10 & 0.68 & 0 & 3.20 & 4 & 0 & 1 & 0\\
17 & 2.01 & 9.78  & 19 &  9 & 60 & 12 & 10 & 0.20 & 0.10 & 11.00 & 4 & 0 & 1 & 0\\
18 & 4.22 & 5.49  &  9 & 11 & 50 & 12 &  6 & 0.18 & 0.18 & 8.90 & 2 & 0 & 1 & 0\\
19 & 2.76 & 8.63  & 12 &  8 & 55 & 13 &  6 & 0.14 & 0 & 12.40 & 2 & 0 & 1 & 0\\
20 & 2.55 & 20.31 & 12 &  7 & 60 & 12 & 10 & 0.05 & 0.99 & 7.80 & 4 & 0 & 1 & 0\\
21 & 1.89 & 40.09 & 15 & 13 & 55 & 12 &  8 & 0.05 & 0.12 & 9.60 & 4 & 0 & 1 & 0\\
22 & 2.34 & 11.81 &  8 &  8 & 60 & 12 & 10 &    0 & 0 & 4.30 & 2 & 0 & 1 & 0\\
23 & 2.83 & 11.39 &  5 &  9 & 50 & 12 &  8 &    0 & 0.09 & 11.10 & 2 & 0 & 1 & 0\\
24 & 1.81 & 22.00 &  5 & 15 & 60 & 12 &  7 &    0 & 0 & 06.80 & 2 & 0 & 1 & 0\\
25 & 9.23 & 3.58  & 23 &  6 & 40 & 12 &  2 & 0.56 & 2.51 & 53.00 & 4 & 0 & 0 & 1\\
26 & 8.60 & 3.23  & 13 &  6 & 45 & 12 &  2 & 0.31 & 0.93 & 17.30 & 2 & 0 & 0 & 1\\
27 & 8.21 & 7.73  &  7 &  8 & 55 & 12 &  8 & 0.13 & 0.52 & 27.30 & 2 & 0 & 0 & 1\\
28 & 2.93 & 14.41 & 10 & 10 & 55 & 12 &  6 &    0 & 0.07 & 18.00 & 2 & 0 & 0 & 1\\
29 & 7.48 & 11.54 & 12 &  7 & 45 & 12 &  3 & 0.09 & 0.09 & 30.20 & 2 & 0 & 0 & 1\\
30 & 2.57 & 11.10 &  9 &  8 & 60 & 12 &  7 &    0 & 0 & 10.30 & 2 & 0 & 0 & 1\\
31 & 5.77 & 22.09 &  4 &  8 & 45 & 11 &  3 &    0 & 0.14 & 18.20 & 2 & 0 & 0 & 1\\
32 & 2.90 & 9.39  &  5 & 10 & 55 & 13 &  1 &    0 & 0 & 12.30 & 2 & 0 & 0 & 1\\
33 & 2.97 & 19.49 &  4 & 13 & 55 & 12 &  4 &    0 & 0 & 7.10 & 2 & 0 & 0 & 1\\
34 & 1.84 & 21.01 &  5 & 12 & 55 & 10 &  8 &    0 & 0.10 & 14.00 & 2 & 0 & 0 & 1\\
35 & 3.78 & 27.16 &  2 & 10 & 55 & 12 &  3 & 0.04 & 0.04 & 11.30 & 2 & 0 & 0 & 1\\
36 & 2.76 & 14.03 &  3 &  8 & 50 & 12 &  4 & 0.07 & 0 & 16.30 & 2 & 0 & 0 & 1\\
37 & 4.27 & 20.63 &  1 & 11 & 55 & 11 &  4 &    0 & 0 & 9.60 & 2 & 0 & 0 & 1\\
38 & 3.05 & 20.06 &  3 & 11 & 60 & 12 &  8 &    0 & 0 & 9.00 & 2 & 0 & 0 & 0\\
39 & 4.12 & 12.91 &  1 & 10 & 55 & 12 &  3 &    0 & 0 & 10.40 & 2 & 0 & 0 & 0\\
\end{tabular}
\caption{Highway Data \label{tab:highway}}
\end{table}
The variables in the table are:

\hfil\begin{tabular}{rrl} 
$y$ & RATE & 1973 accident rate per million vehicle miles \\
$x_1$ & LEN & length of segment in miles \\
$x_2$ & ADT & average daily traffic couunt in thousands (estimated) \\
$x_3$ & TRKS & truck volume as a percent of total volume \\
$x_4$ & SLIM & speed limit (in 1973, before the 55 mph limits) \\
$x_5$ & LWID & lane width in feet \\
$x_6$ & SHLD & width in feet of outer shoulder on the roadway \\
$x_7$ & ITG & number of freeway-type interchanges per mile in the 
segment \\
$x_8$ & SIGS & number of signalized interchanges per mile in the 
segment \\
$x_9$ & ACPT & number of access points per mile in the segment \\
$x_{10}$ & LANE & total number of lanes of traffic in both directions \\
$x_{11}$ & FAI & 1 if Federal aid interstate highway, 0 otherwise \\
$x_{12}$ & FAI & 1 if principal arterial interstate highway, 0 otherwise \\
$x_{13}$ & FAI & 1 if minor arteial interstate highway, 0 otherwise \\
\end{tabular}

Notice that columns 11--13 are dummy variables which are different aspects of
one factor.  We would like to choose a model for $\mb{y}$ from among the
$\mb{x}$ variables (columns of $X$).

Let $X_{39 \times 14}$ be the full design matrix (i.e the table of data
plus a column of ones).  We want a subset of the columns to select as a model.
One way of thinking about this is that we want a matrix 
$X'_{39 \times p}$ which includes all the ``important'' variables in
the problem.  The aim is to make $k$ small---so that the resulting matrix
can be easily interpreted.  Usually we sould want to column of 1's included
in any model, but that still leaves us with $2^{13} = 8192$ possible models.
How does one proceed?

One approach would be to fit all the models.  This is feasible, but it is
not a method we would try with (say) 100 variables.  We need a more organized
approach to the problem.

The first thing we should consider is what subsets belong together and which
ones are distinct and can be split.  In the example the dummy variables from
a logical group---they index the type of road and it would seem reasonable
that either all or none of this set should be included.  By grouping the
variables into sets that we wish to consider together we can reduce the
number of models whcih need to be fit.  There are limits to this 
approach---there are still many hundreds of models which we can consider.

There are two areas of concern:
\begin{enumerate}
\item How do we compare different subsets, i.e., what is a ``measure'' of a
good-fitting model.
\item Given that we didn't want to fit all possible models, how do we find a
procedure for selecting models.
\end{enumerate}

We will first consider the problem of comparing different models.  Some
methods for this are:
\begin{itemize}
\item $\mathbf R^2$\\
$R^2$  can be (and is) used to compare models with the {\em same
number} of parameters.  Thus for the highway data we could compare all
models which include ($x_{11}$--$x_{13}$) and (say) one extra variable by
using the appropriate $R^2$ values.  The model with the largest value of
$R^2$ would be considered as (a candidate for) the ``best''
\item {\bf 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.
\end{itemize}

Neither of the above methods allows us to compare models which are not
nested and have a different number of predictors, e.g., 
$(x_1,x_3,x_5)$ and $(x_2,x_3,x_4,x_7)$

\begin{itemize}
\item {\bf Adjusted $\mathbf R^2$} \\
\[
  \bar{R}^2_p = 1 - \left( \frac{n-1}{n-p} \right ) 
                \left ( 1 - R^2_p \right )
\]
where $p$ is the number of parameters included and $R^2_p$ is the usual $R^2$.
This is an adjustment which is designed to compensate for the fact that
otherwise as $p \nearrow $ then $R^2_p \nearrow$.  Note that there are some
problems with this in that $\bar{R}^2_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 $\bar{R}^2_p$.
\item {\bf Mallows $\mathbf C_p$ (recommended among these methods)} \\
\begin{equation}
  C_p = \frac{RSS_p}{\sigsq} + 2p - n
  \label{mallows}
\end{equation}
% \[
%   C_p = \frac{RSS_p}{\sigsq} + 2p - n
% \]
where $RSS_p$ is computed from the ``current'' model with $p$
parameters, 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.  For the ``correct'' model (i.e., a
model which fits but does not include too many terms), it is possible
to show that $E[C_p] \approx p $.  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$.
%Figure~\ref{fig:cp} gives an example
%of a $C_p$ plot.  The solid line is $C_p = p $. 
The full model always has
$C_{k'} = k'$, 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$.
\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}





\subsection{A Derivation of $C_p$}

% Recall that we had the ``magic'' formula for $C_p$ as
% \begin{equation}
%   C_p = \frac{RSS_p}{\sigsq} + 2p - n
%   \label{mallows}
% \end{equation}
% but that we did not give a good justification for the formula.  We
% remarked that under the {\em correct} model, $ E[C_p]\approx p $, so values
% of $C_p$ near p could be used to indicate models which may be worth
% pursuing.  In order to understand what is going on, we need to explore
% this a little further.  
The following is modeled after Weisberg (1985, pp. 293--295).  A good
historical source, with different notation, is C.~L.~Mallows, {\em
Some Comments on $C_p$}, Technometrics, Vol.~15, p.~661, 1973.

Suppose we have included enough regressors in the ``full model'' 
\begin{equation}
        \mb y_{n\times1} = X_{n\times k'}\mb\beta_{k'\times1} +
        \epsilon_{n\times1} , 
        \label{full}
\end{equation}
to have unbiased estimates of the $\beta$'s and so that $E[\mb y] =
X_{n\times k'}\mb\beta_{k'\times1}$.  Now suppose we partition
\[
        X ~ = ~ [ ~ X_1 ~~ X_2 ~ ]
\]
where $X_1$ is $n\times p$ and $X_2$ is $n\times (k'-p)$, and let us
consider also the smaller model
\begin{equation}
        \mb y_{n\times1} = X_1\mb\beta_1 + \epsilon^*_{n\times1} .
        \label{reduced}
\end{equation}
Denote the fitted values of the full model \eq{full} as $\hat
Y_i$, and let the full model hat matrix be $H = X (X\tp X)^{-1} X\tp$
with diagonal element $h_{ii}$.  Let the fitted values of the reduced
model \eq{reduced} be $\hat y_i$ and let the reduced model hat matrix
be $U = X_1 (X_1\tp X_1)^{-1} X_1\tp$ with diagonal elements $u_{ii}$.

A widely used measure of the performance of a model (the reduced model
in our case) is the (sum of) mean-squared errors of the fitted values,
\[
        J_p = \frac{1}{\sigma^2_\epsilon} \sum_{i=1}^n E[(\hat y_i - E[y_i])^2]
        = \frac{1}{\sigma^2_\epsilon} \sum_{i=1}^n MSE_i
\]
where the expected values are taken under the {\em full} model \eq{full},
and we can readily calculate that the mean-square error can be
decomposed as usual into $($variance$) + ($bias$)^2$:
\begin{eqnarray*}
        MSE_i & = & E[(\hat y_i - E[y_i])^2] 
        ~ = ~ \Var(\hat y_i) + (E[\hat y_i] - E[y_i])^2 \\
        & = & \sigma^2\epsilon u_{ii} + (E[\hat y_i] - E[y_i])^2 \\
        & = & \sigma^2\epsilon u_{ii} + (E[\hat y_i] - E[\hat Y_i])^2 ,
\end{eqnarray*}
where the last line follows since the full model \eq{full} is assumed
to be unbiased.  Weisberg shows that
\[
        (E[\hat y_i] - E[\hat Y_i])^2 = E[(\hat y_i - \hat Y_i)^2]
        - \sigma^2_\epsilon(h_{ii} - u_{ii})
\]
(this result is not obvious, though it seems to be related to our
earlier results for the expected value of a quadratic form).  If we
substitute these identities into the expression for $MSE_i$, estimate
$E[(\hat y_i - \hat Y_i)^2]$ with its observed value $(\hat y_i - \hat
Y_i)^2$, and estimate $\sigma^2\epsilon$ with $\sigsq$ from the full
model, we obtain
\[
        \widehat{MSE}_i = (\hat y_i - \hat Y_i)^2 + \sigsq u_{ii}
        - \sigsq(h_{ii} - u_{ii})
\]
and hence we may estimate $J_p$ with 
\begin{eqnarray*}
        C_p & \equiv & \hat J_p 
        ~ = ~ \frac{1}{\sigsq} \sum_{i=1}^n \widehat{MSE}_i \\
        & = & \sum_{i=1}^n \left[\frac{(\hat y_i - \hat Y_i)^2}{\sigsq} 
                + u_{ii} - (h_{ii} - u_{ii}) \right]
\end{eqnarray*}
The terms in square brackets may be interpreted as how much each
observation contributes to $C_p$, and are sometimes analyzed
individually.  To obtain the form of $C_p$ in Equation \eq{mallows} we
need to substitute the following facts into this expression for $C_p$:
\begin{itemize}
\item $\sum_i (\hat y_i - \hat Y_i)^2 = RSS_p - RSS_{full}$ (draw the
right pythagorean theorem triangles in $n$-space);
\item $\sum_i u_{ii} = \mbox{trace}\,U = p$, for the reduced model;
\item $\sum_i h_{ii} = \mbox{trace}\,H = p$, for the full model.
\end{itemize}
Doing so we obtain
\[
        C_p = \frac{RSS_p - RSS_{full}}{\sigsq} + 2p-n
\]
and we note that $RSS_{full}/\sigsq$ does not vary with $p$, so it is
dropped from $C_p$.

If you have ever heard of AIC (Akaike's Information Criterion), you
will note that $C_p$ looks an awful lot like
``$-2\log(\mbox{lik'hood}) + 2(\mbox{\# param's})$'', the
AIC-penalized log-likelihood criterion for model selection.


\subsection{Stepwise Selection Procedures}

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.}
\begin{enumerate}\tight
\item Stop after a given number of steps or (equivalently) when a subset of
a certain size has been reached
\item Use F-tests to decide when to enter a variable, and stop when no
variable has a big F-test
\item Stop when adding the next variable would result in collinearity 
$ \Longleftrightarrow $ Stop when $R^2$ gets too big or $1 - R^2$ (the
tolerance) gets too small
\end{enumerate}

\paragraph{Backward Elimination.}
\begin{enumerate}\tight
\item Stop after a given number of steps
\item Use F-tests to remove variables, and stop when it is no
longer possible
\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.

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.

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)

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.

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.

Splus has a function {\tt step} that will perform stepwise regression,
though it is based on minimizing $C_p$ rather than using entry and
removal criteria with $F$ statistics.  Since one can show that
\[
        C_p = (k'-p)(F_p-1) + p
\]
where $k$ is the maximum number of predictors in the problem, asking
for $C_p\approx p$ is like asking for $F_p\approx 1$, i.e. asking for
what we expect of $F$ under the null hypothesis that these $p$
predictors are adequate.  Generally speaking, $C_p$ is minimized at
about the same place that $C_p\approx p$; when these criteria are used
in stepwise regression they perform similarly to having F-IN = F-OUT =
2, so they tend to add more new variables and not take them out
later, than the conventional choices, e.g. F-IN = F-OUT = 4.



\subsection{All subsets regressions}


There are algorithms
(e.g., {\tt BMDP9R}, and {\tt leaps} in S) which will consider all the
possible regression (all subsets regression) and report the best ones.
The algorithms are ``smart'' in that they don't actually fit {\em all}
the possible models---only the promising ones.  This is possible
because it is sometimes possible to rule out a entire class a models
based on the information found from one model.  Suppose that $(x_1,
x_3, x_5)$ has a {\em really} bad fit.  It is possible (and even
likely) for the fit to be bad enough that all submodels containing a
subset of $(x_1, x_3, x_5)$ could not possibly be good fits. (Remember
that the submodels can not fit any better than the 3 term model).  The
BMDP and S routines gives choices of selection criterea, either $C_p$,
or adjusted $R^2$ or $R^2$ can be used as the screening method.




\section{Some Examples of Model Selection}

We have read the highway accident data into Splus with the following
variable labels:
{\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
and the data itself is
{\small\begin{verbatim}
Rate   Len   Adt  Trk Sp  lw  shl itg   sigs   acpt lan type
 4.58  4.99   69   8  55  12  10  1.20  0      4.60  8  1
 2.86  16.11  73   8  60  12  10  1.43  0      4.40  4  1
 3.02  9.75   49  10  60  12  10  1.54  0      4.70  4  1
 2.29  10.65  61  13  65  12  10  0.94  0      3.80  6  1
 1.61  20.01  28  12  70  12  10  0.65  0      2.20  4  1
 6.87  5.97   30   6  55  12  10  0.34  1.84  24.80  4  2
 3.85  8.57   46   8  55  12   8  0.47  0.70  11.00  4  2
 6.12  5.24   25   9  55  12  10  0.38  0.38  48.50  4  2
 3.29  15.79  43  12  50  12   4  0.95  1.39   7.50  4  2
 5.88  8.26   23   7  50  12   5  0.12  1.21   8.20  4  2
 4.20  7.03   23   6  60  12  10  0.29  1.85   5.40  4  2
 4.61  13.28  20   9  50  12   2  0.15  1.21  11.20  4  2
 4.80  5.40   18  14  50  12   8  0     0.56  15.20  2  2
 3.85  2.96   21   8  60  12  10  0.34  0      5.40  4  2
 2.69  11.75  27   7  55  12  10  0.26  0.60   7.90  4  2
 1.99  8.86   22   9  60  12  10  0.68  0      3.20  4  2
 2.01  9.78   19   9  60  12  10  0.20  0.10  11.00  4  2
 4.22  5.49    9  11  50  12   6  0.18  0.18   8.90  2  2
 2.76  8.63   12   8  55  13   6  0.14  0     12.40  2  2
 2.55  20.31  12   7  60  12  10  0.05  0.99   7.80  4  2
 1.89  40.09  15  13  55  12   8  0.05  0.12   9.60  4  2
 2.34  11.81   8   8  60  12  10  0     0      4.30  2  2
 2.83  11.39   5   9  50  12   8  0     0.09  11.10  2  2
 1.81  22.00   5  15  60  12   7  0     0      6.80  2  2
 9.23  3.58   23   6  40  12   2  0.56  2.51  53.00  4  3
 8.60  3.23   13   6  45  12   2  0.31  0.93  17.30  2  3
 8.21  7.73    7   8  55  12   8  0.13  0.52  27.30  2  3
 2.93  14.41  10  10  55  12   6  0     0.07  18.00  2  3
 7.48  11.54  12   7  45  12   3  0.09  0.09  30.20  2  3
 2.57  11.10   9   8  60  12   7  0     0     10.30  2  3
 5.77  22.09   4   8  45  11   3  0     0.14  18.20  2  3
 2.90  9.39    5  10  55  13   1  0     0     12.30  2  3
 2.97  19.49   4  13  55  12   4  0     0      7.10  2  3
 1.84  21.01   5  12  55  10   8  0     0.10  14.00  2  3
 3.78  27.16   2  10  55  12   3  0.04  0.04  11.30  2  3
 2.76  14.03   3   8  50  12   4  0.07  0     16.30  2  3
 4.27  20.63   1  11  55  11   4  0     0      9.60  2  3
 3.05  20.06   3  11  60  12   8  0     0      9.00  2  4
 4.12  12.91   1  10  55  12   3  0     0     10.40  2  4
\end{verbatim}}\noindent

\subsection{Variable Selection by Hand}

Lets start by reading the data in S, doing some plots and thinking
about the modeling. 

{\small\begin{verbatim}
> high <- read.table("highway.dat", header=T)
\end{verbatim}}\noindent
I've chosen to split the plots into two groups, because I found one
big scatterplot matrix too much to handle.
The first six variables are first, then the response plus the second
six variables.
\begin{figure}[htbp]
  \centerline{\psfig{figure=high.pairs1.eps,height=3.5in,angle=-90}}
  \caption{First six variables}
  \label{fig:pairs1}
\end{figure}
\begin{figure}[htbp]
  \centerline{\psfig{figure=high.pairs2.eps,height=3.5in,angle=-90}}
  \caption{Response plus second six variables}
  \label{fig:pairs2}
\end{figure}

We gain some insight by looking at these plots, but the general
picture is still confused.

Perhaps the next thing we might try is fitting a model with {\bf all}
of the predictors included to see what happens.

Before we can proceed, lets make sure that our factor (the type of
highway) really is a factor.
{\small\begin{verbatim}
> is.factor(high$type)
[1] F
\end{verbatim}}\noindent
It is not, so we better fix it.
{\small\begin{verbatim}
> high$type <- as.factor(high$type)
> is.factor(high$type)             
[1] T
\end{verbatim}}\noindent

Now, on to the analysis.
{\small\begin{verbatim}
> high.lm<-lm(Rate~Len+Adt+Trk+Sp+lw+shl+itg+sigs+acpt+lan+type,data=high)
> summary(high.lm)

Call: lm(formula = Rate ~ Len + Adt + Trk + Sp + lw + shl + itg + sigs + acpt + lan + 
        type, data = high)
Residuals:
    Min      1Q  Median     3Q   Max 
 -1.992 -0.5853 0.03342 0.5361 2.671

Coefficients:
               Value Std. Error  t value Pr(>|t|) 
(Intercept)  13.6011   6.5909     2.0636   0.0496
        Len  -0.0586   0.0320    -1.8343   0.0785
        Adt  -0.0017   0.0323    -0.0517   0.9591
        Trk  -0.1140   0.1094    -1.0428   0.3070
         Sp  -0.1230   0.0699    -1.7592   0.0908
         lw  -0.1212   0.5462    -0.2220   0.8261
        shl  -0.0026   0.1467    -0.0177   0.9860
        itg  -0.6455   1.2228    -0.5279   0.6022
       sigs   0.8392   0.4378     1.9171   0.0667
       acpt   0.0593   0.0259     2.2869   0.0309
        lan  -0.0313   0.2721    -0.1149   0.9094
      type1  -0.9460   0.6257    -1.5119   0.1431
      type2  -0.1452   0.3611    -0.4021   0.6910
      type3   0.0607   0.2523     0.2407   0.8117
\end{verbatim}}\noindent

Nothing looks very significant! So we may have some problems here.
There are a few variables that might hold some hope.
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}

The standard residual diagnostics aren't too exciting.  The residuals
look a little heavy tailed and perhaps have somewhat unstable
variance, but it is nothing to write home about.  One residual is a
little large (2.7) and there seem to be no highly influential points,
although two points have substantially higher Cooks Distances than the
others (at least one of them clearly has high leverage according to
the top row of plots). 
% \begin{figure}[htbp]
%   \centerline{\psfig{figure=high.rvp.eps,height=3.5in,angle=-90}}
%   \caption{Residuals vs. predict}
%   \label{fig:rvp}
% \end{figure}

A most disappointing aspect of the problem is that the model doesn't
fit very well even with all the predictors.

A standard heuristic method is to begin removing the least
significant variables, viz:
{\small\begin{verbatim}
> summary(high.lmw <- update(high.lm, . ~ . - shl))

Call: lm(formula = Rate ~ Len + Adt + Trk + Sp + lw + itg + sigs + acpt + lan + type, 
        data = high)
Residuals:
    Min      1Q  Median     3Q   Max 
 -1.993 -0.5861 0.03621 0.5378 2.666

Coefficients:
               Value Std. Error  t value Pr(>|t|) 
(Intercept)  13.5644   6.1337     2.2115   0.0360
        Len  -0.0586   0.0310    -1.8886   0.0701
        Adt  -0.0018   0.0313    -0.0565   0.9554
        Trk  -0.1133   0.1000    -1.1339   0.2672
         Sp  -0.1240   0.0455    -2.7263   0.0113
         lw  -0.1158   0.4404    -0.2628   0.7948
        itg  -0.6444   1.1974    -0.5382   0.5950
       sigs   0.8416   0.4080     2.0627   0.0493
       acpt   0.0591   0.0227     2.6049   0.0150
        lan  -0.0310   0.2665    -0.1164   0.9082
      type1  -0.9472   0.6101    -1.5524   0.1327
      type2  -0.1425   0.3213    -0.4435   0.6611
      type3   0.0616   0.2422     0.2545   0.8011

Residual standard error: 1.12 on 26 degrees of freedom
Multiple R-Squared: 0.7826 
F-statistic: 7.798 on 12 and 26 degrees of freedom, the p-value is 6.757e-06 
\end{verbatim}}\noindent

Continuing in this vein we get,

{\small\begin{verbatim}
> summary(high.lmw <- update(high.lmw, . ~ . - Adt))

Call: lm(formula = Rate ~ Len + Trk + Sp + lw + itg + sigs + acpt + lan + type, data
         = high)
Residuals:
    Min      1Q   Median     3Q   Max 
 -1.987 -0.5805 0.005454 0.5402 2.667

Coefficients:
               Value Std. Error  t value Pr(>|t|) 
(Intercept)  13.5445   6.0095     2.2539   0.0325
        Len  -0.0584   0.0303    -1.9297   0.0642
        Trk  -0.1134   0.0981    -1.1565   0.2576
         Sp  -0.1236   0.0441    -2.8011   0.0093
         lw  -0.1160   0.4322    -0.2684   0.7904
        itg  -0.6859   0.9281    -0.7391   0.4662
       sigs   0.8404   0.3998     2.1018   0.0450
       acpt   0.0592   0.0222     2.6623   0.0129
        lan  -0.0384   0.2277    -0.1688   0.8672
      type1  -0.9398   0.5848    -1.6069   0.1197
      type2  -0.1377   0.3040    -0.4529   0.6542
      type3   0.0647   0.2317     0.2791   0.7823

Residual standard error: 1.099 on 27 degrees of freedom
Multiple R-Squared: 0.7825 
F-statistic: 8.833 on 11 and 27 degrees of freedom, the p-value is 2.175e-06 
\end{verbatim}}\noindent

{\small\begin{verbatim}
> summary(high.lmw <- update(high.lmw, . ~ . - lan))
> summary(high.lmw <- update(high.lmw, . ~ . - lw)) 
> summary(high.lmw <- update(high.lmw, . ~ . - itg))
> summary(high.lmw <- update(high.lmw, . ~ . - Trk))
\end{verbatim}}\noindent

At this stage I was not sure if {\tt type} was important or not, so I
looked at the ANOVA table.
{\small\begin{verbatim}
> anova(high.lmw)
Analysis of Variance Table

Response: Rate

Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F) 
Len        1  32.44949 32.44949 28.87014 0.0000074
Sp         1  54.84578 54.84578 48.79600 0.0000001
sigs       1   7.72198  7.72198  6.87021 0.0134579
acpt       1  12.91667 12.91667 11.49189 0.0019217
type       3   7.10872  2.36957  2.10820 0.1193586
Residuals 31  34.84342  1.12398                   
\end{verbatim}}\noindent

It looks like type might not be important, so remove it and see where
we are
{\small\begin{verbatim}
> summary(high.lmw <- update(high.lmw, . ~ . - type))

Call: lm(formula = Rate ~ Len + Sp + sigs + acpt, data = high)
Residuals:
    Min      1Q  Median     3Q   Max 
 -1.673 -0.8402 -0.2407 0.7347 2.906

Coefficients:
              Value Std. Error t value Pr(>|t|) 
(Intercept)  9.9756  2.3544     4.2370  0.0002 
        Len -0.0619  0.0254    -2.4410  0.0200 
         Sp -0.1157  0.0388    -2.9819  0.0053 
       sigs  0.5943  0.3312     1.7944  0.0816 
       acpt  0.0682  0.0211     3.2355  0.0027 

Residual standard error: 1.111 on 34 degrees of freedom
Multiple R-Squared: 0.7201 
F-statistic: 21.87 on 4 and 34 degrees of freedom, the p-value is 5.258e-09 

\end{verbatim}}\noindent

The variable {\tt sigs} still looks marginal, but at this stage I'm
prepared to call it quits.   I {\bf did not} keep any of my data in my
pocket, so I can't really cross-validate.   (We'll look at that
analysis later on \ldots).  

Our ``reasonable'' model comes out as
{\small\begin{verbatim}
 Rate ~ Len + Sp + sigs + acpt
\end{verbatim}}\noindent

\subsection{Some more ``Hand'' Details}

Well, that was a manual way of proceeding.  S has some more
automatic procedures, including {\tt add1} and {\tt drop1}, 
{\tt stepwise} and {\tt step}.   We will start with {\tt drop1}
which is just a better presentation of the summaries we have already
used.  For example, the summary from the full model, 
{\small\begin{verbatim}
> drop1(high.lm)
Single term deletions

Model:
Rate ~ Len + Adt + Trk + Sp + lw + shl + itg + sigs + acpt + lan + type
       Df Sum of Sq      RSS       Cp 
<none>              32.58891 69.08848
   Len  1  4.386136 36.97504 70.86751
   Adt  1  0.003490 32.59240 66.48486
   Trk  1  1.417476 34.00638 67.89885
    Sp  1  4.034260 36.62317 70.51563
    lw  1  0.064248 32.65315 66.54562
   shl  1  0.000408 32.58931 66.48178
   itg  1  0.363299 32.95221 66.84467
  sigs  1  4.790802 37.37971 71.27217
  acpt  1  6.817750 39.40666 73.29912
   lan  1  0.017216 32.60612 66.49859
  type  3  4.785352 37.37426 66.05250

\end{verbatim}}\noindent

would lead us to drop the {\tt shl} variable (as we did).  The columns
in the display are the Sum of Squares for the deleted term (like an
ANOVA with that term last), the Residual Sum of Squares for the model
without that term and the Cp statistic for that model.  This is
clearly a faster way of proceeding than the raw method we chose, but
still a manual process.   

There is also an {\tt add1} function, which takes a given model and
tries adding terms to it.  As an example
{\small\begin{verbatim}
> formula(high.lmw)
Rate ~ Len + Sp + sigs + acpt
> formula(high.lm) 
Rate ~ Len + Adt + Trk + Sp + lw + shl + itg + sigs + acpt + lan + type
> add1(high.lmw, high.lm)
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  3  7.108722 34.84342 54.58560
\end{verbatim}}\noindent

This is pretty clear.  The SS doesn't increase a lot (enough) for
adding any variable, and Cp basically goes up no matter what we do.
(Cp is already too big, we would like it to go down).   Just to finish
up this theme, lets try {\tt drop1} on the model we ended up with.
{\small\begin{verbatim}
> drop1(high.lmw)
Single term deletions

Model:
Rate ~ Len + Sp + sigs + acpt
       Df Sum of Sq      RSS       Cp 
<none>              41.95214 54.29100
   Len  1   7.35220 49.30434 59.17543
    Sp  1  10.97166 52.92379 62.79489
  sigs  1   3.97312 45.92526 55.79635
  acpt  1  12.91667 54.86881 64.73990
\end{verbatim}}\noindent

Again, this is pretty clear.  Dropping anything is a significant
reduction in SS, and Cp goes up.

\subsection{Stepwise variable selection}

Let us now turn to one of the more automatic procedures, {\tt
  step}.
For starters, lets do an example, and then discuss what is going on.  
{\small\begin{verbatim}
> step(high.lm)     
Start:  AIC= 69.0885 
 Rate ~ Len + Adt + Trk + Sp + lw + shl + itg + sigs + acpt + lan + type 

Single term deletions

Model:
Rate ~ Len + Adt + Trk + Sp + lw + shl + itg + sigs + acpt + lan + type

scale:  1.303556 

       Df Sum of Sq      RSS       Cp 
<none>              32.58891 69.08848
   Len  1  4.386136 36.97504 70.86751
   Adt  1  0.003490 32.59240 66.48486
   Trk  1  1.417476 34.00638 67.89885
    Sp  1  4.034260 36.62317 70.51563
    lw  1  0.064248 32.65315 66.54562
   shl  1  0.000408 32.58931 66.48178
   itg  1  0.363299 32.95221 66.84467
  sigs  1  4.790802 37.37971 71.27217
  acpt  1  6.817750 39.40666 73.29912
   lan  1  0.017216 32.60612 66.49859
  type  3  4.785352 37.37426 66.05250

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

Single term deletions

Model:
Rate ~ Len + Adt + Trk + Sp + lw + shl + itg + sigs + acpt + lan

scale:  1.303556 

       Df Sum of Sq      RSS       Cp 
<none>              37.37426 66.05250
   Len  1   4.12668 41.50094 67.57207
   Adt  1   0.05918 37.43344 63.50456
   Trk  1   2.57323 39.94749 66.01862
    Sp  1   1.34287 38.71713 64.78825
    lw  1   1.05079 38.42505 64.49617
   shl  1   1.77280 39.14706 65.21819
   itg  1   0.05093 37.42519 63.49631
  sigs  1   1.99657 39.37083 65.44195
  acpt  1  13.43369 50.80795 76.87908
   lan  1   0.02185 37.39611 63.46723

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

Single term deletions

Model:
Rate ~ Len + Adt + Trk + Sp + lw + shl + itg + sigs + acpt

scale:  1.303556 

       Df Sum of Sq      RSS       Cp 
<none>              37.39611 63.46723
   Len  1   4.27099 41.66710 65.13111
   Adt  1   0.03756 37.43367 60.89768
   Trk  1   2.55160 39.94771 63.41172
    Sp  1   1.38925 38.78535 62.24937
    lw  1   1.05204 38.44814 61.91216
   shl  1   1.84401 39.24012 62.70413
   itg  1   0.06174 37.45784 60.92186
  sigs  1   2.02920 39.42531 62.88932
  acpt  1  13.42385 50.81995 74.28397

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

Single term deletions

Model:
Rate ~ Len + Trk + Sp + lw + shl + itg + sigs + acpt

scale:  1.303556 

       Df Sum of Sq      RSS       Cp 
<none>              37.43367 60.89768
   Len  1   4.27413 41.70780 62.56470
   Trk  1   2.52254 39.95621 60.81311
    Sp  1   1.53385 38.96752 59.82442
    lw  1   1.01633 38.44999 59.30690
   shl  1   1.90372 39.33738 60.19428
   itg  1   0.81888 38.25255 59.10945
  sigs  1   2.28295 39.71662 60.57352
  acpt  1  13.60634 51.04001 71.89691

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

Single term deletions

Model:
Rate ~ Len + Trk + Sp + lw + shl + sigs + acpt

scale:  1.303556 

       Df Sum of Sq      RSS       Cp 
<none>              38.25255 59.10945
   Len  1   5.27342 43.52597 61.77575
   Trk  1   2.37446 40.62701 58.87680
    Sp  1   1.53567 39.78822 58.03801
    lw  1   0.93877 39.19132 57.44111
   shl  1   1.45448 39.70703 57.95682
  sigs  1   2.66569 40.91824 59.16803
  acpt  1  12.93989 51.19244 69.44222

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

Single term deletions

Model:
Rate ~ Len + Trk + Sp + shl + sigs + acpt

scale:  1.303556 

       Df Sum of Sq      RSS       Cp 
<none>              39.19132 57.44111
   Len  1   4.36685 43.55817 59.20085
   Trk  1   2.13786 41.32918 56.97186
    Sp  1   2.58149 41.77282 57.41549
   shl  1   0.91295 40.10427 55.74695
  sigs  1   2.60101 41.79234 57.43501
  acpt  1  13.00508 52.19640 67.83908

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

Single term deletions

Model:
Rate ~ Len + Trk + Sp + sigs + acpt

scale:  1.303556 

       Df Sum of Sq      RSS       Cp 
<none>              40.10427 55.74695
   Len  1   3.74506 43.84933 56.88490
   Trk  1   1.84787 41.95214 54.98770
    Sp  1  10.14561 50.24988 63.28544
  sigs  1   2.35236 42.45664 55.49220
  acpt  1  12.44929 52.55357 65.58913

Step:  AIC= 54.9877 
 Rate ~ Len + Sp + sigs + acpt 

Single term deletions

Model:
Rate ~ Len + Sp + sigs + acpt

scale:  1.303556 

       Df Sum of Sq      RSS       Cp 
<none>              41.95214 54.98770
   Len  1   7.35220 49.30434 59.73279
    Sp  1  10.97166 52.92379 63.35224
  sigs  1   3.97312 45.92526 56.35371
  acpt  1  12.91667 54.86881 65.29726
Call:
lm(formula = Rate ~ Len + Sp + sigs + acpt, data = high)

Coefficients:
 (Intercept)         Len         Sp      sigs       acpt 
    9.975636 -0.06188818 -0.1157323 0.5943472 0.06824634

Degrees of freedom: 39 total; 34 residual
Residual standard error (on weighted scale): 1.110804
\end{verbatim}}\noindent

That was pretty painless (and happened literally faster than I could
watch it scroll by), and we ended up with the same model.  
The procedure seems to have been using something called AIC as the
thing it was optimizing.  What is that all about?

Before we talk about that, lets see the same sort of thing going
forwards (instead of backwards), or better yet, both ways.
{\small\begin{verbatim}
> formula(high.lmb)                                   
Rate ~ Sp
> step(high.lmb,formula(high.lm),direction="both")
Start:  AIC= 89.0676 
 Rate ~ Sp 

Single term deletions

Model:
Rate ~ Sp

scale:  2.172381 

       Df Sum of Sq      RSS       Cp 
<none>               80.3781  89.0676
    Sp  1  69.50797 149.8861 154.2308
Single term additions

Model:
Rate ~ Sp

scale:  2.172381 

       Df Sum of Sq      RSS       Cp 
<none>              80.37810 89.06762
   Len  1  17.78730 62.59079 75.62508
   Adt  1   3.02205 77.35604 90.39033
   Trk  1  15.87336 64.50474 77.53902
    lw  1   0.57414 79.80396 92.83825
   shl  1   1.93260 78.44549 91.47978
   itg  1   3.09570 77.28240 90.31668
  sigs  1  14.65073 65.72736 78.76165
  acpt  1  23.76243 56.61566 69.64995
   lan  1   3.48992 76.88817 89.92246
  type  3   3.42142 76.95668 98.68049

Step:  AIC= 69.6499 
 Rate ~ Sp + acpt 

Single term deletions

Model:
Rate ~ Sp + acpt

scale:  2.172381 

       Df Sum of Sq      RSS       Cp 
<none>              56.61566 69.64995
    Sp  1  15.72195 72.33762 81.02714
  acpt  1  23.76243 80.37810 89.06762
Single term additions

Model:
Rate ~ Sp + acpt

scale:  2.172381 

       Df Sum of Sq      RSS       Cp 
<none>              56.61566 69.64995
   Len  1  10.69041 45.92526 63.30430
   Adt  1   3.66056 52.95510 70.33415
   Trk  1   9.80055 46.81512 64.19417
    lw  1   0.33659 56.27908 73.65813
   shl  1   0.22289 56.39277 73.77182
   itg  1   3.44796 53.16771 70.54676
  sigs  1   7.31132 49.30434 66.68339
   lan  1   3.00228 53.61338 70.99243
  type  3   3.88890 52.72677 78.79534

Step:  AIC= 63.3043 
 Rate ~ Sp + acpt + Len 

Single term deletions

Model:
Rate ~ Sp + acpt + Len

scale:  2.172381 

       Df Sum of Sq      RSS       Cp 
<none>              45.92526 63.30430
    Sp  1  14.84094 60.76620 73.80048
  acpt  1  16.66554 62.59079 75.62508
   Len  1  10.69041 56.61566 69.64995
Single term additions

Model:
Rate ~ Sp + acpt + Len

scale:  2.172381 

       Df Sum of Sq      RSS       Cp 
<none>              45.92526 63.30430
   Adt  1  0.683412 45.24184 66.96565
   Trk  1  3.468619 42.45664 64.18045
    lw  1  0.305279 45.61998 67.34379
   shl  1  0.283523 45.64173 67.36554
   itg  1  0.742574 45.18268 66.90649
  sigs  1  3.973117 41.95214 63.67595
   lan  1  0.802226 45.12303 66.84684
  type  3  3.559542 42.36571 72.77905
Call:
lm(formula = Rate ~ Sp + acpt + Len, data = high)

Coefficients:
 (Intercept)         Sp       acpt         Len 
    11.10415 -0.1312275 0.07591309 -0.07254915

Degrees of freedom: 39 total; 35 residual
Residual standard error (on weighted scale): 1.145491 

\end{verbatim}}\noindent
We have again converged on a model the same as (at least similar to)
what we had before.

So what is the AIC stuff.  For our purposes, AIC is another version of
Cp, and we are aiming to minimize it (rather than get $CP = p$.)  For
models that are not linear, Gaussian, there is a more general
statistical call Akaike's Information Criterion, that reduces to Cp
for our models.  Some details of the connection between AIC and $C_p$
are contained on pp. xxx--xxx [Ch. 6?] of Venables and
Ripley\footnote{There is also something call BIC (for Bayesian
Information Criterion), and there is Splus code to compute it---but we
will ignore that for now.}.

Lets try a raw stepwise, and see where that leads us.
{\small\begin{verbatim}
> formula(mod1)                                       
Rate ~ 1
> formula(mod2)
Rate ~ Len + Adt + Trk + Sp + lw + shl + itg + sigs + acpt + lan + type
( Lots of output edited out below......)
> step(mod1, formula(mod2), direction="both", trace=T)
 Rate ~ 1 
Step:  AIC= 88.1151 
 Rate ~ acpt 

Step:  AIC= 80.2819 
 Rate ~ acpt + Sp 

Step:  AIC= 77.4802 
 Rate ~ acpt + Sp + Len 

Call:
lm(formula = Rate ~ acpt + Sp + Len, data = high)

Coefficients:
 (Intercept)       acpt         Sp         Len 
    11.10415 0.07591309 -0.1312275 -0.07254915

Degrees of freedom: 39 total; 35 residual
Residual standard error (on weighted scale): 1.145491 
\end{verbatim}}\noindent

\subsection{All subsets regression}

That is probably enough on stepwise methods.  It is possible to also
do all possible subsets, using the {\tt leaps} command.

{\small\begin{verbatim}
> attach(high)
> result _ leaps(cbind(Len,Adt,Trk,Sp,lw,shl,itg,sigs,acpt,lan,type),Rate,nbest=2)
> result
$Cp:
 [1] 17.700243 23.557983  8.246302  9.864009  2.457998  3.106291  1.563455
 [8]  1.930997  2.217225  3.109609  3.552115  3.927820  4.868188  5.012033
[15]  6.271606  6.289220  8.070722  8.093827 10.002276 10.063006 11.999997

$size:
 [1]  2  2  3  3  4  4  5  5  6  6  7  7  8  8  9  9 10 10 11 11 12

$label:
 [1] "acpt"                                        
 [2] "Sp"                                          
 [3] "Sp,acpt"                                     
 [4] "Trk,acpt"                                    
 [5] "Len,Sp,acpt"                                 
 [6] "Trk,Sp,acpt"                                 
 [7] "Len,Sp,sigs,acpt"                            
 [8] "Len,Trk,Sp,acpt"                             
 [9] "Len,Trk,Sp,sigs,acpt"                        
[10] "Len,Sp,shl,sigs,acpt"                        
[11] "Len,Trk,Sp,shl,sigs,acpt"                    
[12] "Len,Trk,Sp,lw,sigs,acpt"                     
[13] "Len,Trk,Sp,lw,shl,sigs,acpt"                 
[14] "Len,Trk,Sp,shl,itg,sigs,acpt"                
[15] "Len,Trk,Sp,lw,shl,itg,sigs,acpt"             
[16] "Len,Adt,Trk,Sp,lw,shl,sigs,acpt"             
[17] "Len,Adt,Trk,Sp,lw,shl,sigs,acpt,type"        
[18] "Len,Trk,Sp,lw,shl,itg,sigs,acpt,type"        
[19] "Len,Adt,Trk,Sp,lw,shl,itg,sigs,acpt,type"    
[20] "Len,Adt,Trk,Sp,lw,shl,sigs,acpt,lan,type"    
[21] "Len,Adt,Trk,Sp,lw,shl,itg,sigs,acpt,lan,type"

$which:
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] 
 [1,]    F    F    F    F    F    F    F    F    T     F     F
 [2,]    F    F    F    T    F    F    F    F    F     F     F
 [3,]    F    F    F    T    F    F    F    F    T     F     F
 [4,]    F    F    T    F    F    F    F    F    T     F     F
 [5,]    T    F    F    T    F    F    F    F    T     F     F
 [6,]    F    F    T    T    F    F    F    F    T     F     F
 [7,]    T    F    F    T    F    F    F    T    T     F     F
 [8,]    T    F    T    T    F    F    F    F    T     F     F
 [9,]    T    F    T    T    F    F    F    T    T     F     F
[10,]    T    F    F    T    F    T    F    T    T     F     F
[11,]    T    F    T    T    F    T    F    T    T     F     F
[12,]    T    F    T    T    T    F    F    T    T     F     F
[13,]    T    F    T    T    T    T    F    T    T     F     F
[14,]    T    F    T    T    F    T    T    T    T     F     F
[15,]    T    F    T    T    T    T    T    T    T     F     F
[16,]    T    T    T    T    T    T    F    T    T     F     F
[17,]    T    T    T    T    T    T    F    T    T     F     T
[18,]    T    F    T    T    T    T    T    T    T     F     T
[19,]    T    T    T    T    T    T    T    T    T     F     T
[20,]    T    T    T    T    T    T    F    T    T     T     T
[21,]    T    T    T    T    T    T    T    T    T     T     T

$int:
[1] T
\end{verbatim}}\noindent
{\t leaps} has given us the two best regressions (in the sense of
minimizing $C_p$, which is the same as minimizing $RSS$ for fixed $p$,
for each $p=2, 3, 4, \ldots, 12$.  A convenient way to summarize part
of the results is
{\small\begin{verbatim}
> print(prt_cbind(result$Cp, result$size))  # or: data.frame(label,size,Cp)
           [,1] [,2] 
 [1,] 17.700243    2
 [2,] 23.557983    2
 [3,]  8.246302    3
 [4,]  9.864009    3
 [5,]  2.457998    4
 [6,]  3.106291    4
 [7,]  1.563455    5
 [8,]  1.930997    5
 [9,]  2.217225    6
[10,]  3.109609    6
[11,]  3.552115    7
[12,]  3.927820    7
[13,]  4.868188    8
[14,]  5.012033    8
[15,]  6.271606    9
[16,]  6.289220    9
[17,]  8.070722   10
[18,]  8.093827   10
[19,] 10.002276   11
[20,] 10.063006   11
[21,] 11.999997   12
\end{verbatim}}\noindent
If we take minimizing $C_p$ as the goal, then the best model is model 7:
{\small\begin{verbatim}
> names(high)[-1][result$which[7,]]
[1] "Len"  "Sp"   "sigs" "acpt"
\end{verbatim}}\noindent
If we take finding the least $p$ for which $C_p\approx p$ then model 6
wins:
{\small\begin{verbatim}
> names(high)[-1][result$which[6,]]
[1] "Trk"  "Sp"   "acpt"
\end{verbatim}}\noindent

Let's take a look at the parameter estimates and diagnostic plots of
the three models that came up most frequently with these methods.

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

\centerline{\psfig{figure=plot-lm2.ps,height=3.5in}}

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

\centerline{\psfig{figure=plot-lm3.ps,height=3.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}}

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



