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.
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)
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) )
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.
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
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.