Note: This document was converted to R-Markdown from this page by M. Drew LaMar. You can download the R-Markdown here.

Download the R code on this page as a single file here

New methods

Hover over a function argument for a short description of its meaning. The variable names are plucked from the examples further below.

Fixed effects analysis of variance (ANOVA):

circadianAnova <- lm(shift ~ treatment, data = circadian)
anova(circadianAnova)

Unplanned comparisons (Tukey-Kramer tests) using the multcomp package:

circadianPlanned <- glht(circadianAnova, linfct = mcp(treatment = c(“knee - control = 0”)))
summary(circadianPlanned)

Kruskal-Wallis test:

kruskal.test(shift ~ treatment, data = circadian)

Other new methods:

Example 15.1. Knees who say night

Analysis of variance, comparing phase shift in the circadian rhythm of melatonin production in participants given alternative light treatments. Also, the nonparametric Kruskal-Wallis test. Finally, we use the same data to demonstrate planned comparisons and unplanned comparisons (Tukey-Kramer tests).

Read and inspect the data.

circadian <- read.csv(url("http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter15/chap15e1KneesWhoSayNight.csv"))

Set the preferred ordering of groups in tables and graphs.

circadian$treatment <- factor(circadian$treatment, levels = c("control", "knee", "eyes"))

Table of descriptive statistics by treatment group (Table 15.1-1).

meanShift <- tapply(circadian$shift, circadian$treatment, mean)
sdevShift <- tapply(circadian$shift, circadian$treatment, sd)
n         <- tapply(circadian$shift, circadian$treatment, length)
data.frame(mean = meanShift, std.dev = sdevShift, n = n)
##               mean   std.dev n
## control -0.3087500 0.6175629 8
## knee    -0.3357143 0.7908193 7
## eyes    -1.5514286 0.7063151 7

Strip chart of circadian rhythm shift by treatment, with standard error bars. This makes use of the descriptive statistics calculated in an earlier step. The error bars are added as line segments, offset along the x-axis by the amount adjustAmount. The means are added using points.

stripchart(shift ~ treatment, data = circadian, method = "jitter", vertical = TRUE)
seShift <- sdevShift / sqrt(n)
adjustAmount <- 0.15
segments(c(1,2,3) + adjustAmount, meanShift - seShift, c(1,2,3) + adjustAmount, meanShift + seShift)
points(meanShift ~ c( c(1,2,3) + adjustAmount ))

Commands for a stripchart with more options are shown here.

par(bty="l")
adjustAmount <- 0.15
stripchart(shift ~ treatment, data = circadian, method = "jitter",
    vertical = TRUE, las = 1, pch = 1, xlab = "Light treatment",
    ylab = "Shift in circadian rhythm (h)", col = "firebrick", 
    cex = 1.2, ylim = c(-3, 1))
segments( c(1,2,3) + adjustAmount, meanShift - seShift, 
      c(1,2,3) + adjustAmount, meanShift + seShift )
points(meanShift ~ c( c(1,2,3) + adjustAmount ), pch = 16, col = "firebrick")

Fixed effects ANOVA table (Table 15.1-2). This is done in two steps. The first step involves fitting the ANOVA model to the data using lm (“lm” stands for “linear model”, of which ANOVA is one type). Then we use the command anova to assemble the ANOVA table.

circadianAnova <- lm(shift ~ treatment, data = circadian)
anova(circadianAnova)
## Analysis of Variance Table
## 
## Response: shift
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## treatment  2 7.2245  3.6122  7.2894 0.004472 **
## Residuals 19 9.4153  0.4955                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\(R^{2}\) indicating the fraction of variation in the response variable “explained” by treatment. This is again done in two steps. The first step calculates a bunch of useful quantities from the ANOVA model object previously created with a lm command. The second step shows the \(R^2\) value.

circadianAnovaSummary <- summary(circadianAnova)
circadianAnovaSummary$r.squared
## [1] 0.4341684

Kruskal-Wallis test, a nonparametric method to compare more than two groups. The method is not needed for the circadian rhythm data, because assumptions of ANOVA are met, but we include it here to demonstrate the method. The formula is the same as that used with lm.

kruskal.test(shift ~ treatment, data = circadian)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  shift by treatment
## Kruskal-Wallis chi-squared = 9.4231, df = 2, p-value = 0.008991

Planned and unplanned comparisons between means. A planned comparison is one planned during the design of the study, before the data were collected. To use the method you need a good justification to focus on a specific comparison of two treatments. Only a small number of planned comparisons are allowed. If you don’t have a good prior justification for a specific comparison, use unplanned comparison instead. Unplanned comparisons typically involve testing differences between all pairs of means, and methods provide needed protection against rising Type I errors that would otherwise result from multiple testing.

If you haven’t already done so, you’ll need to install the multicomp package (this needs to be done just once per computer).

if (!require("multcomp")) {install.packages("multcomp", dependencies = TRUE, repos="http://cran.rstudio.com/")}
library(multcomp)

Planned comparison between “control” and “knee” treatments. Begin by loading the multicomp package. The commands shown will give the 95% confidence interval and the planned \(t\)-test of a difference between the treatment means.

circadianPlanned <- glht(circadianAnova, linfct = mcp(treatment = c("knee - control = 0")))
confint(circadianPlanned)
## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: lm(formula = shift ~ treatment, data = circadian)
## 
## Quantile = 2.093
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##                     Estimate lwr      upr     
## knee - control == 0 -0.02696 -0.78951  0.73558
summary(circadianPlanned)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: lm(formula = shift ~ treatment, data = circadian)
## 
## Linear Hypotheses:
##                     Estimate Std. Error t value Pr(>|t|)
## knee - control == 0 -0.02696    0.36433  -0.074    0.942
## (Adjusted p values reported -- single-step method)

Unplanned comparisons (Tukey-Kramer tests) between all pairs of means. The raw data for the “Wood-wide web” example (Example 15.4) are not available, so we have used the circadian rhythm data to demonstrate the R method here instead. The output table will say “t” but it is actually “q” as we describe in the book.

library(multcomp)
circadianTukey <- glht(circadianAnova, linfct = mcp(treatment = "Tukey"))
summary(circadianTukey)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = shift ~ treatment, data = circadian)
## 
## Linear Hypotheses:
##                     Estimate Std. Error t value Pr(>|t|)   
## knee - control == 0 -0.02696    0.36433  -0.074  0.99699   
## eyes - control == 0 -1.24268    0.36433  -3.411  0.00784 **
## eyes - knee == 0    -1.21571    0.37628  -3.231  0.01198 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Example 15.6. Walking stick limbs

Random effects ANOVA to estimate variance components and calculate repeatability of measurements of femur length in walking stick insects.

Read and inspect data. Data are in “long” format. Femur length is one column, with another variable indicating specimen identity. Each specimen was measured twice and therefore takes up two rows.

walkingstick <- read.csv(url("http://www.zoology.ubc.ca/~schluter/WhitlockSchluter/wp-content/data/chapter15/chap15e6WalkingStickFemurs.csv"))
head(walkingstick)
##   specimen femurLength
## 1        1        0.26
## 2        1        0.26
## 3        2        0.23
## 4        2        0.19
## 5        3        0.25
## 6        3        0.23

Descriptive statistics for each specimen, which we need to produce the strip chart. tapply is used to gets the mean, smallest measurement, largest measurement, and specimen id of each specimen.

meanFemur <- tapply(walkingstick$femurLength, walkingstick$specimen, mean)
minFemur <- tapply(walkingstick$femurLength, walkingstick$specimen, min)
maxFemur <- tapply(walkingstick$femurLength, walkingstick$specimen, max)
specimen <- tapply(walkingstick$specimen, walkingstick$specimen, unique)

Stripchart, with line segments connecting the two measurements.

stripchart(femurLength ~ specimen, data = walkingstick, vertical = TRUE, xlab = "Specimen")
segments(specimen, minFemur, specimen, maxFemur)

Here are commands for a prettier figure in which specimens are ordered along the x-axis by mean femur length (Figure 15.6-1).

par(bty="l") 
stripchart(minFemur[order(meanFemur, specimen)] ~ c(1:length(specimen)),
    vertical = TRUE, xaxt = "n", pch = 16, col = "firebrick", las = 1,
    cex = 1.2, ylim = range(c(minFemur, maxFemur)), ylab = "Femur length (mm)",
    xlab = "Individual walking sticks")
stripchart(maxFemur[order(meanFemur, specimen)] ~ c(1:length(specimen)),
    vertical = TRUE, add = TRUE, pch = 16, col = "firebrick", las = 1, cex = 1.2)
segments(c(1:length(specimen)), minFemur[order(meanFemur, specimen)], 
     c(1:length(specimen)), maxFemur[order(meanFemur, specimen)])

Fit the random effects ANOVA using lme. The random effects ANOVA function requires two formulas, rather than just one. The first formula (beginning with fixed =) is for the fixed effect. The walking stick insect example doesn’t include a fixed-effect variable, so we just provide a symbol for a constant in the formula (~ 1), representing the grand mean. The second formula (beginning with random =) is for the random effect. In this example, the individual specimens are the random groups, and the second formula indicates this (the ~ 1 in the formula below indicates that each specimen has its own mean). You will need to load the nlme library to begin.

library(nlme)
walkingstickAnova <- lme(fixed = femurLength ~ 1, random = ~ 1|specimen, data = walkingstick)

Obtain the variance components for the random effects using VarCorr. The output includes the standard deviation and variance for both components of random variation in the random effects model for this example. The first is the variance among the specimen means. This is the variance among groups, and is confusingly labeled “Intercept” in the output. The second component is the variance among measurements made on the same individuals. This is the within group variance, also known as the error mean square, and is labeled “Residual” in the output.

walkingstickVarcomp <- VarCorr(walkingstickAnova)
walkingstickVarcomp
## specimen = pdLogChol(1) 
##             Variance     StdDev    
## (Intercept) 0.0010539182 0.03246411
## Residual    0.0003559996 0.01886795

Calculate the repeatability of the walking stick femur measurements using the estimates of the variance components.

varAmong  <- as.numeric( walkingstickVarcomp[1,1] )
varWithin <- as.numeric( walkingstickVarcomp[2,1] )
repeatability <- varAmong / (varAmong + varWithin)
repeatability
## [1] 0.7475033

Note that lme doesn’t test the significance of the random effects (whether the specimen means are significantly different from one another), since the method basically assumes the presence of variance among random means. As a result, there is no ANOVA table for random effects.