BIOL 5473 Tutorial 6
Starting off
Remember to set your working directory before you start. (See Week 1’s tutorial if you forget how.) You also need to load today’s packages.
library(afex)
library(effectsize)
library(emmeans)
library(GAD)
library(ggplot2)
library(ggthemes)
library(MASS)
library(multcomp)
library(patchwork)
library(predictmeans)Split-unit designs
A split-unit design is an experimental design in which the large units used for one treatment factor are split into smaller units for levels of another treatment factor. Because they originated in agriculture, split-unit designs are commonly known as “split-plot” designs.
As an example of a modern ecological split-plot design, we will be examining data from Kluane National Park in Yukon, Canada. The Kluane Project experimentally investigated the effects of fertilizing and herbivory on vegetation dynamics in the boreal forest ecosystem of Kluane National Park (Krebs, C.J., Boutin, S. & Boonstra, R., eds (2001a) Ecosystem dynamics of the Boreal Forest. The Kluane Project. Oxford University Press, New York).
The data are from a study of the effects of plant resources and herbivory on the defensive chemistry of understory plant species. Each of 16 5×5m plots was randomly assigned one of four treatments:
- Surrounded by a fence exclosure to exclude herbivores
- Fertilized with N-P-K fertilizer
- Fenced and fertilized
- Untreated control
Thus, at the whole-plot level, this is a fully crossed design between two factors: fenced and fertilized.
Each of the 16 plots was divided in two. One (randomly chosen) side of each plot received the treatment continually over the 20 years of study. The other half of each plot received the treatment for the first ten years, after which it was left to revert to its untreated state. Thus, there are 32 subplots in total. Each subplot receives a combined treatment of fenced, fertilized, and permanence.
The response variable is the concentration of phenolics (a crude measure of plant defense compounds) in yarrow (Achillea millefolium), a herb common in the plots. The measurement units are mg tannic acid equivalents per gram of dry mass. This description of the dataset and the dataset itself (stored on D2L as kluane.csv) are modified from: https://www.zoology.ubc.ca/~bio501/R/workshops/lme.html
Let’s read in the data. Because levels of our factors are all uninformatively coded (i.e., 0s and 1s), let’s change the labels to give them informative names we can understand easily. We will also check for balance.
d <- read.csv("Data/kluane.csv")
d$plot <- factor(d$plot)
d$fenced <- factor(d$fenced,
levels = c(0, 1),
labels = c("Unfenced", "Fenced"))
d$fertilized <- factor(d$fertilized,
levels = c(0, 1),
labels = c("Unfertilized", "Fertilized"))
d$permanent <- factor(d$permanent,
levels = c(0, 1),
labels = c("Non-permanent", "Permanent"))
replications(phen.ach ~ fenced * fertilized * permanent + Error(plot), data = d)## fenced fertilized
## 16 16
## permanent fenced:fertilized
## 16 8
## fenced:permanent fertilized:permanent
## 8 8
## fenced:fertilized:permanent
## 4The design is balanced.
Let’s look at the data.
ggplot(data = d, aes(x = fenced, y = phen.ach, col = permanent)) +
facet_wrap(~fertilized) +
geom_point(alpha = 0.2, position = position_dodge(width = 0.5)) +
geom_pointrange(stat = "summary", fun.data = "mean_se",
position = position_dodge(width = 0.5)) +
scale_colour_colorblind(name=NULL) +
labs(x = NULL, y = "Tannic acid equiv. (mg) / dry mass (g)") +
theme_few()
Immediately, I notice that:
- Fertilized plots seem to have much lower phenolics on average than unfertilized plots.
- Fenced plots seem to have slightly higher phenolics on average than unfenced plots.
- In unfertilized plots, permanent vs. non-permanent seems to have no effect. But in fertilized plots, permanent vs. non-permanent seems to have a major effect, with permanent plots having much lower phenolics on average. This suggests an interaction.
Analysis in GAD
Define variables as fixed/random, fit the model with
lm(), and summarize with gad().
fen <- as.fixed(d$fenced)
fer <- as.fixed(d$fertilized)
plo <- as.random(d$plot)
per <- as.fixed(d$perm)
kluane.lm <- lm(phen.ach ~ fen * fer * per + # Main effects & interactions
(plo %in% fen:fer) + # Random plot effects
(plo %in% fen:fer):per, data = d) # Plot:permanence interaction
(kluane.tab <- gad(kluane.lm, quasi.f = TRUE))## $f.versus
## Numerator Denominator Num.Df Den.Df
## fen fen fen:fer:plo 1 12
## fer fer fen:fer:plo 1 12
## per per fen:fer:per:plo 1 12
## fen:fer fen:fer fen:fer:plo 1 12
## fen:per fen:per fen:fer:per:plo 1 12
## fer:per fer:per fen:fer:per:plo 1 12
## fen:fer:per fen:fer:per fen:fer:per:plo 1 12
## fen:fer:plo fen:fer:plo Residuals 12 0
## fen:fer:per:plo fen:fer:per:plo Residuals 12 0
##
## $anova
## Analysis of Variance Table: Quasi F-ratios
##
## All of your terms has an appropriated F-ratio, Quasi F-test is not necessary
##
## Response: phen.ach
## Df Sum Sq Mean Sq F value Pr(>F)
## fen 1 354.31 354.31 5.4865 0.037228 *
## fer 1 3084.27 3084.27 47.7598 1.626e-05 ***
## per 1 477.87 477.87 15.2815 0.002075 **
## fen:fer 1 2.78 2.78 0.0431 0.838977
## fen:per 1 29.53 29.53 0.9443 0.350353
## fer:per 1 375.24 375.24 11.9997 0.004682 **
## fen:fer:per 1 56.66 56.66 1.8118 0.203167
## fen:fer:plo 12 774.95 64.58 NaN NaN
## fen:fer:per:plo 12 375.25 31.27 NaN NaN
## Residuals 0 0.00 NaN
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Because we have only one subplot for each combination of
fen:fer:perm in each plot, subplots are confounded with the
plot:permanence interaction. Therefore, the model is a perfect fit, and
we have no residual degrees of freedom.
To get the table of expected mean squares and the variance
components, we can use GAD’s comp.var() function.1
comp.var(kluane.lm, anova.tab = kluane.tab)## $mse
## Mean square estimates
## fen "1*Residuals + 2*fen:fer:plo + 16*fen"
## fer "1*Residuals + 2*fen:fer:plo + 16*fer"
## per "1*Residuals + 1*fen:fer:per:plo + 16*per"
## fen:fer "1*Residuals + 2*fen:fer:plo + 8*fen:fer"
## fen:per "1*Residuals + 1*fen:fer:per:plo + 8*fen:per"
## fer:per "1*Residuals + 1*fen:fer:per:plo + 8*fer:per"
## fen:fer:per "1*Residuals + 1*fen:fer:per:plo + 4*fen:fer:per"
## fen:fer:plo "1*Residuals + 2*fen:fer:plo"
## fen:fer:per:plo "1*Residuals + 1*fen:fer:per:plo"
## Residuals "1*Residuals"
##
## $comp.var
## Type Estimate Square root Percentage
## fen Fixed 18.108 4.255 13.153
## fer Fixed 188.73 13.738 42.463
## per Fixed 27.912 5.283 16.330
## fen:fer Fixed -7.724 0.000
## fen:per Fixed -0.218 0.000
## fer:per Fixed 42.996 6.557 20.268
## fen:fer:per Fixed 6.347 2.519 7.787
## fen:fer:plo Random 0.000
## fen:fer:per:plo Random NaN 0.000
## Residuals Random NaN 0.000
##
## $note
## [1] "There is one or more negative estimates of components of variation in your model. Perhaps you should consider pooling these terms"Refitting in aov()
Because the model is a perfect fit, all the residuals are 0. Therefore, if we try to conduct any downstream analyses, like producing interaction plots, we will get an error message. Therefore, we will refit the model by removing the plot:permanence interaction.
To refit the model in aov(), we need to specify an
Error() term. Here,
Error(fenced:fertilized:plot) fits a random effect for
plot.
op <- options(contrasts = c("contr.sum", "contr.poly")) # Set sum-to-zero contrasts
kluane.aov <- aov(phen.ach ~ fenced * fertilized * permanent +
Error(fenced:fertilized:plot), data = d) # Fit the model
summary(kluane.aov) # Show the ANOVA table##
## Error: fenced:fertilized:plot
## Df Sum Sq Mean Sq F value Pr(>F)
## fenced 1 354.3 354.3 5.487 0.0372 *
## fertilized 1 3084.3 3084.3 47.760 1.63e-05 ***
## fenced:fertilized 1 2.8 2.8 0.043 0.8390
## Residuals 12 774.9 64.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## permanent 1 477.9 477.9 15.281 0.00208 **
## fenced:permanent 1 29.5 29.5 0.944 0.35035
## fertilized:permanent 1 375.2 375.2 12.000 0.00468 **
## fenced:fertilized:permanent 1 56.7 56.7 1.812 0.20317
## Residuals 12 375.3 31.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The results are the same as produced using GAD above, indicating that
we have specified our Error() term correctly.
To check the residuals, we will use the residplot()
function from the predictmeans package. It will show us the
residuals at both the plot and subplot level:
residplot(kluane.aov)
Residuals at the plot level are shown in the Q-Q plot labelled “Normal Plot for Random Intercept”. These look pretty good, although at least one plot appears to be an outlier. This may be worth further investigation. All other plots look acceptable.
Therefore, we can say that fencing has a significant effect (\(F_{1,12}=5.486, p=0.037\)), and there is a significant interaction between fertilizing and permanence (\(F_{1,12}=12.000, p=0.005\)).
Only functions in the GAD package recognize random effects declared
using as random(). All other functions in R think that
all terms in the GAD model are fixed effects. Therefore, for
downstream analysis (e.g., Tukey’s method) of GAD models containing
random effects, the easiest option is to refit the model using
aov() or lmer().
We do not need to conduct Tukey’s method on fencing because there are only two levels and we already know that they differ significantly. However, we can calculate the estimated marginal means.
emmeans(kluane.aov, ~fenced)## fenced emmean SE df lower.CL upper.CL
## Unfenced 32.8 2.01 12 28.5 37.2
## Fenced 39.5 2.01 12 35.1 43.9
##
## Results are averaged over the levels of: fertilized, permanent
## Warning: EMMs are biased unless design is perfectly balanced
## Confidence level used: 0.95Let’s get the estimated marginal means for each fertilized:permanent combination. We will also calculate Tukey-adjusted pairwise contrasts.
emm <- emmeans(kluane.aov, ~fertilized:permanent)
pairs(emm, adjust = "tukey")## contrast estimate SE df
## (Unfertilized Non-permanent) - (Fertilized Non-permanent) 12.79 3.46 21.4
## (Unfertilized Non-permanent) - Unfertilized Permanent 0.88 2.80 12.0
## (Unfertilized Non-permanent) - Fertilized Permanent 27.36 3.46 21.4
## (Fertilized Non-permanent) - Unfertilized Permanent -11.91 3.46 21.4
## (Fertilized Non-permanent) - Fertilized Permanent 14.58 2.80 12.0
## Unfertilized Permanent - Fertilized Permanent 26.48 3.46 21.4
## t.ratio p.value
## 3.694 0.0067
## 0.315 0.9887
## 7.905 <.0001
## -3.440 0.0119
## 5.214 0.0011
## 7.651 <.0001
##
## Results are averaged over the levels of: fenced
## P value adjustment: tukey method for comparing a family of 4 estimatesNow let’s see the “compact letter display” to show which combinations differ significantly.
(tukey_cld <- cld(emm))## fertilized permanent emmean SE df lower.CL upper.CL .group
## Fertilized Permanent 19.1 2.45 21.4 14.0 24.1 1
## Fertilized Non-permanent 33.6 2.45 21.4 28.6 38.7 2
## Unfertilized Permanent 45.5 2.45 21.4 40.5 50.6 3
## Unfertilized Non-permanent 46.4 2.45 21.4 41.3 51.5 3
##
## Results are averaged over the levels of: fenced
## Warning: EMMs are biased unless design is perfectly balanced
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.To visualize the effects, let’s create two interaction plots
side-by-side. On the left, we will show the significant effect of
fencing, and on the right, the significant interaction between fencing
and permanence. We will force them to share a y-axis and
combine them into a single figure using the patchwork
package.
# Plot fencing
plot1 <- emmip(kluane.aov, ~fenced, CI = TRUE, col = "skyblue") +
labs(x = "Fencing", y = "Tannic acid equiv. (mg) / dry mass (g)", colour = NULL) +
ylim(c(10, 55)) +
theme_few()
# Plot fertilizing:permanence interaction
plot2 <- emmip(kluane.aov, permanent ~ fertilized, CI = TRUE) +
labs(x = "Fertilizing", y = NULL, colour = NULL) +
scale_colour_colorblind() +
ylim(c(10, 55)) +
theme_few() +
theme(legend.position = "top")
# Combine plots
plot1 + plot2
To add “a”, “b”, “ab” (etc.) indicators to the right-hand plot (corresponding to “1”, “2”, “12” (etc.) in the “compact letter display”), the easiest way is to edit the figure in Adobe Illustrator, Inkscape, etc.
Reset contrasts.
options(op)Refitting in lmer()
If you wanted to refit the model using REML in lmer(),
here is how you would do so. (1|fenced:fertilized:plot)
fits a random effect for plot. Because the design is balanced, it
produces the same results as aov() above.
library(car) # Load packages
library(lme4)
op <- options(contrasts = c("contr.sum", "contr.poly")) # Set sum-to-zero contrasts
kluane.lmer <- lmer(phen.ach ~ fenced * fertilized * permanent +
(1|fenced:fertilized:plot), data = d) # Fit model
(lmer.tab <- Anova(kluane.lmer, test = "F", type = 3)) # Show ANOVA table
omega_squared(lmer.tab) # Show effect sizes
S(kluane.lmer) # Show variance components
confint(kluane.lmer, oldNames = FALSE) # Show confidence intervals
residplot(kluane.lmer) # Check residuals
emmeans(kluane.lmer, ~fenced) # Show EMM for fencing
(emm <- emmeans(kluane.lmer, ~fertilized:permanent)) # Show EMM for interaction
pairs(emm, adjust = "tukey") # Tukey's method for interaction
(tukey_cld <- cld(emm)) # Show CLD
# Plot fencing
plot1 <- emmip(kluane.lmer, ~fenced, CI = TRUE, col = "skyblue") +
labs(x = "Fencing", y = "Tannic acid equiv. (mg) / dry mass (g)", colour = NULL) +
ylim(c(10, 55)) +
theme_few()
# Plot fertilizing:permanence interaction
plot2 <- emmip(kluane.lmer, permanent ~ fertilized, CI = TRUE) +
labs(x = "Fertilizing", y = NULL, colour = NULL) +
scale_colour_colorblind() +
ylim(c(10, 55)) +
theme_few() +
theme(legend.position = "top")
# Combine plots
plot1 + plot2
options(op) # Reset contrastsRepeated-measures designs
Although it is possible to conduct a repeated-measures ANOVA like a standard split-plot ANOVA as above, this is generally not recommended because the assumption of sphericity is likely to be violated in repeated-measures experiments. Therefore, corrections are often needed.
In afex
The afex package contains tests and corrections for
repeated measures. Note that the terminology used in the
afex package comes from Psychology; the word
subject is used to refer to an experimental unit. Levels of a
factor may vary between subjects (e.g. a different drug given
to each patient) or within subjects (e.g., the times at which
measurements are taken for each subject.)
As a demonstration of repeated measures, we will use a dataset containing information about the breathing rates of cane toads. In this dataset, each toad was initially either a buccal breather or a lung breather, so breathing type is a between-subjects factor. O2 concentration (in %) was varied over time for each toad. Therefore, O2 level is a within-subjects factor. The response variable was the rate of buccal breathing, which was square-root transformed to reduce positive skewness.
Read in the data, convert to factors, and check for balance.
d <- read.csv("Data/mullens.csv")
d$TOAD <- as.factor(d$TOAD)
d$BRTH.TYP <- as.factor(d$BRTH.TYP)
d$O2LEVEL <- as.ordered(d$O2LEVEL)
replications(SFREQBUC ~ O2LEVEL * BRTH.TYP + Error(TOAD), data = d)## $O2LEVEL
## [1] 21
##
## $BRTH.TYP
## BRTH.TYP
## buccal lung
## 104 64
##
## $`O2LEVEL:BRTH.TYP`
## BRTH.TYP
## O2LEVEL buccal lung
## 0 13 8
## 5 13 8
## 10 13 8
## 15 13 8
## 20 13 8
## 30 13 8
## 40 13 8
## 50 13 8The experiment is unbalanced. We have 21 replicates for each level of
O2: 13 buccal breathers and 8 lung breathers. Therefore, we
need to make sure that we use Type II or Type III sums of squares.
Functions in the afex package will do this for us.
The format for the aov_ez() and afex_plot()
functions is quite different than for aov() or
gad():
idspecifies the subject ID. What do you have repeated measures on? In our case, we have repeated measures on each toad, and their IDs are given in theTOADcolumn.dvspecifies the dependent variable, in this caseSFREQBUC.betweenspecifies the between-subject factors. In our case, each toad was either a buccal or a lung breather (BRTH.TYP).withinspecifies within-subject factors. In our case,O2LEVELwas varied over time.
Let’s now run aov_ez() to conduct the analysis. This
experiment is unbalanced, so we will use Type III sums of squares.
aov_ez() will automatically set contr.sum
contrasts and use Type III sums of squares for us. We will also use
afex_plot() to take a look at the data, which is not
possible until we fit the model.
set_sum_contrasts()
toad.anova <- aov_ez(id = "TOAD", dv = "SFREQBUC",
between = "BRTH.TYP", within = "O2LEVEL", data = d)
afex_plot(toad.anova, x = "O2LEVEL", trace = "BRTH.TYP", legend_title = "Breathing type",
error = "none") +
labs (x = "O2 concentration (%)", y = "Rate of buccal breathing (square-root transformed)") +
scale_x_discrete(labels = c("0", "5", "10", "15", "20", "30", "40", "50")) +
theme_few()
From this plot, it appears clear that there is an interaction between breathing type and O2 level: as O2 increases, breathing rate increases for lung breathers but decreases for buccal breathers.
Let’s plot the residuals. It is difficult to generate all the usual
residual plots from afex results, but inside the
toad.anova object there is a feature called lm
that contains the residuals. We can plot the most important ones:
toad.lm <- toad.anova$lm
opar <- par(mfrow = c(1, 2), pty = "s")
plot(fitted(toad.lm), studres(toad.lm), main = "Residuals vs. fitted")
abline(h = 0, lty = 2)
qqnorm(residuals(toad.lm))
qqline(residuals(toad.lm))
par(opar)
The residuals look pretty good. Although there is perhaps one outlier observation visible in the QQ plot that might warrant further investigation.
Let’s summarize the results. Given the choice, the MANOVA results are probably the more sensible results to use.
summary(toad.anova)##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 1695.15 1 131.63 19 244.6759 2.627e-12 ***
## BRTH.TYP 39.92 1 131.63 19 5.7622 0.02678 *
## O2LEVEL 25.75 7 100.17 133 4.8841 6.258e-05 ***
## BRTH.TYP:O2LEVEL 56.37 7 100.17 133 10.6928 1.228e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## O2LEVEL 0.013754 1.3415e-05
## BRTH.TYP:O2LEVEL 0.013754 1.3415e-05
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## O2LEVEL 0.42818 0.004333 **
## BRTH.TYP:O2LEVEL 0.42818 1.132e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## O2LEVEL 0.5172755 2.208358e-03
## BRTH.TYP:O2LEVEL 0.5172755 1.869971e-06print(toad.anova$Anova)##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.92794 244.676 1 19 2.627e-12 ***
## BRTH.TYP 1 0.23270 5.762 1 19 0.02678 *
## O2LEVEL 1 0.88489 14.277 7 13 3.542e-05 ***
## BRTH.TYP:O2LEVEL 1 0.67477 3.853 7 13 0.01730 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set_default_contrasts()As expected from looking at the graph, there is a significant interaction between breathing type and O2 level. We may wish to split up our data by breathing type to examine the effects more closely. This is left as an exercise for the reader.
This week’s task
Conduct an appropriate ANOVA for the Kluane dataset using the
aov_ez() function. Do you get the same results as with
aov() and gad()?
In GAD 2.0, there is a bug in this function. It calculates percentages using standard deviations instead of variances (and it includes the contribution of fixed effects). To obtain more accurate estimates of variance components (with confidence intervals), functions in the
VCApackage can be used.↩︎