\section{Computational Techniques}
This section presents some background on the computations shown in
previous sections, including numeric and statistical results
to justify the computations.
The purpose is to make clearer the computations available and how
they might be further extended.
The topic is a rich one and our discussion here is necessarily
brief and incomplete.
However, it should help to clarify what is going on and why.

\subsection{Basic Computational Theory}

There are many ways to approach computations for analysis of variance,
using various algebraic and/or geometric arguments.
Most older treatments do not give a general discussion, perhaps
feeling that the mathematical level would have to be too high for
the intended audience.
In the context of this book, however, we do need some general
results to delimit the validity of the computations and suggest
directions for extension.
Fortunately, the numerical algorithms themselves provide a natural
way to derive many of the results.
In particular, we will begin by laying out the basic algebra of the
computations, using the corresponding results from section \LinearChapter.4.

We start by rephrasing the statistical model appropriately for this chapter.
Suppose there are $k$ terms in the formula, with
expressions involving \Co{*} or \Co{/} expanded, the result simplified
to eliminate duplicate terms, and with the intercept if any included
as the first term.
The $n$ observations on the response $y$ can be written:
\begin{equation}\label{model.eqn}
y = t_1 + t_2 + \cdots + t_k + \epsilon
\end{equation}
where each $t_j$ represents the expected contribution from the $j$th term and
$\epsilon$ is the error, conventionally assumed independently normal with
zero mean and variance $\sigma^2$.
The model in (\ref{model.eqn}) is equivalent to a linear model with a model
matrix, say $X$.
The term $t_j$ uses
contrast definitions for all the factors appearing in this term, computed
to include all combinations of levels in the term except those
identically included in earlier terms, as discussed in section \ref{contrasts}.
The contrasts, applied to the dummy variables for the factors, create some columns of coded variables for each of the factors.
These columns are then multiplied together in all possible combinations
to create $X_j$, the submatrix of the model matrix corresponding to $t_j$. 
We can write the model matrix as:
\begin{equation}\label{X.eqn}
X = [ X_1  X_2  \cdots  X_k ]
\end{equation}
and let $X_j$ have $d_j$ columns.
Keep in mind that $d_j$ is the maximum number of degrees of freedom
for $t_j$, as determined from the form of the model (\ref{model.eqn}),
before the actual data is even examined.
The computable number of degrees of freedom may be less than $d_j$,
but if the design is complete or chosen so that all terms are estimable,
then $X$ will be of full rank.
Since the notation is a bit simpler, let us assume to begin with that $X$ has full rank,
$d = \sum d_j$,
and then note later the changes necessary in the general case.

The numerical properties of the computations are easy to see
in terms of the orthogonal decomposition discussed in section \LinearChapter.4.
This expresses $X$ in terms of an
$n$ by $n$ orthogonal matrix $Q$ and an $n$ by $d$ upper-triangular
matrix $R$ such that
\begin{eqnarray}\label{decomposition}
Q^t \cdot X & = & R \nonumber \\
Q^t \cdot y & = & c \nonumber \\
 & = & [ c_1 c_2 \cdots c_k c_{k+1} ]
\end{eqnarray}
where $c_j$ is the set of effects ($d_j$ of them by the assumption of full
rank) associated with term $t_j$, and $c_{k+1}$ is the set of residual
effects, if any.
These  effects are the key computational tool for standard anova summaries.
Their essential properties follow directly from the construction of the
orthogonal decomposition:
\begin{enumerate}
\item
The elements of $c$ are all uncorrelated, both within and between
terms, regardless of whether the design is balanced or not.
\item
The distribution of $c_j$ is unaffected by any of the preceding terms $t_{j'}$, $j' < j$;
in particular, $t_{j'}$ contributes
nothing to the expected value of $c_j$.  Again, this is true regardless of
balance.
\item
In  addition, {\em if} a following term
is balanced with respect to the $j$-th term, the same is true of it; that is,
the distribution of $c_j$ is unaffected by $t_{j'}$, $j' > j$.
\end{enumerate}
Properties  2 and 3 are worth showing algebraically.
To do so, write $Q$ in columns
corresponding to the terms,
\[
Q = [ Q_1  Q_2  \cdots  Q_k Q_{k+1} ]
\]
In the Householder algorithm,
$Q$ is formed from the columns of $X$ and not stored explicitly as an $n$ by $n$
matrix, but the matrix form is nevertheless fully defined.
We also
write out the linear model for $y$ in terms of coefficients corresponding
to the chosen contrasts,
\begin{equation}\label{beta.eqn}
y = X_1 \beta_1 + X_2 \beta_2 + \cdots + X_k \beta_k + \epsilon
\end{equation}
and partition the upper-triangular matrix $R$ into rows and columns
corresponding to the terms; that is, $R_{1 1}$ is the first $d_1$ rows
and $d_1$ columns, $R_{2 1}$ the next $d_2$ rows of the same columns,
and so on.
Since $R$ is upper-triangular,
\[
R_{j j'} = 0, j' < j
\]
In addition, if the $j'$-th term is balanced with respect to the $j$-th
term,
\[
R_{j j'} = 0, j' > j
\]
as well.
The second and third properties then follow from writing  out $c_j$:
\begin{eqnarray}\label{orth.effects}
c_j  & = & Q_j^t \cdot y \nonumber \\
 & = & Q_j^t \cdot ( X_1 \beta_1 + X_2 \beta_2 + \cdots + X_k \beta_k + \epsilon) \nonumber \\
 & = & \sum_{j'=1}^k R_{j j'} \cdot \beta_{j'} + Q_j^t \cdot \epsilon \nonumber \\
 & = & R_{j j} \cdot \beta_j + \cdots + R_{j k} \cdot \beta_k  + Q_j^t \cdot \epsilon
\end{eqnarray}
and $c_j = R_{j j} \cdot \beta_j$, if all the terms are balanced.

\subsection{Aliasing; Rank-deficiency; Confounding}

When the computed rank $r$ of $X$ is less than $d = \sum d_j$, the decomposition pivots
columns so that the first $r$ columns are linearly independent.
In this case, the numerical decomposition is written:
\begin{eqnarray}\label{alias.decomposition}
Q^t \cdot X \cdot P & = & R \nonumber \\
 & = & \left[
	\begin{array}{l l}
		R_{+} & R_a \\
		0 & 0
	\end{array}
	\right]
\end{eqnarray}
Here $P$ is a permutation matrix, representing the permutation of columns
described by the component \Co{pivot} returned from the decomposition and
$R_{+}$ is an $r$ by $r$ upper-triangular, nonsingular matrix.
The computation has decided that the diagonal elements in the triangular
matrix corresponding to the pivoted columns are all effectively
zero.

The matrix $R_a$ contains the aliasing information. To see this, let
$X_0$ be the columns of $X$ found to be linearly dependent; i.e., the
last $d-r$ columns of $X \cdot P$ in (\ref{alias.decomposition}).
Similarly, let $X_{+}$ be the first $r$ columns of $X \cdot P$.
and $Q_0$, $Q_{+}$ the corresponding columns of $Q$.
Then
\[
Q_{+}^t \cdot X_0 = R_a
\]
These are the effects from fitting $X_0$ to $X_{+}$.
In the sense developed in section \LinearChapter.4, the elements in $R_a$
define the contribution of the contrasts in $X_{+}$ to predicting the
contrasts in $X_0$.
But since the fit is perfect, the effects actually define $X_0$ in
terms of $X_+$.
An alternative way to express the alias would be in terms of the coefficients
from fitting $X_0$ to $X_{+}$, rather than the effects; namely,
\[
R_{+}^{-1} \cdot R_a
\]
Neither version is inherently right, but we have followed a general preference
for the effects as being more stable computationally in returning
$R_a^t$ as the definition of the aliasing in section \ref{more.aliasing}.

This discussion of aliasing, or complete confounding, makes clear the connection with partial confounding.
When aliasing occurs, the contrasts in $X_0$ are regressed on $X_{+}$ with effectively
zero residuals.
When the $j$-th contrast is not aliased, we can still consider fitting
it to the preceding $j-1$ columns of $X$, and it remains true that
the $j$-th column of $R$ represents the effects of this fit.
If the first $j-1$ elements
of this column are effectively zero, then the $j$-th column
is orthogonal to (unconfounded with) the preceding columns.
Otherwise,  partial confounding exists, and can be summarized by the effects, or by
some rescaled version of those effects, if we wish to remove the
scale of the $j$-th column of $X$.
Once again, there is no single definitive numerical summary.
The columns of $R$, whether or not rescaled, have the nice property that
if the $j'$-th element is
zero, this means that there is no effect of the $j'$-th
contrast in fitting the $j$-th contrast, which seems to match intuitively the idea
that there is no confounding between the two.

\subsection{Error Terms}
This subsection explains the computations used
for models with an explicit error term.
The method is nicely simple and general, and uses the
basic linear model computations in a neat recursive way.
The general problem is to fit an analysis of variance model
when the error is assumed to be the sum of
errors due to the levels of certain factors occuring in the
error model, as well as a common residual error. 

To begin, we need to write such a model explicitly.
Let's submerge the specifics of the terms
and just let $\mu$ stand for the
sum of all the terms in (\ref{model.eqn}) on page \pageref{model.eqn}, {\em except} for the error.
This is often referred to as the treatment model.
Suppose the error model has $s$ terms.
Then (\ref{model.eqn}) can be rewritten, in terms of individual elements,
\begin{equation}\label{error.eqn}
y_{i} = \mu_{i} + \sum_{j=1}^s \epsilon_{l(i,j)} + \epsilon_{i}
\end{equation}
where the notation $l(i, j)$ just means that
the $j$-th error term is at level $l(i, j)$
for the $i$-th observation.

The key to the computations is to write the element-wise model (\ref{error.eqn}) 
in terms of the equivalent model matrix, $B$.
In matrix form:
\[
y = \mu + B \cdot \epsilon_B + \epsilon
\]
But now we can compute an orthogonal matrix $Q_B$ from the decomposition of $B$.
As in (\ref{decomposition}), $Q_B$ applied to $y$ produces orthogonal effects, which we write $y^*$:
\begin{eqnarray}\label{error.effects}
Q_B^t \cdot y & = & y^* \nonumber \\
 & = & \mu^* + R_B \cdot \epsilon_B + Q_B^t \cdot \epsilon \nonumber \\
 & = & ( y_1^* y_2^* \cdots  y_s^* y_{s+1}^* )
\end{eqnarray}
where
\[
\mu^* = Q_B^t \cdot \mu
\]
Since $Q_B$ is orthogonal, $Q_B^t \cdot \epsilon$
has the same distribution as $\epsilon$ under the normal error assumptions.
The analogy with (\ref{decomposition}) extends to  equation (\ref{orth.effects}).
In particular, the $j'$-th error term contributes nothing to the
variance of $y_j^*$, for $j' < j$.
Statistically, this says that the analysis of the error model
produces a set of $b_j$ effects constituting the components
affected by the $j$-th error term (but not by previous
error terms), and a set of $n- \sum b_j$
components affected only by the intra-block errors.
The appropriate analyis then applies the treatment model separately
to these sets of effects, producing analyses for each that can be
related to the corresponding error estimates.

Computationally, this works out quite simply.
The quantity we have written as $\mu$ is really the sum of the
components from each of the terms in (\ref{model.eqn}) and can be written in
terms of the model matrix (\ref{X.eqn}):
\[
\mu =  X \cdot \beta
\]
and therefore
\[
\mu^* = X^* \cdot \beta
\]
where $X^* = Q_B^t \cdot X$.
Equation (\ref{error.effects}) fits $y_1^*$ to the first $b_1$ rows of $X^*$,
$y_2^*$ to the next $b_2$ rows, and so on.
Each fit may produce estimates for any of the terms in the treatment model,
which are to be compared to the relevant error estimates for
that fit, assuming there are degrees of freedom left for residuals in that
fit.
It is these separate fits that correspond to the error {\em strata}.

To summarize the computations:
\begin{enumerate}
\item
Compute the orthogonal transformation, $Q_B$, for the error model.
\item
Apply this to the response
and to the model matrix for the treatment model, producing
$y^*$ and $X^*$
\item
For each term in the error model, fit the elements of $y^*$ assigned
to that term
to the corresponding rows of $X^*$.
\item
The object representing the complete analysis of variance is the list of these
models.
\end{enumerate}

\subsection{Computations for Random Effects\label{random.effects.theory}}

The random effects model adds to the standard anova model the assumption
that the parameters $\beta_j$ for some or all of the terms $t_j$ are
themselves random variables; specifically, that $\beta_j$ has a normal
distribution with mean $0$ and variance $\sigma_j^2$, to be estimated.
The usual distribution for $y$, as in (\ref{beta.eqn}) on
page \pageref{model.eqn} becomes the {\em conditional} distribution
given $\beta$,
and $\beta$ is no longer a parameter to be estimated, but an intermediate
quantity on the way to estimating the $\sigma_j$.
In the sense that this reduces the number of parameters, the model is now
simpler.
The catch is that it is no longer a linear model.

If we assume that the error $\epsilon$ in (\ref{beta.eqn}) has a normal
distribution with mean $0$ and variance $\sigma_0^2$, and write the
residuals as
\[
z = y - \sum_{j=1}^k X_j \cdot \beta_j
\]
then the probability density for $z$ is
\[
\exp ( \sum_i z^2 / \sigma_0^2 ) / ( \sqrt{2 \pi } \sigma_0 )
\]
The distribution of $z$ from the random effects model, however, has a
density which is the product of this conditional density multiplied by
the density of $\beta$,
\[
\prod_j \exp ( { \parallel \beta_j \parallel }^2 / \sigma_j^2 ) / ( \sqrt{2 \pi } \sigma_j )
\]
where ${ \parallel \beta_j \parallel }^2$ is the  sum of the squares of the
elements of $\beta_j$.
For estimation purposes, let $\ell(\sigma)$, where
$\sigma = (\sigma_0, \sigma_1, \cdots, \sigma_k)$,
be the negative log-likelihood:
\begin{equation}\label{reffects.eqn}
\ell(\sigma ) = \sum_{j=0}^k \log \sigma_j + ( 
	\sum_{i=1}^n z_i^2/\sigma_0^2 + \sum_{j=1}^k { \parallel \beta_j \parallel }^2 / \sigma_j^2 ) /2
\end{equation}
In spite of the nonlinearity of the model, one could go ahead to estimate
$\sigma$ computationally.
A rather neat way to do this is to rewrite (\ref{reffects.eqn}) in terms of
$d = \sum d_j$ new pseudo-observations.
The first $d_1$ of them have values of $X_1$ equal to
\[
\sigma_0 {\bf I}_{d_1} / \sigma_1
\]
where $\bf I$ is the identity matrix.
The values of the other $X_j$ and of $y$ are all zero.
Similarly, the next $d_2$ pseudo-observations have nonzero values for $X_2$
and so on.
Letting $y_*$, $X_*$ and $z_*$ be the data with these $d$ pseudo-observations
appended, then (\ref{reffects.eqn}) can be written:
\[
\sum_{j=0}^k \log \sigma_j + ( 
	\sum_{i=1}^{n+d} z_i^{* 2}/\sigma_0^2
\]
For fixed $\sigma$ this is just a linear model in $\beta$ with data  $y_*$ and $X_*$.
The solution can be computed by the standard linear model methods.
Of course, the result depends on the values assigned to $\sigma$, so the
computational procedure would be to compute maximum-likelihood estimates
of $\sigma$, say by some of the techniques of Chapter \NonlinChapter.
The function to be optimized would not be (\ref{reffects.eqn}) but
\[
\begin{array}[t]{c}{\rm max } \\
\beta 
\end{array} \ell(\sigma )
\]
with the conditional maximization for $\beta$ done by a linear computation.
The linear step might take advantage of the special form of $X_*$; form example, by noting
that only the last $d$ observations change from one step to another.

Back when random effects models were introduced, however, this
was beyond the available computational (and perhaps, mathematical) capabilities.
The approach instead was to ignore the nonlinear nature of the model, to assume that
the experiment was balanced, and to equate the standard sums-of-squares terms
to their expected values under the random effects model.
The general idea can be expressed by looking at the variance of the effects, $c = Q^t \cdot y$.
In matrix form:
\[
\mbox{ var\ } c = \sigma_0^2 {\bf I}_n + R \cdot \Sigma \cdot R^t
\]
where $\Sigma$ is the variance matrix of $\beta$ from the random effects
model and $R$ is the $n$ by $p$ upper-triangular matrix from the decomposition\footnote{
See, for example, C. R. Rao, (1973, section 4a.11) for the variance of $y$, from which
the variance of $c$ follows}.
Under the random effects model $c$ has zero mean and the mean squares
of the effects for a particular term can be equated to its expectation.
In the case of nested, balanced models, this reduces to a formula in terms
of the replications of the various factors appearing in the model.
Our function \Co{raov()} computes these coefficients, and the summary method
for \Co{aov} objects solves for the estimates of $\sigma_j^2$.
See the functions \Co{replication()}, \Co{raov()}, and \Co{summary.aov()}
for details.
