[R] Confused about Tukey mult. comp. after ANCOVA

[Chabot Denis]

Hi,

I am reposting this as I fear my original post (on Oct. 4th) got  
buried by all the excitement of the R 2.6 release...

I had a first occasion to try multiple comparisons (of intercepts, I  
suppose) following a significant result in an ANCOVA. As until now I  
was doing this with JMP, I compared my results and the post-hoc  
comparisons were different between R and JMP.

I chose to use an example data set from JMP because it was small, so  
I can show it here. It is not the best example for an ANCOVA because  
the factor "Drug" does not have a significant effect, but it will do.

 >drug$x
  [1] 11  8  5 14 19  6 10  6 11  3  6  6  7  8 18  8 19  8  5 15 16  
13 11  9 21 16 12
[28] 12  7 12
 >
 > drug$y
  [1]  6  0  2  8 11  4 13  1  8  0  0  2  3  1 18  4 14  9  1  9 13  
10 18  5 23 12  5
[28] 16  1 20
 > drug$Drug
  [1] a a a a a a a a a a d d d d d d d d d d f f f f f f f f f f
Levels: a d f

I did not manage to get TukeyHSD to work if I fitted the ANCOVA with  
lm, so I used aov:

my.anc <- aov(y~x+Drug, data=drug)

 > summary(my.anc)
             Df Sum Sq Mean Sq F value    Pr(>F)
x            1 802.94  802.94 50.0393 1.639e-07 ***
Drug         2  68.55   34.28  2.1361    0.1384
Residuals   26 417.20   16.05
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

I tried this to compare the Drugs, correcting for the effect of x.

 > TukeyHSD(my.anc, "Drug")
   Tukey multiple comparisons of means
     95% family-wise confidence level

Fit: aov(formula = y ~ x + Drug, data = drug)

$Drug
           diff       lwr      upr     p adj
d-a 0.03131758 -4.420216 4.482851 0.9998315
f-a 3.04677613 -1.404758 7.498310 0.2239746
f-d 3.01545855 -1.436075 7.466992 0.2305187

Warning message:
non-factors ignored: x in: replications(paste("~", xx), data = mf)

I am not sure about the Warning, maybe it is the reason the  
differences shown here are different from those shown in JMP for the  
same analysis. Maybe TukeyHSD is not meant to be used with non- 
factors (i.e. not valid for ANCOVAs)?

I just found the package multcomp and am not sure I understand it  
well yet, but its Tukey comparisons gave the same results as JMP.

 > summary(glht(m3, linfct=mcp(Drug="Tukey")))

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = y ~ x + Drug, data = drug)

Linear Hypotheses:
            Estimate Std. Error t value p value
d - a == 0    0.109      1.795   0.061   0.998
f - a == 0    3.446      1.887   1.826   0.181
f - d == 0    3.337      1.854   1.800   0.189
(Adjusted p values reported)


I would very much like to understand why these two "Tukey" tests gave  
different results in R.

Thanks in advance,

Denis