Factorial ANOVA (without interaction)

Unlike the one-way (or one factor) ANOVA we have worked with thus far, factorial ANOVA refers to designs with 2 or more factors. For example, consider the design below:

Two researchers setup a field experiment in the shallows of a small lake on Vancouver Island (British Columbia) to measure how fish abundance affects the diversity and abundance of prey species. To minimize noise from background variation in prey availability at different locations in the lake, they randomly selected five locations, and within each location they setup three fish abundance treatments – low, high, and control.

In the low treatment, 30 small fish were added to a mesh cage, while in the high treatment 90 fish were added to a second cage nearby. An adjacent space was left as is to serve as the control treatment. After 13 days diversity of zooplankton prey was measured using an index (Levins’ D) that takes both the number of species and their rarity into account. Does treatment affect zooplankton diversity?

\(H_0\): Treatment does not impact zooplankton diversity
\(H_A\): Treatment does impact zooplankton diversity

\(H_0\): Location does not impact zooplankton diversity
\(H_A\): Location does impact zooplankton diversity

zoo = read.csv("http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter18/chap18e2ZooplanktonDepredation.csv")
names(zoo); summary(zoo)
## [1] "treatment" "diversity" "block"
##    treatment   diversity         block  
##  control:5   Min.   :1.000   Min.   :1  
##  high   :5   1st Qu.:1.400   1st Qu.:2  
##  low    :5   Median :2.200   Median :3  
##              Mean   :2.133   Mean   :3  
##              3rd Qu.:2.550   3rd Qu.:4  
##              Max.   :4.100   Max.   :5

Location is recorded as block so let us convert it into a factor. We should also make control the reference group, followed by low, and ending with high.

zoo$location = factor(zoo$block, levels = c(1:5), labels = c("Location 1", "Location 2", "Location 3", "Location 4", "Location 5"))
zoo$groups = factor(zoo$treatment, levels = c("control", "low", "high"))
zoo$groups = relevel(zoo$groups, ref = "control")

Now we fit the two-way ANOVA model.

lm.1 = lm(diversity ~ groups + location, data = zoo)
anova(lm.1)
## Analysis of Variance Table
## 
## Response: diversity
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## groups     2 6.8573  3.4287 16.3660 0.001488 **
## location   4 2.3400  0.5850  2.7924 0.101031   
## Residuals  8 1.6760  0.2095                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Evidently Location has no impact on diversity but treatment does.

Pairwise Comparisons

Once again, we could look at pairwise comparisons and see which treatments really differ from one another. Before we do so, however, let us take a look at some options here. (1) We could stick with Tukey’s HSD if we wanted all unplanned comparisons, or opt for (2) Dunnett Contrasts if we are only interested in comparisons with the control group as the reference group. Both are shown below:

library(multcomp)
uc1 = glht(lm.1, linfct = mcp(groups = "Tukey"))
summary(uc1)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = diversity ~ groups + location, data = zoo)
## 
## Linear Hypotheses:
##                     Estimate Std. Error t value Pr(>|t|)   
## low - control == 0   -1.0200     0.2895  -3.524  0.01872 * 
## high - control == 0  -1.6400     0.2895  -5.665  0.00104 **
## high - low == 0      -0.6200     0.2895  -2.142  0.14231   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
uc2 = glht(lm.1, linfct = mcp(groups = "Dunnett"))
summary(uc2)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Dunnett Contrasts
## 
## 
## Fit: lm(formula = diversity ~ groups + location, data = zoo)
## 
## Linear Hypotheses:
##                     Estimate Std. Error t value Pr(>|t|)    
## low - control == 0   -1.0200     0.2895  -3.524 0.014135 *  
## high - control == 0  -1.6400     0.2895  -5.665 0.000876 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Factorial ANOVA (with interaction)

Here we have a third null hypothesis – that of no interaction between the two factors. Let us set this up with the inter-tidal algae example.

Harley (2003) investigated how herbivores affect the abundance of plants living in the intertidal habitat of coastal Washngton using field transplants of a red alga, Mazaella parksii. The experiment also investigated whether the effects of herbivores on the algae depended upon where in the intertidal zone the plants were growing. Thirty-two study plots were establshed just above the low-tide mark, and another 32 plots were set at the mid-point between the low- and high-tide marks. Plots were cleared of all organisms, and a constant amount of new algae was glued onto the rock surface of each plot. Using copper fencing, herbivores (mainly limpets and snails) were excluded from a randomly chosen half of the plots at each height. The remaining plots were left accessible to herbivores. The design included every combination of the height and herbivore treatments. At the end of the experiment the surface area covered by algae (in \(cm^2\)) was measured in each plot. The data were square-root transformed to meet the assumptions of normality and equal variances. The data are given below:

algae = read.csv("http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter18/chap18e3IntertidalAlgae.csv")
names(algae); summary(algae)
## [1] "height"     "herbivores" "sqrtArea"
##  height   herbivores    sqrtArea      
##  low:32   minus:32   Min.   : 0.7071  
##  mid:32   plus :32   1st Qu.: 5.9204  
##                      Median :22.2549  
##                      Mean   :22.8382  
##                      3rd Qu.:37.2306  
##                      Max.   :61.0960

\(H_0\): Herbivory treatments have no impact on algae cover
\(H_A\): Herbivory treatments have an impact on algae cover

\(H_0\): Height has no impact on algae cover
\(H_A\): Height has an impact on algae cover

\(H_0\): There is no interaction between Herbivory treatments and Height
\(H_A\): There is an interaction between Herbivory treatments and Height

lm.2 = lm(sqrtArea ~ herbivores * height, data = algae)
anova(lm.2)
## Analysis of Variance Table
## 
## Response: sqrtArea
##                   Df  Sum Sq Mean Sq F value   Pr(>F)   
## herbivores         1  1512.2 1512.18  6.3579 0.014360 * 
## height             1    89.0   88.97  0.3741 0.543096   
## herbivores:height  1  2617.0 2616.96 11.0029 0.001549 **
## Residuals         60 14270.5  237.84                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
with(algae, interaction.plot(herbivores, height, sqrtArea) )

with(algae, interaction.plot(height, herbivores, sqrtArea) )

Adjusting for a Covariate via ANCOVA

In observational studies one invariably has to control for the effects of a confounding variable – variables that, unless included in the model, could distort the true effect of the treatment of interest. For example, a group of researchers wanted to assess physiological differences between two types of worker model-rats – frequent workers and infrequent workers. They compared daily energy expenditures of wild mole-rats during a dry season, and because energy expenditure varies with body mass but frequent workers are heavier than infrequent workers they had to include a measure not only of body mass but also which worker caste(frequent/infrequent) a mole-rat belonged to.

Trying to control for potential confounding variables is a typical concern in observational studies. In the language of ANOVA, analyses that involve these designs are referred to as ANCOVA models.

In these cases, one usually proceeds by testing for an interaction between the variables. If the interaction is not significant, that term is dropped from the model. In the case of these mole rates, the model looks like the following:

\[lnEnergy = lnMass + caste + lnMass \times caste\]

If the interaction term, \(lnMass \times caste\) is not significant and is dropped the model reduces to:

\[lnEnergy = lnMass + caste\]

Let us see how the modeling proceeds.

rats = read.csv("http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter18/chap18e4MoleRatLayabouts.csv")
summary(rats)
##     caste        lnMass         lnEnergy    
##  lazy  :14   Min.   :3.850   Min.   :3.555  
##  worker:21   1st Qu.:4.248   1st Qu.:3.902  
##              Median :4.511   Median :4.190  
##              Mean   :4.540   Mean   :4.193  
##              3rd Qu.:4.844   3rd Qu.:4.489  
##              Max.   :5.263   Max.   :5.043

We begin with a model that includes the interaction term.

lm3 = lm(lnEnergy ~ lnMass + caste + lnMass:caste, data = rats)
anova(lm3)
## Analysis of Variance Table
## 
## Response: lnEnergy
##              Df  Sum Sq Mean Sq F value    Pr(>F)    
## lnMass        1 1.31061 1.31061 14.9100 0.0005362 ***
## caste         1 0.63747 0.63747  7.2521 0.0113217 *  
## lnMass:caste  1 0.08956 0.08956  1.0188 0.3206094    
## Residuals    31 2.72494 0.08790                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The results indicated that there was no interaction between body mass and caste, and hence the reduced model can then be estimated as:

lm4 = lm(lnEnergy ~ lnMass + caste, data = rats)
anova(lm4)
## Analysis of Variance Table
## 
## Response: lnEnergy
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## lnMass     1 1.31061 1.31061 14.9013 0.0005178 ***
## caste      1 0.63747 0.63747  7.2478 0.0111984 *  
## Residuals 32 2.81450 0.08795                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Clearly both body mass and caste matter for daily energy expenditure. Caste is of primary interest here so normally that would be the variable one would focus on.

Random-effects ANOVA

Repeatability ranges from 0 to 1 and expresses the proportion of variation in some measurement that is due to differences among the units not due to differences within the units. One way to measure this is via a ratio, called repeatability, where repeatability is measured as:

\[r = \dfrac{s^2_A}{s^2 + s^2_A}\]

where \(s^{2}_{A}\) is the variance among the units, and \(s^2\) is the variance within the units. These are, in turn, calculated as follows:

\[s^2 = {MS}_{within}\]
\[s^{2}_{A} = \dfrac{MS_{among} - MS_{within}}{n}\]
\[n = \text{number of measurements per unit}\]

If \(r\) tends to 1 then we have high repeatability; most of the variation is due to the variance among the units. If \(r\) tends to 0, then repeatability is low since most of the variation in the measurements seems to be coming from within the units themselves.

The walking stick Timema cristinae is a wingless herbivorous insect that lives on plants in chaparral habitats in California. In a study of the insect’s adaptations to different plant species, researchers measured a variety of traits using digital photographs of specimens collected from a study site. They used a computer to measure various traits in the photographs. Because they were concerned about measurement error they took two separate photographs of each specimen after some gap in time. Very often the second set of measurements were different from the first set, indicating measurement error. How large was the measurement error compared with the real variation among individuals in the trait?

Let us load the data to begin with.

sticks = read.csv("http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter15/chap15e6WalkingStickFemurs.csv")
summary(sticks)
##     specimen   femurLength    
##  Min.   : 1   Min.   :0.1400  
##  1st Qu.: 7   1st Qu.:0.2000  
##  Median :13   Median :0.2300  
##  Mean   :13   Mean   :0.2244  
##  3rd Qu.:19   3rd Qu.:0.2475  
##  Max.   :25   Max.   :0.3100

Now we fit the random-effects ANOVA.

library(nlme)
lm5 = lme(fixed = femurLength ~ 1, random = ~ 1|specimen, data = sticks)
anova(lm5)
##             numDF denDF  F-value p-value
## (Intercept)     1    25 1021.889  <.0001

The variance components can be recovered via the VarCorr() command. You will see two variances reported – one for the (Intercept) and the the other for the Residual.

lm5.Varcomp = VarCorr(lm5)
lm5.Varcomp
## specimen = pdLogChol(1) 
##             Variance     StdDev    
## (Intercept) 0.0010539182 0.03246411
## Residual    0.0003559996 0.01886795
  • Variance among the insects is the (Intercept) and = 0.001053
  • Variance within the insects is the Residual and = 0.000355

Repeatability can then be calculated very easily by extracting, from lm5.Varcomp, the variance in the first row, first column, and then the variance in the second row, first column.

varAmong  = as.numeric( lm5.Varcomp[1,1] )
varWithin = as.numeric( lm5.Varcomp[2,1] )
repeatability = varAmong / (varAmong + varWithin)
repeatability
## [1] 0.7475033

For the walking sticks, then, repeatability is about 0.7475, and so not bad at all. In plain words, this repeatability measure indicates that almost 75% of the variation in femur length is coming from differences between the specimens rather than from measurement error per se.