Hi, the function TukeyHSD gives incorrect results for balanced incomplete block designs, as the example below shows, but I can only half fix it. There are two problems, 1. It uses model.tables to estimate treatment means, 2. It uses the wrong standard error The first problem can be fixed using dummy.coef, if the lines> TukeyHSD.aovfunction (x, which = seq(along = tabs), ordered = FALSE, conf.level = 0.95, ...) { mm <- model.tables(x, "means") tabs <- mm$tables[-1] tabs <- tabs[which] ... are altered to> TukeyHSD.aovfunction (x, which = seq(along = tabs), ordered = FALSE, conf.level = 0.95, ...) { mm <- model.tables(x, "means") tabs <- dummy.coef(x) ## This line changed tabs <- tabs[which] ... it calculates the right treatment differences (provided the block term preceeds the treatment in the model formula). To solve the second problem, I can manually get the correct se from summary.lm, but I can't figure out how to extract the right se from summary.lm programmatically. I don't think these changes would make TukeyHSD fail on other cases, but I haven't tested so I'm not sure. I hope these few observations can help someone a bit more knowledgable than I fix this problem. Simon. ---
Sorry, Forgot to include the example> dConsumer Pillows Comfort 1 1 A 59 2 1 B 26 3 1 C 38 4 2 D 85 5 2 E 92 6 2 F 69 7 3 G 74 8 3 H 52 9 3 I 27 10 4 A 63 11 4 D 70 12 4 G 68 13 5 B 26 14 5 E 98 15 5 H 59 16 6 C 31 17 6 F 60 18 6 I 35 19 7 A 62 20 7 E 85 21 7 I 30 22 8 B 23 23 8 F 73 24 8 G 75 25 9 C 49 26 9 D 74 27 9 H 51 28 10 A 52 29 10 F 76 30 10 H 43 31 11 B 18 32 11 D 79 33 11 I 41 34 12 C 42 35 12 E 84 36 12 G 81 TukeyHSD gives> fit <- aov(Comfort~ Consumer+Pillow,data=d) > TukeyHSD(fit,"Pillow")Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = Comfort ~ Consumer + Pillow, data = d) $Pillow diff lwr upr B-A -31.00 -45.8641937 -16.1358063 C-A -15.50 -30.3641937 -0.6358063 ... But simint from the multcomp library gives> library(multcomp) > simint(Comfort~ Consumer+Pillow,data=d,whichf="Pillow",type="Tukey")Simultaneous confidence intervals: Tukey contrasts Call: simint.formula(formula = Comfort ~ Consumer + Pillow, data = d, whichf = "Pillow", type = "Tukey") 95 % confidence intervals Estimate lower CI upper CI PillowB-PillowA -41.333 -58.498 -24.168 PillowC-PillowA -20.667 -37.832 -3.502 We get the right diffs from dummy.coef> pillow <- dummy.coef(fit)$Pillow > outer(pillow,pillow,"-")A B C A 0.000000 41.33333 20.666667 ... B -41.333333 0.00000 -20.666667 ... C -20.666667 20.66667 0.000000 ... and the right se from summary.lm> summary.lm(fit)Call: aov(formula = Comfort ~ Consumer + Pillow, data = d) Residuals: Min 1Q Median 3Q Max -6.00000 -3.36111 0.05556 2.94444 9.00000 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 61.667 4.404 14.001 2.14e-10 *** ... Consumer12 1.111 5.334 0.208 0.83761 PillowB -41.333 4.825 -8.567 2.26e-07 *** PillowC -20.667 4.825 -4.284 0.00057 *** PillowD 14.333 4.825 2.971 0.00901 ** PillowE 25.333 4.825 5.251 7.92e-05 *** PillowF 9.333 4.825 1.934 0.07094 . PillowG 14.000 4.825 2.902 0.01040 * PillowH -11.333 4.825 -2.349 0.03200 * PillowI -25.667 4.825 -5.320 6.91e-05 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Residual standard error: 5.909 on 16 degrees of freedom Multiple R-Squared: 0.9669, Adjusted R-squared: 0.9275 F-statistic: 24.57 on 19 and 16 DF, p-value: 2.007e-08> qtukey(0.95,9,16)/sqrt(2)*4.825[1] 17.16474> -41.333333 + 17.16474[1] -24.16859> -41.333333 - 17.16474[1] -58.49807 ---