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

\input info.tex

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

%\slides

\begin{document}

\begin{center}\bf
Notes for Tuesday February 3, 1998
\end{center}

These notes borrow from Neter, Wasserman and Kutner, Chapters 18 and
20.

\section{Multifactor Studies}

The data for an ANOVA model involving a single factor is sometimes
called a {\em ``one-way layout''}, since a table of the data would
usually be laid out one way along the factor as a row or column.

When there are two factors, the data is usually laid out as a two
dimensional table with the levels of the first factor in rows (say)
and levels of the second factor in columns; this is called a {\em
``two-way layout''}, and the corresponding model is called a {\em
two-way ANOVA model}.

In general there may be many factors that a researcher can control
(experimental factors) or must take account of (classification
factors) in a study.

\subsection{Some examples}

\begin{quote}
{\bf Example 1.}

A steel company studied the effects of carbon content and tempering on
the strength of steel. Carbon content was investigated at a {\bf high}
level and a {\bf low} level.  Tempering was also studied at a {\bf
  high} temperature and a {\bf low} temperature (the precise
definitions of carbon and temperature levels are not important here).
Thus there were $2\times2=4$ possible treatments:
\begin{center}
\begin{tabular}{cc}
\hline
\em Treatment & \em Description \\
\hline
1 & High carbon, high temperature \\
2 & High carbon, low  temperature \\
3 & Low  carbon, high temperature \\
4 & Low  carbon, low  temperature \\
\hline
\end{tabular}
\end{center}
The four treatments were then randomly assigned to twelve production
batches and the strength of steel in each batch was measured.
$\Box$
\end{quote}

This is an experiment, indeed it is a {\em completely randomized design.}
If the random assignment was done so that exactly three batches were
assigned to each treatment, then it is also a {\em balanced design}
(when the sample sizes in the treatments are unequal, the design is
called {\em unbalanced}).  

Laid out in a two-way table, the data might look like this:
\[
\begin{array}{|c|r|r|}
\hline
 & \multicolumn{2}{c|}{\em Temperature} \\
\cline{2-3}
\em Carbon & \rm High & \rm Low \\
\hline
\rm High &  60.35 & -16.82 \\
         & -45.63 & -24.02 \\
         &  76.88 & -32.09 \\
\hline
\rm Low  & -17.00 &  11.58 \\
         &  33.48 & -96.78 \\
         &  51.44 &  63.11 \\
\hline
\end{array}
\]

If you think of the design as four individual treatments (as they were
listed, 1, 2, 3, 4, above) you could analyze the data using one-way
ANOVA.  However, when there are two (or more) factors we usually want
to analyze the data in special ways that allow us to assess the
effects of the two factors, individually and jointly, on the outcome.
The individual effects of the factors are called {\em main effects}
and their joint effects are called {\em interactions}.  The separate
study of main effects and interactions is why the two-way layout (the
second table) and its associated analysis, two-way ANOVA, is
important.

\begin{quote}
{\bf Example 2.} 

An analyst studied the effects of family income (under \$10,000,
\$10,000--\$29,999, \$30,000--\$49,999, \$50,000 and more) and stage
in the life cycle of the family (1, 2 or 3) on appliance purchases.
Thus $4\times3=12$ treatments are defined:
\begin{center}
\begin{tabular}{cc}
\hline
\em Treatment & \em Description \\
\hline
1 & Under \$10,000, Family stage 1 \\
2 & Under \$10,000, Family stage 2 \\
3 & Under \$10,000, Family stage 3 \\
4 & \$10,000--\$29,999, Family stage 1 \\
\vdots & \vdots \\
12 & \$50,000 or more, Family stage 3 \\
\hline
\end{tabular}
\end{center}
The analyst selected 20 families whose income and family stage matched
each ``treatment'' above, yielding 240 famlies for the study.
$\Box$
\end{quote}

This is an observational study, because the data were obtained without
the analyst assigning ``treatments'' to each family.  It is a balanced
design (equal sample sizes for each treatment), and there are two
factors, so the data could be displayed in a two-way layout as in the
previous example.

\begin{quote}
{\bf Example 3.}

A psychologist interested in human-computer interaction wants to know
what makes audio menus (like the recorded messages ``Press 1 for
blah-blah, press 2 for blah-blah'' in bank-by-phone and voice mail
systems) easy or difficult for people to use.  Before beginning the
study, he measures subjects' reading span (a measure of the capacity
of their short-term memory) and categorizes them as having {\bf high},
{\bf medium}, or {\bf low} reading spans.  He then gives each subject
a series of tasks to perform with a bank-by-phone system, so that each
subject sees each combination of the following task features:
\begin{itemize}
\item A {\bf broad menu} system that requires two selections from among
  nine options each to reach the goal in the menu; versus a {\bf deep
    menu} system that requires four selections from among three
  options each to reach the goal;
\item A task that is expressible in terms of {\bf numerical} information,
  like ``deposit \$1000'', versus one that is expressible in {\bf
    verbal} terms, such as ``change mailing address''.
\item A task that is {\bf first} in a sequence of two tasks, and a
  task that is {\bf second} in a sequence of two tasks.
\end{itemize}
The number of ``treatments'' in this study is
$3\times2\times2\times2=24$.  Within each reading span group, the
psychologist selects four people for each of the eight combinations of
task features, for a total of 96 subjects.  He then assigns 10 tasks
to each person and records the proportion of tasks correctly performed.  
$\Box$
\end{quote}

This design is a four-factor, or four-way, design.  Sometimes it would
be referred to as a $3\times 2^3$ design.  It is also balanced---the
same number of subjects are assigned to each ``treatment.''

In this study there are three {\em experimental factors}---the types
of menu and task, that the experimenter explicity assigns to
subjects---and one {\em classification} factor---the reading span of
the subjects, which is not determined by the experimenter.  

\begin{quote}

{\em {\bf Aside}, about the definition of units and the definition of the
outcome measure:} From the psychologist's point of view, this may
not be the best design for understanding the influences of the various
factors on {\em tasks}, since people, not tasks, are the fundamental
units being studied.
\begin{itemize}
\item One alternative design would treat the tasks as the units, so
  that (perhaps) the outcome measure would be the proportion of
  persons getting each task correct.
  
\item Another alternative design would recognize that the psychologist
  might want to learn how each human subject responds to tasks with
  differing menu depths, numerical or verbal descriptors, etc.  This
  leads to a more complicated analysis, since now the {\em same} human
  subjects would get used in {\em different} treatment cells of the
  study, which introduces {\em dependence} between the ANOVA error
  terms $\epsilon_{ij}$ in different treatment cells.

\end{itemize}

\end{quote}

\subsection{Complete and Fractional Factorial Designs}

The three examples above are {\em complete factorial designs} : this
means, all possible combinations of the factors were studied.

In some studies it may be infeasable or undesirable to study all
possible combinations of factors.  In {\bf Example 1}, if the steel
company wished to study six temperatures, five carbon levels and four
methods of cooling the steel, the complete factorial design would have
$6\times5\times4=120$ different treatments.  It may be too expensive
or too time consuming to study three batches of steel under each of
these treatments.  

A possible alternative is a {\em fractional factorial design}, in
which only a fraction of the 120 possible treatments are studied.  If
carefully designed, a fractional factorial study can still give us
good information about the main effects of the three factors, as well
as some of their interactions.

\subsection{Advantages of Multifactor (Factorial) Studies}

An older alternative to factorial and fractional factorial designs is
to study two or more factors in separate one-way experiments, each of
which manipulates only one factor at a time.

\begin{itemize}
  
\item {\bf Information about interactions.}  At best the method of
  studying one factor at a time can give you information about the
  main effects of the factors, but it cannot give you information
  about the interactions.
  
  In {\bf Example 1}, suppose that the influence of carbon content on
  the strength of the steel is different, depending on whether the
  steel is tempered at a high or low temperature; this is the essence
  of an {\em interaction}.  This cannot be discovered by a pair of
  studies, the first of which holds the carbon level fixed while it
  varies the temperature, and the second of which holds the
  temperature fixed while varying the carbon level. 

  Moreover, if there is an interaction between carbon and temperature
  levels, then the fixed value of one factor that one chooses strongly
  influences the outcome as one varies the other factor, so that one's
  estimates of the main effects are biased.  Factorial studies can
  separate out interactive effects, since they look at all possible
  combinations of the factors.

\item {\bf Efficiency.}  Even if there is no interaction, factorial
  studies use the data more efficiently.  Suppose there is no
  interaction in the steel example, {\bf Example 1}, and suppose we
  can only afford to run 12 batches of steel.  In the factorial
  experiment the comparison of temperature levels can be based on 6
  batches at each level, {\em and} the comparison of carbon levels can
  be based on 6 batches at each level (the batches are effectively
  re-used).   However the sequence of two one-factor studies would be
  forced to base its temperature comparisons on three batches per
  level, and its carbon comparisons on three batches per level: once
  the steel is made, the batches cannot be reused.

\item {\bf Validity of findings.}  In the one-factor-at-a-time study,
  an assumtion must be made---that there are no interactions---to
  interpret the results, and there is no direct way to check the
  validity of the assumption (and therefore the validity of the
  study's conclusions).  In a complete factorial study, a hypothesis
  test can be performed to check whether an interaction is present or
  not.

\end{itemize}

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

\section{Two-Way Layouts} 

\subsection{Example: Balanced two-way layout}

A consumer organization studied the effects of age and gender of the
seller (car owner) on the the size of a cash offer for a used car.  A
single, medium-priced, six year old car was offered for sale to 36 car
dealers selected at random from all dealers in the region.  The 36
dealers were randomly ``assigned'' in groups of six, to each
``owner.''  The offers, in 100's of dollars, were as follows:
\begin{verbatim}
   "Owners'"         Dealers'
Age     Gender   Offers (x $100)
------- ------- -----------------
YOUNG   MALE    21 23 19 22 22 23
YOUNG   FEMALE  21 22 20 21 19 25
MIDDLE  MALE    30 29 26 28 27 27
MIDDLE  FEMALE  26 29 27 28 27 29
ELDERLY MALE    25 22 23 21 22 21
ELDERLY FEMALE  23 19 20 21 20 20
\end{verbatim}
To move the data into Splus and format it in a way that is useful with
{\tt aov()} we do the following (there is more than one way to skin
this cat, but\ldots):
\begin{verbatim}
402 > raw.offers _ read.table("car-offers.txt")
402 > raw.offers
       V1     V2 V3 V4 V5 V6 V7 V8 
1   YOUNG   MALE 21 23 19 22 22 23
2   YOUNG FEMALE 21 22 20 21 19 25
3  MIDDLE   MALE 30 29 26 28 27 27
4  MIDDLE FEMALE 26 29 27 28 27 29
5 ELDERLY   MALE 25 22 23 21 22 21
6 ELDERLY FEMALE 23 19 20 21 20 20
402 > offer.vector _ c(as.matrix(raw.offers[,3:8]))
402 > age _ factor(raw.offers[,1])
402 > gender _ factor(raw.offers[,2])
402 > car.offers _ data.frame(age=rep(age,6),gender=rep(gender,6),
+ offer=offer.vector)
402 > rm(age,gender,offer.vector) 
402 > car.offers
       age gender offer 
 1   YOUNG   MALE    21
 2   YOUNG FEMALE    21
 3  MIDDLE   MALE    30
 4  MIDDLE FEMALE    26
 5 ELDERLY   MALE    25
 6 ELDERLY FEMALE    23
 7   YOUNG   MALE    23
 8   YOUNG FEMALE    22
 9  MIDDLE   MALE    29
10  MIDDLE FEMALE    29
11 ELDERLY   MALE    22
12 ELDERLY FEMALE    19
.
.
.
34  MIDDLE FEMALE    29
35 ELDERLY   MALE    21
36 ELDERLY FEMALE    20
\end{verbatim}
We can get a general look at the data with {\tt plot.design} and {\tt
  plot.factor}:
\begin{verbatim}
402 > par(mfrow=c(2,2))
402 > plot.design(car.offers,fun=mean)
402 > plot.design(car.offers,fun=median)    
402 > plot.factor(car.offers)
\end{verbatim}

\psfig{figure=car-offer-1.ps,width=6in,angle=270}

Some things are already clear from these plots: middle-aged sellers
were offered a better price than the young or the old, and although
the mean and median offers to men were better than for women, this
difference is likely to be swamped by variability in the data.  We can
also see that the mean and median offers in each treatment were
approximately equal, so there is no strong skewing in the data; this
impression is confirmed by the boxplots.

Another plot that is interesting to look at with two-way data is the
{\tt interaction.plot}.  

\begin{verbatim}
402 > attach(car.offers)
402 > interaction.plot(age,gender,offer)
402 > interaction.plot(gender,age,offer)
402 > detach()
\end{verbatim}

\psfig{figure=car-offer-2.ps,width=6in,height=3in}

\begin{itemize}

\item 
The interaction plot on the left shows male-female differences in the
cell means, at each level of the age factor.  If gender and age
operate independently of one another to influence the offer made on
the car, these differences should be the same at each level of age,
i.e. the ``curves'' in the plot should be parallel.

\item
The interaction plot on the right shows differences between offers
made to each age group, for males and for females.  Again, if gender
and age operate independently to influence the offer made, these lines
should be parallel (the age differences in offer should be the same
for male as for female sellers).

\end{itemize}

Clearly, the two interaction plots contain the same information.  But
sometimes one plot makes it easier to ``see'' whether there is an
interaction than the other.  In this case, the potential for an
interaction is hard to see in the plot on the left, but it is easy to
see in the plot on the right.  

\subsubsection{The cell-means model}

By itself the cell-means model is not very interesting, but it gives
us some useful notation to use later.  The cell-means model for
two-way data is
\[
	y_{ijk} = \mu_{ij} + \epsilon_{ijk}
\]
where  $i$ is the level of the first (row) factor, $j$ is the level of
the second (column) factor, and $k$ indexes observation in treatment
cell $ij$.  $\mu_{ij}$ is the mean of the observations in treatment
cell $ij$, and $\epsilon_{ijk}$ are independent $N(0,\sigma^2)$ errors
for each observation.  It is useful to think of the data as laid out
in rows and columns, like so
\[
\begin{array}{lcccc}
\hline
               & \multicolumn{3}{c}{Age}                & Row \\
\cline{2-4}
Gender         & Young~(j=1) & Middle~(j=2) & Old~(j=3) & Average \\
\hline
Male~(i=1)     & \mu_{11}    & \mu_{12}     & \mu_{13}  & \mu_{1.} \\
Female~(i=2)   & \mu_{21}    & \mu_{22}     & \mu_{23}  & \mu_{2.} \\
\hline
Column~Average & \mu_{.1}    & \mu_{.2}     & \mu_{.3}  & \mu_{..}  
\end{array}
\]

The extra parameters in the table are defined as one might expect:
\begin{eqnarray*}
	\mu_{i.} & = & \frac1J \sum_{j=1}^J \mu_{ij} \\
	\mu_{.j} & = & \frac1I \sum_{i=1}^I \mu_{ij} \\
	\mu_{..} & = & \frac{1}{IJ} \sum_{i=1}^I \sum_{j=1}^J \mu_{ij}
\end{eqnarray*}
where $I$ is the number of levels of the row factor and $J$ is the
number of levels of the column factor.  The cell means are estimated
in the natural way, as 
\[
    \hat\mu_{ij} = \bar Y_{ij} = \frac{1}{n_{ij}} \sum_{k=1}^{n_{ij}} y_{ijk}
\]
where the $n_{ij}$ are the number of observations in each cell (when
all the $n_{ij}$ are equal, it is a {\em balanced design}).  Note that
the interaction plots above are just plots of the rows or columns of
the $\bar Y_{ij}$'s in the two-way layout.  Marginal totals
$\mu_{i.}$, $\mu_{.j}$ and $\mu_{..}$ would be estimated by the
corresponding averages of $bar Y$'s, i.e. $\bar Y_{i.}$, $\bar Y_{.j}$ and
$\bar Y_{..}$. 

The fitted values are $\hat y_{ij} = \bar Y_{ij}$, and the residuals
are $r_{ijk} = y_{ijk} - \hat y_{ij}$.  We could partition the total
sum of squares, as for one-way ANOVA, into between- and within- sums
of squares, 
\begin{eqnarray*}
	\sum_{i=1}^I\sum_{j=1}^J \sum_{k=1}^{n_{ij}} (\bar Y_{..} - y_{ijk})^2
	& = & 
	\sum_{i=1}^I\sum_{j=1}^J n_{ij} (\bar Y_{..} - \hat y_{ij})^2
	+ 
	\sum_{i=1}^I\sum_{j=1}^J \sum_{k=1}^{n_{ij}} (\hat y_{ij} - y_{ijk})^2
	\mbox{, or}\\
	SST & = & SSTR + SSE
\end{eqnarray*}
and this gives rise to the usual ANOVA table, as before.
\begin{center}
\begin{tabular}{lllr@{=}l}
\em Source & \em SS & \em df & \multicolumn{2}{c}{\em MS} \\
\hline
\\
Factors & SSTR & $IJ-1$ & MSTR & $SSTR/(IJ-1)$ \\
Error   & SSE  & $n-IJ$ & MSE & $SSE/(n-IJ)$ \\
\\
\hline
\\
Total   & SST  & $n-1$ & MST & $SST/(n-1)$
\end{tabular}
\end{center}
which we can examine immediately in SPLUS
\begin{verbatim}
402 > car.aov _ aov( offer ~ gender:age,data=car.offers)
402 > anova(car.aov)
Analysis of Variance Table

Response: offer

Terms added sequentially (first to last)
           Df Sum of Sq  Mean Sq  F Value         Pr(F) 
gender:age  5  327.2222 65.44444 27.39535 2.451114e-010
 Residuals 30   71.6667  2.38889            
\end{verbatim}
Note the use of ``:'' to indicate the row and column variables (more
generally ``:'' is used to construct interactions, but more about that
later).  Unsurprisingly (given the plots we did before), something
about the factors matters, but this analysis does not reveal clearly
what it is.

\subsubsection{Main Effects and Interactions}

To learn more about which of the two factors are important and how the
factors may or may not interact, we now define and examine the {\em
main effects} and {\em interactions} of the two factors for this
problem.

\begin{itemize}
\item Main row (first factor) effects: 
\[
	\alpha_i = \mu_{i.} - \mu_{..};
\]

\item Main column (second factor) effects: 
\[
	\beta_j = \mu_{.j} - \mu_{..};
\]
\end{itemize}

It is easy to verify that $\sum_i \alpha_i = 0$ and $\sum_j \beta_j =
0$, and it is easy to see how to estimate $\alpha_i$ and $\beta_i$ in
terms of the $\bar Y$'s.

To see whether the first factor ``matters'' we want to construct a
test of $H_0: \alpha_i \equiv 0$ vs $H_1:~not~all~0$; and similarly to
see whether the second factor ``matters'' we should test $H_0: \beta_j
\equiv 0$ vs $H_1:~not~all~0$. 

In one-way ANOVA we say that the cell means model and the grand means
model were related by the equation $\mu_i = \mu_{.} + \alpha_i$.  The
corresponding equation here,
\begin{eqnarray*}
	\mu_{ij} & \stackrel{?}{=} & \mu_{..} + \alpha_i + \beta_j \\
	& = & \mu_{i.} + \mu_{.j} - \mu_{..}
\end{eqnarray*}
is not true in general.  The differences between the left hand sides of
this ``equation'' and the right hand sides are called the {\em interactions},
\begin{itemize}
\item Interactions: 
\begin{eqnarray*}
	(\alpha\beta)_{ij} & = & \mu_{ij} - (\alpha_i + \beta_i + \mu_{..}) \\
	& = & \mu_{ij} - \mu_{i.} - \mu_{.j} + \mu_{..}
\end{eqnarray*}
\end{itemize}
It is also easy to see that $\sum_i\sum_j (\alpha\beta)_{ij} = 0$.
Now we have equality:
\begin{eqnarray*}
	\mu_{ij} & = & \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} \\
	& = & \mu + (\mu_{i.} - \mu) + (\mu_{.j} - \mu) + 
		(\mu_{ij} - \mu_{i.} - \mu_{.j} + \mu) ,
\end{eqnarray*}
where $\mu=\mu_{..}$ is the grand mean; and this gives rise to the
``grand mean plus effects'' model for two-way ANOVA, namely
\[
    y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \epsilon_{ijk}
\]
subject to the conditions that the $\alpha_i$'s sum to zero, the
$\beta_j$'s sum to zero, and the $(\alpha\beta)_{ij}$'s sum to zero.

In any two-way ANOVA, it turns out that the group of main effects for
the first factor is orthogonal to the group of main effects for the
second factor, and both of these are orthogonal to the interactions.
This means we can partition SSR in the ANOVA table above into
contributions from each set of main effects and the interactions:
\begin{eqnarray*}
	\sum_{i=1}^I\sum_{j=1}^J n_{ij} (\bar Y_{..} - \hat y_{ij})^2 
	& = & \\
	\sum_{i=1}^I\sum_{j=1}^J n_{ij} (\bar Y_{..} - \bar Y_{ij})^2 
	& = & 
	\sum_{i=1}^I n_{i.} (\bar Y_{..} - \bar Y_{i.})^2
	+
	\sum_{j=1}^J n_{.j} (\bar Y_{..} - \bar Y_{.j})^2
	+
	n\sum_{i=1}^I\sum_{j=1}^J (\bar Y_{..} - \bar Y_{i.}
		- \bar Y_{.j} + \bar Y_{ij})^2 \\
	\mbox{or}~~~~~~~~~~~~~~\\
	SSR & = & SS_{\alpha} + SS_{\beta} + SS_{(\alpha\beta)}
\end{eqnarray*}
and this allows us to expand the ANOVA table.  

In the case of a balanced design, say $n_{ij} \equiv m$ for all cells,
the expanded ANOVA table looks like this:
\begin{center}
\begin{tabular}{lllr@{=}l}
\em Source & \em SS & \em df & \multicolumn{2}{c}{\em MS} \\
\hline
\\
Factors & SSTR & $IJ-1$ & MSTR & $SSTR/(IJ-1)$ \\
\\
Main effects ($\alpha$) & SS$_{\alpha}$        & $I-1$ &
			  MS$_{\alpha}$ & SS$_{\alpha}/(I-1)$ \\
Main effects ($\beta$)  & SS$_{\beta}$         & $J-1$ &
			  MS$_{\beta}$ & SS$_{\beta}/(J-1)$ \\
Interactions            & SS$_{(\alpha\beta)}$ & $(I-1)(J-1)$ &
			  MS$_{(\alpha\beta)}$ & SS$_{(\alpha\beta)}/((I-1)(J-1))$ \\
\\
\hline
\\
Error   & SSE  & $IJ(m-1)$ & MSE & $SSE/(IJ(m-1))$ \\
\\
\hline
\\
Total   & SST  & $mIJ-1$ & MST & $SST/(mIJ-1)$
\end{tabular}
\end{center}
This is also easy to obtain in SPLUS:
\begin{verbatim}
402 > car.aov _ aov( offer ~ gender*age,data=car.offers)
402 > anova(car.aov)
Analysis of Variance Table

Response: offer

Terms added sequentially (first to last)
           Df Sum of Sq  Mean Sq  F Value    Pr(F) 
    gender  1    5.4444   5.4444  2.27907 0.141592
       age  2  316.7222 158.3611 66.29070 0.000000
gender:age  2    5.0556   2.5278  1.05814 0.359700
 Residuals 30   71.6667   2.3889                  
\end{verbatim}
Note the use of the ``*'' to indicate that main effects {\tt gender},
{\tt age} as well as the interactions {\tt gender:age} should be
included in the model.  Scanning the $F$-tests here, it seems clear
that age matters and gender does not (since the gender and age effects
are orthogonal it does not matter in which order we examine them, as
can be seen from the alternate order for the factors below).
\begin{verbatim}
402 > car.aov _ aov( offer ~ age*gender,data=car.offers)
402 > summary(car.aov)
           Df Sum of Sq  Mean Sq  F Value    Pr(F) 
       age  2  316.7222 158.3611 66.29070 0.000000
    gender  1    5.4444   5.4444  2.27907 0.141592
age:gender  2    5.0556   2.5278  1.05814 0.359700
 Residuals 30   71.6667   2.3889                  
\end{verbatim}

Moreover, the interaction doesn't seem to matter.  When we can safely
set the interaction terms to zero in a two-way ANOVA, we arrive at
what is called the {\em additive model},
\[
    y_{ijk} = \mu + \alpha_i + \beta_j + \epsilon_{ijk}
\]
In the additive model it is true that $\mu_{ij} = \mu + \alpha_i +
\beta_j$, so differences in cell means in one factor are the same, no
matter the value of the second factor:
\[
	\mu_{ij} - \mu_{i'j} = \mu_{ij'} - \mu_{i'j'}
\]
This is why the curves in an interaction plot should be parallel when
there is no interaction between the variables.

\subsubsection{Model Selection and Diagnostics}

Usually, it is only sensible to include or exclude entire factors or
interactions in an ANOVA model, rather than looking at individual
factor levels.  In the car offer problem, we probably would not include
the interaction {\tt gender:age}, and we would probably also exclude
the gender effects.

To illustrate diagnostics we will examine the diagnostic plots of all
three models under consideration.

\begin{verbatim}
402 > car.aov _ aov( offer ~ gender*age,data=car.offers)
402 > par(mfrow=c(2,3)); plot(car.aov)
402 > par(mfrow=c(1,2))
402 > hist(residuals(car.aov)); plot(residuals(car.aov))
\end{verbatim}

\psfig{figure=car-offer-3.ps,width=6in,height=3in}

\hfil\psfig{figure=car-offer-4.ps,width=6in,height=3in}

\begin{verbatim}
402 > car.aov _ aov( offer ~ gender + age,data=car.offers)
402 > anova(car.aov)
          Df Sum of Sq  Mean Sq  F Value     Pr(F) 
   gender  1    5.4444   5.4444  2.27082 0.1416397
      age  2  316.7222 158.3611 66.05069 0.0000000
Residuals 32   76.7222   2.3976             
402 > par(mfrow=c(2,3))
402 > plot(car.aov)
402 > par(mfrow=c(1,2))
402 > hist(residuals(car.aov))
402 > plot(residuals(car.aov))
\end{verbatim}

\psfig{figure=car-offer-5.ps,width=6in,height=3in}

\psfig{figure=car-offer-6.ps,width=6in,height=3in}

\clearpage

\begin{verbatim}
402 > car.aov _ aov( offer ~ age,data=car.offers)
402 > anova(car.aov)
          Df Sum of Sq  Mean Sq  F Value         Pr(F) 
      age  2  316.7222 158.3611 63.60142 4.768963e-012
Residuals 33   82.1667   2.4899       
402 > par(mfrow=c(2,3))
402 > plot(car.aov)
402 > par(mfrow=c(1,2))
402 > hist(residuals(car.aov))
402 > plot(residuals(car.aov))
\end{verbatim}

\psfig{figure=car-offer-7.ps,width=6in,height=3in}

\psfig{figure=car-offer-8.ps,width=6in,height=3in}

Other diagnostics, for example examining the relationship between
within-cell variances and estimated cell means to find a
variance-stabilizing transformation, proceed in much the same way as
one-way ANOVA.

\subsubsection{Fitted Values, Parameter Estimates}

Fitted values are easy to obtain in SPLUS from the {\tt aov} object,
using {\tt fitted()}, or by computing cell means by hand.  
\begin{itemize}
\item If interactions are {\em absent}, it is usually most interesting
to obtain estimates of the marginal means $\mu_{i.}$ and $\mu_{.j}$.
\item If interactions are {\em present}, then one usually must focus on the
cell means $\mu_{ij}$ themselves.
\end{itemize}

Parameter estimates are a little tricky, again because SPLUS's default
ANOVA parameterization gives estimates of the Helmert contrasts
instead of something sensible like the factor effects and
interactions.

If one wants estimates of the factor effects and interactions, one can
obtain them much as with one-way ANOVA.  Here are is how to do it with
the {\tt C( , )} function; for example {\tt C(age,sum)} says to
parameterize the model so that age is estimated with factor effects
that sum to zero.

\begin{verbatim}
402 > car.aov _ aov(offer ~ C(gender,sum) + C(age,sum), data=car.offers)
402 > summary.lm(car.aov)

Call: aov(formula = offer ~ C(gender, sum) + C(age, sum), data = car.offers)
Residuals:
    Min     1Q  Median     3Q   Max 
 -2.889 -1.028 -0.1111 0.9444 3.889

Coefficients:
                  Value Std. Error  t value Pr(>|t|) 
   (Intercept)  23.5556   0.2581    91.2765   0.0000
C(gender, sum)  -0.3889   0.2581    -1.5069   0.1416
  C(age, sum)1  -2.1389   0.3650    -5.8606   0.0000
  C(age, sum)2   4.1944   0.3650    11.4928   0.0000

Residual standard error: 1.548 on 32 degrees of freedom
Multiple R-Squared: 0.8077 
F-statistic: 44.79 on 3 and 32 degrees of freedom, the p-value is 1.466e-011 

Correlation of Coefficients:
               (Intercept) C(gender, sum) C(age, sum)1 
C(gender, sum)  0.0                                   
  C(age, sum)1  0.0         0.0                       
  C(age, sum)2  0.0         0.0           -0.5     
\end{verbatim}

\subsection{Example: One observation per cell}

The data below (used by Box, Hunter, and Hunter (1978)) were collected
to determine the effect of treatments A, B, C, and D on the yield of
penicillin in a penicillin manufacturing process.  The values of the
response variable {\tt yield} are the numbers in the table, and the
columns of the table correspond to the levels A, B, C, and D of the
treatment factor. 

There was a second factor, namely the blend factor, since a separate
blend of the corn-steep liquor had to be made for each application of
the treatments.  Your main interest is in determining whether the
treatment factor affects yield.  The blend factor is of only secondary
interest; it is a {\em blocking} variable introduced to increase the
sensitivity of the inference for treatments.  The order of the
treatments within blocks was chosen at random.  

This is a {\em randomized blocks} experiment: four penicillin cultures
within each blend were randomly assigned to the four treatments.

\begin{verbatim}
Effect of four treatments on penicillin yield

         Treatment
Block    A  B  C  D
Blend-1 89 88 97 94
Blend-2 84 77 92 79
Blend-3 81 87 87 85
Blend-4 87 92 89 84
Blend-5 79 81 80 88
\end{verbatim}

Here is one way to convert this to a form SPLUS will like
\begin{verbatim}
402 > pen.raw _ read.table("penicillin.txt",header=T)
402 > pen.raw _ c(as.matrix(pen.raw))
402 > pen.raw
 [1] 89 84 81 87 79 88 77 87 92 81 97 92 87 89 80 94 79 85 84 88
402 > pendata _ data.frame(blend=factor(rep(1:5,4)),
+ tx=rep(LETTERS[1:4],rep(5,4)), yield=pen.raw)
402 > pendata
   blend tx yield 
 1     1  A    89
 2     2  A    84
 3     3  A    81
 4     4  A    87
 5     5  A    79
 .
 .
 .
20     5  D    88
\end{verbatim}
We begin with the usual plots
\begin{verbatim}
402 > par(mfrow=c(2,2))
402 > plot.design(pendata)
402 > plot.design(pendata,fun=median)
402 > plot.factor(pendata)
\end{verbatim}

\psfig{figure=pendata-1.ps,width=6in,height=3in}


It turns out that you can't fit an ANOVA model with full,
unconstrained interactions to a data sets that has only one
observation per cell, if you want any degrees of freedom left over for
the error variance.

To see this, recall the ANOVA table for balanced designs ($n_{ij}
\equiv m$ for all cells),
\begin{center}
\begin{tabular}{lllr@{=}l}
\em Source & \em SS & \em df & \multicolumn{2}{c}{\em MS} \\
\hline
\\
Factors & SSTR & $IJ-1$ & MSTR & $SSTR/(IJ-1)$ \\
\\
Main effects ($\alpha$) & SS$_{\alpha}$        & $I-1$ &
			  MS$_{\alpha}$ & SS$_{\alpha}/(I-1)$ \\
Main effects ($\beta$)  & SS$_{\beta}$         & $J-1$ &
			  MS$_{\beta}$ & SS$_{\beta}/(J-1)$ \\
Interactions            & SS$_{(\alpha\beta)}$ & $(I-1)(J-1)$ &
			  MS$_{(\alpha\beta)}$ & SS$_{(\alpha\beta)}/((I-1)(J-1))$ \\
\\
\hline
\\
Error   & SSE  & $IJ(m-1)$ & MSE & $SSE/(IJ(m-1))$ \\
\\
\hline
\\
Total   & SST  & $mIJ-1$ & MST & $SST/(mIJ-1)$
\end{tabular}
\end{center}
we see that, when $m=1$, if we include the interaction terms, there
will be no degrees of freedom left for calculating the SSE and MSE.
We can see this directly in SPLUS:
\begin{verbatim}
402 > pen.aov _ aov(yield ~ blend*tx,data=pendata)
402 > summary(pen.aov)
         Df Sum of Sq  Mean Sq 
   blend  4       264 66.00000
      tx  3        70 23.33333
blend:tx 12       226 18.83333
\end{verbatim}
(and no error line is given in the ANOVA table).  However, if we omit
the interaction term, then we can use the interaction sum of squares
as the error sum of squares:
\begin{verbatim}
402 > pen.aov _ aov(yield ~ blend + tx, data=pendata)
402 > summary(pen.aov)
          Df Sum of Sq  Mean Sq  F Value     Pr(F) 
    blend  4       264 66.00000 3.504425 0.0407462
       tx  3        70 23.33333 1.238938 0.3386581
Residuals 12       226 18.83333                   
\end{verbatim}
Unfortunately, the treatment does not seem to be having much effect on
the penicillin yield. On the other hand, we seem to have accidentally
discovered that the blocking factor, the blend of the corn-steep
liquor used, is marginally significant!

The diagnostic plots do not yield any substantial difficulties with
the regression.  A time series plot of the residuals, not reproduced
here, also does not reveal any substantial problems.

\psfig{figure=pendata-3.ps,width=6in,height=3in}

One might also ask whether the interaction plots can be used to assess
whether there ``might be'' an interaction that the ANOVA simply can't
pick up.  We can maket the plots, but they are very difficult to
interpret.  For example the interaction plots for this data set are
quite jumbled and nonparallel, yet there is no strong indication of
interaction in them.

\psfig{figure=pendata-2.ps,width=6in,height=2in}

Chapter 13 of the SPLUS Guide to Statistics has more information about
testing and estimating specially-constrained interaction terms;
however for many purposes these interaction terms would be hard to
interpret.

Finally, it might be interesting to compare these results with the
results one would get if one ignored the blocking factor and did a
one-way ANOVA on the four treatments:
\begin{verbatim}
402 > pen.aov _ aov(yield ~ tx, data = pendata)
402 > summary(pen.aov)
          Df Sum of Sq  Mean Sq   F Value     Pr(F) 
       tx  3        70 23.33333 0.7619048 0.5317827
Residuals 16       490 30.62500                
\end{verbatim}
Note that the $p$-value is even higher for {\tt tx}, the treatment,
than it was under the two-way, randomized blocks analysis.  A more
careful inspection of the two ANOVA tables explains why: the error SS
in the one-way analysis is the sum of the error variance in the
two-way analysis plus the SS due to the blocks: in other words, the
variation due to differing blocks (differing corn-liquor blends) has
been lumped into the error variance of the one-way analysis.  This is
why blocking is such an effective idea in experimental design. 

\subsection{Example: Unbalanced two-way layout}

Human growth hormone was administered at a clinical research center to
growth-hormone-deficient, short children who had not yet reached
puberty.  The investigator was interested in the effects of gender
(factor A: male or female) and bone development (factor B: severely (1),
moderately (2), or mildly (3), depressed) on the rate of growth
induced by the hormone.  
\begin{verbatim}
M 1 1.4
M 1 2.4
M 1 2.2
M 2 2.1
M 2 1.7
M 3 0.7
M 3 1.1
F 1 2.4
F 2 2.5
F 2 1.8
F 2 2.0
F 3 0.5
F 3 0.9
F 3 1.3
\end{verbatim}
Reading the data in this format is straightforward, so we do this and
then produce what we hope are useful initial plots.
\begin{verbatim}
402 > hormone.data _ read.table("hormone.txt")
402 > names(hormone.data) 
[1] "V1" "V2" "V3"
402 > names(hormone.data)  _ c("gender","severity","growthrate")
402 > hormone.data$severity _ factor(hormone.data$severity)
402 > par(mfrow=c(2,2))
402 > plot.design(hormone.data)
402 > plot.design(hormone.data, fun=median)
402 > plot.factor(hormone.data)
402 > attach(hormone.data)
402 > par(mfrow=c(1,2))
402 > interaction.plot(gender,severity,growthrate)
402 > interaction.plot(severity,gender,growthrate)
402 > detach()
\end{verbatim}

\psfig{figure=hormone-1.ps,width=6in,height=3in}

\psfig{figure=hormone-2.ps,width=6in,height=2in}

It appears that severity will matter, but gender will probably not.
There does not appear to be much evidence of interactions in the
interaction plots.

Clearly the data are unbalanced; the cell sample sizes are
\begin{verbatim}
Cell  n_{ij}
M 1     3
M 2     2
M 3     2
F 1     1
F 2     3
F 3     3
\end{verbatim}
In principle there is little difference between doing an ANOVA with
unbalanced designs and doing ANOVA with balanced designs.  Calculation
of sums of squares, and presentation of the ANOVA table, is more
difficult, because the sums of squares in the formulas
\begin{eqnarray*}
		SST & \stackrel{?}{=} & SSTR + SSE \mbox{, and} \\
	 	SSTR & \stackrel{?}{=} & SS_{\alpha} + SS_{\beta}
 			+ SS_{(\alpha\beta)} 
\end{eqnarray*}
are no longer orthogonoal---these equalities no longer hold.  Instead,
sequential sums of squares as in regression must be calculated.  In
fact, in many computer packages, the simplest way to get correct
parameter estimates and sums of squares for unbalanced designs is to
analyze the data using dummy variables in a regression package.

These problems can be illustrated in SPLUS without resorting to
regression, because {\tt aov()} does properly handle unbalanced ANOVA
problems.  
\begin{verbatim}
402 > hormone.aov _ aov(growthrate ~ gender*severity,data=hormone.data)
402 > summary(hormone.aov)
                Df Sum of Sq  Mean Sq  F Value     Pr(F) 
         gender  1  0.002857 0.002857  0.01758 0.8977853
       severity  2  4.396000 2.198000 13.52615 0.0027133
gender:severity  2  0.075429 0.037714  0.23209 0.7980337
      Residuals  8  1.300000 0.162500              


402 > hormone.aov _ aov(growthrate ~ severity*gender,data=hormone.data)
402 > summary(hormone.aov)
                Df Sum of Sq  Mean Sq  F Value     Pr(F) 
       severity  2  4.306286 2.153143 13.25011 0.0028912
         gender  1  0.092571 0.092571  0.56967 0.4720219
severity:gender  2  0.075429 0.037714  0.23209 0.7980337
      Residuals  8  1.300000 0.162500               

402 > par(mfrow=c(2,3))
402 > plot(hormone.aov)
\end{verbatim}

Either way, it appears that only severity is an important influence on
the growth rate, as we suspected from the boxplots.

\psfig{figure=hormone-3.ps,width=6in,height=3in}

The diagnostic residual plots look reasonably good here too.  One
thing to watch out for when using these plots with unbalanced data is
that if you are judging the ``equal variance'' assumption by eye, be
careful to note that even if the variances in all cells {\em are} the
same, cells with fewer observations will appear to look less variable.

Defining the grand mean and the marginal means in an unbalanced
two-way layout can be done several different ways, just as for
unbalanced one-way problems.  The two main choices are
\begin{itemize}
\item {\em Unweighted means.}  The definitions of $\mu$ and of
$\mu_{i.}$ and $\mu_{.j}$ are exactly as for the balanced case.
\item {\em Sample size weighting.}  The definitions are changed to be
\begin{eqnarray*}
	\mu_{i.} & = & \frac{1}{n_{i.}}\sum_{j=1}^J n_{ij} \mu_{ij} \\
	\mu_{.j} & = & \frac{1}{n_{.j}}\sum_{i=1}^I n_{ij} \mu_{ij} \\
	\mu_{..} & = & \frac{1}{n_{..}}\sum_i\sum_j n_{ij} \mu_{ij}
\end{eqnarray*}
and corresponding changes in the ``sum to zero'' constraints on the
$\alpha$'s, $\beta$'s, etc.
\end{itemize}

In SPLUS, the {\tt dummy.coef()} function standardizes factor and
treatment effects so that they sum to zero, corresponding to
unweighted means.

\begin{verbatim}
402 > dummy.coef(hormone.aov)
$"(Intercept)":
 (Intercept) 
         1.7

$severity:
   1   2    3 
 0.5 0.3 -0.8

$gender:
   F    M 
 0.1 -0.1

$"severity:gender":
  1F             2F   3F   1M            2M  3M 
 0.1 -2.081668e-017 -0.1 -0.1 2.081668e-017 0.1
\end{verbatim}

The {\tt model.tables()} function produces factor effect estimates
that correspond to sample-size weighting.

\begin{verbatim}
402 > model.tables(hormone.aov,se=T)
Refitting model to allow projection
Standard error information not returned as design is unbalanced. 
Standard errors can be obtained through se.contrast.
Tables of effects

 severity 
         1      2       3 
    0.4571 0.3771 -0.7429
rep 4.0000 5.0000  5.0000

 gender 
          F        M 
    0.07714 -0.07714
rep 7.00000  7.00000

 severity:gender 
Dim 1 : severity
Dim 2 : gender 
         F      M 
  1  0.171 -0.057
rep  1.000  3.000
  2  0.011 -0.017
rep  3.000  2.000
  3 -0.069  0.103
rep  3.000  2.000
\end{verbatim}

Note that {\tt model.tables()} would not produce standard errors.  In
the unbalanced case, the standard errors for comparing two means
depend on the replication factors, hence it could be very complex to
tabulate all combinations. Instead, they can be computed directly with
{\em se.contrast}.  For instance, to compare the first and third bone
severity levels for boys in our hormone study:

\

\begin{verbatim}
402 > se.contrast(hormone.aov,contrast=list(gender=="M",severity==1,
+ gender=="M",severity==3),data=hormone.data)
Refitting model to allow projection
[1] 0.1703463
\end{verbatim}
By default, the standard error of the difference of the means specified by
contrast is computed.  The contrast itself can be computed using {\tt
model.tables()} again:
\begin{verbatim}
402 > model.tables(hormone.aov,type="means")
Refitting model to allow projection
Tables of means
Grand mean
       
 1.643

 severity 
      1    2   3 
    2.1 2.02 0.9
rep 4.0 5.00 5.0

 gender 
       F     M 
    1.72 1.566
rep 7.00 7.000

 severity:gender 
Dim 1 : severity
Dim 2 : gender 
      F   M 
  1 2.4 2.0
rep 1.0 3.0
  2 2.1 1.9
rep 3.0 2.0
  3 0.9 0.9
rep 3.0 2.0
\end{verbatim}
Thus, an approximate 95\% confidence interval for the difference in
growthrate for the first and third bone severity levels among boys is
$(2.0-0.9) \pm  qt(0.975,3)*(0.1703463)$, which runs from 
0.557882 to 1.642118 ({\tt qt(0.975,3) = 3.182446}).

Obtaining SE's for other contrasts, using {\tt se.contrast()}, is
explained in Chapter 14 of the SPLUS Guide to Statistics.

\end{document}

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


