formulas.Rnw

\documentclass[11pt,letterpaper]{article}
\usepackage[letterpaper,left=0.8in,top=1.2in]{geometry}
\usepackage{amsmath,amsfonts,amssymb,amsthm,bm,Sweave,paralist,hyperref,fancyvrb}
\setkeys{Gin}{width=\textwidth}
<<preliminaries,echo=FALSE,results='hide'>>=
library(knitr)
library(grid)
library(ggplot2)
opts_chunk$set(fig.width=8)
options(width=85, show.signif.stars = FALSE, str = strOptions(strict.width = "cut"))
@
\newcommand{\R}{\textsf{R}}
\newcommand{\code}[1]{\texttt{\small #1}}
\newcommand{\trans}{\ensuremath{^\prime}}
\begin{document}
\section{Linear model formulas}
\label{sec:Formula}

In \R{} linear models are specified using a \emph{model formula},
which is an expression that contains a tilde (the $\sim$ character).
The response is on the left-hand side of the tilde, typically as the
name of a variable, e.g.~\code{optden}, but it can also be a function
of a variable, e.g.~\code{log(BrainWt)}.

The right-hand side of the formula is composed of \emph{model terms}
separated by plus signs.  In the formulas below we write the response
as \code{y}, continuous covariates as \code{x}, \code{z}, \code{u},
$\dots$ and categorical covariates as \code{f} and \code{g}.  Note
that the categorical covariates are assumed to be stored as factors
(which includes ordered factors).

Some of the formulas for typical models are:
\begin{description}
\item[Simple linear regression] The formula
<<slr,echo=FALSE>>=
y ~ x
@
denotes the simple linear regression model
\begin{displaymath}
  y_i=\beta_0+\beta_1 x_i+\epsilon_i,\quad i=1,\dots,n
\end{displaymath}
In the formula shown above, the intercept term is implicit.  If you
prefer to make it explicit you can write the formula as
<<slr1,echo=FALSE>>=
y ~ 1 + x
@

\item[Regression through the origin] If you do not want the intercept
  term in the model, you must suppress it using the formula
<<regthr0,echo=FALSE>>=
y ~ 0 + x
@
or, alternatively,
<<regthr0alt,echo=FALSE>>=
y ~ x - 1
@
  The model specified in this way can be written as
  \begin{displaymath}
    y_i=\beta x_i+\epsilon_i,\quad i=1,\dots,n
  \end{displaymath}
  Notice that you can remove terms, even the implicit intercept, with
  a negative sign.

\item[Multiple linear regression] Multiple covariates can be listed on the right hand side, as in
<<mlr,echo=FALSE>>=
y ~ 1 + x + z + u
@
corresponding to the model
\begin{displaymath}
  y_i=\beta_0+\beta_1 x_i+\beta_2 z_i+\beta_3 u_i+\epsilon_i,\quad i=1,\dots,n
\end{displaymath}
\item[Polynomial regression] To include polynomial terms in the model
  you must protect the circumflex operator by surrounding the term
  with \code{I()}, which is the identity operator.  It implies that
  the expression inside is to be taken literally in terms of the
  arithmetic operators, not as the formula language operators.  The model
  \begin{displaymath}
    y_i=\beta_0+\beta_1 x_i+\beta_2 x_i^2 + \beta_3 x_i^3+\epsilon_i,\quad i=1,\dots,n
  \end{displaymath}
  is written
<<polyreg,echo=FALSE>>=
y ~ x + I(x^2) + I(x^3)
@
  Another specification for a polynomial regression model uses the
  \code{poly()} function which generates \emph{orthogonal polynomial}
  terms.  The fitted responses will be the same from the model shown
  above and from
<<polyfunc,echo=FALSE>>=
y ~ poly(x, 3)
@
  but the coefficients will be different because they are defined with
  respect to the orthogonal polynomials. These have some advantages if
  you are doing the calculations by hand but in practice you don't
  expect to be doing so.
\item[One categorical covariate] The model described as a one-way
  analysis of variance for the levels of factor, \code{f}, corresponds
  to the formula
<<one-way,echo=FALSE>>=
y ~ f
@
  Often we use the function \code{aov()} instead of \code{lm()} to fit
  such models.  \code{aov()} is the same as \code{lm()} except that it
  puts an extra tag on the fitted model that designates it as only
  having categorical covariates.  This changes, for example, the
  \code{summary()} method, which produces an analysis of variance
  table instead of a summary of the estimated coefficients.  The model
  that is fit is sometimes written as
  \begin{displaymath}
    y_{ij}=\mu+\alpha_i+\epsilon_{ij},\quad i=1,\dots,I\;j=1,\dots,n_i
  \end{displaymath}
  although it is not fit in that form.
\item[Two categorical covariates, additive] The formula for an
  additive two-factor analysis of variance model,
  \begin{displaymath}
    y_{ijk}=\mu+\alpha_i+\beta_j+\epsilon_{ijk},\quad i=1,\dots,I\;j=1,\dots,J\;k=1,\dots,n_{ij} ,
  \end{displaymath}
  is
<<two-way-additive,echo=FALSE>>=
y ~ f + g
@
  This produces a two-factor analysis of variance table.  In the
  balanced case the analysis of variance table for this model is
  equivalent to that for the model
<<two-way-additive-rev,echo=FALSE>>=
y ~ g + f
@
in the sense that, although the rows of the table will be in a
different order, they are otherwise the same.  For unbalanced data the
order of the factors is important.  The sums of squares in the table
are \emph{sequential} sums of squares corresponding to the
contribution of the first factor, given the intercept, then the
contribution of the second factor, given the first factor and the
intercept, and so on. In particular, \emph{blocking factors}, which
represent uncontrolable sources of variability, should be listed
before experimental factors.

\item[Two categorical covariates, allowing for interactions] If the
  data include replicate observations (more than one observation at
  the same combination of covariate values) we can fit and analyze a
  model with interaction terms with a formula like
<<two-way-interaction,echo=FALSE>>=
y ~ f + g + f:g
@
  where an expression like \code{f:g} is a two-factor interaction.
  Similar expressions are used for higher-order interactions.  This
  model can also be expressed as
<<two-way-interaction-alt,echo=FALSE>>=
y ~ f * g
@
  In general the asterisk operator, (\code{*}), generates the main
  effects plus interactions.  A three-factor model with all the main
  effects, two-factor interactions and the three-factor interaction
  can be written as
<<three-factor-interaction,echo=FALSE>>=
y ~ f * g * h
@
\item[Combination of continuous and categorical covariates] What is
  sometimes called an \emph{analysis of covariance} model incorporates
  both categorical and numeric covariates.  If there is only one
  numeric covariate, \code{x}, then the model can be described in
  terms of the lines formed by the fitted values on the \code{y}
  versus \code{x} plot.  The most common models are the parallel lines
  (different intercepts, same slope) generated by
<<ancova,echo=FALSE>>=
y ~ f + x
@
  and the model in which slopes and intercepts both vary according to
  the levels of \code{f}
<<ancovafull,echo=FALSE>>=
y ~ f * x
@
  which is equivalent to
<<ancovafullalt,echo=FALSE>>=
y ~ f + x + f:x
@
  Occasionally we incorporate an interaction term without a
  main-effect for \code{f}.
<<zerodose,echo=FALSE>>=
y ~ x + f:x
@
  I call this the ``zero-dose'' model because it is used in the case
  that \code{x} represents something like a dose and the levels of
  \code{f} corresponds to different treatments.  We don't have a main
  effect for \code{f} in such a model because a zero dose of treatment
  1 is the same as a zero dose of treatment 2.  Thus the lines for the
  different levels of the factors should coincide at \code{x=0}.
\end{description}

\section{Examples}
\label{sec:examples}

The \code{datasets} package contains several sample datasets that have
been used in different texts.  By default, this package is attached in
an \R{} session.

\paragraph{Simple linear regression}

<<carbplot,echo=FALSE,fig.height=5,fig.width=6,fig.align='center',fig.pos='tb',fig.cap='Observations of optical density versus carbohydrate amount from the calibration of a Formaldehyde assay.'>>=
print(qplot(carb,optden,data=Formaldehyde,xlab="Carbohydrate amount (ml)",
            ylab="Optical density")+geom_smooth(method="lm",se=FALSE))
@
The \code{Formaldehyde} data are a simple example from a calibration
study consisting of 6 observations of the carbohydrate content (ml.)
(variable \code{carb}) and the corresponding optical density (variable
\code{optden}).  Figure~\ref{fig:carb} is a data plot with the fitted simple linear regression line.  This model is fit as
<<lm1>>=
summary(lm1 <- lm(optden ~ 1 + carb, Formaldehyde))
@
\noindent(In what follows we will often skip the full summary output and
concentrate on the coefficients table, produced by
\code{coef(summary())}, or the analysis of variance table, produced by
\code{anova()}.)

\paragraph{Regression through the origin}

To constrain the line to pass through the origin (that is, to suppress
the \code{(Intercept)} term) we fit the model as
<<lm1a>>=
coef(summary(lm1a <- lm(optden ~ 0 + carb, Formaldehyde)))
@
A comparative analysis of variance of these two models
<<lm1lm1acomp>>=
anova(lm1a,lm1)
@
\noindent produces the same p-value as the t-test on the intercept coefficient
in model \code{lm1}, which is as it should be, because these are two
versions of the same test.

\paragraph{Polynomial regression}

Alternatively, we could fit \code{optden} as a quadratic function of
\code{carb} using
<<lm1bshow,eval=FALSE>>=
coef(summary(lm1b <- lm(optden ~ 1 + carb + I(carb^2), Formaldehyde)))
@
<<lm1b,echo=FALSE>>=
printCoefmat(coef(summary(lm1b <- lm(optden ~ 1 + carb + I(carb^2), Formaldehyde))))
@
\noindent Notice that the quadratic term is significant at the 5\% level and
generally we would retain it in the model.  The reason we don't see
much curvature in the data plot (Fig.~\ref{fig:carb}) is because there
is such a strong linear trend that it masks any nonlinear behaviour.

An alternative specification is
<<lm1cshow,eval=FALSE>>=
coef(summary(lm1c <- lm(optden ~ poly(carb, 2), Formaldehyde)))
@
<<lm1c,echo=FALSE>>=
printCoefmat(coef(summary(lm1c <- lm(optden ~ poly(carb, 2), Formaldehyde))))
@

\paragraph{Multiple linear regression}

The \code{trees} data are measurements of the volume of usable lumber
(variable \code{Volume}) from a sample of 31 black cherry trees.
Covariates are a measurement of the girth (\code{Girth}), which is
comparatively easy to measure (you just walk up to the tree and loop a
tape measure around it), and the height (\code{Height}), which is
somewhat more difficult to measure.  (There is some confusion in the
description of the data regarding whether the girth has been converted
to an equivalent diameter - we'll assume it is the girth.)  If we
consider the tree to have the shape of as a cylinder or a cone we
would expect that the volume would be related to the square of the
girth times the height.  In Fig.~\ref{fig:trees}
\begin{figure}[tb]
  \centering
<<treesplot,echo=FALSE,fig.height=3.5,fig.width=8.5>>=
print(qplot(Girth,Volume,data=trees,xlab="Girth of tree (inches)",
            ylab=expression("Volume of usable timber (ft"^3*")")) +
            geom_smooth(method="lm",se=FALSE),
      TRUE, vp=viewport(0.25,0.5,0.5,1.0))
print(qplot(Girth,Volume,data=trees,log="xy",xlab="Girth of tree (inches)",
            ylab=expression("Volume of usable timber (ft"^3*")")) +
            geom_smooth(method="lm",se=FALSE),
      FALSE, vp=viewport(0.75,0.5,0.5,1.0))
@
\caption{Scatterplot of the volume of usable lumber versus the girth
  of the tree for 31 black cherry trees.  The left panel is on the
  original scale.  The right panel is on the log-log scale.}
  \label{fig:trees}
\end{figure}
we show the volume versus the girth on the original scale and on a
log-log scale.  There is not a tremendous difference in the patterns
but careful examination shows better linear behavior in the log-log
scale.

Our initial model is
<<lm2show,eval=FALSE>>=
coef(summary(lm2 <- lm(log(Volume) ~ log(Girth), trees)))
@
<<lm2,echo=FALSE>>=
printCoefmat(coef(summary(lm2 <- lm(log(Volume) ~ log(Girth), trees))))
@
To fit a model corresponding to a conical or cylindrical shape
we add a term in \code{log(Height)} (recall that we are on the log-log scale)
<<lm2ashow,eval=FALSE>>=
coef(summary(lm2a <- lm(log(Volume) ~ log(Girth) + log(Height), trees)))
@
<<lm2a,echo=FALSE>>=
printCoefmat(coef(summary(lm2a <- lm(log(Volume) ~ log(Girth) + log(Height), trees))))
@

\paragraph{Testing specific combinations of parameters}

At this point we may want to check if a version of the formula for the
volume of a cylinder or of a cone, both of which have the form
\begin{displaymath}
  V=k\,d^2h
\end{displaymath}
where $k$ is a constant, $d$ is the diameter (or, equivalently, the
girth or circumference at the base) and $h$ is the height.  Such an
expression would correspond to a value of 2 for the \code{log(Girth)}
coefficient and 1 for the \code{log(Height)} term.  The
\code{log(Height)} term is highly significant and the coefficients of
\code{log(Girth)} and \code{log(Height)} are reasonably close to 2 and
1.  In particular, confidence intervals on these coefficients include
2 and 1
<<lm2aconf>>=
confint(lm2a)
@

The confidence intervals do not, by themselves, answer the question of
whether a model of the form
\begin{displaymath}
  \log(\text{Volume}_i)=\beta_0+2\log(\text{Girth}_i)+\log(\text{Height}_i)+\epsilon_i,
  \quad i=1,\dots,31
\end{displaymath}
is a reasonable fit.  To fit this model we use an \code{offset}
expression in the model formula.
<<lm2c>>=
lm2c <- lm(log(Volume) ~ 1 + offset(2*log(Girth) + log(Height)), trees)
@
\noindent and perform a comparative analysis of variance
<<anovalm2>>=
anova(lm2c,lm2a)
@
The large p-value indicates that the more complex model (general
values of the coefficients for \code{log(Girth)} and
\code{log(Height)}) does not fit significantly better than the simpler
model (assuming \code{2*log(Girth) + log(Height)}), thus we prefer the
simpler model.

\paragraph{One-way analysis of variance}

Next consider the \code{InsectSprays} data with a response, \code{count},
related to a categorical covariate, \code{spray}.  Comparative
boxplots (Fig.~\ref{fig:sprays})
\begin{figure}[tb]
  \centering
<<spraysplot,echo=FALSE,fig.height=4.5,fig.width=8.5>>=
print(qplot(reorder(spray,count), count, data=InsectSprays, geom="boxplot",
            xlab="Spray", ylab="Insect count") + coord_flip(),
      TRUE, vp=viewport(0.5,0.75,1,0.5))
print(qplot(reorder(spray,count), sqrt(count), data=InsectSprays, geom="boxplot",
            xlab="Spray", ylab=expression(sqrt("count"))) + coord_flip(),
      FALSE, vp=viewport(0.5,0.25,1,0.5))
@
\caption{Comparative boxplots of the insect count by spray type in the
  \code{InsectSprays} data.  The sprays have been reordered according
  to increasing mean response.  In the lower panel the response is the
  square root of the count.}
  \label{fig:sprays}
\end{figure}
show that the square root of the count is a more reasonable scale for
the response and that there is considerable differences in the
response according to the spray type.

Althought we can fit a model with categorical covariates using the
\code{lm()} function, there is an advantage in using the \code{aov()}
function instead, because it allows us to extract some additional
information that applies only to categorical factors.  Also,
\code{summary()} applied to an \code{aov()} model produces the
analysis of variance table, which for such models, is more interesting
than the coefficients table.
<<av1>>=
summary(av1 <- aov(sqrt(count) ~ spray, InsectSprays))
@
\noindent If we want to express the model in terms of ``effects'' we can obtain these as
<<av1eff>>=
model.tables(av1)
@
\noindent or, if we are interested in the estimates of the means for each group,
<<av1eff2>>=
model.tables(av1, type="means")
@
\noindent Various types of ``multiple comparisons'' methods are also available.
We will discuss these later.

\paragraph{Multi-factor analysis of variance}

When we have more than one categorical covariate, as in the
\code{OrchardSprays} data,
<<strOrchard>>=
str(OrchardSprays)
@
\noindent we simply include them in the model formula.  It happens that for
this experiment there are two blocking factors, \code{rowpos} and
\code{colpos}, and one experimental factor, \code{treatment}, so we
put the blocking factors first.
<<av2>>=
summary(av2 <- aov(decrease ~ factor(rowpos) + factor(colpos) + treatment, OrchardSprays))
@
\noindent These data are arranged in what is called a ``Latin square'' design,
which is a special type of fractional replication. There are 64
observations on three factors, each at 8 levels, so not only are there
no replications, we don't even have an observation in each of the
possible $8\times8\times8=512$ combinations, and cannot try to fit a
model with interaction terms.

\paragraph{Multi-factor anova with replications}

The \code{warpbreaks} data, shown in Fig.~\ref{fig:warpbreaks},
<<warpbreaks,echo=FALSE,fig.height=4.5,fig.width=8.5,fig.align='center',fig.pos='tb',fig.cap='Comparative boxplots of the number of warp breaks per loom by tension setting for the warpbreaks data.  Panels are determined by wool type.  The upper two panels are on the original scale of number of breaks.  The lower two panels are on the reciprocal scale (i.e. number of looms per break).'>>=
warpbreaks <- within(warpbreaks,invbreaks <- 1/breaks)
print(qplot(tension, breaks, data=warpbreaks, geom="boxplot",
            xlab="Tension setting", ylab="Breaks per loom") +
      coord_flip() + facet_grid(wool~.),
      TRUE, vp=viewport(0.5,0.75,1,0.5))
print(qplot(tension, invbreaks, data=warpbreaks, geom="boxplot",
            xlab="Tension setting", ylab="Inverse breaks per loom") +
      coord_flip() + facet_grid(wool~.),
      FALSE, vp=viewport(0.5,0.25,1,0.5))
@
\noindent are counts of the number of warp breaks per loom (a length of wool)
according to the tension setting for the wool and the type of wool.
We see that on the original scale of the number of breaks per loom
there is increasing variance with an increasing level of the response,
whereas on the reciprocal scale (number of looms per break) the
variability is much closer to being constant.

Because there are 9 replications at each of the wool/tension combinations
<<xtabsbreaks>>=
xtabs(~ wool + tension, warpbreaks)
@
\noindent we can fit a model with main effects for \code{wool} and for
\code{tension} and the \code{wool:tension} interaction.
<<av3>>=
summary(av3 <- aov(breaks ~ wool * tension, warpbreaks))
@
\noindent In this model the interaction is significant.  When an interaction is
significant we typically retain both of the main effects in the model.

However, if we fit the model on the reciprocal scale
<<av3a>>=
summary(av3a <- aov(1/breaks ~ wool * tension, warpbreaks))
@
\noindent we no longer have a significant interaction and could reduce the model
to the main effects only
<<av3b>>=
summary(av3b <- aov(1/breaks ~ tension +  wool, warpbreaks))
@

Here we have reordered the factors \code{tension} and \code{wool} so
that \code{wool} is the last term and thus the second row of the
analysis of variance table corresponds to a test of the main effect of
the wool given that tension had been taken into account.  (If you look
closely at the sums of squares, degrees of freedom and mean squares
you will see that they are consistent in models \code{av3b} and
\code{av3a} but that is a consequence of the data being completely
balanced with respect to these factors.  To be safe, always make the
factor you are going to test be the last one in the model formula.)
The \code{wool} factor is not significant and we can reduce the model
to a single factor model
<<av3c>>=
summary(av3c <- aov(1/breaks ~ tension, warpbreaks))
@
\noindent corresponding to Fig.~\ref{fig:warpbreaks2}
<<breaksplot2,echo=FALSE,fig.height=1.5,fig.width=8.5,fig.align='center',fig.pos='tb',fig.cap='Comparative boxplots of the reciprocal of the number of warp breaks per loom by tension setting for the warpbreaks data.'>>=
print(qplot(tension, invbreaks, data=warpbreaks, geom="boxplot",
            xlab="Tension setting", ylab="Reciprocal of breaks per loom") +
       coord_flip())
@
\noindent in which we can see a trend across the three ordered levels of
tension; low tension gives a low reciprocal number of breaks
(corresponding to a higher frequency of breaks), medium tension gives
an intermediate reciprocal number and high tension gives the highest
reciprocal number.

This is a common situation with a factor like \code{tension} whose
levels are in a natural ordering, \code{L < M < H}.  Details will be
given later but, for now, it is enough to see that if we convert the
factor to an ordered factor
<<orderedtension>>=
str(warpbreaks <- within(warpbreaks, tension <- ordered(tension)))
@
\noindent and fit the model as before,
<<av3d>>=
summary(av3d <- aov(1/breaks ~ tension, warpbreaks))
@
\noindent we get the same analysis of variance table but now the two degrees of
freedom for \code{tension} are divided into a linear trend and a
quadratic relationship in addition to the linear trend
<<av3dcoefshow,eval=FALSE>>=
coef(summary.lm(av3d))
@
<<av3dcoef,echo=FALSE>>=
printCoefmat(coef(summary.lm(av3d)))
@
\noindent With a p-value of 98.4\%, the quadratic term is not at all
significant, indicating that we could reduce to only the linear trend.

\paragraph{Analysis of covariance models}

The data shown in Fig.~\ref{fig:housing}
\begin{figure}[tb]
  \centering
  \setkeys{Gin}{width=0.6\textwidth}
<<housing,echo=FALSE,fig.height=6,fig.width=8>>=
DataURL <- "http://www.stat.wisc.edu/~bates/data/"
assessed <- read.table(paste(DataURL, "assessedval.txt", sep=''),
                       header=FALSE)
names(assessed) <- c("Selling", "Assessed", "Corner")
assessed <- within(assessed, Corner <- factor(Corner, labels=c("N","Y")))
print(qplot(Assessed, Selling, data=assessed, color=Corner,
            linetype=Corner, shape=Corner)+geom_smooth(method="lm",se=FALSE))
@
  \caption{Selling price versus assessed value for a random sample of 16 single-family homes located on corner lots and 48 single-family homes not on corner lots.}
  \label{fig:housing}
\end{figure}
\noindent are from a study (conducted several
years ago) described as
\begin{quote}
  A tax consultant studied the current relation between
  selling price and assessed valuation of single-family residential
  dwellings in a large tax district by obtaining data for a random
  sample of 16 "arm's length" sales transactions of single-family
  dwellings located on corner lots and for a random sample of 48
  recent sales of single-family dwellings not located on corner lots.
\end{quote}
Assuming a linear relationship between the selling price and the
assessed value, there are three different models we would typically consider:
<<eval=FALSE>>=
coef(summary(lm3 <- lm(Selling ~ 1 + Assessed, assessed)))
@
<<lm3,echo=FALSE>>=
printCoefmat(coef(summary(lm3 <- lm(Selling ~ Assessed, assessed))))
@
\noindent in which there is no distinction between corner lots and non-corner lots,
<<eval=FALSE>>=
coef(summary(lm3a <- lm(Selling ~ 1 + Corner + Assessed, assessed)))
@
<<lm3a,echo=FALSE>>=
printCoefmat(coef(summary(lm3a <- lm(Selling ~ 1 + Corner + Assessed, assessed))))
@
\noindent in which there is a common slope but a different intercept for the corner lots, and
<<eval=FALSE>>=
coef(summary(lm3b <- lm(Selling ~ 1 + Corner * Assessed, assessed)))
@
<<lm3b,echo=FALSE>>=
printCoefmat(coef(summary(lm3b <- lm(Selling ~ 1 + Corner * Assessed, assessed))))
@
\noindent in which there are different intercepts and slopes for the corner lots.

We could use a comparative analysis of variance to compare a simpler
model ($H_0$) versus a more complex model ($H_a$).
<<anovaab>>=
anova(lm3, lm3a)
anova(lm3a, lm3b)
@

Because these F tests have one numerator degree of freedom there will
be a corresponding t-test for a coefficient.  This will be the t-test
on the coefficient in $H_a$ that distinguishes it
from $H_0$.

\section{Classes and methods for linear models}
\label{sec:classes}

A model fit with \code{lm()} has class \code{"lm"}
<<classlm2>>=
class(lm2)
@
\noindent for which there are several methods defined
<<methodslm>>=
methods(class="lm")
@
\noindent We have already seen several of these in use:
\begin{description}
\item[anova] Return the (sequential) analysis of variance table for a
  single fitted model or a comparative analysis of variance for
  multiple fitted models.
\item[confint] Return confidence intervals on the coefficients
\item[deviance] Return the residual sum of squares (RSS) for the
  model. (This is a misnomer because the RSS is related to but not
  exactly the same as the deviance.)
\item[formula] Return the model formula.
\item[kappa] Return the condition number of the model matrix or an
  upper bound on its condition number.
\item[logLik] Return the value of the log-likelihood at the estimated
  parameter values.
\item[model.frame] Return the model frame to which the model was
  actually fit.
\item[model.matrix] Return the model matrix.
\item[plot] Produce some common residual plots for evaluating the
  model fit.
\item[predict] Returns evaluations of the fitted model, and optionally
  their standard errors, at the observed or newly specified values of
  the covariates.
\item[residuals] Returns the residuals from the fit.
\item[rstandard] Returns the ``standardized residuals'' (to be
  described later).
\item[rstandard] Returns the ``Studentized residuals'' (to be
  described later).
\item[simulate] Return a matrix of simulated response vectors
  according to the model assuming that the fitted values of the
  parameters are the true parameter values.
\item[summary] Return a summary of the fitted model
\item[vcov] Return the (estimated) variance-covariance matrix of
  $\widehat{\bm\beta}$ (i.e. the matrix that could be expressed as
  $s^2(\bm X\trans\bm X)^{-1}$).
\end{description}

Other extractor functions such as \code{coef} and \code{fitted} do not
have specific methods for class \code{"lm"} but instead apply the
default method to objects of this class.

A model fit by \code{aov} has class
<<classaov>>=
class(av3)
@
\noindent \code{"aov"} and also class \code{"lm"}.  This means that methods
for class \code{"aov"} will be chosen, if they exist, otherwise
methods for class \code{"lm"} and, finally, the default method.

Specific methods for class \code{"aov"} are
<<methodsaov>>=
methods(class="aov")
@
\noindent from which we can see that specific methods for the \code{coef},
\code{extractAIC}, \code{print}, \code{proj} and \code{summary}
generics are available for this class and will be chosen in preference
to the \code{"lm"} or default method.  Furthermore there are specific
methods for the \code{model.tables}, \code{se.contrast} and
\code{TukeyHSD} generics.

\section{Simulating linear model fits}
\label{sec:simulating}

It is possible to simulate a large number of replications of linear
model fits quite quickly, if you do it carefully.  The trick is to
realize that, if the expression on the left-hand size of the model
formula in \code{lm()} or \code{aov()} is a matrix with $n$ rows and
$N$ columns then the model is fit to all $N$ of the response vectors
simultaneously.

As shown above, one of the methods that can be applied to a fitted
model is called \code{simulate} and it generates such a matrix using
the parameter estimates (both the coefficients and the variance,
$\sigma^2$ of the ``random noise'') for the simulation.  Consider
again the model \code{lm1} with coefficient table
<<lm1revisited,echo=FALSE>>=
printCoefmat(coef(summary(lm1)))
@
\noindent To simulate 10,000 response vectors simulated from this model is
fast
<<lm1sim>>=
set.seed(1234321)
system.time(Ylst <- simulate(lm1, 10000))
str(Ylst, 0)
@
\noindent but, unfortunately, it is not a matrix, which is what we want.  To get
a matrix from a data frame we use
<<Ymat>>=
str(Ymat <- data.matrix(unname(Ylst)))
@
\noindent (the reason for using \code{unname()} is to get rid of the 10,000 column names).

Finally, we fit the model to all 10,000 response vectors simultaneously.
<<lm1sim2>>=
system.time(lm1sim <- lm(Ymat ~ 1 + carb, Formaldehyde))
@
\noindent which is much faster than you could ever hope to do with a loop.

Now the coefficients are in the form of a 2 by 10,000 matrix
<<simcoef>>=
str(simcoef <- coef(lm1sim))
@
\noindent for which we could produce density plots, etc.

To obtain a density plot of the estimated intercepts we use
<<intercepts>>=
intercepts <- simcoef[1,]
slopes <- simcoef[2,]
@
\noindent from which we produce the plots in Fig.~\ref{fig:dens}
\begin{figure}[htb]
  \centering
<<coefdens,echo=FALSE,fig.height=4.5>>=
print(qplot(intercepts,geom="density"), TRUE, vp=viewport(0.25,0.5,0.5,1))
print(qplot(slopes, geom="density"), FALSE, vp=viewport(0.75,0.5,0.5,1))
@
  \caption{Estimated intercepts and slopes from data simulated according to model \code{lm1}}
  \label{fig:dens}
\end{figure}

We can also check that the mean of the estimated intercept is close to
the value used in the simulation (i.e. the estimated intercept in
model \code{lm1}) and standard deviation of the simulated intercepts
is close to the standard error of the that estimate.
<<meanstd>>=
c(mean=mean(intercepts), sd=sd(intercepts))
c(mean=mean(slopes), sd=sd(slopes))
printCoefmat(coef(summary(lm1)))
@

\end{document}