Factorial
Nancy
26/10/2021
Two treatment factors, completely randomized design
data(poisons, package = "SMPracticals")
pmod <- lm(time ~ poison + treat, data = poisons)
anova(pmod)## Analysis of Variance Table
##
## Response: time
## Df Sum Sq Mean Sq F value Pr(>F)
## poison 2 1.03301 0.51651 20.643 5.704e-07 ***
## treat 3 0.92121 0.30707 12.273 6.697e-06 ***
## Residuals 42 1.05086 0.02502
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1with(poisons, interaction.plot(treat, poison, time))
with(poisons, interaction.plot(poison,treat,time))
Randomized block, 1 treatment factor and 1 blocking factor
data(oatvar, package = "faraway")
xtabs(yield ~ variety + block, data = oatvar)## block
## variety I II III IV V
## 1 296 357 340 331 348
## 2 402 390 431 340 320
## 3 437 334 426 320 296
## 4 303 319 310 260 242
## 5 469 405 442 487 394
## 6 345 342 358 300 308
## 7 324 339 357 352 220
## 8 488 374 401 338 320with(oatvar, interaction.plot(variety,block,yield))
with(oatvar, interaction.plot(block, variety, yield))
lmod <- lm(yield ~ block + variety, data = oatvar)
anova(lmod)## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## block 4 33396 8348.9 6.2449 0.001008 **
## variety 7 77524 11074.8 8.2839 1.804e-05 ***
## Residuals 28 37433 1336.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmod2 <- lm(yield ~ variety + block, data = oatvar)
anova(lmod2)## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## variety 7 77524 11074.8 8.2839 1.804e-05 ***
## block 4 33396 8348.9 6.2449 0.001008 **
## Residuals 28 37433 1336.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Faraway writes “An examination of which varieties give the highest yields and which are significantly better than others can now follow”, but he doesn’t do it. Here’s a start:
means <- with(oatvar, tapply(yield, variety, mean))
sort(means)## 4 7 6 1 3 2 8 5
## 286.8 318.4 330.6 334.4 362.6 376.6 384.2 439.4lmod3 <- aov(yield ~ variety + block, data = oatvar)
TukeyHSD(lmod3, "variety") # this is better than the code in the text, we don't want to compare block means.## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = yield ~ variety + block, data = oatvar)
##
## $variety
## diff lwr upr p adj
## 2-1 42.2 -33.42354 117.8235403 0.6095159
## 3-1 28.2 -47.42354 103.8235403 0.9192042
## 4-1 -47.6 -123.22354 28.0235403 0.4638341
## 5-1 105.0 29.37646 180.6235403 0.0021695
## 6-1 -3.8 -79.42354 71.8235403 0.9999998
## 7-1 -16.0 -91.62354 59.6235403 0.9965700
## 8-1 49.8 -25.82354 125.4235403 0.4078261
## 3-2 -14.0 -89.62354 61.6235403 0.9985165
## 4-2 -89.8 -165.42354 -14.1764597 0.0116440
## 5-2 62.8 -12.82354 138.4235403 0.1594614
## 6-2 -46.0 -121.62354 29.6235403 0.5061924
## 7-2 -58.2 -133.82354 17.4235403 0.2295334
## 8-2 7.6 -68.02354 83.2235403 0.9999741
## 4-3 -75.8 -151.42354 -0.1764597 0.0491482
## 5-3 76.8 1.17646 152.4235403 0.0445643
## 6-3 -32.0 -107.62354 43.6235403 0.8569155
## 7-3 -44.2 -119.82354 31.4235403 0.5549007
## 8-3 21.6 -54.02354 97.2235403 0.9798036
## 5-4 152.6 76.97646 228.2235403 0.0000094
## 6-4 43.8 -31.82354 119.4235403 0.5658118
## 7-4 31.6 -44.02354 107.2235403 0.8644352
## 8-4 97.4 21.77646 173.0235403 0.0050799
## 6-5 -108.8 -184.42354 -33.1764597 0.0014099
## 7-5 -121.0 -196.62354 -45.3764597 0.0003484
## 8-5 -55.2 -130.82354 20.4235403 0.2858981
## 7-6 -12.2 -87.82354 63.4235403 0.9993886
## 8-6 53.6 -22.02354 129.2235403 0.3194030
## 8-7 65.8 -9.82354 141.4235403 0.1237199