In this lesson, we will examine three new ANOVA models (Models 1, 2, 3), as well as the effects model (Model 4) from the previous lesson, defined as follows:

Each of these four models can be written as a general linear model (GLM): \(\mathbf{Y}=\mathbf{X}\boldsymbol{\beta}+\boldsymbol{\mathcal{E}}\) simply by changing the design matrix \(\mathbf{X}\). Thus to perform the data analysis, in terms of the computer coding instructions, the appropriate numerical values for the \(\mathbf{X}\) matrix elements will need to be inputted.

A note on how ANOVA is calculated: In the past lessons, we carried out the ANOVA computations conceptually in terms of deviations from means. For the calculation of total variance, we used the deviations of the individual observations from the overall mean, while the treatment SS was calculated using the deviations of treatment level means from the overall mean, and the residual or error SS was calculated using the deviations of individual observations from treatment level means. In practice, however, to achieve higher computational efficiency, SS for ANOVA is computed utilizing the following mathematical identity:

\(SS=\Sigma(Y_i-\bar{Y})^2=\Sigma Y_{i}^{2}-\dfrac{(\Sigma Y_i)^2}{N}\)

This identity is commonly called the "working formula" or "machine formula". The second term on the right-hand side is often referred to as the correction factor (CF).

When computing the SS for the total variance of the responses, the formula above can be used as it is, but modifications need be made for the others. For example, to compute the treatment SS, the above equation has to be modified as:

\(SS_{treatment} = \sum_{i=1}^{T}\dfrac{\left( \sum_{j=1}^{n_i}Y_{ij} \right)^2}{n_i}-\dfrac{(\sum_{i=1}^{T} \sum_{j=1}^{n_i} Y_{ij})^2}{N}\)

These working formulas will be utilized throughout the remaining course material. 

To understand how various facades of the model relate to each other, let us look at a toy example with 3 treatments (or factor levels) and 2 replicates of each treatment. 

We have 6 observations, which means that \(\mathbf{Y}\) is a column vector of dimension 6 and so is the error vector \(\boldsymbol{\mathcal{E}}\) where its elements are the random error values associated with the 6 observations. In the GLM model, \(\mathbf{Y}=\mathbf{X}\boldsymbol{\beta}+\boldsymbol{\mathcal{E}}\), the design matrix \(\mathbf{X}\) for the overall mean model turns out to be a 6-dimensional column vector of ones. The parameter vector, \(\boldsymbol{\beta}\) is a scalar equal to \(\mu\), the overall population mean.

That is,

\(\mathbf{Y}=\begin{bmatrix}
2\\
1\\
3\\
4\\
6\\
5
\end{bmatrix} \text{,}\;\;\;\;
\mathbf{X}=\begin{bmatrix}
1\\
1\\
1\\
1\\
1\\
1
\end{bmatrix} \text{,}\;\;\;\;
\boldsymbol{\beta}=[\mu]  \text{,}\;\;\;\; \text{ and } 
\boldsymbol{\epsilon}=\begin{bmatrix}
\epsilon_1\\
\epsilon_2\\
\epsilon_3\\
\epsilon_4\\
\epsilon_5\\
\epsilon_6
\end{bmatrix}\)

Using the method of least squares, the estimates of the parameters in \(\boldsymbol{\beta}\) are obtained as:

\(\hat{\boldsymbol \beta}= (\mathbf{X}^{'} \mathbf{X} ) ^{-1} \mathbf{X}^{' } \mathbf{Y}\)

Using the estimate \(\hat{\boldsymbol{\beta}}\), the \(i^{th}\) predicted response \(\mathbf{\hat{y}}_i\) can be computed as \(\mathbf{\hat{y}}_i\)=\(\mathbf{x}_i^{'} \hat{\boldsymbol{\beta}} \) where \(\mathbf{x}_i^{'}\) denotes the \(i^{th}\) row vector of the design matrix.

In this simplest of cases, we can see how the matrix algebra works. The term \(\mathbf{X}^{'} \mathbf{X}\) would be:

\(\begin{bmatrix}
1 & 1 & 1 & 1 &1 & 1
\end{bmatrix}\ast \begin{bmatrix}
1\\
1\\
1\\
1\\
1\\
1
\end{bmatrix}=1+1+1+1+1+1=6 = n\)

The term \(\mathbf{X}^{'} \mathbf{Y}\) would be:

 \(\begin{bmatrix}
1 & 1 & 1 & 1 &1 & 1
\end{bmatrix}\ast \begin{bmatrix}
2\\
1\\
3\\
4\\
5\\
6
\end{bmatrix}=2+1+3+4+5+6=21 = \Sigma Y_i\)

So in this case, the estimate \(\hat{\boldsymbol{\beta}}\) as expected is simply the overall mean = \(\hat{\boldsymbol{\mu}} = \bar{y}_{\cdot\cdot} = 21/6 = 3.5\)

Note that the exponent of \(\mathbf{X}^{'} \mathbf{X}\) in the formula above indicates arithmetic division as \(\mathbf{X}^{'} \mathbf{X}\)  is a scalar in this case. In the more general setting, the superscript of '-1 ' indicates the inverse operation in matrix algebra.

To perform these matrix operations in SAS IML, we will open a regular SAS editor window, and then copy and paste three components from the file (IML Design Matrices):

The cell means model does not fit an overall mean, but instead fits an individual mean for each of the treatment levels. Let us examine this model for the same toy dataset assuming that a pair of observations arise from \(T = 3\) treatment levels. Each column will represent a specific treatment level using indicator coding, a ‘1’ for the rows corresponding to the observations receiving the specified treatment level, and a ‘0’ for the other rows.

To write the cell means model as a GLM, let \(\mathbf{X} = \begin{bmatrix} 1 & 0 & 0\\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} \mathbf{x}_1^{'} \\ \mathbf{x}_2^{'} \\ \mathbf{x}_3^{'} \\ \mathbf{x}_4^{'} \\ \mathbf{x}_5^{'} \\  \mathbf{x}_6^{'} \end{bmatrix}\), where \(\mathbf{x}_i^{'}\) is the \(i^{th}\) row vector of the design matrix. It can be seen that \(r=2\) is the number of replicates for each treatment level. Observe that column 1 generates the mean for treatment level 1, column 2 for treatment level 2, and column 3 for treatment level 3.

The parameter vector \(\boldsymbol{\beta}\) is a 3-dimensional column vector and is defined by,

\(\boldsymbol{\beta}= \begin{bmatrix}
\beta_0\\
\beta_1\\
\beta_2
\end{bmatrix} = \begin{bmatrix}
\mu_1\\
\mu_2\\
\mu_3
\end{bmatrix}.\)

The parameter estimates \(\hat{\boldsymbol{\beta}}\) can again be found using the least squares method. One can verify that \(\mu_{i} = \bar{y}_i\), the \(i^{th}\) treatment mean, for \(i=1,2,3\). Using this estimate, the resulting estimated regression equation for the cell means model is,

\(\hat{\mathbf{Y}} = \mathbf{X}\hat{\boldsymbol{\beta}}\)

which produces \(\hat{\mathbf{y}}_i = \mathbf{x}_i^{'}\begin{bmatrix} \hat{\mu}_1 \\ \hat{\mu}_2 \\ \hat{\mu}_3 \end{bmatrix}\).

The GLM applied to data with categorical predictors can be viewed from a regression modeling perspective as an ordinary multiple linear regression (MLR) with ‘dummy’ coding, also known as indicator coding, for the categorical treatment levels. Typically software performing the MLR will automatically include an intercept, which corresponds to the first column of the design matrix and is a column of 1's. This automatic inclusion of the intercept can lead to complications when interpreting the regression coefficients (discussed below).

Consider again the toy example with a pair of observations from 3 different treatment levels. The design matrix for this dummy variable model is as follows:

\(\bf{X} = \begin{bmatrix}
1 & 1 & 0 \\
1 & 1 & 0 \\
1 & 0 & 1 \\
1 & 0 & 1 \\
1 & 0 & 0 \\
1 & 0 & 0
\end{bmatrix}  \)
 
Notice that in the above design matrix, there are only two indicator columns even though there are three treatment levels in the study. If we were to have a design matrix with another indicator column representing the third treatment level (as seen below), the resulting 4 columns would form a set of linearly dependent columns, a mathematical condition which hinders the computation process any further. This over-parameterized design matrix would look as follows: 

\(\bf{X} = \begin{bmatrix}
1 & 1 & 0 & 0 \\
1 & 1 & 0 & 0\\
1 & 0 & 1 & 0\\
1 & 0 & 1 & 0\\
1 & 0 & 0 & 1\\
1 & 0 & 0 & 1
\end{bmatrix}  \)

This matrix containing all 4 columns has the property that the sum of columns 2, 3, and 4 will equal the first column representing the intercept. As a result, the matrix is what is known as 'singular' and the matrix computations will not run. To remedy this, one of the treatment levels is omitted from the coding in the correct design matrix above and the eliminated level is called the ‘reference’ level. In SAS, typically, the treatment level with the highest label is defined as the reference level (treatment level 3 in our example).

Note that the parameter vector for the dummy variable regression model is, \(\boldsymbol{\beta} = \begin{bmatrix}
\mu\\
\mu_1\\
\mu_2
\end{bmatrix}\)

In the effects model that we discussed in lesson 3, the treatment means were not estimated, but instead, the deviations of treatment means from the overall mean were estimated (the \(\tau_i\)'s). The model must include the overall mean, thus we need to include \(\mu_., \tau_1 , ... , \tau_T\) as elements in the parameter vector, \(\boldsymbol{\beta}\) of the GLM model. Note that, in the case of equal replication at each factor level, the deviations satisfy the following constraint:

\(\sum_{i=1}^{T} \tau_i =0\)

This implies that one of the \(\tau_i\)  parameters is not needed since it can be expressed in terms of the other T - 1 parameters and need not be included in the \(\boldsymbol{\beta}\) parameter vector. We shall drop the parameter \(\tau_T\) from the regression equation as it can be expressed in terms of the other T - 1 parameters \(\tau_i\) as follows:

\(\tau_T\ = -\tau_1 - \tau_2 - \cdots \tau_{T-1}\)

Thus, the \(\boldsymbol{\beta}\) vector of the GLM is a \(T \times 1\) vector containing only the parameters \(\mu, \tau_1 , ... , \tau_{T-1}\) for the linear model.

To illustrate how a linear model is developed with this approach, consider again a single-factor study with T = 3 factor levels with \(n_1=n_2=n_3=2\). The \(\mathbf{Y}\), \(\mathbf{X}\), \(\boldsymbol{\beta}\) and \(\boldsymbol{\epsilon}\) matrices for this case are as follows:

\(\mathbf{Y} = \begin{bmatrix}Y_{11}\\ Y_{12}\\ Y_{21}\\ Y_{22}\\ Y_{31}\\ Y_{32}\end{bmatrix}\text{, }\mathbf{X} = \begin{bmatrix}1 & 1 & 0\\ 1 & 1 & 0\\ 1 & 0 & 1\\ 1 & 0 & 1\\ 1 & -1 & -1\\ 1 & -1 & -1\end{bmatrix}\text{, }\boldsymbol{\beta} = \begin{bmatrix}\beta_{0}\\ \beta_{1}\\ \beta_{2}\end{bmatrix}\text{, and } \boldsymbol{\epsilon} = \begin{bmatrix}\epsilon_{11}\\ \epsilon_{12}\\ \epsilon_{21}\\ \epsilon_{22}\\ \epsilon_{31}\\ \epsilon_{32}\end{bmatrix}\)

where \(\beta_0, \beta_1 \text{ and } \beta_2\), corresponds to \(\mu, \tau_1 \text{ and } \tau_2\) respectively. 

Note that the vector of expected values \(\mathbf{E}\{\mathbf{Y}\} =\mathbf{X}\boldsymbol{\beta}\), yields the following:

\begin{aligned} \mathbf{E}\{\mathbf{Y}\} &= \mathbf{X}\boldsymbol{\beta} \\  \begin{bmatrix}E\{Y_{11}\}\\ E\{Y_{12}\}\\ E\{Y_{21}\}\\ E\{Y_{22}\}\\E\{Y_{31}\}\\E\{Y_{32}\}\end{bmatrix}  &= \begin{bmatrix}1 & 1 & 0\\ 1 & 1 & 0\\ 1 & 0 & 1\\ 1 & 0 & 1\\ 1 & -1 & -1\\ 1 & -1 & -1\end{bmatrix}\begin{bmatrix}\beta_{0}\\ \beta_{1}\\ \beta_{2}\end{bmatrix} \\ &= \begin{bmatrix}\mu_{.}+\tau_{1}\\ \mu+\tau_{1}\\ \mu+\tau_{2}\\ \mu+\tau_{2}\\ \mu-\tau_{1}-\tau_{2}\\ \mu-\tau_{1}-\tau_{2}\end{bmatrix} \end{aligned}

Since \(\tau_3=-\tau_{1}-\tau_{2}\), see above, we see that \(E\{Y_{31}\} =E\{Y_{32}\}=\mu+\tau_{3}\). Thus, the above \(\mathbf{X}\) matrix and \(\boldsymbol{\beta}\) vector representation provides the appropriate expected values for all factor levels as expressed below:

\(E\{Y_{ij}\} =\mu+\tau_{i}\)

In Lesson 1.1 we encountered a brief description of an experiment. The description of an experiment provides a context for understanding how to build an appropriate statistical model. All too often mistakes are made in statistical analyses because of a lack of understanding of the setting and procedures in which a designed experiment is conducted. Creating a study diagram is one of the best and intuitive ways to address this. A study diagram is a schematic diagram that captures the essential features of the experimental design. Here, as we explore the computations for a single factor ANOVA in a simple experimental setting, the study diagram may seem trivial. However, in practice and in lessons to follow, the ability to create an accurate study diagram usually makes a substantial difference in getting the model right.

In our example in Lesson 1.1, a plant biologist thinks that plant height may be affected by the fertilizer type and three types of fertilizer were chosen to investigate this claim.  Next, 24 plants were randomly chosen and 4 batches of 6 plants each were assigned to the 3 fertilizer types and as well as a control (untreated) group. The individual plants were randomly assigned the treatment levels to produce 6 replications within each level. The researchers kept all the plants under equivalent conditions in a single greenhouse throughout the experiment. Here is the resulting data:

We have a description of the treatment levels, how they were assigned to individual experimental units (the potted plant), and we can see the data organized in a table... What are we missing? A key question is: how was the experiment conducted? This question is a practical one and is answered with a study diagram. These are usually hand-drawn depictions of a real setting, indicating the treatments, levels of treatments, and how the experiment was laid out. For this example, we need to draw a greenhouse bench, capable of holding the 4 × 6 = 24 experimental units:

The diagram identified the response variable, listed the treatment levels, and indicated the random assignment of treatment levels to these 24 experimental units on the greenhouse bench.

This randomization and the subsequent experimental layout we would identify as a Completely Randomized Design (CRD). We know from this schematic diagram that we need a statistical model that is appropriate for a one-way ANOVA in a Completely Randomized Design (CRD).

Furthermore, once the plant heights are recorded at the end of the study, the experimenter may observe that the variability in the growth may be influenced by additional factors besides the fertilizer level (e.g. position on the bench). A careful examination of the layout of the plants in the study diagram may perhaps reveal this additional factor. For example, if the growth is higher in the plants placed on the row nearer to the windows, it is reasonable to assume that sunlight also plays a role. In this case, an experimenter may redesign the experiment as a randomized completely block design (RCBD) with rows as a blocking factor. Note that design aspects of experiments (such as CRD and RCBD) are covered in Lessons 7 and 8. 

Being able to draw and reproduce a study diagram is very useful in identifying the components of the ANOVA models.

This lesson covered four different versions of single-factor ANOVA models. They are: Overall Mean, Cell Means, Dummy Variable Regression, and Effects Coded Regression models. This lesson also provided the coding compatible with the SAS IML procedure, which facilitates the ANOVA computations using Matrix Algebra in a GLM setting. The method of least squares was used to estimate model parameters yielding a prediction equation for the response in terms of the treatment level. This prediction tool will show to be more useful in ANCOVA settings where model predictors are both categorical and numerical (more details on ANCOVA in Lessons 9  and 10). The prediction process can be utilized effectively only with a sound knowledge of the parameterization process for each ANOVA model. To facilitate this understanding, the models were run with IML code and the resulting parameter vector was used for interpreting the prediction (regression) equations.

Finally, using the greenhouse example, the concept of a study diagram was discussed. Though a simple visual tool, a study diagram may play an important role in identifying new predictors so that a pre-determined ANOVA model can be extended to include additional factors creating a multi-factor model (discussed in Lessons 5 and 6). In addition to identifying the treatment design, the study diagram also helps in choosing an appropriate randomization design (a topic discussed in Lessons 7 and 8).