\section[Advanced Use]{The S Functions: Advanced Use}
In this section, we present some more advanced techniques for
using the S functions to analyse data from experiments with
categorical variables and to adapt those functions for specialized
applications.

\subsection{Parametrization; Contrasts}\label{contrasts}

If we think of the anova model as implying coefficients
for each level of each factor or interaction, then the model
is inherently over-parametrized.
The sum of all the coefficients for a main effect is trying to
estimate the same  thing as the intercept term,
and the sum of the coefficients over one factor
in a two-factor interaction is estimating the same thing as a single level
of the other factor.
No matter how many observations we take, these parameters are functionally
aliased.
In addition, if the experiment does not include all the possible
combinations of levels for the factors in the model,
even parameters that are not functionally aliased may be aliased in
the design.

The resolution of these ambiguities is done in two steps.
First, a parametrization of the model is constructed according
to a choice of {\em contrasts} for each factor appearing in the model,
to eliminate functional aliasing.
Second, the numerical fitting of the model checks for
and identifies any further design-dependent
aliasing.

For the first step, either the user or the software
needs to choose a set of linearly independent contrasts; that is,
for a factor with $k$ levels, $k-1$
linear combinations of the levels that are independent of each
other and of the vector of all \Co{1}'s.
These define $k-1$ new parameters which can all be estimated,
if the design contains data on all levels of the factor.
For statistical purposes, it is preferrable that the contrasts
be orthogonal to each other and to the vector of \Co{1}'s;
in particular, the sum of the elements of each contrast should be zero.
If the contrasts are not orthogonal, there will be correlations among the
coefficients for a particular term.
If the contrasts do not sum to zero,
the coefficients will also be correlated with earlier terms
in the model.
However, the software does not enforce either of these requirements
and some standard choices of contrasts do not satisfy them.

In some sense, all valid parametrizations are equivalent, in that
they produce the same fitted values, the same residuals, and
the same sum of squares
for each term in the model.
Effects are computed to be orthogonal regardless of the
parametrization.
However, the experimenter often wants to study particular combinations
of factor levels, in which case the 
choice of contrasts is important,
and interpretation of individual effects as shown in a q-q plot of
effects will be easier if the individual contrasts are meaningful.

Contrasts are an attribute of the
individual factors.
The function \Co{contrasts()} extracts or sets this attribute.
If the attribute is not \Co{NULL}, it should be a matrix whose columns
are the chosen contrasts; that is, for a factor with $k$ levels,
the contrasts are a $k$ by $k-1$ matrix.
A number of functions generate common choices for contrasts:
\begin{Example}
	contr.sum(n)
	contr.poly(n)
	contr.helmert(n)
	contr.treatment(n)
\end{Example}
The argument to these functions is either
the levels attribute or just the number of levels.
As an example, let's take a category, \Co{region}, that specifies the geographic
regions (north-east, south, north-central, or west) to which the
states of the United States belong.
\begin{Example}
> l <- levels(region);  l
> [1] "NE" "S"  "NC" "W" 
\end{Example}
Contrasts for this factor will be 4 by 3 matrices:
\begin{Example}
> contr.sum(l)
   NE  S NC 
NE  1  0  0
 S  0  1  0
NC  0  0  1
 W -1 -1 -1
> contr.poly(l)
           .L   .Q         .C 
NE -0.6708204  0.5 -0.2236068
 S -0.2236068 -0.5  0.6708204
NC  0.2236068 -0.5 -0.6708204
 W  0.6708204  0.5  0.2236068
> contr.helmert(l)
    S NC  W 
NE -1 -1 -1
 S  1 -1 -1
NC  0  2 -1
 W  0  0  3
> contr.treatment(l)
   S NC W 
NE 0  0 0
 S 1  0 0
NC 0  1 0
 W 0  0 1
\end{Example}
The interpretation of the various contrast functions is as follows:
\begin{itemize}
\item
\Co{contr.sum} produces a parametrization that contrasts each of the
first $k-1$ levels with the last.
This corresponds to a traditional anova parametrization got by
applying the constraint that the sum of the coefficients be zero.
\item
\Co{contr.poly} chooses orthogonal polynomials (linear, quadratic, $\cdots$)
up to degree $k-1$ over the points given as argument to \Co{contr.poly},
with equally-spaced points chosen as default.
\item
\Co{contr.helmert} is a set of contrasts that take the difference
of level $j, j = 2, \cdots, k$ versus the average of the previous
$j-1$ levels.
\item
\Co{contr.treatment} just uses coefficients
levels $2, \cdots, k$.
This corresponds to models in which the first level of all factors
is regarded as the control, with a coefficient defined to be zero.
\end{itemize}
The first and second are the default choice for contrasts for
unordered and ordered factors, respectively.
Only the second and third are orthogonal contrasts, and the
last is not even orthogonal to the mean.

Contrasts can be set by assigning the contrast attribute as the
value of one of these functions, or of any suitable contrast matrix:
\begin{Example}
contrasts(region) <- contr.helmert(levels(region))
\end{Example}
It is possible, also, to assign fewer than $k-1$ contrasts if the
remaining contrasts are of no particular interest.
There are two distinct ways to do this:
\begin{Example}
contrasts(x) <- value
attr(x,"contrasts") <- value
\end{Example}
The first version uses a special assignment function that checks the
value and fills it out to give a complete parametrization.
The second version assigns {\em only} the supplied value; it must be used
more carefully, but allows including only a partial term
in the analysis.

Suppose we wanted to make sure the first effect for a four-level factor
contrasted the average of the first and third against the average
of the second and fourth levels:
\begin{Example}
> contrasts(region) <- c(1,-1,1,-1)
> contrasts(region)
   C1    2    3 
NE  1 -0.7 -0.1
 S -1  0.1 -0.7
NC  1  0.7  0.1
 W -1 -0.1  0.7
\end{Example}
The assignment function for contrasts has
appended two additional, orthogonal columns to the supplied constrast.
The assignment also ensures that the contrast matrix will have
column names, for use in forming coefficient names.
The supplied contrasts have either user-supplied names or \Co{"C1"}, etc.
The remaining columns are distinguished by being labelled with the column
number only.

For an example of a partial term, we look again at the wafer data in Table \ref{wafer.table}.
The factor called \Co{visc.tem} was a three-level factor concocted from
three of the four possible levels got from combining two levels each
for viscosity (204 and 206) and baking temperature (90 and 105).
In analyzing this factor, one needs to choose contrasts specially.
If one chose to assume that
the temperature factor had no effect, then the quadratic effect from three-level
orthogonal polynomial contrasts turns out to be equivalent to the single
contrast for viscosity.
One way to install this assumption in the anova fitting would then be:
\begin{Example}
> levs <- levels(wafer$visc.tem)
> visc.contrast <- contr.poly(levs)[, 2, drop=F]
> visc.contrast
                .Q 
 204,90  0.4082483
 206,90 -0.8164966
204,105  0.4082483
\end{Example}
Contrast attributes must be matrices with as many rows as levels in
the factor, so we extracted the single column of the contrasts with
\Co{drop=F} to keep it a matrix.
Now we assign the contrasts attribute and fit the \Co{post.mean}
response:
\begin{Example}
> wterms <- allterms(wafer, response= post.mean \Twiddle .)
> attr(wafer$visc.tem,"contrasts") <- visc.contrast
> wvaov <- aov(wterms, wafer)
> summary(wvaov)
          Df Sum of Sq Mean Sq     F 
  maskdim  1    0.6767  0.6767 24.24
 visc.tem  1    2.5122  2.5122 89.99
   spinsp  2    1.4244  0.7122 25.51
 baketime  2    0.1477  0.0739  2.65
 aperture  2    0.0254  0.0127  0.46
  exptime  2    1.5577  0.7789 27.90
  devtime  2    0.9970  0.4985 17.86
 etchtime  2    0.0960  0.0480  1.72
Residuals  3    0.0838  0.0279  1.00
\end{Example}
Notice that this time only one degree of freedom goes to \Co{visc.tem},
the other  being included in the residuals.


\subsection{Terms Objects}\label{terms}
In the process of fitting models, the formula given to \Co{aov()} is expanded to
an object, of class \Co{terms}, containing information about the individual terms that
will appear in the model.
The terms object was introduced in Chapter \ModelChapter, but it comes into its own
with anova models, where more structured terms arise in crossed, nested, and
multi-strata models.
This subsection illustrates the information in the terms object that can be used
to handle special modeling situations.
Terms objects are {\em not} needed to specify models that can be expressed
easily by a formula, so this subsection will be useful mainly in situations
where computations need to depend on attributes of individual
terms or where repeated fitting needs to modify the model in special ways.
A \Co{terms} object can always be provided as the \Co{formula} argument to 
\Co{aov()}.
Doing so saves some computation and allows finer control
of the model in some situations.

A \Co{terms} object is a vector of mode \Co{"expression"} whose elements
are the individual terms, expanded and possibly simplified from the
formula.
The function \Co{terms()} generates the object; e.g.,
\begin{Example}
myterm <- terms(Yield \Twiddle Temp * Conc * Cat)
\end{Example}
The object has several special attributes, the more frequently useful
being:
\begin{description}
\item[\Co{factors}:] 
A matrix with one column per term, nonzero elements in
a column saying which variables
appear in that term.
Elements of the matrix are 1 if the corresponding factor
(with $k$ levels, say) should enter
as $k-1$ contrasts and 2 if it should enter as $k$ dummy variables.
\item[\Co{variables}:]
All the variables appearing in the formula
returned as a vector of mode \Co{"expression"}.
The response is included.
These variables define the rows of
the \Co{factors} attribute.
\item[\Co{order}:]
The order of the terms, 1 for main effects, 2 for two-way interactions, etc.
\item[\Co{intercept}:]
1 if the model  includes an intercept, 0 otherwise.
\end{description}
There are two optional arguments to \Co{terms()}: \Co{abb} allows
explicit specification of permitted abbreviations for variable names
in the formula, and \Co{specials} names some special functions that
are to be noted in the \Co{terms} object.
When \Co{abb} is included in the call to \Co{terms()}, it should
be an object like that returned by \Co{abbr()}, giving the abbreviations
as a character vector with the corresponding full names as the \Co{"names"}
attribute.
When the \Co{specials} argument is supplied and any of the special functions
appear in the formula, a corresponding \Co{specials} attribute will be
included in the terms object, identifying each of the occurrences.
See the detailed documentation of \Co{terms()} for an explanation of the
form of this component.

\subsection{Aliasing; Confounding; Balance}\label{more.aliasing}

These three concepts are used in anova to refer to aspects of linear dependency
and
correlation among the predictor variables in the models.
Discussions in the anova literature tend to be written from a restricted view
of possible designs.
We try to be more general, which means that the functions corresponding
to the concepts sometimes return the results in a more general form
than is  typical.
This subsection describes the
general form of the results, and supplements this with some examples of
how the results can be summarized in a simpler, perhaps more familiar way.

Aliasing refers to the property that some terms in the model may
be unestimable, wholly or in part, because they are {\em aliased} to terms that appeared
earlier in the model.
Numerically, this says that an effect or coefficient
that would have been estimated is found to be equal to some
linear combination of previously estimated parts of the model.
To be precise, a column in the model matrix representing a
contrast for one of the terms
is found to be linearly dependent on previous columns of that matrix.
Since we are not directly interested in the model matrix, we will express
this by saying that the  contrast itself is aliased.

Aliasing then is represented by a matrix of dependencies; each
row equates an unestimable contrast to a linear combination
of previously included contrasts.
The expression \Co{alias(fit)} returns the aliasing pattern
appropriate to the object \Co{fit}: for a linear model or an anova
with a single  error stratum the result is a matrix or
\Co{NULL}; for a model
with multiple error strata the result is a list of the non-null alias
matrices for the strata.

The numerical definition of aliasing is needed for further computations,
but for seeing the pattern visually, some simplified form of display
may be preferrable.
The \Co{summary()} method for anova, as shown on page \pageref{half.aov},
just prints an array with the character \Co{*} where some dependencies
exist.
It does this  by using a function \Co{prpattern()} that prints such a pattern
for any numeric object or list of numeric objects:
\begin{Example}
> prpattern(alias(half.aov))
       (Int.) Te Co Ca 
  TeCo               *
  TeCa            *   
  CoCa         *      
TeCoCa      *         
\end{Example}
In this example the interpretation of the alias pattern is very simple,
since each unestimable contrast is exactly aliased with one earlier
contrast.
For example, one can say that the two-factor interaction \Co{TeCo}
is {\em aliased} with or {\em confounded} with the main effect \Co{Ca},
hence the terminology.
When several nonzero coefficients appear in single rows of the alias
matrix, interpretation may be more difficult.

Confounding is used more generally to refer also to the situation
that two contrasts are correlated, but not exactly linearly dependent.
The terms partial and complete confounding are used to distinguish
the cases.
The numerical treatment of partial confounding and aliasing are
quite different, although the intuitive idea of aliasing as a limiting
case of confounding is quite reasonable numerically.
For partial confounding, some numerical measure of the partial dependence
of each coefficient on previous coefficients would be the
obvious analogue to the summary of aliasing.
A number of possible measures could be used; the one we use gives
in each row the effects obtained from fitting the corresponding single
variable in the model matrix to the previous terms.
The key point, qualitatively, is that the measure of counfounding
should be zero when the two variables in question are orthogonal
to one another.

To examine partial confounding, we return to the wafer-processing example in
Table \ref{wafer.table} on page \pageref{wafer.table}.
The analysis in the referenced paper, and our own analysis so far, have ignored
the fact that three of the 18 runs were done on one wafer, meaning that they
had 5 instead of 10 repeated observations to compute the mean and standard deviation
of the line widths.
Under these circumstances, the usual linear model computations would weight the
observations proportionally, since the variance of the mean line width would
be inversely proportional to the number of repeated observations.
This is easy to do with the \Co{aov()} function.
Any arguments not recognized by \Co{aov()} are passed on to the linear model function.
In particular, we can provide weights to be used in the fitting.
\begin{Example}
wwaov <- aov( wterms, wafer, weight = N)
\end{Example}
The weighted linear model is not balanced, and we can examine the
pattern of confounding.
For example, the following function extracts the triangular part of
the decomposition, removes the portion dealing with the intercept,
sets the diagonal to zero and uses \Co{prpattern()} to show the
result:
\begin{Example}
prconfounding <- function(fit) \{
    r <- (fit$R)[-1, -1]    # without the intercept
    p <- dim(r)[2]
    ii <- seq(1, p * p, by = p + 1)    # diagonal elements
    r[ii] <- 0
    prpattern(r, round.to = .Options$digits - 2)
    invisible(r)
\}
\end{Example}
If we apply it to the previous fit, the only confounding printed
is between the two levels of the aperture factor, because of non-orthogonal
contrasts.
With the unbalanced weighting, nearly all the effects are confounded:
\begin{Example}
> prconfounding(wwaov)
    A B C.L C.Q D.L D.Q E1 E2 F.L F.Q G.L G.Q H.L H.Q 
  A   *   *   *   *   *  *  *   *   *   *   *       *
  B       *   *   *   *  *  *   *   *   *   *       *
C.L           *   *   *  *  *   *   *   *   *       *
C.Q               *   *  *  *   *   *   *   *       *
D.L                   *  *  *   *   *   *   *       *
D.Q                      *  *   *   *   *   *       *
 E1                         *   *   *   *   *       *
 E2                             *   *   *   *       *
F.L                                 *   *   *       *
F.Q                                     *   *       *
G.L                                         *       *
G.Q                                                 *
H.L                                                  
H.Q                                                  
\end{Example}
The general pattern of counfounding is typical of missing or unequal
weighting.
Some nested designs tend to produce less global patterns of confounding.

Balance is a term used originally to describe a relation between pairs
of factors in
which each pair of levels is replicated the same number of times
in the design,
in particular the situation in which each level of a treatment factor
appeared the same number of times in each block (that is, level of
a blocking factor).
If all the pairs of factors are balanced in this sense, the
variables in the model matrix will be orthogonal.
We provide a function, \Co{replications()}, that presents
the relevant information about balance as part of a more general result.
The arguments to  \Co{replications()} are a formula and an optional
data frame.
The value of the function is in general a list of tables.
For each term in the formula,
the function computes a table of the number of replications of
each distinct level defined by the term
(a one-way table for main effects, a two-way table for two-factor interactions,
and so on).
If the term is balanced in the design, all the numbers in this table
will be equal.
In this case \Co{replications()} replaces the table by this single
number.
And if all terms are balanced, \Co{replications()} replaces the list
by a numeric vector containing for each term the number of replications.
The response in the formula given to \Co{replications()} is irrelevant
and can be omitted.

