\section{S Functions for Categorical Experiments}
\sectionmark{S Functions}
This sections describes the basic use of S functions to do analysis
of variance, to display the results of experiments graphically, and
to generate designs.

The modeling techniques of the analysis of variance
apply the ideas of formulas and data from previous chapters,
expanding on the range of formulas
and introducing some summaries
of particular relevance for categorical experiments.
Subsection \ref{aov.graphics} presents some summary and diagnostic techniques,
chiefly graphical.
These are particularly valuable and we encourage you to use the plots
as preliminary views of the data and as communication tools.
Subsection \ref{design.functions} presents functions
for generating designs and running experiments.
This subsection will not be relevant if you are only analyzing data
after its collection.
The first two subsections, however, should be read by anyone interested
in analysing categorical experiments.

\subsection{Analysis of Variance Models\label{aov.models}}

The basic expression for fitting an analysis of variance is
\begin{Example}
aov(formula, data)
\end{Example}
where \Co{formula} is a model formula, relating the response
to appropriate factors,
and \Co{data} is an optional data frame containing the data from
the experiment.
Suppose the catalyst experiment in Table \ref{bhh.table} is contained
in data frame \Co{catalyst}.
To fit a simple additive model in the three factors,
\begin{Example}\label{cat.example1}
aovcat <- aov(Yield \Twiddle Temp + Conc + Cat, catalyst)
\end{Example}
The object returned is of class \Co{"aov"}.
It is very similar to a linear model object, but has some special methods
defined for it.
In particular, it prints quite differently:
\begin{Example}
> aovcat
Call:
aov(formula = Yield \Twiddle Temp + Conc + Cat, data = catalyst)

Terms:
                  Temp   Conc    Cat Residuals 
 Sum of Squares 1058.0   50.0    4.5     205.0
Deg. of Freedom      1      1      1         4

Degrees of freedom: 7 total; 4 residual
Residual standard error: 51.25 
\end{Example}
Where a printed linear model shows individual coefficients,
printed analysis of variance objects show terms, which
may correspond to several coefficients.
In this, and in other ways we will demonstrate, our focus will shift
somewhat in this chapter.
Nothing has been lost, however, because methods for linear models can
always be applied explicitly to \Co{aov} objects.
To print \Co{aovcat} like a linear model:
\begin{Example}
> print.lm(aovcat)
Call:
aov(formula = Yield \Twiddle Temp + Conc + Cat)

Coefficients:
 (Int.)  Temp Conc   Cat 
  64.25 -11.5  2.5 -0.75

Degrees of freedom: 8 total; 4 residual
Residual standard error: 51.25 
\end{Example}
A summary table for the analysis of variance can be produced by the \Co{summary()}
function:
\begin{Example}
> summary(aovcat)
          Df Sum of Sq Mean Sq     F 
     Temp  1    1058.0  1058.0 20.64
     Conc  1      50.0    50.0  0.98
      Cat  1       4.5     4.5  0.09
Residuals  4     205.0    51.2  1.00
\end{Example}
We could also have attached the data frame:
\begin{Example}
attach(catalyst)
\end{Example}
to allow us to refer to the individual variables (\Co{Temp},  etc.) in arbitrary
S expressions.
In this case, the second argument to \Co{aov()} is omitted.

\subsection*{Crossed and Nested Terms}

In models from categorical experiments, interactions
between the factors are often of importance.
When possible, good data analysis suggests that one should ask whether the effect
of one factor depends on the levels of one or more other terms.
This leads us to a richer use of formulas than was typical of Chapter \LinearChapter.

The expression \Co{A\ *\ B} in a model formula puts into the model the terms for
\Co{A}, \Co{B}, and the interaction of \Co{A} and \Co{B} (represented as \Co{A:B}).
The notation can be used in any general way that makes sense: the operands of
\Co{*} can be any expression that evalutates to a category,
and crossed terms can be combined with other operators,
and parenthesized to indicate grouping.
With a small factorial experiment, we may want  fit the full model
including all possible effects.
This model is specified by combining all the terms with \Co{*} (an alternative
specification avoids writing an explicit formula at all, as we will show on page \pageref{allterms}).
In our previous example:
\begin{Example}
> aovall <- aov(Yield \Twiddle Temp * Conc * Cat, catalyst)
> aovall
Call:
aov(formula = Yield \Twiddle Temp * Conc * Cat, data = catalyst)

Terms:
                  Temp   Conc    Cat Temp:Conc Temp:Cat 
 Sum of Squares 1058.0   50.0    4.5       4.5    200.0
Deg. of Freedom      1      1      1         1        1

                Conc:Cat Temp:Conc:Cat 
 Sum of Squares      0.0           0.5
Deg. of Freedom        1             1

Degrees of freedom: 7 total; 0 residual
\end{Example}
The specified model fits three main effects, three two-way interactions
and one three-way interaction.
Notice that the term \Co{"Residuals"} did not appear.
When all interactions in this experiment are included in the model, no degrees
of freedom are left for residuals.

When the levels of a factor \Co{B} are meaningful only given the level of some
other factor \Co{A}, \Co{B} is said to be nested in \Co{A}.
The notation in a formula is \Co{A\ /\ B}, and the result is to fit first \Co{A}
and then the effects of \Co{B} in \Co{A}.
Main effects for \Co{B} are not meaningful in this case.
The operands to \Co{/} can also be arbitrary expressions.

An experiment on methods for firing naval guns was reported by 
Hicks (1973, page 194).
Two methods were tested by gunners corresponding to three different
physiques (slight, average, and heavy).
Nine gunners of each physique were divided into three teams, and each
team tested the two loading methods twice,
for a total of 36 runs.
The response was the number of rounds fired per minute.
The data is on data frame \Co{gun}; a sample of 6 runs is as follows:
\begin{Example}
   Method Physique Team Rounds 
 3     M1        A   T1   22.0
 4     M2        A   T1   14.1
14     M2        S   T3   12.5
25     M1        S   T2   26.9
29     M1        H   T2   23.7
34     M2        A   T3   16.0
\end{Example}
For this experiment, the \Co{Team} factor is only meaningful within
each physique, so one
possible model is:
\begin{Example}
Rounds \Twiddle Method + Physique/Team
\end{Example}
which will fit effects for \Co{Method}, \Co{Physique} and \Co{Team}
within \Co{Physique}:
\begin{Example}
> gunaov <- aov(R \Twiddle M + P/T, gun)
> summary(gunaov)
                   Df Sum of Sq Mean Sq    F 
            Method  1       652  651.95 1.25
          Physique  2        16    8.03 0.02
Team %in% Physique  6        39    6.54 0.01
         Residuals 26     13509  519.60 1.00
\end{Example}
Nested terms are often related to random effects models, as discussed
on page \pageref{random.effects} below.

\subsection*{Models without Formulas}\label{allterms}

When an experiment involves a substantial number of factors, writing a formula
to contain them all is tedious.
An alternative is to generate a model specification including all the possible terms,
up to a specified order of interaction.

The function \Co{allterms()} constructs a model specification from a data frame
or from specifications of the factors and response to be used.
As an illustration, suppose we want to fit
a model to the wafer etching experiment in Table \ref{wafer.table} on page \pageref{wafer.table}, assigned as the dataframe \Co{wafer}.
\begin{Example}
waov1 <- aov( allterms(wafer), wafer)
\end{Example}
The call to \Co{allterms()} replaces the formula specification.
When \Co{allterms()} is given a data frame as its argument, it constructs
a model using all factorsin the data frame as predictors and the first numeric
variable as the response.
In this case the result is to fit the response \Co{pre.mean} to the 8 factors
shown in Table \ref{wafer.table}.
By default, \Co{allterms()} generates a model with only main effects; by supplying
the argument \Co{order=2}, for example, we could have fit a model with all possible
interactions up to order 2.

To obtain more flexibility in the model, one can subset columns of the data
frame to change either the set of factors or the variable chosen as the
response:
\begin{Example}
waov2 <- aov( allterms(wafer[,-c(9,10)]), wafer)
\end{Example}
drops the first two responses, causing \Co{post.mean} to be chosen as the response.
The response can also be supplied as a separate argument, allowing transformed variables
to be used.
The \Co{response} argument is usually best supplied as a formula, only the left
side being used.
For example, to fit the log of the pre-etch standard deviation:
\begin{Example}
waov3 <- aov( allterms(wafer,
    response = log(post.dev) \Twiddle .), wafer)
\end{Example}
Model updating, discussed in section \ref{updating}, can also be used
to refine the specification, by adding or dropping individual terms.
The object returned by \Co{allterms()} is a \Co{terms} object.
Section \ref{terms} will discuss such objects;
in simple use as illustrated here, you need not worry about the object itself.

\subsection*{Multiple-Response Models}\label{multiresponse}

The examples described so far have dealt with single response variables.
It is also possible to supply a numeric matrix as the response variable.
Columns of the matrix are interpreted as the individual responses.
The fit carries through as before, with the distinctions that
the class of the object is now \Co{"maov"} and that the
effects, coefficients, residuals, and fitted values will be matrices
having as many columns as the response.

Multi-response models can be described in all the same ways as single-response
models.
As an example, we use the \Co{allterms()} function discussed in the previous
paragraphs, and fit two responses, \Co{pre.mean} and \Co{post.mean}, from
the \Co{wafer} data.
\begin{Example}
> wmaov <- aov( allterms(wafer,
+   response = cbind(pre.mean, post.mean) \Twiddle .), wafer)
\end{Example}
Summaries for multi-response models usually just repeat the univariate
summary for each response:
\begin{Example}
> summary(wmaov)

Response: Y1 
          Df Sum of Sq Mean Sq    F 
  maskdim  1      0.65   0.652 0.02
 visc.tem  1      1.33   1.326 0.03
   spinsp  2      0.77   0.383 0.01
 baketime  2      0.00   0.001 0.00
 aperture  2      0.03   0.016 0.00
  exptime  2      0.54   0.272 0.01
  devtime  2      0.28   0.140 0.00
 etchtime  2      0.10   0.052 0.00
Residuals  3    129.11  43.036 1.00

Response: Y2 
          Df Sum of Sq Mean Sq    F 
  maskdim  1      0.68   0.677 0.01
 visc.tem  1      2.51   2.512 0.04
   spinsp  2      1.42   0.712 0.01
 baketime  2      0.15   0.074 0.00
 aperture  2      0.03   0.013 0.00
  exptime  2      1.56   0.779 0.01
  devtime  2      1.00   0.499 0.01
 etchtime  2      0.10   0.048 0.00
Residuals  3    187.68  62.561 1.00
\end{Example}


\subsection*{Aliasing; Over-specified Models\label{aliasing}}

The specification of a model for the analysis of variance may
include more terms than can actually be estimated from the design.
In this case the fit goes through, but some of the terms or coefficients
specified in the formula may not appear in the fit.
In printing a summary of such models, a table of {\em aliasing} information
will be produced.
Let's consider a small example.  Suppose that instead of the complete
8 observations in the \Co{catalyst} data shown in Table \ref{bhh.table}
on page \pageref{bhh.table}, we could only afford 4 observations.
By asking for a fractional factorial design, and assuming the yields
were the same as in Table \ref{bhh.table}, we would obtain the following
data, in data frame \Co{catalyst2}, say:
\begin{Example}
> catalyst2
  Temp Conc Cat Yield 
1 180    20    A      60
2 160    40    A      54
3 160    20    B      52
4 180    40    B      80
\end{Example}
What happens if we fit the complete model in
the three factors?
\begin{Example}\label{half.aov}
> half.aov <- aov(Yield \Twiddle Temp*Conc*Cat, catalyst2)
> half.aov
Call:
aov(formula = Yield \Twiddle Temp * Conc * Cat, data = catalyst2)

Terms:
                Temp Conc Cat 
 Sum of Squares   289   121   81
Deg. of Freedom     1     1    1

Degrees of freedom: 3 total; 0 residual
\end{Example}
The fit includes three terms plus the intercept, using up all 4 degrees of
freedom available.
Looking at the summary gives some more information about what happened:
\begin{Example}
> summary(half.aov)
      Df Sum of Sq Mean Sq 
Temp  1       289     289
Conc  1       121     121
 Cat  1        81      81

Alias pattern:
       (Int.) Te Co Ca 
  TeCo               *
  TeCa            *   
  CoCa         *      
TeCoCa      *         
\end{Example}
The second portion of the summary is new, and is only printed
when the model asked for more estimates than the data could produce.
The alias pattern shows, in the rows of the table,
 what single degrees of freedom could not
be estimated because they were linearly dependent on things
already included in the model, shown as the columns of the table.
Each row expresses the corresponding linear dependency.
For example, the first row shows that the estimate for the \Co{Temp:Conc}
interaction was aliased to the \Co{Cat} main effect, already included.
Statistically, the implication is that the design cannot distinguish
the two terms.
Only if one is willing a priori to treat the \Co{Temp:Conc} interaction
as known to be zero can the sum of squares assigned to \Co{Cat} be
attributed unambiguously to that term.
For more discussion of aliasing, see page \pageref{more.aliasing}.

Incomplete designs are widely used and useful when experiments are
expensive to run or when many factors need to be studied simultaneously.
It's perhaps preferable in these cases to understand beforehand what
can be estimated, and to phrase the model to make the assumptions
plain.
In this case, for example, if we had specified only the main effects,
as in the initial example on page \pageref{cat.example1}, we would
have got the same fit, but without aliasing.
On the other hand, if you are unsure just what can be estimated
from a particular design, fitting a complete model and then examining
the alias pattern is one way to find out.


\subsection*{Error Terms\label{error.terms}}
The analysis of experiments with categorical variables departs
from standard linear models in
an important way when the model includes multiple error terms.
An example, described in Federer (1955, page 274) will
suggest the idea.
In an experiment on 8 varieties of guayule (a Mexican plant yielding rubber),
four different treatments were to be applied to the seeds.
The relevant response was the number of plants germinating per 100
seeds planted on an individual greenhouse flat.
The experimenter reasoned that plants grown on different flats were likely
to have differences due to the flats.
These differences are not of particular interest and can be modelled
as a random quantity depending on the individual flat.

In order to allow some relevant terms in the model to be more accurately
estimated, the experimenter divides each flat ({\em plot} in the classic
terminology) into subplots and assigns the other factors (variety and
treatment in this example) so as to allow estimation that is
orthogonal to the plots.
For the moment we are not concerned with the choice of this design (that
will come up in subsection \ref{design.functions}), but with
how the error model will be specified.
The data in the example consisted of 96 observations, with 8 levels
for varieties, 4 for treatments and 24 for flats; that is, there were
24 flats each supplying 4 observations.
The data for the first 4 flats is shown below:
\begin{Example}
> guayule[1:16,]
   variety treatment flats plants 
 1      V1        T1     1     66
 2      V1        T2     1     12
 3      V1        T3     1     13
 4      V1        T4     1      6
 5      V2        T1     2     77
 6      V2        T2     2     26
 7      V2        T3     2     27
 8      V2        T4     2     15
 9      V3        T1     3     51
10      V3        T2     3      8
11      V3        T3     3     20
12      V3        T4     3     10
13      V4        T1     4     52
14      V4        T2     4      4
15      V4        T3     4     19
16      V4        T4     4     13
\end{Example}
In the analysis of variance, we want to specify an error term  corresponding
to the levels of \Co{flats}, which we do by including the expression
\Co{Error(flats)} in the model formula:
\begin{Example}
> attach(guayule)
> gaov <- aov(plants \Twiddle variety * treatment + Error(flats))
\end{Example}
This model will produce a separate fit for the portion of the data
corresponding to the effects of the variable \Co{flats} and for the residuals
from these effects -- in classical anova terminology these are the
whole plot and subplot error strata.

As you might expect, since there is a separate fit for each error stratum,
the overall fit is represented as a list of \Co{aov} objects.
Its class is \Co{aovlist}, and methods exist for the usual functions
like \Co{print} and \Co{summary}:
\begin{Example}
> gaov

Call:
aov(formula = plants \Twiddle variety * treatment + Error(flats))

Stratum 1: flats

Terms:
                 variety Residuals 
 Sum of Squares  763.156  1415.833
Deg. of Freedom        7        16

Degrees of freedom: 23 total; 16 residual
Residual standard error: 88.48958 

Stratum 2: Within

Terms:
                treatment variety:treatment Residuals 
 Sum of Squares  30774.28           2620.14   1162.83
Deg. of Freedom         3                21        48

Degrees of freedom: 72 total; 48 residual
Residual standard error: 24.22569 
\end{Example}
Notice that the experimenter chose a design such that \Co{treatment}
effects, which the experiment was particularly anxious to estimate
accurately, were entirely orthogonal to \Co{flats}.
They do not appear in the first (less accurate) stratum at all.
The experimenter's assumption that flat-to-flat variation would
be substantial is supported by the residual standard error being
over three times smaller in the \Co{Within} stratum.

The error model specified can include more than one term; the
argument to \Co{Error} can be anything that would go on the
right side of a formula.
The example we used was, in fact, a little more complicated than
disclosed so far.
The 96 flats were themselves organized into 3 replicates of
a design with 24 observations divided into 6 plots.
Depending on how the experiment was run, one might wonder whether
flats within one replication would tend to be more alike than
those in different replications.
This reasoning suggests having two special error terms: \Co{reps}
and \Co{flats}.
The formula
\begin{Example}
plants \Twiddle variety * treatment + Error(reps + flats)
\end{Example}
Fitting with this error model produces an analysis divided into three
strata: \Co{reps}, \Co{flats} after removing \Co{reps}, and residuals
from both.
In this case,  \Co{flats} are actually defined within \Co{reps}, but nothing
in the computations requires such an arrangement.
\begin{Example}
> gaov2

Call:
aov(formula = plants \Twiddle variety * treatment +
     Error(reps + flats) )

Intercept: 25.30208 

Stratum 1: reps

Terms:
                Residuals 
 Sum of Squares  38.58333
Deg. of Freedom         2

Degrees of freedom: 2 total; 2 residual
Residual standard error: 19.29167 

Stratum 2: flats

Terms:
                 variety Residuals 
 Sum of Squares  763.156  1377.250
Deg. of Freedom        7        14

Degrees of freedom: 21 total; 14 residual
Residual standard error: 98.375 

Stratum 3: Within
      {\em As before }

\end{Example}
This experiment was designed so that
treatments, varieties, and their interaction were balanced within
each replication.
As a result, there are no effects from \Co{varieties*treatment} at all in the
\Co{reps} stratum.

\subsection*{Abbreviations in Formulas}\label{abbr.formula}

In anova models, and generally in the models in this book, variable
names can be abbreviated to a specified set of short forms.
The function \Co{abbr()} will generate the vector of shortest abbreviations
that distinguish a given set of names.
The argument to \Co{abbr()} can be a character vector or a data frame.
In the latter case the abbreviations are for the variable names in
the data frame.
For the \Co{catalyst} data used at the beginning of the section:
\begin{Example}
> abbr(catalyst)
 Temp Conc  Cat Yield 
 "T"  "Co" "Ca" "Y"  
\end{Example}
These are the abbreviations recognized by default in formulas.
An \Co{"abbreviations"} attribute may be assigned to the data
frame to provide any abbreviations you like.
Thus, the formula specified before as
\begin{Example}
Yield \Twiddle Temp * Conc * Cat
\end{Example}
could equally well have been given as
\begin{Example}
Y \Twiddle T * Co * Ca
\end{Example}
(or some of the names could have been written out fully and others
abbreviated).
Each name must be either exactly the abbreviation or exactly
the full form; for example \Co{Con} would not be recognized.
The abbreviation  algorithm can be used in
other applications as well; see the documentation for \Co{abbreviate()}.

The abbreviations attribute can be any (unique) set of character strings, and need not have
any relationship to the names.
If \Co{abbr()} is given a second argument, this is taken to be the
desired abbreviations.
For example, for those who like to write formulas using \Co{"A"}, \Co{"B"},
and so on, the following suggests a way:
\begin{Example}
> attr(catalyst,"abbreviations") <- abbr(catalyst,LETTERS)
> aov(D \Twiddle A + B + C, catalyst)
Call:
aov(formula = Yield \Twiddle Temp + Conc + Cat, data = catalyst)

Terms:
                  Temp   Conc    Cat Residuals 
 Sum of Squares 1058.0   50.0    4.5     205.0
Deg. of Freedom      1      1      1         4

Degrees of freedom: 7 total; 4 residual
Residual standard error: 51.25 
\end{Example}
Notice that the full names are used in the fitted model, regardless
of the abbreviation.

\subsection*{Random Effects Models\label{random.effects}}

Some factors in experiments have the property that the levels chosen
are not so much interesting in themselves as examples of ``typical''
values for the underlying variable.
Samples of counties within a state, households on a block, or
animal or human subjects for testing may have been chosen in the
hope that conclusions made from the experiment can be generalized
to the overall population from which the samples are chosen.
This assumption can be built into the model explicitly by treating
the coefficients or effects for each level of the factor as
a random variable, sampled from some population distribution
of such coefficients.
The term {\em random effects} is often used for such factors, emphasizing
that the effects are not now themselves regarded as interesting
parameters in the model.

The focus of the analysis now moves back from the specific
coefficients for the sampled levels of the factor to some
parametrization of the assumed underlying distribution for those
coefficients.
In a sense, this just carries further the general emphasis in anova
on attributing variation to terms in the model, rather than
on looking at coefficients or predicted values as was more typical
in Chapter \LinearChapter.
The standard distributional assumption is that the coefficients
for individual levels are distributed with zero mean and some
unknown standard deviation.
Summarizing the fitted model then includes
estimating the variances for those factors whose levels are
assumed chosen at random.
The fitting of the model is essentially unchanged, except that
some additional information is computed to allow
the \Co{summary()} method to estimate the standard deviations.

Random effects models are one case in which we will require
the design to be balanced in the relevant factors.
Computations for a general, unbalanced situation could be
developed, but they would be considerably more complicated, and the
underlying statistical analysis is always phrased in terms of
a balanced situation.
Older computations, in fact, are usually expressed in a way that
{\em assumes} balance; indeed, assumes essentially a complete
factorial design.
This we don't do.
Instead our software tests for balance, in the process of computing
the number of replications for each of the terms involved.
Subsection \ref{random.effects.theory} will outline what happens.

To obtain a fitted analysis of variance including random effects
information, use
\begin{Example}
raov(formula, data)
\end{Example}
instead of \Co{aov()}.
In addition, the relevant factors in the experiment have to be marked
as random.
The function \Co{is.random()} tests or sets the attribute
\Co{random} in a category or in variables of a design.
Let's do an example.
Box, Hunter and Hunter (1978, page 572) report an experiment in which
15 batches of a pigment were each sampled analytically twice and two
repeated analyses were performed from each sample, to measure the moisture
content.
Suppose we have this experiment in a data frame, \Co{pigment}
\begin{Example}
> pigment
   Batch Sample Test Moisture 
 1 B1      S1   T1         40
 2 B2      S1   T1         26
 3 B3      S1   T1         29
 4 B4      S1   T1         30
 5 B5      S1   T1         19
\vspace{.5em}
          \(\cdots\)
\end{Example}
By using \Co{is.random()} on the left of an assignment, we can set
the \Co{random} attribute of specific categories, or of all categories:
\begin{Example}
is.random(pigment) <- TRUE
\end{Example}
(Non-categories will not be changed by the assignment.)
Now we can fit an anova model, with random effects assumed for
all the factors.
\begin{Example}
> praov <- raov(Moisture \Twiddle Batch/Sample, pigment)
> summary(praov)
                  Df Sum of Sq Mean Sq 
            Batch 14    1210.9  86.495
Sample %in% Batch 15     869.8  57.983
        Residuals 30      27.5   0.917

Random Effects -- estimated variances:
 Batch Sample %in% Batch Residuals 
 7.128            28.533     0.917
\end{Example}
We said that the necessary information about replications was computed
automatically, with a built-in check for balance.
However, you can alternatively supply the replications information.
For very large problems this might save some time, but the more
likely use of this argument is to cheat about balance.
The random-effects summary uses a component \Co{ems.coef} attached
to the \Co{aov} object by \Co{raov()}.
You can make the summary behave as if the experiment was a different
(balanced) design, so far as the random effects are concerned, simply
by setting this component to the corresponding component from
a balanced design:
\begin{Example}
> unbalaov <- aov(formula, data)
> unbalaov$ems.coef <- balancedaov$ems.coef
\end{Example}
For small unbalance, this may be justified.



\subsection{Graphical Methods and Diagnostics\label{aov.graphics}}
There are a number of useful graphical methods for displaying highly structured
categorical data.
We have already seen an example of one of them in 
Figure \ref{bhh.figure},
which shows a plot of the means at each level of each factor
for the data given in Table \ref{bhh.table}.
It was obtained by
\begin{Example}
design.plot(catalyst)
\end{Example}
We now show further examples on the data from the Dallas experiment on solderability
described in Chapter \DataChapter.
The experiment, in data frame \Co{dallas.balanced}, has factors \Co{Opening},
\Co{Samt},
\Co{Stypeth},
\Co{Padtyp},
replication factor \Co{Panel},
and response \Co{skips}.

\Co{design.plot} can be used to show summary statistics 
other than means, and additional arguments can be supplied to
the summary function used.
Figure \ref{Ddp} is a summary plot of the 5 percent trimmed mean of the response for 
each level of each category in the data frame obtained by
\begin{Example}
design.plot(dallas.balanced, FUN = mean, 
  ARGS = list(trim = .05))
\end{Example}
\begin{figure}
\psfig{figure=dp.dallas.p,width=4in,silent=}
\caption{ \label{Ddp} 5 percent trimmed average number of solder skips at each level of each figure}
\end{figure}
Any function which returns a single value may be used.
If the response is not contained in the data frame, it may be specified by the
parameter \Co{RESPONSE}.

Besides summary statistics, we may want to look at the distributions of the
response structured by the design factors. 
Boxplots, shaded bar plots, and binary plots can be produced by 
\Co{factor.plot}.
The choice of representation depends on the number of values that the response can 
have.
Figure \ref{Dfpbox} shows four boxplots corresponding to the four factors of
\Co{dallas.balanced}: 

\begin{Example}
factor.plot(dallas.balanced[ , -5])
\end{Example}
\begin{figure}
\psfig{figure=bp.dallas.p,width=4in,silent=}
\caption{Box plots of solder skips by factors}
\label{Dfpbox}
\end{figure}

We remember from the previous description that this is a hierarchical
experiment with the factors 
\Co{Opening},
\Co{Samt}, and
\Co{Stypeth}
varying over boards, and 
\Co{Padtyp} and
\Co{Panel} varying within a board.
To look at the distribution of defects for each board 
we can derive a factor which denotes boards by combining the three factors
which vary over boards with the function \Co{interaction}.
\begin{Example}
Boards <- interaction(Opening, Samt, Stypeth)
\end{Example}
Then one way to evaluate the boards might be to display the
responses to emphasize 0 defects and many defects.
We construct a derived response, \Co{skips2}, which is 0 if \Co{skips} = 0, 
1 if \Co{skips} = 1 or 2, 2 if $3<= \Co{skips} <= 10$ and 3 if $\Co{skips} >10$.
Figure \ref{Dshbar} is a shaded bar plot that shows good boards (lots of white space) 
and bad boards (lots of black space).
\begin{Example}
factor.plot(Boards, RESPONSE = skips2, style = "shbar"")
\end{Example}
\begin{figure}
\psfig{figure=shbar.dallas.p,width=4in,silent=}
\caption{Shaded bar plot of derived function of solder skips}
\label{Dshbar}
\end{figure}
We can see immediately that the first board is the best with no defects, and
the last board is the worst with more than 10 defects everywhere.
Intermediately we can also see that within each \Co{Stypeth}, the combination
of \Co{S Opening} and \Co{Thin Samt} gives the most defects, and the combination of
\Co{L Opening} and \Co{Thick Samt} gives the fewest defects.
Rearranging the order of bars in the plot to group by \Co{Opening} and \Co{Samt}
\begin{Example}
Boards <- interaction(Stypeth, Samt, Opening )
\end{Example}
would emphasize this pattern.

Plotting parameters may be given as arguments to \Co{design.plot} and
\Co{factor.plot}.

A final type of display which may be useful is that of two-way plots of 
functions of the response.  These plots are related to the \Co{matplot()}
function (\s, page 67),
but compute a summary function and provide more appropriate labeling.
Which plots to compute for data sets with many factors may be determined
by looking at an analysis of variance table where the two-factor
interactions have been included in the model, or chosen from previous
knowledge of the factors.
Figure \ref{D2way} shows the mean plot for \Co{Panel} by \Co{Stypeth} and
the 95th-percent-point plot for \Co{Padtype} by \Co{Opening}.
\begin{Example}
int.plot(skips, x.factor = Panel, trace.factor = Stype.th)
int.plot(skips, quantile, Padtype, Opening, 
  ARGS = list(probs = .95) )
\end{Example}
\begin{figure}
\psfig{figure=int.dallas.p,width=4in,silent=}
\caption{ \label{D2way} Two-way plots of solder skips}
\end{figure}

\subsection*{Quantile-Quantile Plots of Effects}

Under the standard assumptions for linear models (independently
normal errors with constant variance), the effects computed in fitting
the model are also independently normal.
If true value of the coefficient for a particular contrast
were zero, the corresponding effect would have mean zero.
Thus, large values on a normal quantile-quantile plot of the effects
suggest important effects.
Since the signs of the effects are of secondary interest, the
tradition is to plot the absolute values of the effects (usually referred
to as a ``half-normal'' plot).
If there are a number of insignificant effects, these should form a roughly
linear pattern in the lower left of the half-normal plot; the slope
of this pattern estimates the residual standard errror in the fit.
The expression
\begin{Example}
qqpe(fit)
\end{Example}
produces one or more such plots, with \Co{fit} being a fitted model.
One plot will be produced for each response in a multiple-response model,
and separate plots of effects will be produced for each stratum in a
multiple-stratum analysis of variance.
A number of optional arguments are available to \Co{qqpe()}; one
that is frequently useful is \Co{label}.
If \Co{label=T}, the user will be prompted to locate interactively
(e.g., with a mouse) interesting points on the plot.
The corresponding effect labels will be plotted at the located position.


As an example, we show the half-normal plot of effects from the
analysis of post-etch mean in the \Co{wafer}  data discussed
in the previous subsection, with the five largest effects identified
interactively:
\begin{Example}
qqpe(wpostaov, identify=T)
\end{Example}

\begin{figure}
\psfig{figure=qqpe.p,width=3in,silent=}
\caption{ \label{qqpeplot} Half-normal plots of post-etch mean effects
in the wafer experiment.}
\end{figure}

\subsection{Generating Designs\label{design.functions}}


Two functions are provided to generate data frames describing
experimental designs.
In the simplest call, the functions take an argument \Co{levels},
the number of levels for each of the factors in the design:
\begin{Example}
fac.design(levels)
oa.design(levels)
\end{Example}
\Co{fac.design()} produces a factorial design, with a row of the
data frame for each possible combination of the factors in the
design.
To generate a $2^3$ design as in Table \ref{bhh.table} at the beginning
of the Chapter, we ask for a design with three factors, each at 2 levels:
\begin{Example}
> fac.design(c(2,2,2))
     A    B    C 
1 A1   B1   C1  
2 A2   B1   C1  
3 A1   B2   C1  
4 A2   B2   C1  
5 A1   B1   C2  
6 A2   B1   C2  
7 A1   B2   C2  
8 A2   B2   C2  
\end{Example}
\Co{oa.design()} produces one of the {\em orthogonal array} designs,
specially selected designs that aim to study a large number of factors
in as few runs as possbile, on the assumption that only main effects
and possibly some low-order interactions need to be included in the
model.
The wafer experiment in Table \ref{wafer.table}  is an example of
such a design.
The underlying design had 9 factors,
three with two levels and 6 with three levels.
As run, however, two of the two-level factors  were combined
in a special way to define a three-level factor, creating
an experiment with one two-level factor and 7 three-level factors.
One could ask for an orthogonal-array design of this form by:
\begin{Example}\label{orth.array}
> oa.design(c(2,rep(3,7)))
      A    B    C    D    E    F    G    H 
 1   A1   B1   C1   D1   E1   F1   G1   H1
 2   A1   B1   C2   D2   E2   F2   G2   H2
 3   A1   B1   C3   D3   E3   F3   G3   H3
 4   A1   B2   C1   D1   E2   F2   G3   H3
 5   A1   B2   C2   D2   E3   F3   G1   H1
 6   A1   B2   C3   D3   E1   F1   G2   H2
 7   A1   B3   C1   D2   E1   F3   G2   H3
 8   A1   B3   C2   D3   E2   F1   G3   H1
 9   A1   B3   C3   D1   E3   F2   G1   H2
10   A2   B1   C1   D3   E3   F2   G2   H1
11   A2   B1   C2   D1   E1   F3   G3   H2
12   A2   B1   C3   D2   E2   F1   G1   H3
13   A2   B2   C1   D2   E3   F1   G3   H2
14   A2   B2   C2   D3   E1   F2   G1   H3
15   A2   B2   C3   D1   E2   F3   G2   H1
16   A2   B3   C1   D3   E2   F3   G1   H2
17   A2   B3   C2   D1   E3   F1   G2   H3
18   A2   B3   C3   D2   E1   F2   G3   H1
\end{Example}
Orthogonal array designs are usually defined for factors with
two or three levels only.
The \Co{oa.design()} function uses a table of known designs
to handle a limited range of requests.
Outside of that range, it gives up:
\begin{Example}
> oa.design(rep(2,100))
Error in call to "oa.design": Can't generate an oa design
   that needs  > 32 runs
\end{Example}
Only two- and three-level factors are allowed, but other factors
can sometimes be accomodated by adjusting the design; for example,
a 4-level factor can be formed out of two two-level
factors, provided the interaction of those factors is estimable.

The names of the factors and the names of the levels of the factors
were generated in a standard form in the examples above.
These names may be provided by the \Co{names} argument.
If supplied, \Co{names} can either be a character vector or a list;
in either case it should have length equal to the number of factors.
A character vector provides names for the factors.
The names for the levels are then constructed from the abbreviated
form of the factor names:
\begin{Example}
> fac.design(c(2,2,2),names = c("Temp","Conc","Cat"))
  Temp  Conc   Cat 
1   T1   Co1   Ca1
2   T2   Co1   Ca1
3   T1   Co2   Ca1
4   T2   Co2   Ca1
5   T1   Co1   Ca2
6   T2   Co1   Ca2
7   T1   Co2   Ca2
8   T2   Co2   Ca2
\end{Example}
If \Co{names} is a list, its elements provide names for the levels
of the corresponding factors, and the names attribute of the list gives the
factor names.
\begin{Example}
> nlist <- list(Temp = c(160, 180), Conc = c(20, 40),
+  Cat = LETTERS[1:2])
> fac.design(c(2,2,2), names = nlist)
   Temp Conc Cat 
1   160   20   A
2   180   20   A
3   160   40   A
4   180   40   A
5   160   20   B
6   180   20   B
7   160   40   B
8   180   40   B
\end{Example}
The expression \Co{factor.names(design)} returns a similar list
giving the factor and level names of \Co{design}.
Factor names can be assigned to an existing design by
using \Co{factor.names()}  on the left of an assignment.
The same interpretation is given to the object on the right
of the assignment as was outlined for the \Co{names} argument:
\begin{Example}
> cdesign <- fac.design(rep(2,2,2))
> factor.names(cdesign) <- nlist
\end{Example}
This produces the same design as in the previous example.
A \Co{names} argument can be used in calling \Co{oa.design()} as well,
with the same interpretation.
Both functions also take an optional argument \Co{row.names} to provide
names for the rows of the resulting data frame:
\begin{Example}
> cdesign <- fac.design(rep(2,2,2),
+     row.names = paste("Run",1:8))
\end{Example}
The row names are the attribute \Co{"row.names"} of the resulting data
frame, and can be extracted or assigned as this attribute.

Factorial designs can be generated with multiple replications of
each combination of factor levels, or as fractional designs that include
only some of the possible combinations.
The optional argument \Co{replications}
generates the specified number of replications of the entire design.
Fractional factorials can be specified for designs where all factors
have two levels.
The argument \Co{fraction} defines what fractional design is wanted.
The argument can just be a numeric fraction, in which case the function
picks from a standard set of fractions:
\begin{Example}
> fac.design(rep(2,5),fraction=1/4)
     A    B    C    D    E 
1   A1   B2   C1   D1   E1
2   A2   B1   C2   D1   E1
3   A2   B2   C1   D2   E1
4   A1   B1   C2   D2   E1
5   A2   B1   C1   D1   E2
6   A1   B2   C2   D1   E2
7   A1   B1   C1   D2   E2
8   A2   B2   C2   D2   E2

Fraction: y \Twiddle A:B:D + B:C:E
\end{Example}
Alternatively, \Co{fraction} can specify
one or more \Co{defining contrasts}.
In terms of our model formulas, a defining  contrast is an interaction
of factors, such as \Co{A:B:C}.
In a two-level design, each such interaction defines one constrast
with two possible values.
Therefore each such contrast implicitly divides the design into two halves,
one of which will be chosen in forming the fractional design.
If one contrast is supplied, a half replicate is produced; if two
contrasts are supplied, a quarter replicate, and so on.
For each contrast, the function chooses the half for which the contrast
has a positive value, by default.
The optional argument \Co{which.half} can override this (see the documentation
for \Co{fac.design()}).
The printing method for designs will show the defining contrasts, whether
chosen automatically or specified in the call.

Defining contrasts should be supplied as formulas, which need not have
any lefthand side:
\begin{Example}
> fac.design(rep(2,5), fraction = \Twiddle A:B:C:D + B:C:D:E)
     A    B    C    D    E 
1   A1   B1   C1   D1   E1
2   A1   B2   C2   D1   E1
3   A1   B2   C1   D2   E1
4   A1   B1   C2   D2   E1
5   A2   B2   C1   D1   E2
6   A2   B1   C2   D1   E2
7   A2   B1   C1   D2   E2
8   A2   B2   C2   D2   E2

Fraction:  \Twiddle A:B:C:D + B:C:D:E
\end{Example}
This produces a quarter replicate, the portion of the full design
for which the two interactions specified both have positive sign.
It's important to notice that we used \Co{":"}, not \Co{"*"} in
defining the contrasts.
We wanted {\em only} the four-factor interactions, not the expanded
model with all the main effects, etc.
The latter defines a fractional factorial also, but not likely the
one we want.

\subsection*{Randomization}\label{randomization}

In carrying out an experiment, one may want to randomize the order
in which the runs are to take place.
A randomized order could be generated by
\begin{Example}
mydesign[ sample(1:nrow(mydesign), ]
\end{Example}
which permutes the rows of the design in a random order.
Unfortunately, things are not quite so simple.
In practice, some factors may be difficult to vary.
We want to restrict randomization to leave these factors alone,
so they can be varied as infrequently as possible.
The function \Co{randomize} takes a design and the names of some
factors to be restricted, and returns
an
ordering, first by levels of the restricted factors,
and then randomized within those levels.
Running the experiment in the order given by this permutation would
provide the restricted randomization requested.
For example, suppose we wanted to randomize the design in Table \ref{bhh.table},
but not over the \Co{Cat} factor:
\begin{Example}
> perm <- randomize(catalyst,"Cat")
> perm
[1] 1 3 4 2 8 6 5 7
> catalyst[perm,]
      Temp Conc Cat Yield 
1      160   20   A    60
3      160   40   A    54
4      180   40   A    68
2      180   20   A    72
8      180   40   B    80
6      180   20   B    83
5      160   20   B    52
7      160   40   B    45
\end{Example}
The result is to do all the \Co{"A"} catalyst runs first, then all
the \Co{"B"} catalyst, while randomly permuting within each level of
the catalyst factor.
The \Co{restrict} argument is essentially a subscript argument on
the factors and can have any form suitable for such a subscript; e.g.,
numeric or logical.
