Aliases: se.contrast se.contrast.aov se.contrast.aovlist
Keywords: models
### ** Examples ## From Venables and Ripley (2002) p.165. N <- c(0,1,0,1,1,1,0,0,0,1,1,0,1,1,0,0,1,0,1,0,1,1,0,0) P <- c(1,1,0,0,0,1,0,1,1,1,0,0,0,1,0,1,1,0,0,1,0,1,1,0) K <- c(1,0,0,1,0,1,1,0,0,1,0,1,0,1,1,0,0,0,1,1,1,0,1,0) yield <- c(49.5,62.8,46.8,57.0,59.8,58.5,55.5,56.0,62.8,55.8,69.5, 55.0, 62.0,48.8,45.5,44.2,52.0,51.5,49.8,48.8,57.2,59.0,53.2,56.0) npk <- data.frame(block = gl(6,4), N = factor(N), P = factor(P), K = factor(K), yield = yield) ## Set suitable contrasts. options(contrasts = c("contr.helmert", "contr.poly")) npk.aov1 <- aov(yield ~ block + N + K, data = npk) se.contrast(npk.aov1, list(N == "0", N == "1"), data = npk)
[1] 1.609175
# or via a matrix cont <- matrix(c(-1,1), 2, 1, dimnames = list(NULL, "N")) se.contrast(npk.aov1, cont[N, , drop = FALSE]/12, data = npk)
N 1.609175
## test a multi-stratum model npk.aov2 <- aov(yield ~ N + K + Error(block/(N + K)), data = npk) se.contrast(npk.aov2, list(N == "0", N == "1"))
[1] 1.812166
## an example looking at an interaction contrast ## Dataset from R.E. Kirk (1995) ## 'Experimental Design: procedures for the behavioral sciences' score <- c(12, 8,10, 6, 8, 4,10,12, 8, 6,10,14, 9, 7, 9, 5,11,12, 7,13, 9, 9, 5,11, 8, 7, 3, 8,12,10,13,14,19, 9,16,14) A <- gl(2, 18, labels = c("a1", "a2")) B <- rep(gl(3, 6, labels = c("b1", "b2", "b3")), 2) fit <- aov(score ~ A*B) cont <- c(1, -1)[A] * c(1, -1, 0)[B] sum(cont) # 0
[1] 0
sum(cont*score) # value of the contrast
[1] -18
se.contrast(fit, as.matrix(cont))
Contrast 1 14.24547
(t.stat <- sum(cont*score)/se.contrast(fit, as.matrix(cont)))
Contrast 1 -1.26356
summary(fit, split = list(B = 1:2), expand.split = TRUE)
Df Sum Sq Mean Sq F value Pr(>F) A 1 18.78 18.78 2.221 0.14661 B 2 62.00 31.00 3.666 0.03763 * B: C1 1 1.50 1.50 0.177 0.67662 B: C2 1 60.50 60.50 7.155 0.01199 * A:B 2 81.56 40.78 4.823 0.01527 * A:B: C1 1 13.50 13.50 1.597 0.21612 A:B: C2 1 68.06 68.06 8.049 0.00809 ** Residuals 30 253.67 8.46 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## t.stat^2 is the F value on the A:B: C1 line (with Helmert contrasts) ## Now look at all three interaction contrasts cont <- c(1, -1)[A] * cbind(c(1, -1, 0), c(1, 0, -1), c(0, 1, -1))[B,] se.contrast(fit, cont) # same, due to balance.
Contrast 1 Contrast 2 Contrast 3 14.24547 14.24547 14.24547
rm(A, B, score) ## multi-stratum example where efficiencies play a role ## An example from Yates (1932), ## a 2^3 design in 2 blocks replicated 4 times Block <- gl(8, 4) A <- factor(c(0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1, 0,1,0,1,0,1,0,1,0,1,0,1)) B <- factor(c(0,0,1,1,0,0,1,1,0,1,0,1,1,0,1,0,0,0,1,1, 0,0,1,1,0,0,1,1,0,0,1,1)) C <- factor(c(0,1,1,0,1,0,0,1,0,0,1,1,0,0,1,1,0,1,0,1, 1,0,1,0,0,0,1,1,1,1,0,0)) Yield <- c(101, 373, 398, 291, 312, 106, 265, 450, 106, 306, 324, 449, 272, 89, 407, 338, 87, 324, 279, 471, 323, 128, 423, 334, 131, 103, 445, 437, 324, 361, 302, 272) aovdat <- data.frame(Block, A, B, C, Yield) fit <- aov(Yield ~ A + B * C + Error(Block), data = aovdat) cont1 <- c(-1, 1)[A]/32 # Helmert contrasts cont2 <- c(-1, 1)[B] * c(-1, 1)[C]/32 cont <- cbind(A = cont1, BC = cont2) colSums(cont*Yield) # values of the contrasts
A BC 10.40625 -20.90625
se.contrast(fit, as.matrix(cont))
A BC 3.377196 3.899650
## No test: # comparison with lme library(nlme) fit2 <- lme(Yield ~ A + B*C, random = ~1 | Block, data = aovdat) summary(fit2)$tTable # same estimates, similar (but smaller) se's.
Value Std.Error DF t-value p-value (Intercept) 291.59375 3.287945 20 88.685710 1.944005e-27 A1 10.40625 3.287945 20 3.164971 4.870019e-03 B1 70.96875 3.287945 20 21.584530 2.507178e-15 C1 93.34375 3.287945 20 28.389692 1.234427e-17 B1:C1 -20.90625 3.287945 20 -6.358455 3.324598e-06
## End(No test)