This report was automatically generated with the R package knitr (version 1.5).
library(faraway)
data(oatvar, package = "faraway")
xtabs(yield ~ variety + block, 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 320
plot(yield ~ variety, oatvar)
plot(yield ~ block, oatvar)
with(oatvar, interaction.plot(variety, block, yield))
with(oatvar, interaction.plot(block, variety, yield))
require(ggplot2)
ggplot(oatvar, aes(x = variety, y = yield, group = block, linetype = block)) +
geom_line() + theme(legend.position = "top", legend.direction = "horizontal")
ggplot(oatvar, aes(x = block, y = yield, group = block, linetype = block)) +
geom_line()
lmod <- lm(yield ~ block + variety, oatvar)
anova(lmod)
Analysis of Variance Table
Response: yield
Df Sum Sq Mean Sq F value Pr(>F)
block 4 33396 8349 6.24 0.001
variety 7 77524 11075 8.28 1.8e-05
Residuals 28 37433 1337
anova(lm(yield ~ block, oatvar))
Analysis of Variance Table
Response: yield
Df Sum Sq Mean Sq F value Pr(>F)
block 4 33396 8349 2.54 0.057
Residuals 35 114957 3284
anova(lm(yield ~ variety + block, oatvar))
Analysis of Variance Table
Response: yield
Df Sum Sq Mean Sq F value Pr(>F)
variety 7 77524 11075 8.28 1.8e-05
block 4 33396 8349 6.24 0.001
Residuals 28 37433 1337
anova(lm(yield ~ block + variety, subset = -1, oatvar))
Analysis of Variance Table
Response: yield
Df Sum Sq Mean Sq F value Pr(>F)
block 4 38581 9645 8.41 0.00015
variety 7 75339 10763 9.38 7.3e-06
Residuals 27 30968 1147
anova(lm(yield ~ variety + block, subset = -1, oatvar))
Analysis of Variance Table
Response: yield
Df Sum Sq Mean Sq F value Pr(>F)
variety 7 75902 10843 9.45 6.8e-06
block 4 38018 9504 8.29 0.00017
Residuals 27 30968 1147
drop1(lm(yield ~ variety + block, subset = -1, oatvar), test = "F")
Single term deletions
Model:
yield ~ variety + block
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 30968 284
variety 7 75339 106307 319 9.38 7.3e-06
block 4 38018 68986 308 8.29 0.00017
plot(lmod, which = 1:2, caption = NULL)
varcoefs <- c(0, coef(lmod)[6:12])
blockcoefs <- c(0, coef(lmod)[2:5])
ab <- rep(varcoefs, each = 5) * rep(blockcoefs, 8)
h <- update(lmod, . ~ . + ab)
anova(h)
Analysis of Variance Table
Response: yield
Df Sum Sq Mean Sq F value Pr(>F)
block 4 33396 8349 6.06 0.0013
variety 7 77524 11075 8.03 2.8e-05
ab 1 213 213 0.15 0.6974
Residuals 27 37220 1379
lmcrd <- lm(yield ~ variety, oatvar)
summary(lmcrd)$sig
[1] 47.05
summary(lmod)$sig
[1] 36.56
(47.047/36.564)^2
[1] 1.656
data(abrasion, package = "faraway")
xtabs(wear ~ run + position, abrasion)
position
run 1 2 3 4
1 235 236 218 268
2 251 241 227 229
3 234 273 274 226
4 195 270 230 225
matrix(abrasion$material, 4, 4, byrow = TRUE)
[,1] [,2] [,3] [,4]
[1,] "C" "D" "B" "A"
[2,] "A" "B" "D" "C"
[3,] "D" "C" "A" "B"
[4,] "B" "A" "C" "D"
ggplot(abrasion, aes(x = run, y = wear, shape = position, group = material)) +
geom_point() + geom_line(aes(linetype = material))
ggplot(abrasion, aes(x = position, y = wear, shape = run, group = material)) +
geom_point() + geom_line(aes(linetype = material))
lmod <- lm(wear ~ material + run + position, abrasion)
drop1(lmod, test = "F")
Single term deletions
Model:
wear ~ material + run + position
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 368 70.1
material 3 4622 4989 105.9 25.15 0.00085
run 3 987 1354 85.0 5.37 0.03901
position 3 1468 1836 89.9 7.99 0.01617
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 254.75 6.19 41.17 1.4e-08
materialB -45.75 5.53 -8.27 0.00017
materialC -24.00 5.53 -4.34 0.00489
materialD -35.25 5.53 -6.37 0.00070
run2 -2.25 5.53 -0.41 0.69842
run3 12.50 5.53 2.26 0.06466
run4 -9.25 5.53 -1.67 0.14566
position2 26.25 5.53 4.74 0.00318
position3 8.50 5.53 1.54 0.17545
position4 8.25 5.53 1.49 0.18661
n = 16, p = 10, Residual SE = 7.83, R-Squared = 0.95
qtukey(0.95, 4, 6) * 5.53/sqrt(2)
[1] 19.14
scoefs <- c(0, coef(lmod)[2:4])
outer(scoefs, scoefs, "-")
materialB materialC materialD
0.00 45.75 24.00 35.25
materialB -45.75 0.00 -21.75 -10.50
materialC -24.00 21.75 0.00 11.25
materialD -35.25 10.50 -11.25 0.00
lmodr <- lm(wear ~ material, abrasion)
(summary(lmodr)$sig/summary(lmod)$sig)^2
[1] 3.84
data(rabbit, package = "faraway")
xtabs(gain ~ treat + block, rabbit)
block
treat b1 b10 b2 b3 b4 b5 b6 b7 b8 b9
a 0.0 37.3 40.1 0.0 44.9 0.0 0.0 45.2 44.0 0.0
b 32.6 0.0 38.1 0.0 0.0 0.0 37.3 40.6 0.0 30.6
c 35.2 0.0 40.9 34.6 43.9 40.9 0.0 0.0 0.0 0.0
d 0.0 42.3 0.0 37.5 0.0 37.3 0.0 37.9 0.0 27.5
e 0.0 0.0 0.0 0.0 40.8 32.0 40.5 0.0 38.5 20.6
f 42.2 41.7 0.0 34.3 0.0 0.0 42.8 0.0 51.9 0.0
ggplot(rabbit, aes(y = gain, x = block, shape = treat)) + geom_point() + theme(legend.position = "top",
legend.direction = "horizontal")
ggplot(rabbit, aes(y = gain, x = treat)) + geom_point(position = position_jitter(width = 0.1))
lmod <- lm(gain ~ block + treat, rabbit)
drop1(lmod, test = "F")
Single term deletions
Model:
gain ~ block + treat
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 151 78.4
block 9 596 747 108.4 6.59 0.00076
treat 5 159 309 90.0 3.16 0.03817
plot(lmod, which = 1:2)
sumary(lmod)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 36.0139 2.5886 13.91 5.6e-10
blockb10 3.2972 2.7960 1.18 0.2567
blockb2 4.1333 2.6943 1.53 0.1458
blockb3 -1.8028 2.6943 -0.67 0.5136
blockb4 8.7944 2.7960 3.15 0.0067
blockb5 2.3056 2.7960 0.82 0.4225
blockb6 5.4083 2.6943 2.01 0.0631
blockb7 5.7778 2.7960 2.07 0.0565
blockb8 9.4278 2.7960 3.37 0.0042
blockb9 -7.4806 2.7960 -2.68 0.0173
treatb -1.7417 2.2418 -0.78 0.4493
treatc 0.4000 2.2418 0.18 0.8608
treatd 0.0667 2.2418 0.03 0.9767
treate -5.2250 2.2418 -2.33 0.0341
treatf 3.3000 2.2418 1.47 0.1617
n = 30, p = 15, Residual SE = 3.17, R-Squared = 0.86
qtukey(0.95, 6, 15)
[1] 4.595
4.59 * 2.24/sqrt(2)
[1] 7.27
tcoefs <- c(0, coef(lmod)[11:15])
abs(outer(tcoefs, tcoefs, "-")) > 7.27
treatb treatc treatd treate treatf
FALSE FALSE FALSE FALSE FALSE FALSE
treatb FALSE FALSE FALSE FALSE FALSE FALSE
treatc FALSE FALSE FALSE FALSE FALSE FALSE
treatd FALSE FALSE FALSE FALSE FALSE FALSE
treate FALSE FALSE FALSE FALSE FALSE TRUE
treatf FALSE FALSE FALSE FALSE TRUE FALSE
lmodt <- lm(gain ~ treat, rabbit)
(summary(lmodt)$sig/summary(lmod)$sig)^2
[1] 3.095
The R session information (including the OS info, R version and all packages used):
sessionInfo()
R version 3.1.0 (2014-04-10)
Platform: x86_64-apple-darwin13.1.0 (64-bit)
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] graphics grDevices utils datasets methods stats base
other attached packages:
[1] faraway_1.0.6 knitr_1.5 ggplot2_0.9.3.1
loaded via a namespace (and not attached):
[1] colorspace_1.2-4 dichromat_2.0-0 digest_0.6.4
[4] evaluate_0.5.3 formatR_0.10 grid_3.1.0
[7] gtable_0.1.2 labeling_0.2 MASS_7.3-31
[10] munsell_0.4.2 plyr_1.8.1 proto_0.3-10
[13] RColorBrewer_1.0-5 Rcpp_0.11.1 reshape2_1.2.2
[16] scales_0.2.3 stringr_0.6.2 tools_3.1.0
Sys.time()
[1] "2014-06-16 14:02:37 BST"