[R] TukeyHSD and BIBD

[Simon Wotherspoon]

Sorry,
	Forgot to include the example

> d
   Consumer 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
---