Notes for Tuesday February 3, 1998

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

1 Multifactor Studies

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

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

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

1.1 Some examples

Example 1.

A steel company studied the effects of carbon content and tempering on the strength of steel. Carbon content was investigated at a high level and a low level. Tempering was also studied at a high temperature and a low temperature (the precise definitions of carbon and temperature levels are not important here). Thus there were 2 × 2 = 4 possible treatments:

Treatment  Description
1  High carbon, high temperature
2  High carbon, low temperature
3  Low carbon, high temperature
4  Low carbon, low temperature

The four treatments were then randomly assigned to twelve production batches and the strength of steel in each batch was measured. □

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

Laid out in a two-way table, the data might look like this:

Temperature

Carbon High  Low
High  60.35 −16.82
−45.63 −24.02
76.88 −32.09
Low −17.00 11.58
33.48 −96.78
51.44 63.11

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

Example 2.

An analyst studied the effects of family income (under $10,000, $10,000–$29,999, $30,000– $49,999, $50,000 and more) and stage in the life cycle of the family (1, 2 or 3) on appliance purchases. Thus 4 × 3 = 12 treatments are defined:

Treatment  Description
1  Under $10,000, Family stage 1
2  Under $10,000, Family stage 2
3  Under $10,000, Family stage 3
4  $10,000–$29,999, Family stage 1
⋮  ⋮
12  $50,000 or more, Family stage 3

The analyst selected 20 families whose income and family stage matched each “treatment” above, yielding 240 famlies for the study. □

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

Example 3.

A psychologist interested in human-computer interaction wants to know what makes audio menus (like the recorded messages “Press 1 for blah-blah, press 2 for blah-blah” in bank-by-phone and voice mail systems) easy or difficult for people to use. Before beginning the study, he measures subjects’ reading span (a measure of the capacity of their short-term memory) and categorizes them as having high, medium, or low reading spans. He then gives each subject a series of tasks to perform with a bank-by-phone system, so that each subject sees each combination of the following task features:

• A broad menu system that requires two selections from among nine options each to reach the goal in the menu; versus a deep menu system that requires four selections from among three options each to reach the goal;

• A task that is expressible in terms of numerical information, like “deposit $1000”, versus one that is expressible in verbal terms, such as “change mailing address”.

• A task that is first in a sequence of two tasks, and a task that is second in a sequence of two tasks.

The number of “treatments” in this study is 3 × 2 × 2 × 2 = 24. Within each reading span group, the psychologist selects four people for each of the eight combinations of task features, for a total of 96 subjects. He then assigns 10 tasks to each person and records the proportion of tasks correctly performed. □

This design is a four-factor, or four-way, design. Sometimes it would be referred to as a 3 × 2^{3} design. It is also balanced—the same number of subjects are assigned to each “treatment.”

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

Aside, about the definition of units and the definition of the outcome measure: From the psychologist’s point of view, this may not be the best design for understanding the influences of the various factors on tasks, since people, not tasks, are the fundamental units being studied.

• One alternative design would treat the tasks as the units, so that (perhaps) the outcome measure would be the proportion of persons getting each task correct.

• Another alternative design would recognize that the psychologist might want to learn how each human subject responds to tasks with differing menu depths, numerical or verbal descriptors, etc. This leads to a more complicated analysis, since now the same human subjects would get used in different treatment cells of the study, which introduces dependence between the ANOVA error terms ϵ_{ij} in different treatment cells.

1.2 Complete and Fractional Factorial Designs

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

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

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

1.3 Advantages of Multifactor (Factorial) Studies

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

• Information about interactions. At best the method of studying one factor at a time can give you information about the main effects of the factors, but it cannot give you information about the interactions.

In Example 1, suppose that the influence of carbon content on the strength of the steel is different, depending on whether the steel is tempered at a high or low temperature; this is the essence of an interaction. This cannot be discovered by a pair of studies, the first of which holds the carbon level fixed while it varies the temperature, and the second of which holds the temperature fixed while varying the carbon level.

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

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

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

2 Two-Way Layouts

2.1 Example: Balanced two-way layout

A consumer organization studied the effects of age and gender of the seller (car owner) on the the size of a cash offer for a used car. A single, medium-priced, six year old car was offered for sale to 36 car dealers selected at random from all dealers in the region. The 36 dealers were randomly “assigned” in groups of six, to each “owner.” The offers, in 100’s of dollars, were as follows:

"Owners'"         Dealers'
Age     Gender   Offers (x $100)
------- ------- -----------------
YOUNG   MALE    21 23 19 22 22 23
YOUNG   FEMALE  21 22 20 21 19 25
MIDDLE  MALE    30 29 26 28 27 27
MIDDLE  FEMALE  26 29 27 28 27 29
ELDERLY MALE    25 22 23 21 22 21
ELDERLY FEMALE  23 19 20 21 20 20

To move the data into Splus and format it in a way that is useful with aov() we do the following (there is more than one way to skin this cat, but…):

402 > raw.offers _read.table("car-offers.txt")
402 > raw.offers
       V1     V2 V3 V4 V5 V6 V7 V8
1   YOUNG   MALE 21 23 19 22 22 23
2   YOUNG FEMALE 21 22 20 21 19 25
3  MIDDLE   MALE 30 29 26 28 27 27
4  MIDDLE FEMALE 26 29 27 28 27 29
5 ELDERLY   MALE 25 22 23 21 22 21
6 ELDERLY FEMALE 23 19 20 21 20 20
402 > offer.vector _c(as.matrix(raw.offers[,3:8]))
402 > age _factor(raw.offers[,1])
402 > gender _factor(raw.offers[,2])
402 > car.offers _data.frame(age=rep(age,6),gender=rep(gender,6),
+ offer=offer.vector)
402 > rm(age,gender,offer.vector)
402 > car.offers
       age gender offer
 1   YOUNG   MALE    21
 2   YOUNG FEMALE    21
 3  MIDDLE   MALE    30
 4  MIDDLE FEMALE    26
 5 ELDERLY   MALE    25
 6 ELDERLY FEMALE    23
 7   YOUNG   MALE    23
 8   YOUNG FEMALE    22
 9  MIDDLE   MALE    29
10  MIDDLE FEMALE    29
11 ELDERLY   MALE    22
12 ELDERLY FEMALE    19
.
.
.
34  MIDDLE FEMALE    29
35 ELDERLY   MALE    21
36 ELDERLY FEMALE    20

We can get a general look at the data with plot.design and plot.factor:

402 > par(mfrow=c(2,2))
402 > plot.design(car.offers,fun=mean)
402 > plot.design(car.offers,fun=median)
402 > plot.factor(car.offers)

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

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

402 > attach(car.offers)
402 > interaction.plot(age,gender,offer)
402 > interaction.plot(gender,age,offer)
402 > detach()

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

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

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

2.1.1 The cell-means model

By itself the cell-means model is not very interesting, but it gives us some useful notation to use later. The cell-means model for two-way data is

y_{ijk} = μ_{ij} + ϵ_{ijk}

where i is the level of the first (row) factor, j is the level of the second (column) factor, and k indexes observation in treatment cell ij. μ_{ij} is the mean of the observations in treatment cell ij, and ϵ_{ijk} are independent N(0, σ^{2}) errors for each observation. It is useful to think of the data as laid out in rows and columns, like so

Age  Row
Gender  Y oung (j = 1) Middle (j = 2) Old (j = 3) Average
Male (i = 1)  μ_{11}  μ_{12}  μ_{13}  μ_{1.}
F emale (i = 2)  μ_{21}  μ_{22}  μ_{23}  μ_{2.}
Column Average  μ_{.1}  μ_{.2}  μ_{.3}  μ_{..}

The extra parameters in the table are defined as one might expect:

μ_{i.} = \frac{1}{J} ∑_{j=1}^{J} μ_{ij}

μ_{.j} = \frac{1}{I} ∑_{i=1}^{I} μ_{ij}

μ_{..} = \frac{1}{IJ} ∑_{i=1}^{I} ∑_{j=1}^{J} μ_{ij}

where I is the number of levels of the row factor and J is the number of levels of the column factor. The cell means are estimated in the natural way, as

μ̂_{ij} = Ȳ_{ij} = \frac{1}{n_{ij}} ∑_{k=1}^{n_{ij}} y_{ijk}

where the n_{ij} are the number of observations in each cell (when all the n_{ij} are equal, it is a balanced design). Note that the interaction plots above are just plots of the rows or columns of the Ȳ_{ij}’s in the two-way layout. Marginal totals μ_{i.}, μ_{.j} and μ_{..} would be estimated by the corresponding averages of barY ’s, i.e. Ȳ_{i.}, Ȳ_{.j} and Ȳ_{..}.

The fitted values are ŷ_{ij} = Ȳ_{ij}, and the residuals are r_{ijk} = y_{ijk} − ŷ_{ij}. We could partition the total sum of squares, as for one-way ANOVA, into between- and within- sums of squares,

∑_{i=1}^{I} ∑_{j=1}^{J} ∑_{k=1}^{n_{ij}}(Ȳ_{..} − y_{ijk})^{2} = ∑_{i=1}^{I} ∑_{j=1}^{J} n_{ij}(Ȳ_{..} − ŷ_{ij})^{2} + ∑_{i=1}^{I} ∑_{j=1}^{J} ∑_{k=1}^{n_{ij}}(ŷ_{ij} − y_{ijk})^{2}, or
SST = SST R + SSE

and this gives rise to the usual ANOVA table, as before.

Source SS  df  MS
Factors SSTR IJ − 1 MSTR=SST R/(IJ − 1)
Error SSE n − IJ MSE=SSE/(n − IJ)

Total SST n − 1  MST=SST/(n − 1)

which we can examine immediately in SPLUS

402 > car.aov _aov( offer ~ gender:age,data=car.offers)
402 > anova(car.aov)
Analysis of Variance Table
Response: offer
Terms added sequentially (first to last)
           Df Sum of Sq  Mean Sq  F Value         Pr(F)
gender:age  5  327.2222 65.44444 27.39535 2.451114e-010
 Residuals 30   71.6667  2.38889

Note the use of “:” to indicate the row and column variables (more generally “:” is used to construct interactions, but more about that later). Unsurprisingly (given the plots we did before), something about the factors matters, but this analysis does not reveal clearly what it is.

2.1.2 Main Effects and Interactions

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

• Main row (first factor) effects:

α_{i} = μ_{i.} − μ_{..};

• Main column (second factor) effects:

β_{j} = μ_{.j} − μ_{..};

It is easy to verify that ∑_{i} α_{i} = 0 and ∑_{j} β_{j} = 0, and it is easy to see how to estimate α_{i} and β_{i} in terms of the Ȳ ’s.

To see whether the first factor “matters” we want to construct a test of H_{0} : α_{i} ≡ 0 vs H_{1} : not all 0; and similarly to see whether the second factor “matters” we should test H_{0} : β_{j} ≡ 0 vs H_{1} : not all 0.

In one-way ANOVA we say that the cell means model and the grand means model were related by the equation μ_{i} = μ_{.} + α_{i}. The corresponding equation here,

μ_{ij} =^{?} μ_{..} + α_{i} + β_{j}
= μ_{i.} + μ_{.j} − μ_{..}

is not true in general. The differences between the left hand sides of this “equation” and the right hand sides are called the interactions,

• Interactions:

(αβ)_{ij} = μ_{ij} − (α_{i} + β_{i} + μ_{..})
= μ_{ij} − μ_{i.} − μ_{.j} + μ_{..}

It is also easy to see that ∑_{i} ∑_{j}(αβ)_{ij} = 0. Now we have equality:

μ_{ij} = μ + α_{i} + β_{j} + (αβ)_{ij}
= μ + (μ_{i.} − μ) + (μ_{.j} − μ) + (μ_{ij} − μ_{i.} − μ_{.j} + μ),

where μ = μ_{..} is the grand mean; and this gives rise to the “grand mean plus effects” model for two-way ANOVA, namely

y_{ijk} = μ + α_{i} + β_{j} + (αβ)_{ij} + ϵ_{ijk}

subject to the conditions that the α_{i}’s sum to zero, the β_{j}’s sum to zero, and the (αβ)_{ij}’s sum to zero.

In any two-way ANOVA, it turns out that the group of main effects for the first factor is orthogonal to the group of main effects for the second factor, and both of these are orthogonal to the interactions. This means we can partition SSR in the ANOVA table above into contributions from each set of main effects and the interactions:

∑_{i=1}^{I} ∑_{j=1}^{J} n_{ij}(Ȳ_{..} − ŷ_{ij})^{2} =

∑_{i=1}^{I} ∑_{j=1}^{J} n_{ij}(Ȳ_{..} − Ȳ_{ij})^{2} = ∑_{i=1}^{I} n_{i.}(Ȳ_{..} − Ȳ_{i.})^{2} + ∑_{j=1}^{J} n_{.j}(Ȳ_{..} − Ȳ_{.j})^{2} + n ∑_{i=1}^{I} ∑_{j=1}^{J}(Ȳ_{..} − Ȳ_{i.} − Ȳ_{.j} + Ȳ_{ij})^{2}

or

SSR = SS_{α} + SS_{β} + SS_{(αβ)}

and this allows us to expand the ANOVA table.

In the case of a balanced design, say n_{ij} ≡ m for all cells, the expanded ANOVA table looks like this:

Source  SS  df  MS
Factors  SSTR IJ − 1  MSTR=SST R/(IJ − 1)
Main effects (α) SS_{α}  I − 1  MS_{α}=SS_{α}/(I − 1)
Main effects (β) SS_{β}  J − 1  MS_{β}=SS_{β}/(J − 1)
Interactions  SS_{(αβ)} (I − 1)(J − 1) MS_{(αβ)}=SS_{(αβ)}/((I − 1)(J − 1))

Error  SSE  IJ(m − 1)  MSE=SSE/(IJ(m − 1))

Total  SST mIJ − 1  MST=SST/(mIJ − 1)

This is also easy to obtain in SPLUS:

402 > car.aov _aov( offer ~ gender*age,data=car.offers)
402 > anova(car.aov)
Analysis of Variance Table
Response: offer
Terms added sequentially (first to last)
           Df Sum of Sq  Mean Sq  F Value    Pr(F)
    gender  1    5.4444   5.4444  2.27907 0.141592
       age  2  316.7222 158.3611 66.29070 0.000000
gender:age  2    5.0556   2.5278  1.05814 0.359700
 Residuals 30   71.6667   2.3889

Note the use of the “*” to indicate that main effects gender, age as well as the interactions gender:age should be included in the model. Scanning the F -tests here, it seems clear that age matters and gender does not (since the gender and age effects are orthogonal it does not matter in which order we examine them, as can be seen from the alternate order for the factors below).

402 > car.aov _aov( offer ~ age*gender,data=car.offers)
402 > summary(car.aov)
           Df Sum of Sq  Mean Sq  F Value    Pr(F)
       age  2  316.7222 158.3611 66.29070 0.000000
    gender  1    5.4444   5.4444  2.27907 0.141592
age:gender  2    5.0556   2.5278  1.05814 0.359700
 Residuals 30   71.6667   2.3889

Moreover, the interaction doesn’t seem to matter. When we can safely set the interaction terms to zero in a two-way ANOVA, we arrive at what is called the additive model,

y_{ijk} = μ + α_{i} + β_{j} + ϵ_{ijk}

In the additive model it is true that μ_{ij} = μ + α_{i} + β_{j}, so differences in cell means in one factor are the same, no matter the value of the second factor:

μ_{ij} − μ_{i′j} = μ_{ij′} − μ_{i′j′}

This is why the curves in an interaction plot should be parallel when there is no interaction between the variables.

2.1.3 Model Selection and Diagnostics

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

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

402 > car.aov _aov( offer ~ gender*age,data=car.offers)
402 > par(mfrow=c(2,3)); plot(car.aov)
402 > par(mfrow=c(1,2))
402 > hist(residuals(car.aov)); plot(residuals(car.aov))

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

36-402: Advanced Data Analysis II Spring 1998 13

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

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

2.1.4 Fitted Values, Parameter Estimates

Fitted values are easy to obtain in SPLUS from the aov object, using fitted(), or by computing cell means by hand.

• If interactions are absent, it is usually most interesting to obtain estimates of the marginal means μ_{i.} and μ_{.j}.

• If interactions are present, then one usually must focus on the cell means μ_{ij} themselves.

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

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

402 > car.aov _aov(offer ~ C(gender,sum) + C(age,sum), data=car.offers)
402 > summary.lm(car.aov)
Call: aov(formula = offer ~ C(gender, sum) + C(age, sum), data = car.offers)
Residuals:
    Min     1Q  Median     3Q   Max
 -2.889 -1.028 -0.1111 0.9444 3.889
Coefficients:
                  Value Std. Error  t value Pr(>|t|)
   (Intercept)  23.5556   0.2581    91.2765   0.0000
C(gender, sum)  -0.3889   0.2581    -1.5069   0.1416
  C(age, sum)1  -2.1389   0.3650    -5.8606   0.0000
  C(age, sum)2   4.1944   0.3650    11.4928   0.0000
Residual standard error: 1.548 on 32 degrees of freedom
Multiple R-Squared: 0.8077
F-statistic: 44.79 on 3 and 32 degrees of freedom, the p-value is 1.466e-011
Correlation of Coefficients:
               (Intercept) C(gender, sum) C(age, sum)1
C(gender, sum)  0.0
  C(age, sum)1  0.0         0.0
  C(age, sum)2  0.0         0.0           -0.5

2.2 Example: One observation per cell

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

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

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

Effect of four treatments on penicillin yield
         Treatment
Block    A  B  C  D
Blend-1 89 88 97 94
Blend-2 84 77 92 79
Blend-3 81 87 87 85
Blend-4 87 92 89 84
Blend-5 79 81 80 88

Here is one way to convert this to a form SPLUS will like

402 > pen.raw _read.table("penicillin.txt",header=T)
402 > pen.raw _c(as.matrix(pen.raw))
402 > pen.raw
 [1] 89 84 81 87 79 88 77 87 92 81 97 92 87 89 80 94 79 85 84 88
402 > pendata _data.frame(blend=factor(rep(1:5,4)),
+ tx=rep(LETTERS[1:4],rep(5,4)), yield=pen.raw)
402 > pendata
   blend tx yield
 1     1  A    89
 2     2  A    84
 3     3  A    81
 4     4  A    87
 5     5  A    79
 .
 .
 .
20     5  D    88

We begin with the usual plots

402 > par(mfrow=c(2,2))
402 > plot.design(pendata)
402 > plot.design(pendata,fun=median)
402 > plot.factor(pendata)

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

To see this, recall the ANOVA table for balanced designs (n_{ij} ≡ m for all cells),

Source  SS  df  MS
Factors  SSTR IJ − 1  MSTR=SST R/(IJ − 1)
Main effects (α) SS_{α}  I − 1  MS_{α}=SS_{α}/(I − 1)
Main effects (β) SS_{β}  J − 1  MS_{β}=SS_{β}/(J − 1)
Interactions  SS_{(αβ)} (I − 1)(J − 1) MS_{(αβ)}=SS_{(αβ)}/((I − 1)(J − 1))

Error  SSE  IJ(m − 1)  MSE=SSE/(IJ(m − 1))

Total  SST mIJ − 1  MST=SST/(mIJ − 1)

we see that, when m = 1, if we include the interaction terms, there will be no degrees of freedom left for calculating the SSE and MSE. We can see this directly in SPLUS:

402 > pen.aov _aov(yield ~ blend*tx,data=pendata)
402 > summary(pen.aov)
         Df Sum of Sq  Mean Sq
   blend  4       264 66.00000
      tx  3        70 23.33333
blend:tx 12       226 18.83333

(and no error line is given in the ANOVA table). However, if we omit the interaction term, then we can use the interaction sum of squares as the error sum of squares:

402 > pen.aov _aov(yield ~ blend + tx, data=pendata)
402 > summary(pen.aov)
          Df Sum of Sq  Mean Sq  F Value     Pr(F)
    blend  4       264 66.00000 3.504425 0.0407462
       tx  3        70 23.33333 1.238938 0.3386581
Residuals 12       226 18.83333

Unfortunately, the treatment does not seem to be having much effect on the penicillin yield. On the other hand, we seem to have accidentally discovered that the blocking factor, the blend of the corn-steep liquor used, is marginally significant!

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

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

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

Finally, it might be interesting to compare these results with the results one would get if one ignored the blocking factor and did a one-way ANOVA on the four treatments:

402 > pen.aov _aov(yield ~ tx, data = pendata)
402 > summary(pen.aov)
          Df Sum of Sq  Mean Sq   F Value     Pr(F)
       tx  3        70 23.33333 0.7619048 0.5317827
Residuals 16       490 30.62500

Note that the p-value is even higher for tx, the treatment, than it was under the two-way, randomized blocks analysis. A more careful inspection of the two ANOVA tables explains why: the error SS in the one-way analysis is the sum of the error variance in the two-way analysis plus the SS due to the blocks: in other words, the variation due to differing blocks (differing corn-liquor blends) has been lumped into the error variance of the one-way analysis. This is why blocking is such an effective idea in experimental design.

2.3 Example: Unbalanced two-way layout

Human growth hormone was administered at a clinical research center to growth-hormone-deficient, short children who had not yet reached puberty. The investigator was interested in the effects of gender (factor A: male or female) and bone development (factor B: severely (1), moderately (2), or mildly (3), depressed) on the rate of growth induced by the hormone.

M 1 1.4
M 1 2.4
M 1 2.2
M 2 2.1
M 2 1.7
M 3 0.7
M 3 1.1
F 1 2.4
F 2 2.5
F 2 1.8
F 2 2.0
F 3 0.5
F 3 0.9
F 3 1.3

Reading the data in this format is straightforward, so we do this and then produce what we hope are useful initial plots.

402 > hormone.data _read.table("hormone.txt")
402 > names(hormone.data)
[1] "V1" "V2" "V3"
402 > names(hormone.data)  _c("gender","severity","growthrate")
402 > hormone.data$severity _factor(hormone.data$severity)
402 > par(mfrow=c(2,2))

402 > plot.design(hormone.data)
402 > plot.design(hormone.data, fun=median)
402 > plot.factor(hormone.data)
402 > attach(hormone.data)
402 > par(mfrow=c(1,2))
402 > interaction.plot(gender,severity,growthrate)
402 > interaction.plot(severity,gender,growthrate)
402 > detach()

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

Clearly the data are unbalanced; the cell sample sizes are

Cell  n_{ij}
M 1     3
M 2     2
M 3     2
F 1     1
F 2     3
F 3     3

In principle there is little difference between doing an ANOVA with unbalanced designs and doing ANOVA with balanced designs. Calculation of sums of squares, and presentation of the ANOVA table, is more difficult, because the sums of squares in the formulas

SST =^{?} SST R + SSE, and
SST R =^{?} SS_{α} + SS_{β} + SS_{(αβ)}

are no longer orthogonoal—these equalities no longer hold. Instead, sequential sums of squares as in regression must be calculated. In fact, in many computer packages, the simplest way to get correct parameter estimates and sums of squares for unbalanced designs is to analyze the data using dummy variables in a regression package.

These problems can be illustrated in SPLUS without resorting to regression, because aov() does properly handle unbalanced ANOVA problems.

402 > hormone.aov _aov(growthrate ~ gender*severity,data=hormone.data)
402 > summary(hormone.aov)
                Df Sum of Sq  Mean Sq  F Value     Pr(F)
         gender  1  0.002857 0.002857  0.01758 0.8977853
       severity  2  4.396000 2.198000 13.52615 0.0027133
gender:severity  2  0.075429 0.037714  0.23209 0.7980337
      Residuals  8  1.300000 0.162500

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

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

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

Defining the grand mean and the marginal means in an unbalanced two-way layout can be done several different ways, just as for unbalanced one-way problems. The two main choices are

• Unweighted means. The definitions of μ and of μ_{i.} and μ_{.j} are exactly as for the balanced case.

• Sample size weighting. The definitions are changed to be

μ_{i.} = \frac{1}{n_{i.}} ∑_{j=1}^{J} n_{ij}μ_{ij}

μ_{.j} = \frac{1}{n_{.j}} ∑_{i=1}^{I} n_{ij}μ_{ij}
μ_{..} = \frac{1}{n_{..}} ∑_{i} ∑_{j} n_{ij}μ_{ij}

and corresponding changes in the “sum to zero” constraints on the α’s, β’s, etc.

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

402 > dummy.coef(hormone.aov)
$"(Intercept)":
 (Intercept)
         1.7
$severity:

   1   2    3
 0.5 0.3 -0.8
$gender:
   F    M
 0.1 -0.1
$"severity:gender":
  1F             2F   3F   1M            2M  3M
 0.1 -2.081668e-017 -0.1 -0.1 2.081668e-017 0.1

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

402 > model.tables(hormone.aov,se=T)
Refitting model to allow projection
Standard error information not returned as design is unbalanced.
Standard errors can be obtained through se.contrast.
Tables of effects
 severity
         1      2       3
    0.4571 0.3771 -0.7429
rep 4.0000 5.0000  5.0000
 gender
          F        M
    0.07714 -0.07714
rep 7.00000  7.00000
 severity:gender
Dim 1 : severity
Dim 2 : gender
         F      M
  1  0.171 -0.057
rep  1.000  3.000
  2  0.011 -0.017
rep  3.000  2.000
  3 -0.069  0.103
rep  3.000  2.000

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

402 > se.contrast(hormone.aov,contrast=list(gender=="M",severity==1,
+ gender=="M",severity==3),data=hormone.data)
Refitting model to allow projection
[1] 0.1703463

By default, the standard error of the difference of the means specified by contrast is computed. The contrast itself can be computed using model.tables() again:

402 > model.tables(hormone.aov,type="means")
Refitting model to allow projection
Tables of means
Grand mean
 1.643
 severity
      1    2   3
    2.1 2.02 0.9
rep 4.0 5.00 5.0
 gender
       F     M
    1.72 1.566
rep 7.00 7.000
 severity:gender
Dim 1 : severity
Dim 2 : gender
      F   M
  1 2.4 2.0
rep 1.0 3.0
  2 2.1 1.9
rep 3.0 2.0
  3 0.9 0.9
rep 3.0 2.0

Thus, an approximate 95% confidence interval for the difference in growthrate for the first and third bone severity levels among boys is (2.0 − 0.9) ± qt(0.975, 3) ∗ (0.1703463), which runs from 0.557882 to 1.642118 (qt(0.975,3) = 3.182446).

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