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

\let\ul\underline
\newcommand{\eq}[1]{(\ref{#1})}
\long\def\comment#1{}

\begin{document}

\begin{center}\Large\bf
One-way Analysis of Variance Models
\end{center}

\section{ANOVA Tables and ANOVA Models}

In this part of the course we will learn about ANOVA models, and their
application to designed experiments.  ``ANOVA'' is one of the most
confusing terms in applied statistics.  The term ``Analysis of
variance (ANOVA)'' is used to describe both
\begin{itemize}
\item the ANOVA table, which summarizes sums-of-squares and
degree-of-freedom information for any linear regression; and  
\item ANOVA models, which are a stylized way of computing and
interpreting regressions with dummy variables to learn about
differences in group means for grouped data.
\end{itemize}

\subsection{A bit more about ANOVA Tables}
You already know about the ANOVA table used to describe a multiple
regression. 
\begin{equation}
        Y = \alpha + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_3 +
        \epsilon
        \label{generic-regression}
\end{equation}
To simplify the writing a bit, we will write
\begin{eqnarray*}
        SS(X_1,X_2,X_3) & = & RegSS, \mbox{ for the regression {\tt Y \~{}
        X$_1$ + X$_2$ + X$_3$} in \eq{generic-regression} above; } \\
        RSS(X_1,X_2,X_3) & = & RSS, \mbox{ for the same regression;}
\end{eqnarray*}
and similarly for other expressions such as $SS(X_1)=RegSS$ for the
regression {\tt Y \~{} X$_1$}, $SS(X_1,X_2) = RegSS$,  for the
regression {\tt Y \~{} X$_1$ + X$_2$}, etc.  We will also write
\begin{eqnarray*}
        SS(X_3|X_2,X_1) & = & SS(X_1,X_2,X_3) - SS(X_1,X_2) \\  
        SS(X_2|X_1)     & = & SS(X_1,X_2) - SS(X_1)
\end{eqnarray*}
and so forth (for other subscript combinations).  Since 
\begin{equation}
        TSS = SS(X_1, \ldots, X_k) + RSS(X_1, \ldots, X_k) 
        \label{TSS}
\end{equation}
for least-squares regression, and we know for example that
$RSS(X_1,X_2,X_3) \leq RSS(X_1,X_2)$, we know that $SS(X_1,X_2,X_3)
\geq SS(X_1,X_2)$; thus expressions like $SS(X_3|X_1,X_2)$ above are
guaranteed to be nonnegative.  By definition, the quantity
$SS(X_3|X_1,X_2)$ represents the incremental variation (in $Y$)
explained by regressing on $X_3$ after having already regressed on
$X_1$ and $X_2$.  Because of \eq{TSS}, of course, we can also express
each $SS(\cdot|\cdots)$ expression as a difference in RSS's.

Sometimes it is more convenient to refer to the regression
coefficients $\beta_1$, $\beta_2$, etc. than the regressors $X_1$,
$X_2$, etc., so we will also denote
\begin{eqnarray*}
        SS(\beta_1,\beta_2,\beta_3) & = & SS(X_1,X_2,X_3) \\
        SS(\beta_3|\beta_2,\beta_1) & = & SS(X_3|X_2,X_1)
\end{eqnarray*}
and so on.

The standard ANOVA table that summarizes sums of squares for the
regression \eq{generic-regression} can then be written as follows:
\[
\begin{array}{lr@{=}rrr}
\em Source & \multicolumn{2}{c}{\em SS} & \em df \hfil\hbox{}& \em MS\hfil\hbox{} \\
\hline
\\
X_1            & SS(X_1)         & 55443.5 & 1 & 55443.5 \\
X_2|X_1        & SS(X_2|X_1)     & 167927.2 & 1 & 167927.2 \\
X_3|X_1,X_2    & SS(X_3|X_1,X_2) &  78792.9 & 1 & 78792.9 \\
\\
\hline
\\
\rm Regression & SS(X_1,X_2,X_3) & 302163.4 & 3 & 100721.1 \\ 
\\
\hline
\\
\rm Error      & RSS(X_1,X_2,X_3) & 339258.2 & n-4 = 46 & 7375.178 \\
\\
\hline
\\
\rm Total      & TSS              & 641421.6 & n-1 = 49 & 
\end{array}
\]
which we can read off almost directly from the SPLUS output
\begin{verbatim}
> fuel _ read.table("../week04/fuel.dat")
> anova(lm(Fuel ~ Tax + Dlic + Inc,data=fuel))
Analysis of Variance Table

Response: Fuel

Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value       Pr(F) 
      Tax  1   55443.5  55443.5  7.51758 0.008674482
     Dlic  1  167927.2 167927.2 22.76924 0.000018830
      Inc  1   78792.9  78792.9 10.68352 0.002049453
Residuals 46  339258.2   7375.2                     
\end{verbatim}
You can verify that the definitions of $SS(Inc|Tax,Dlic)$, and so
forth, work as indicated above, by separately fitting and examining
the RSS's from the models
\begin{verbatim}
Fuel ~ Tax
Fuel ~ Tax + Dlic
Fuel ~ Tax + Dlic + Inc
\end{verbatim}

This sort of ANOVA table is useful in all regression analysis, and in
particualr it will be useful as we build and assess the stylized
regression models called ``ANOVA models''.

\subsection{ANOVA vs Linear Regression}

Someone who prefers thinking about regression models (like me) will
often say, {\em ``oh well, ANOVA, that's just multiple regression
where all the predictor variables are dummy variables!''.}  This is
mathematically correct (and we will exploit this idea often!) but
ANOVA modeling is more than just running {\tt lm()} with a bunch
of dummy variables for $X's$.  ANOVA models involve:
\begin{itemize}
\item Particular parameterizations of the linear regression model that
help us to understand relationships between the categories of discrete
regressors; and
\item Special arrangements of computing the least-squares estimates
and the sums-of-squares in the ANOVA table.  These arrangements
emphasize the idea that when we are comparing group means in an ANOVA
model we are comparing ``between'' and ``within'' variance estimates;
and they can be much more efficient (for large data sets, e.g.) than
the standard linear regression calculations.
\end{itemize}

\subsection{ANOVA}

Analysis of variance models (ANOVA models) are used to study the
relationship between predictor variables that {\em are not}
necessarily quantitative and response variables that {\em are}
quantitative.  They are structured like linear regression models, and
have the same basic assumptions, but they do not assume a linear
relationship between the predictor and response variables.  We've also
seen some examples of regressing on dummy variables for categorical
predictors.  Here I want to introduce the idea that such regressions
can also be used to explore nonlinear relationships between {\em
continuous} $X$'s and $Y$'s.

{\bf Example.}

Suppose we want to study the relationship between $X$ $=$ {\em the
price of a product}, and $Y$ $=$ {\em the sales volume in one month}.

As a linear regression problem, this is simple.  We want to fit the
model
\[
        Y = \alpha + \beta_1 X + \epsilon
\]
But this explicitly assumes that the relationship between price and
sales is linear, and this need not be the case, as shown in
Figure~\ref{salesdata} figure (for about 70 months' worth of data).
\comment{
> salesdata _ read.table("salesdata.dat")
> dim(salesdata)
[1] 70  2
> salesdata[1:4,]
  price    sales 
1    10 7826.493
2    20 7094.763
3    30 6669.467
4    40 6285.748
> attach(salesdata)
> plot(price,sales)
> points(20,mean(sales[price==20]),pch=16,col=2,cex=2)
> points(50,mean(sales[price==50]),pch=16,col=2,cex=2)
> points(90,mean(sales[price==90]),pch=16,col=2,cex=2)
}

\begin{figure}
\hfil\psfig{figure=nonlinear.ps,width=5.5in,height=3.5in}
\caption{Some fictional sales data.  We could fit the linear
regression $y = \alpha + \beta X + \epsilon$ to all the data, but
instead we will focus on the particular price points \$20.00, \$50.00,
and \$90.00.  (In real life, if one knew one were only going to focus on these
prices, of course one would not bother to collect data at the other
price points.)}\label{salesdata}
\end{figure}

\begin{figure}
\hfil\psfig{figure=plotfactorsales.ps,width=5.5in,height=3.5in}
\caption{The result of {\tt plotfactor()} is always some sort of
parallel boxplot display.}
\end{figure}

An alternative to fitting the complex nonparametric regression models
we talked about at the beginning of the semester is to choose a few
price points, such as \$20.00, \$50.00, and \$90.00, and define dummy
variables
\[
X_1 = \left\{\begin{array}{rr}1 & \mbox{if price$=$\$20.00} \\
                              0 & \mbox{else} \end{array}\right.
~~~~~~~~~~~~
X_2 = \left\{\begin{array}{rr}1 & \mbox{if price$=$\$50.00} \\
                              0 & \mbox{else} \end{array}\right.
~~~~~~~~~~~~
X_3 = \left\{\begin{array}{rr}1 & \mbox{if price$=$\$90.00} \\
                              0 & \mbox{else} \end{array}\right.
\]
and then fit the regression model
\begin{equation}
        Y = \alpha + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_3 +
        \epsilon
        \label{naive}
\end{equation}
As we have seen before, this is an example of a {\em ``one-way ANOVA''
model}.

There are three things to note about this approach:
\begin{itemize}
\item Even though the model is still the familiar linear regression
model, it is possible to detect and explore nonlinearity in the
relationship between $X$ and $Y$, by looking at just these three price
levels, or perhaps by adding two or three more price levels to get a
better look at the curve.  This can be much more efficient than
fitting a polynomial or nonlinear regression model.
\item It may also be more efficient to collect data for the one-way
ANOVA model than for the general regression model; maybe we can get 5
months' worth of data at each of three price points, and have the data
collected in a year or so. It is much harder to say how to collect
data for this problem if all possible prices are allowed.
\item\comment{ In matrix form this (naive form of the) ANOVA model is
\[
\ul y = \left[
\begin{array}{c}
y_{11} \\
y_{21} \\
\vdots \\
y_{n_11} \\
y_{12} \\
y_{22} \\
\vdots \\
y_{n_22} \\
y_{13} \\
y_{23} \\
\vdots \\
y_{n_33} 
\end{array}
\right]
= X\ul\beta + \ul \epsilon= 
\left[
\begin{array}{cccc}
1 & 1 & 0 & 0 \\
1 & 1 & 0 & 0 \\
  & \vdots & & \\
1 & 1 & 0 & 0 \\
1 & 0 & 1 & 0 \\
1 & 0 & 1 & 0 \\
  & & \vdots & \\
1 & 0 & 1 & 0 \\
1 & 0 & 0 & 1 \\
1 & 0 & 0 & 1 \\
  & & & \vdots  \\
1 & 0 & 0 & 1 
\end{array}
\right]
\left[
\begin{array}{c}
\alpha \\
\beta_1 \\
\beta_2 \\
\beta_3
\end{array}
\right]
+
\left[
\begin{array}{c}
\epsilon_{11} \\
\epsilon_{21} \\
\vdots \\
\epsilon_{n_11} \\
\epsilon_{12} \\
\epsilon_{22} \\
\vdots \\
\epsilon_{n_22} \\
\epsilon_{13} \\
\epsilon_{23} \\
\vdots \\
\epsilon_{n_33}
\end{array}
\right]
\]
showing that the model is overparameterized (linearly dependent
columns).  
}
%{\em 
All ANOVA models are overparametrized in their
``natural'' form, and we often must make a choice about what a
sensible ``reduced'' parametrization is.
%}
~$\Box$
\end{itemize}

\section{One-way ANOVA analysis in SPLUS}

We can illustrate the one-way analysis of variance using the sales data:
\begin{verbatim}
> attach(salesdata)
> sales205090 _ salesdata[((price==20)|(price==50)|(price==90)),]
> sales205090$price _ factor(sales205090$price)
> table(sales205090$price)  # There are 11 months of sales at each price
 20 50 90 
 11 11 11

> view(sales205090$price)

      structure(.Data = c(1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3,
      1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3), .Label =
      c("20", "50", "90"), class = "factor")

> plot.factor(sales205090,style.bxp="old")
\end{verbatim}


There is a special function, {\tt aov()} for fitting analysis of
variance models that Splus can apply to data where the regressors are
discrete factors.  {\tt aov()} can calculate a regression more
efficiently in some cases than {\tt lm()}, and {\tt aov()} retains
information that is useful for ANOVA analyses, that {\tt lm()} does
not.  Let's compare the {\tt aov()} output with the usual {\tt lm()}
output; see Table~\ref{comparison}.

\begin{table}
\begin{tabular}{@{}p{3.25in}@{}|p{3.25in}@{}}
\small
\begin{verbatim}
> sales.aov _ aov(sales ~ price,
+ data=sales205090)
> sales.aov
Call:
   aov(formula = sales ~ price,
      data = sales205090)

Terms:
                   price Residuals 
 Sum of Squares 13332276    223371
Deg. of Freedom        2        30

Residual standard error: 86.28854 
Estimated effects are balanced
> summary(sales.aov)
         Df Sum of Sq Mean Sq  F Value Pr(F) 
   price   2 13332276 6666138 895.2989     0
Residuals 30   223371    7446               
 
> anova(sales.aov)
Analysis of Variance Table

Response: sales

Terms added sequentially (first to last)
         Df Sum of Sq Mean Sq  F Value Pr(F) 
   price  2  13332276 6666138 895.2989     0
Residuals 30   223371    7446  
>  13332276/( 13332276 +  223371 )  
+ # calculate R^2
[1] 0.9835219
\end{verbatim}
&
\small
\begin{verbatim}
>  sales.lm _ lm(sales ~ price,data=sales205090)
> sales.lm
Call:
lm(formula = sales ~ price, data = sales205090)

Coefficients:
 (Intercept)    price1    price2 
    6184.607 -466.4094 -359.8494

Degrees of freedom: 33 total; 30 residual
Residual standard error: 86.28854 
> summary(sales.lm)

Call: lm(formula = sales ~ price, data = sales205090)
Residuals:
    Min     1Q Median    3Q   Max 
 -194.8 -55.67 -4.224 78.53 113.8

Coefficients:
                 Value Std. Error    t value   Pr(>|t|) 
(Intercept)  6184.6074    15.0209   411.7333     0.0000
     price1  -466.4094    18.3968   -25.3528     0.0000
     price2  -359.8494    10.6214   -33.8797     0.0000

Residual standard error: 86.29 on 30 degrees of freedom
Multiple R-Squared: 0.9835 
F-statistic: 895.3 on 2 and 30 degrees of 
freedom, the p-value is 0 

Correlation of Coefficients:
       (Intercept) price1 
price1 0                 
price2 0           0     
> anova(sales.lm)
Analysis of Variance Table

Response: sales

Terms added sequentially (first to last)
          Df Sum of Sq Mean Sq  F Value Pr(F) 
    price  2  13332276 6666138 895.2989     0
Residuals 30    223371    7446               
\end{verbatim}
\\
\end{tabular} 
\caption{Comparison of {\tt aov()} and {\tt lm()} output for the
fictional sales data.}\label{comparison}
\end{table}

This tells us that the ``price'' variable explains much of the
variability in the ``sales'' data, but isn't very informative about
details, like what the mean sales was at each price level, etc. 

\subsection{The cell-means parametrization}

As we have seen before, to calculate the mean sales at each price we
want to think of a model like this:
\[
        y_{ij} = \mu_j + \epsilon_{ij}
\]
where $j=1, 2$ or 3, depending on whether the price was 20, 50 or 90,
and $i$ indexes the 11 observations at each price level\footnote{In
general in ANOVA models we write multiple subscripts like this to help
organize our view of the data into observations ($j$) within cells
($i$) In more complicated models we may have to change from two to
three, or more, subscripts.}.  This is called the {\em cell means}
ANOVA parametrization; and it is equivalent to the regression model
\[
        y_i = \mu_1 X_1 + \mu_2 X_2 + \mu_3 X_3 + \epsilon_i
\]
where now $i$ indexes all the data, and the $X_j$'s indicate whether
the price was 20, 50 or 90.

There is a function, {\tt model.tables()}, in Splus, that will
calculate the $\mu_i$'s directly for us:

\vspace*{\fill}

\clearpage

\begin{verbatim}
> model.tables(sales.aov,type="means")
Refitting model to allow projection
Tables of means
Grand mean
      
 6185

 price 
   20   50   90 
 7011 6078 5465
\end{verbatim}
Let's compare these results with fitting the model ``by hand'' using
{\tt lm()}.  
\begin{verbatim}
> lmsales _ sales205090$sales
> lmprice _ sales205090$price
> p20 _ ifelse(lmprice==20,1,0)
> p50 _ ifelse(lmprice==50,1,0)
> p90 _ ifelse(lmprice==90,1,0)

> sales.lm _ lm(lmsales ~ p20 + p50 + p90 - 1)
> summary(sales.lm)

Call: lm(formula = lmsales ~ p20 + p50 + p90 - 1)
Residuals:
    Min     1Q Median    3Q   Max 
 -194.8 -55.67 -4.224 78.53 113.8

Coefficients:
        Value Std. Error   t value  Pr(>|t|) 
p20 7010.8662   26.0170   269.4728    0.0000
p50 6078.0475   26.0170   233.6185    0.0000
p90 5464.9085   26.0170   210.0517    0.0000

Residual standard error: 86.29 on 30 degrees of freedom
Multiple R-Squared: 0.9998 
F-statistic: 57100 on 3 and 30 degrees of freedom, the p-value is 0 

Correlation of Coefficients:
    p20 p50 
p50 0      
p90 0   0  
\end{verbatim}
Actually, as we know, when {\tt lm()} or {\tt aov()} sees a {\tt
factor} variable in a regression, it will construct dummy variables
automatically, but by default it produces a coding that is difficult
for humans to interpret (this coding also makes the regressors {\em
orthogonal}---so like our by-hand analysis the correlations between
the regression coefficients will be zero---but there is more than one
way to do this).  We will learn more about Splus's coding of dummy
variables later.

\subsection{The grand-mean-plus-effects parametrization}

Often we are not directly interested in the mean of each cell in an
ANOVA model, but rather in whether there are substantial differences
between the cells.  An alternative parametrization focuses the
attention away from the overall mean of the data and toward deviations
from the overall mean in each cell.  
In ANOVA notation, we would write the model \eq{naive} as
\[
        y_{ij} = \mu + \alpha_j + \epsilon_{ij}
\]
where again $j=1, 2$ or 3, depending on whether the price was 20, 50
or 90, and $i$ indexes the 11 observations at each price level.

{\em Nomenclature:}
\begin{itemize}
\item $\mu$ is called the {\em overall mean} or {\em grand mean}
\item $\alpha_1$, $\alpha_2$ and $\alpha_3$ are called ``effects''.
Just to be confusing, they are also called ``contrasts.''
\end{itemize}
Note that there are four parameters but only three groups of
data---cells---from which to estimate them.

In the past we have fixed this problem by dropping a regressor, which
is equivalent to adding the linear constraint ``(that coefficient) $=$
0''.  More general linear constraints, of the form
\[
        w_0\mu + \sum_{j=1}^{\#cells} w_j \alpha_j = 0 ,
\]
where $w_0, w_1, \ldots$ are fixed, known weights, will also fix the
problem.  The usual linear constraint for the ``grand means plus
effects'' model is to make the effects sum to zero:
\[
        \sum_{j=1}^{\#cells} \alpha_j = 0
\]
(so $w_0=0$, and the other $w_j$'s $=$ 1).  The $\alpha's$ can again
be estimated with the {\tt model.tables()} function:

\begin{verbatim}
> model.tables(sales.aov,type="effects")
Refitting model to allow projection
Tables of effects

 price 
    20     50     90 
 826.3 -106.6 -719.7
> sum(scan())        # verify that the alpha's sum to zero
1:  826.3 -106.6 -719.7
4: 
[1] -1.136868e-013
\end{verbatim}

How can we think about what {\tt model.tables()} is doing for the
``grand mean plus effects'' parametrization?  

It is generally quite useful to start with the ``cell means''
parametrization expressed as a {\em linear regression} model, impose
the new parameterization on it, and see what linear regression is
required to estimate the parameters in the ANOVA model.

We begin by computing the cell means in the grand mean plis effects
model, and noting the linear restriction:
\begin{eqnarray*}
        \mu_1 & = & \mu + \alpha_1 \\
        \mu_2 & = & \mu + \alpha_2 \\
        \mu_3 & = & \mu + \alpha_3 \mbox{, and}\\
        0     & = & \alpha_1 + \alpha_2 + \alpha_3
\end{eqnarray*}
Now we write out the regression model for cell means, substitute the
values for the cell means above, and use the linear constraint:
\begin{eqnarray*}
        y       & = & \mu_1 X_1 + \mu_2 X_2 + \mu_3 X_3 \\
                & = & (\mu + \alpha_1)X_1 + (\mu + \alpha_2) X_2 + 
                        (\mu + \alpha_3) X_3 + \epsilon \\
                & = & (\mu + \alpha_1)X_1 + (\mu + \alpha_2) X_2 + 
                        (\mu - \alpha_1 - \alpha_2) X_3 + \epsilon \\
                & =&  \mu + \alpha_1(X_1-X_3) + \alpha_2(X_2-X_3) 
                        + \epsilon
\end{eqnarray*}
This suggests that to estimate the effects in the grand mean plus
effects model we should construct new dummies
\begin{eqnarray*}
        D_1 & = & X_1-X_3 ~ = ~
                \left\{\begin{array}{rl}
                1, & \mbox{if price $=$ 20} \\
                0, & \mbox{if price $=$ 50} \\
                -1, & \mbox{if price $=$ 90}
                \end{array}\right. 
        \\
        D_2 & = & X_2-X_3 ~ = ~
                \left\{\begin{array}{rl}
                0, & \mbox{if price $=$ 20} \\
                1, & \mbox{if price $=$ 50} \\
                -1, & \mbox{if price $=$ 90}
                \end{array}\right. 
\end{eqnarray*}
and fit the regression
\[
        y = \mu + \alpha_1 D_1 + \alpha_2 D_2 + \epsilon
\]
Let's try it:
\begin{verbatim}
> d1 _ p20-p90
> d2 _ p50-p90
>  sales.lm _ lm(lmsales ~ d1 + d2)
> summary(sales.lm)

Call: lm(formula = lmsales ~ d1 + d2)
Residuals:
    Min     1Q Median    3Q   Max 
 -194.8 -55.67 -4.224 78.53 113.8

Coefficients:
                 Value Std. Error    t value   Pr(>|t|) 
(Intercept)  6184.6074    15.0209   411.7333     0.0000
         d1   826.2588    21.2428    38.8960     0.0000
         d2  -106.5599    21.2428    -5.0163     0.0000

Residual standard error: 86.29 on 30 degrees of freedom
Multiple R-Squared: 0.9835 
F-statistic: 895.3 on 2 and 30 degrees of freedom, the p-value is 0 

Correlation of Coefficients:
   (Intercept)   d1 
d1  0.0            
d2  0.0        -0.5
\end{verbatim}
what about $\alpha_3$?  Since $\alpha_1 + \alpha_2 + \alpha_3 = 0$
(and since the mle of a function is the function of the mle's), we
know it can be calculated as
\[
        \hat\alpha_3 = - \hat\alpha_1 - \hat\alpha_2 
        = - 826.2588 + 106.5599 = -719.6998
\]
The SE of $\hat\alpha_3$ can also be found from the variance-covariance
matrix of $\hat\alpha_1$ and $\hat\alpha_2$ using a bit of matrix algebra, but
we will leave that for another day.





\section{A more realistic one-way ANOVA example}

\subsection{Initial Analysis and ANOVA table: cell means model}

The following data is from Box, Hunter and Hunter (1978) and is also
analyzed in Chapter 13 of the SPLUS Guide to Statistics.  It gives
blood coagulation times for each of four diets.

\begin{center}
\begin{tabular}{r|rrrr}
Diet  & A & B & C & D \\
\hline
\hline
Coag  & 62 & 63 & 68 & 56 \\
Times & 60 & 67 & 66 & 62 \\
      & 63 & 71 & 71 & 60 \\
      & 59 & 64 & 67 & 61 \\
      &    & 65 & 68 & 63 \\
      &    & 66 & 68 & 64 \\
      &    &    &    & 63 \\
      &    &    &    & 59 \\
\hline
\\
\end{tabular}

Blood Coagulation Times
\end{center}

There are many ways to enter this data.  In our case I will use NA
(Splus's code for ``missing data'') as a placeholder in order to start
with a square array, I'll attach the diet labels, and then delete the
NA's since they are not part of this data set (the data isn't really
missing, it was never intended to be collected).
\begin{verbatim}
> coag.times _ scan()
1:  62 63 68 56
5:  60 67 66 62
9:  63 71 71 60
13:  59 64 67 61
17:  NA 65 68 63
21:  NA 66 68 64
25:  NA NA NA 63
29:  NA NA NA 59
33: 
> diet _ factor(rep(c("A","B","C","D"),8))
> split(coag.times,diet)  # check that the factor labels are right
$A: [1] 62 60 63 59 NA NA NA NA
$B: [1] 63 67 71 64 65 66 NA NA
$C: [1] 68 66 71 67 68 68 NA NA
$D: [1] 56 62 60 61 63 64 63 59
> coag _ data.frame(coag=coag.times,diet=diet) # make the data frame
> coag _ coag[!is.na(coag$coag),] # delete rows with NA's
> split(coag$coag,coag$diet) # recheck that factor labels are still OK
$A: [1] 62 60 63 59
$B: [1] 63 67 71 64 65 66
$C: [1] 68 66 71 67 68 68
$D: [1] 56 62 60 61 63 64 63 59
> sapply(.Last.value,mean)
  A  B  C  D 
 61 66 68 61
> plot.factor(coag,style.bxp="old")

\end{verbatim}

\begin{figure}
\hfil\psfig{figure=coagplotfactor.ps,width=5.5in,height=3.5in}
\caption{Boxplots of coagulation times by diet for the blood
coagulation data.}
\end{figure}

Now we fit the model, look at some diagnostic plots, and consider the
analysis of variance table.  More details on the SPLUS parts of the
problem can be found in Chapter 13 of the SPLUS Guide to Statistics.

\begin{verbatim}
> coag.aov _ aov(coag ~ diet, data=coag)
> par(mfrow=c(2,3))
> plot(coag.aov)
\end{verbatim}

\begin{figure}
\hfil\psfig{figure=coagdiag.ps,width=6in}
\caption{Standard diagnostic plots for one-way ANOVA for the blood
coagulation data.}
\end{figure}

\begin{verbatim}
> model.tables(coag.aov,type="means")
Refitting model to allow projection
Tables of means
Grand mean
    
 64

 diet 
     A  B  C  D 
    61 66 68 61
rep  4  6  6  8
\end{verbatim}

The cell means $Y_{i.}$ and the grand mean $Y_{..}$ are illustrated in
the figure below; the cell means $\bar Y_{i.}$ are the estimates of
the $\mu_i$'s in the {\em cell means model}
\[
        Y_{ij} = \mu_i + \epsilon_{ij}, 
        ~~ i = 1,\ldots,\#cells,
        ~~ J = 1,\ldots, n_i
\]
where $k$ is the number of cells, and there are $n_i$ observations in
the $i^{th}$ cell.  As usual in regression, the error terms 
($\epsilon_{ij}$ in this case) are distributed $N(0,\sigma^2)$ for
some unknown error variance $\sigma^2$. 

\begin{center}
\begin{tabular}{|p{5in}|}
\hline
\sf Note: The roles of $i$ are $j$ are reversed in this example: $i$
now stands for the group or cell, and $j$ stands for individual
observations within a group or cell.
\\ \hline
\end{tabular}
\end{center}

\comment{
attach(coag,pos=1)
# search()
par(mfrow=c(1,1))
length(coag)
plot(as.numeric(diet)+runif(24,-.025,+.025),coag,xaxt="n",
        xlab="diet",ylab="coag times",xlim=c(0.5,4.5))
axis(1, at=1:4,labels=c("A","B","C","D"))
lines(c(0.5,4.5),c(64,64))
lines(c(0.5,1.5),c(61,61))
lines(c(1.5,2.5),c(66,66))
lines(c(2.5,3.5),c(68,68))
lines(c(3.5,4.5),c(61,61))
detach(pos=1)
}

\begin{figure}
{\unitlength=1in
\begin{picture}(6,4)
\put(0,0){\psfig{figure=cellmeans.ps,width=6in}}
\put(5,2.6){\mbox{$\bar Y_{..} = 64$}}
\put(1.5,1.9){\mbox{$\bar Y_{1.} = 61$}}
\put(2.75,2.95){\mbox{$\bar Y_{2.} = 66$}}
\put(3.85,3.35){\mbox{$\bar Y_{3.} = 68$}}
\put(4,1.9){\mbox{$\bar Y_{4.} = 61$}}
\end{picture}
}
\caption{Illustration of the fitted cell means model for the blood
coagulation data.}
\end{figure}

Now recall that we can write the deviation of a single observation
from the grand mean as
\[
        (Y_{ij} - Y_{..}) = (Y_{ij} - Y_{i.}) + (Y_{i.} - Y_{..})
\]
If we square and sum these terms, the magic of orthogonal sums of
squares tells us
\[
        \sum_i\sum_j (Y_{ij} - Y_{..})^2 = 
        \sum_i\sum_j (Y_{ij} - Y_{i.})^2 + 
        \sum_i n_i\cdot(Y_{i.} - Y_{..})^2
\]
with degrees of freedom $n-1$, $n-k$ and $k-1$ respectively.  This
gives rise to the ANOVA table
\begin{center}
\begin{tabular}{lr@{=}lrr@{=}l}
\em Source & \multicolumn{2}{c}{\em SS} & \em df & \multicolumn{2}{c}{\em MS} \\
\hline
\\
Factor & RegSS & $\sum_i n_i\cdot(Y_{i.} - Y_{..})^2$ & $k-1$ & RegMS & $RegSS/(k-1)$ \\
Error  & RSS & $\sum_i\sum_j (Y_{ij} - Y_{i.})^2$ & $n-k$ & RMS & $RSS/(n-k)$ \\
\\
\hline
\\
Total  & TSS & $\sum_i\sum_j (Y_{ij} - Y_{..})^2$ & $n-1$ & TMS & $TSS/(n-1)$
\end{tabular}
\end{center}
Here is SPLUS's ANOVA table for this ANOVA model.
\begin{verbatim}
> anova(coag.aov)     # some output omitted below
Analysis of Variance Table

Response: coag

Terms added sequentially (first to last)
          Df Sum of Sq Mean Sq  F Value         Pr(F) 
     diet  3       228    76.0 13.57143 0.00004658471
Residuals 20       112     5.6             
>  228/(228+ 112) # R^2
[1] 0.6705882
\end{verbatim}
The F statistic for testing whether the factor explains the variation
in $Y$ is
\[
        F^* = \frac{RegMS}{RMS} = \frac{(\mbox{between SS})/(k-1)}
        {(\mbox{within SS})/(n-k)}
\]
Under the null hypothesis
\[
        H_0: \mu_1 = \mu_2 = \cdots = \mu_k
\]
$F^*$ is distributed as $F(k-1,n-k)$.  Large values of $F^*$ argue in
favor of 
\[
        H_a: \mbox{not all $\mu_i$ are equal}
\]

We can also look at the effects (deviations from means), either by
direct calculation or with {\tt model.tables()}:
\begin{verbatim}
> model.tables(coag.aov,type="effects")
Refitting model to allow projection
Tables of effects

 diet 
     A B C  D 
    -3 2 4 -3
rep  4 6 6  8


> x1 _ ifelse(coag$diet=="A",1,0)
> x2 _ ifelse(coag$diet=="B",1,0)
> x3 _ ifelse(coag$diet=="C",1,0)
> x4 _ ifelse(coag$diet=="D",1,0)
> d1 _ x1-x4
> d2 _ x2-x4
> d3 _ x3-x4
> summary(lm(coag$coag ~ x1+x2+x3+x4 - 1))

Call: lm(formula = coag$coag ~ x1 + x2 + x3 + x4 - 1)
Residuals:
 Min    1Q     Median   3Q Max 
  -5 -1.25 -4.36e-016 1.25   5

Coefficients:
     Value Std. Error t value Pr(>|t|) 
x1 61.0000  1.1832    51.5544  0.0000 
x2 66.0000  0.9661    68.3165  0.0000 
x3 68.0000  0.9661    70.3867  0.0000 
x4 61.0000  0.8367    72.9089  0.0000 

Residual standard error: 2.366 on 20 degrees of freedom
Multiple R-Squared: 0.9989 
F-statistic: 4399 on 4 and 20 degrees of freedom, the p-value is 0 

Correlation of Coefficients:
   x1 x2 x3 
x2 0       
x3 0  0    
x4 0  0  0 
> summary(lm(coag$coag ~ d1+d2+d3))

Call: lm(formula = coag$coag ~ d1 + d2 + d3)
Residuals:
 Min    1Q     Median   3Q Max 
  -5 -1.25 1.249e-016 1.25   5

Coefficients:
                Value Std. Error   t value  Pr(>|t|) 
(Intercept)   64.0000    0.4979   128.5367    0.0000
         d1   -3.0000    0.9736    -3.0813    0.0059
         d2    2.0000    0.8453     2.3659    0.0282
         d3    4.0000    0.8453     4.7319    0.0001

Residual standard error: 2.366 on 20 degrees of freedom
Multiple R-Squared: 0.6706 
F-statistic: 13.57 on 3 and 20 degrees of freedom, the p-value is 0.00004658 

Correlation of Coefficients:
   (Intercept)      d1      d2 
d1  0.2106                    
d2 -0.0346     -0.4075        
d3 -0.0346     -0.4075 -0.3061

> -( - 3 + 2 + 4 )
[1] -3

\end{verbatim}

\end{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% FOLLOWING MAY BE USEFUL LATER??

\section{Uses of ANOVA models}

ANOVA models are used to study the effects of the independent
(predictor) variable(s) on the dependent (response variable), without
making assumptions about the ``shape'' of the relationship.  They are
used in randomized experiments, where the assumption of random
uncorrelated (and perhaps normal) errors in the regression model is
carefully built into the data by the experimenter, and in
observational studies, where the assumption of random uncorrelated
(perhaps normal) errors must be argued into the problem.

Many of the ideas in this section were taken from Chapter 14 of
\begin{description}
\item Neter, J., Wasserman, W., and Kutner, M. H. (1990).  {\em
Applied linear statistical models.} Homewood IL: Irwin.
\end{description}

\subsection{Experimental and Observational Studies}

{\bf Example 1.}

\begin{quote}

The effect of four cooking temperatures on the fluffiness of omelets
prepared from a mix was studied by randomly assigning five packages of
mix to each of four cooking temperatures.  This is an {\em experiment}
because the independent variable of interest (cooking temperature) was
controlled by the researcher.  The experimental design was a {\em
completely randomized design}, since the 20 packages of mix were
randomly distributed among the cooking temperatures; this (together
with some care in performing each of the 20 ``runs'') guarantees that
the errors in the ANOVA/regression model are random and independent;
it may {\em not} guarantee that they are normal, or have the same
variability in each of the four cells (cooking temperatures) of the
design.

In a real industrial setting (say, to decide what instructions to
write on boxes of omelet mix), probably all four cooking temperatures
produce somewhat fluffy omelets.  What is of interest is which one is
better, not what the overall average fluffiness is.  Since this is a
common situation in ANOVA applications, ANOVA software often doesn't
bother to report the grand mean, even in the ``grand mean plus
effects'' parameterization. 

If the packages were also chosen randomly from the company's output,
it may be reasonable to assume that these results will be
generalizable to the hypothetical population of all boxes of omelet
mix of this type produced by this company.  ~$\Box$

\end{quote}

{\bf Example 2.}

\begin{quote}

A study of the effects of education and type of past employment
experience of salespeople on their sales volumes was made by selecting
a random sample of salespeople currently employed at the company and
obtaining information on highest degree obtained, type of experience,
and sales volume for each selected employee.  This was an {\em
observational study} because the independent variables (education,
type of experience) was not controlled by the experimenter.  Random
sampling does guarantee that the errors will be random and independent
(but again normality and equal variances across cells is harder).

Because of the random sampling the results of this survey may
generalize to a hypothetical population of employees like those
working in sales for this company right now.  However, it seems less
likely that this is a population of interest---and even if it is, that
it will remain stable over time so that the results can be used with
future employees---than in {\bf Example 1}.

\end{quote}

{\bf Example 3.}

\begin{quote}

An appliance manufacturer operates three training centers in the
United States for training mechanics to service the company's
products.  At each center, two different training programs were
studied, with trainees randomly assigned to each program.  The outcome
is a score on a standardized repair exam.
~$\Box$

\end{quote}

\subsection{Factors}

The word {\em factor} is just another word for independent variable.
Usually factors are qualitative variables; the different values that a
factor can have (such as sales prices of \$20, \$50 or \$90) are
called {\em levels} of the factor.
\begin{itemize}

\item In {\bf Example 1} there is just one factor: the cooking
temperature.  It has four, ordered levels.

\item In {\bf Example 2} there are two factors: Education, and type of
experience.  Education could be a {\em quantitative factor} whose
values are just the number of years of education the salesperson has
had.  More likely it is an {\em ordered qualitative factor} with
levels like {\em High School}, {\em College Degree}, {\em Master's
Degree}, and {\em PhD}.  Type of past employment experience is
probably an {\em unordered qualitative factor}, with levels like {\em
laborer}, {\em office work}, {\em sales}, {\em consulting}, etc.

\item In {\bf Example 3} again there are two qualitative factors (both
unordered): which center the employee is at, and which training
program the employee was assigned to.

\end{itemize}

ANOVA applied to single-factor studies is called {\em one-way ANOVA}.
ANOVA applied to multi-factor studies is called {\em multi-way
ANOVA}. 

Many studies have a mix of ``experimental'' and ``observational''
components.  In a study, the factors that were controlled and randomly
assigned by the researcher are called {\em experimental factors}.
Those that merely classify subjects, but are not under the control of
the researcher, are called {\em classification factors}.  {\bf Example
3} is an example with both an {\bf experimental factor} (which of the
two programs the examinee was assigned to) and a {\bf classification
factor} (which training center the employee was at).

If we tried to analyze {\bf Example 3} as a one-way ANOVA using only
the training program, we might mask some important differences among
the centers (suppose at Center A, method 2 is much better than method
1, but at centers B and C, method 1 is moderately better than method
2; if the data is ``pooled across centers'' then we might conclude
that the two methods are equally effective.

If we analyze {\bf Example 3} as a two-way ANOVA, we still must be
careful: if method 1 is better than method 2 at all three sites, the
evidence in favor of method 1 is strong.  On the other hand, if method
2 does better at some centers and worse at others, these differences
are hard to interpret without getting more information: perhaps one
center has better teachers for method 2, or some centers are in
geographic locations where better education is provided.  At least the
differences among levels of the classification factor (training site)
suggest what sort of information is needed to understand the results.

\subsection{Randomization vs. Blocking}

As illustrated in the discussion of the site location in {\bf Example
3} above, a potential difficulty with factors not under the
researcher's control is that such factors may be hiding a confounding
variable (trainer skills, community education levels, etc.) that make
interpretation of the results difficult.  Still, it is better to
include such factors in one's analysis, for otherwise important
effects might cancel out.  This leads to the following maxim: 

\begin{quote}\em
Randomize what you can, block the rest.
\end{quote}

In other words, in an experiement, some factors will be under your
control and others will not.  {\em Randomize what you can:} i.e.,
randomly assign subjects to levels of all factors that you control;
{\em block the rest:} if some factors are not in your control, they
should still be included in the model.

This maxim makes less sense in observational studies: nothing is under
your control.  In such cases the best the researcher can do is to
include all the blocking variables that (s)he thinks will be
important.  However, standard errors tend to grow when more regressors
are included in the model, so we don't want to include so many that
the ANOVA analysis loses power to tell us anything interesting.  This
tradeoff is one of the most difficult aspects of careful analysis of
observational data (and there is no neat, clean solution!).

\subsection{Treatments or ``Cells''}

In a one-factor study, the treatments are the same as the levels of
the factor.  In a multifactor study, each unique combination of
factors is a treatment.  For example, a two-factor study of three
kinds of deodorant and four levels of exercise has twelve treatments.

When we are thinking about the numerical data more than about the
experiment that generated it we often use the synonym ``cell'' for
``treatment''.  They mean essentially the same thing though.

Choosing which treatments (which factors, how many levels, etc.) to
study is not easy.  This is part of experimental design, which we will
discuss more in a couple of weeks.  For now, here are some main points
to keep in mind:
\begin{itemize}
\item There are generally two goals of a study or series of studies:
one is to understand a process, and the other is to make a prediction
or comparison.

Studies that are intended to help in understanding a process may start
out small: the initial investigation may involve only a few factors,
or a few levels of the factors; and then after the main factors of
interest have been identified a more complete study may be undertaken,
to provide a more detailed picture of the process.

Studies that are intended to help answer a specific practical question
tend to have fewer factors, but even here including some extra levels
to help provide some understanding of the findings is useful.  For
example if your company is considering changing to a new manufacturing
process that is expected to allow you to increase the number of units
per day that you produce, and you want to study defects in the
manufactured products, you may want to try include not only the
``old'' method under the current production schedule, the ``new''
method under the current production schedule, and the ``new'' method
under the accellerated production schedule.  This way you can
disentangle the effect of the change in process from the effect of
speeding up the production schedule, on the number of defects you see.

\item The definition of ``What is a treatment'' is not always obvious.
This is especially clear in educational studies.  Suppose two teaching
methods were being compared, and different teachers were using the
different methods.  Since the teacher has an effect on how well a
method is implemented, it's may not be clear whether the treatment is
just the teaching method, or a combination of teacher and method.

\item A {\em control treatment} consists of withholding all active
treatments but still applying the same experiemental procedures as
would be needed for the active treatments.  Controls are needed, for
example,
\begin{itemize}
\item when the general effectiveness of all the active treatments is
not known or not consistent; or
\item in studies where ``halo'', ``placebo'', or ``experimenter bias''
effects could make an ineffective active treatment look effective.
\end{itemize}
Some studies, such as comparisons of a new treatment for a disease
vs. a current, already effective treatment for the disease, do not
need a control since there is an effective treatment standard, the
current treatment, to which we may compare.  In such studies it may
also be unethical to include a control treatment (withholding
treatment from those in need).

\end{itemize} 

\subsection{Units of Study}

Finally, in addition to ambiguities in what constitutes a
``treatment'' and which treatments should be administered, there may
be substantial ambiguity about what constitues a ``unit of study'' in
an experiment or observational study.  

The ``units of study'' are the rows in the $Y$ vector and $X$ matrix
in the matrix form of the regression equations.

{\bf Example.} 

\begin{quote}

In a study of a new method of teaching mathematics, school district
officials were able to get 60 students in each of nine schools to
participate.  In each school, 30 children were randomly assigned to
the ``current'' method, and 30 were assigned to the ``new'' method.
Three teachers taught three classes each of the ``current'' method
classes , and three more taught three classes each of the ``new''
method classes. ~$\Box$

\end{quote}

In the example, $60\times9=540$ children participated in the study,
with 270 children experiencing each teaching method.  We might think
of these children as the experimental units.   However,
\begin{itemize}
\item Because the children interact with one another in the classoom
(in both methods), their responses to the end-of-semester test will be
dependent on one another within each classroom (even if no one cheats
on the exam!).  Thus perhaps we should think of the classrooms as the
experimental units; now there are only 9 units in the ``current''
method and 9 units in the ``new'' method. 
\item Finally, perhaps the teacher also has an influence on the
outcome of the students; then all the children who have the same
teacher will have dependent responses, so perhaps classes within
teacher should be the units of study; now there are only 3 units in
the ``current'' method, and 3 units in the ``new'' method.
\end{itemize}
This example is frustrating because each of the possible units:
children, classrooms, and teachers, contributes something to the
outcome.  Mixed effects linear models, which we will discuss later,
provide one way to ``hedge our bets'' across the three possible
definitions of units in this problem.

Nevertheless, identifying and clearly stating what the units of study
are is an important point in both communication and in understanding
what the data are telling us.

\end{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



