Statistics for Laboratory Scientists ( 140.615 )
Analysis of Variance - Two-Way Models
Example 1 from class: rats and lard
The data.
lard.dat <- data.frame(rsp = c(709,592,679,538,699,476,657,508,594,505,677,539),
sex = factor(rep(c("M","F"), rep(6,2))),
lard = factor(rep(c("fresh","rancid"), 6)))
lard.dat## rsp sex lard
## 1 709 M fresh
## 2 592 M rancid
## 3 679 M fresh
## 4 538 M rancid
## 5 699 M fresh
## 6 476 M rancid
## 7 657 F fresh
## 8 508 F rancid
## 9 594 F fresh
## 10 505 F rancid
## 11 677 F fresh
## 12 539 F rancidstr(lard.dat)## 'data.frame': 12 obs. of 3 variables:
## $ rsp : num 709 592 679 538 699 476 657 508 594 505 ...
## $ sex : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 1 1 1 1 ...
## $ lard: Factor w/ 2 levels "fresh","rancid": 1 2 1 2 1 2 1 2 1 2 ...summary(lard.dat)## rsp sex lard
## Min. :476.0 F:6 fresh :6
## 1st Qu.:530.5 M:6 rancid:6
## Median :593.0
## Mean :597.8
## 3rd Qu.:677.5
## Max. :709.0Plot the data, one factor at a time.
set.seed(1)
par(las=1)
stripchart(jitter(rsp,factor=2) ~ sex, data=lard.dat, ylab="response", pch=1, vertical=TRUE, xlim=c(0.5,2.5))
stripchart(jitter(rsp,factor=2) ~ lard, data=lard.dat, ylab="response", pch=1, vertical=TRUE, xlim=c(0.5,2.5))
Plot the data, both factors.
stripchart(split(jitter(lard.dat$rsp,factor=2), list(lard.dat$sex, lard.dat$lard)), ylab="response", pch=1,
vertical=TRUE, xlim=c(0.5,4.5))
The one-way ANOVA. The “:” symbol pastes the two factors together.
summary(aov(rsp ~ sex:lard, data=lard.dat))## Df Sum Sq Mean Sq F value Pr(>F)
## sex:lard 3 65904 21968 15.06 0.00118 **
## Residuals 8 11667 1458
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The two-way ANOVA. The “*” symbol includes the interaction.
summary(aov(rsp ~ sex * lard, data=lard.dat))## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 3781 3781 2.593 0.146036
## lard 1 61204 61204 41.969 0.000192 ***
## sex:lard 1 919 919 0.630 0.450255
## Residuals 8 11667 1458
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The model with interaction can also be written this way.
summary(aov(rsp ~ sex + lard + sex:lard, data=lard.dat))## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 3781 3781 2.593 0.146036
## lard 1 61204 61204 41.969 0.000192 ***
## sex:lard 1 919 919 0.630 0.450255
## Residuals 8 11667 1458
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A two-way ANOVA with an additive model. The interaction SS are then included in the “residual” SS.
summary(aov(rsp ~ sex + lard, data=lard.dat))## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 3781 3781 2.704 0.135
## lard 1 61204 61204 43.768 9.75e-05 ***
## Residuals 9 12585 1398
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ANOVA using only the factor lard.
summary(aov(rsp ~ lard, data=lard.dat))## Df Sum Sq Mean Sq F value Pr(>F)
## lard 1 61204 61204 37.4 0.000113 ***
## Residuals 10 16366 1637
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The same inference using a 2-sample t-test.
subset(lard.dat, lard=="fresh")## rsp sex lard
## 1 709 M fresh
## 3 679 M fresh
## 5 699 M fresh
## 7 657 F fresh
## 9 594 F fresh
## 11 677 F freshsubset(lard.dat, lard=="rancid")## rsp sex lard
## 2 592 M rancid
## 4 538 M rancid
## 6 476 M rancid
## 8 508 F rancid
## 10 505 F rancid
## 12 539 F rancidx <- subset(lard.dat, lard=="fresh")$rsp
x## [1] 709 679 699 657 594 677y <- subset(lard.dat, lard=="rancid")$rsp
y## [1] 592 538 476 508 505 539t.test(x, y, var.equal=TRUE)##
## Two Sample t-test
##
## data: x and y
## t = 6.1153, df = 10, p-value = 0.0001134
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 90.7912 194.8755
## sample estimates:
## mean of x mean of y
## 669.1667 526.3333t.test(x, y, var.equal=TRUE)$p.val## [1] 0.0001133845t.test(x, y, var.equal=TRUE)$stat## t
## 6.115285t.test(x, y, var.equal=TRUE)$stat^2## t
## 37.39671Within-group means.
means <- tapply(lard.dat$rsp, list(lard.dat$sex,lard.dat$lard), mean)
means## fresh rancid
## F 642.6667 517.3333
## M 695.6667 535.3333Means by sex.
sex.means <- tapply(lard.dat$rsp, lard.dat$sex, mean)
sex.means## F M
## 580.0 615.5Means by treatment.
lard.means <- tapply(lard.dat$rsp, lard.dat$lard, mean)
lard.means## fresh rancid
## 669.1667 526.3333An interaction plot.
interaction.plot(lard.dat$sex, lard.dat$lard, lard.dat$rsp, lty=1, col=c("blue","red"), lwd=2)
The same plot, but reversing the role of “sex” and “lard”.
interaction.plot(lard.dat$lard, lard.dat$sex, lard.dat$rsp, lty=1, col=c("blue","red"), lwd=2)
Example 2 from class: mouse, food and temperature
The data.
library(SPH.140.615)
example(mouse)##
## mouse> mouse
## rsp food temp
## 1 62.69 restricted 8deg
## 2 54.07 restricted 8deg
## 3 65.73 restricted 8deg
## 4 62.98 restricted 8deg
## 5 72.60 restricted 18deg
## 6 70.97 restricted 18deg
## 7 74.32 restricted 18deg
## 8 53.02 restricted 18deg
## 9 46.22 restricted 18deg
## 10 59.10 restricted 18deg
## 11 61.79 restricted 18deg
## 12 61.89 restricted 18deg
## 13 95.73 unrestricted 8deg
## 14 63.95 unrestricted 8deg
## 15 144.30 unrestricted 8deg
## 16 144.30 unrestricted 8deg
## 17 101.19 unrestricted 18deg
## 18 76.88 unrestricted 18deg
## 19 74.08 unrestricted 18deg
## 20 81.40 unrestricted 18deg
## 21 66.58 unrestricted 18deg
## 22 84.38 unrestricted 18deg
## 23 118.95 unrestricted 18deg
## 24 118.95 unrestricted 18deg
##
## mouse> str(mouse)
## 'data.frame': 24 obs. of 3 variables:
## $ rsp : num 62.7 54.1 65.7 63 72.6 ...
## $ food: Factor w/ 2 levels "restricted","unrestricted": 1 1 1 1 1 1 1 1 1 1 ...
## $ temp: Factor w/ 2 levels "18deg","8deg": 2 2 2 2 1 1 1 1 1 1 ...
##
## mouse> summary(mouse)
## rsp food temp
## Min. : 46.22 restricted :12 18deg:16
## 1st Qu.: 62.49 unrestricted:12 8deg : 8
## Median : 71.78
## Mean : 79.84
## 3rd Qu.: 87.22
## Max. :144.30
##
## mouse> tapply(mouse$rsp,list(mouse$food,mouse$temp),mean)
## 18deg 8deg
## restricted 62.48875 61.3675
## unrestricted 90.30125 112.0700Plot the data.
set.seed(1)
par(las=1)
stripchart(split(jitter(mouse$rsp,amount=3),list(mouse$food, mouse$temp)),
ylab="", pch=1, vertical=TRUE, xlim=c(0.5,4.5))
stripchart(split(jitter(mouse$rsp,amount=3),list(mouse$food, mouse$temp)),
ylab="", pch=1, vertical=TRUE, xlim=c(0.5,4.5), log="y")
stripchart(split(jitter(mouse$rsp,amount=3),list(mouse$temp, mouse$food)),
ylab="", pch=1, vertical=TRUE, xlim=c(0.5,4.5), log="y")
Numbers of individuals per condition,
tapply(mouse$rsp, list(mouse$food, mouse$temp), length)## 18deg 8deg
## restricted 8 4
## unrestricted 8 4The two-way ANOVA, with log-transformed response.
aov.out <- aov(log(rsp) ~ food * temp, data=mouse)
summary(aov.out)## Df Sum Sq Mean Sq F value Pr(>F)
## food 1 1.0725 1.0725 21.421 0.000162 ***
## temp 1 0.0407 0.0407 0.813 0.378025
## food:temp 1 0.0498 0.0498 0.995 0.330479
## Residuals 20 1.0014 0.0501
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The interaction plot
interaction.plot(mouse$food, mouse$temp, log(mouse$rsp), lwd=2, col=c("blue","red"), lty=1)
Reverse role of food and temperature.
interaction.plot(mouse$temp, mouse$food, log(mouse$rsp), lwd=2, col=c("blue","red"), lty=1)
Two-way ANOVA without interaction.
aov.out <- aov(log(rsp) ~ food + temp, data=mouse)
summary(aov.out)## Df Sum Sq Mean Sq F value Pr(>F)
## food 1 1.0725 1.0725 21.426 0.000145 ***
## temp 1 0.0407 0.0407 0.813 0.377463
## Residuals 21 1.0512 0.0501
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ANOVA only using the factor food.
aov.out <- aov(log(rsp) ~ food, data=mouse)
anova(aov.out)## Analysis of Variance Table
##
## Response: log(rsp)
## Df Sum Sq Mean Sq F value Pr(>F)
## food 1 1.0725 1.07250 21.61 0.0001238 ***
## Residuals 22 1.0919 0.04963
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Example 3 from class: flies, density, strain, no replicates
The data.
flies <- data.frame(rsp = c( 9.6, 9.3, 9.3,
10.6, 9.1, 9.2,
9.8, 9.3, 9.5,
10.7, 9.1,10.0,
11.1,11.1,10.4,
10.9,11.8,10.8,
12.8,10.6,10.7),
strain = factor(rep(c("OL","BELL","bwb"),7)),
density = factor(rep(c(60,80,160,320,640,1280,2560),rep(3,7))))
summary(flies)## rsp strain density
## Min. : 9.10 BELL:7 60 :3
## 1st Qu.: 9.30 bwb :7 80 :3
## Median :10.40 OL :7 160 :3
## Mean :10.27 320 :3
## 3rd Qu.:10.80 640 :3
## Max. :12.80 1280:3
## 2560:3Plot the data.
par(las=1)
interaction.plot(flies$density, flies$strain, flies$rsp, lwd=2, col=c("blue","red","green3"), lty=1, type="b",
xlab="density", ylab="response")
Note: we can’t fit an interaction. If we try, we have no ability to get an error SS, so we don’t get p-values either.
flies.aovA <- aov(rsp ~ strain * density, data=flies)
summary(flies.aovA) ## Df Sum Sq Mean Sq
## strain 2 2.789 1.3943
## density 6 12.543 2.0905
## strain:density 12 4.111 0.3426Assume an additive model.
flies.aovB = aov(rsp ~ strain + density, data=flies)
summary(flies.aovB) ## Df Sum Sq Mean Sq F value Pr(>F)
## strain 2 2.789 1.3943 4.069 0.04476 *
## density 6 12.543 2.0905 6.101 0.00395 **
## Residuals 12 4.111 0.3426
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1