Many thanks
 





 


At 2013-03-14 00:23:43,"Richard M. Heiberger" <rmh@temple.edu> wrote:
Meng,


What seems to be going on is that the covariates are handled very differently in TukeyHSD and in glht.


Please see the interaction_average and covariate_average arguments to glht.


I ran your example twice, first as you did, with the covariates after the factor.
x2 is not significant if the factor comes first.


The second time I placed the covariates before the factor.




> result_aov <- aov(y ~ method + x1 + x2, data=meng)
> anova(result_aov)
Analysis of Variance Table


Response: y
          Df Sum Sq Mean Sq F value    Pr(>F)    
method     2 4.4705 2.23525 41.3822 1.170e-08 ***
x1         1 2.8352 2.83519 52.4892 1.363e-07 ***
x2         1 0.0747 0.07469  1.3827    0.2507    
Residuals 25 1.3504 0.05401                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> result2_aov <- aov(y ~ x1 + x2 + method, data=meng)
> anova(result2_aov)
Analysis of Variance Table


Response: y
          Df Sum Sq Mean Sq  F value    Pr(>F)    
x1         1 5.4399  5.4399 100.7113 2.985e-10 ***
x2         1 0.3134  0.3134   5.8017   0.02371 *  
method     2 1.6271  0.8135  15.0616 5.100e-05 ***
Residuals 25 1.3504  0.0540                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 




Looking at the plot of your data shows a very interesting pattern.


> xyplot(y ~ x1 + x2 | method, outer=TRUE, data=meng)


b and c both give higher values of y than does a.  You will also see
that in the table of means.


Although I leave the interpretation of x1 and x2 to you, my inclination is
to drop x2 and look at the ancova of method and x1.


ancova is in the HH package.
## install.packages("HH")  ## if you don't have it yet.
library(HH)


> result3_aov <- ancova(y ~ method + x1, data=meng)
> result3_aov
Analysis of Variance Table


Response: y
          Df Sum Sq Mean Sq F value    Pr(>F)    
method     2 4.4705 2.23525  40.782 9.616e-09 ***
x1         1 2.8352 2.83519  51.728 1.226e-07 ***
Residuals 26 1.4251 0.05481                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> update(attr(result3_aov, "trellis"), ylim=c(3.5, 6.7))
> 


Now it looks like conditional on x1, all three methods differ.

Rich


On Wed, Mar 13, 2013 at 7:15 AM, meng <laomeng_3@163.com> wrote:
Hi all:
I have a question about multi-comparison.

The data is in the attachment.

My purpose:
Compare the predicted means of the 3 methods(a,b,c) pairwisely.

I have 3 ideas:

#idea1
result_aov<-aov(y~ method + x1 + x2)
TukeyHSD(result_aov)
     diff        lwr       upr     p adj
b-a  0.845  0.5861098 1.1038902 0.0000001
c-a  0.790  0.5311098 1.0488902 0.0000002
c-b -0.055 -0.3138902 0.2038902 0.8578386

#idea2
library(multcomp)
summary(glht(result_aov,linfct=mcp(method="Tukey")))
         Estimate Std. Error t value Pr(>|t|)
b - a == 0   0.3239     0.1402   2.309   0.0683 .
c - a == 0  -0.3332     0.1937  -1.720   0.2069
c - b == 0  -0.6570     0.1325  -4.960   <0.001 ***

#idea3
#ref=a
dat$method <- relevel(dat$method, ref="a")
lm_ref_a<-lm(y~method + x1 + x2)
summary(lm_ref_a)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.92202    0.64418   1.431   0.1647
methodb      0.32389    0.14025   2.309   0.0295 *
methodc     -0.33316    0.19372  -1.720   0.0978 .
x1           0.57935    0.09356   6.192 1.78e-06 ***
x2           0.13596    0.11563   1.176   0.2507

#ref=b
dat$method <- relevel(dat$method, ref="b")
lm_ref_b<-lm(y~method + x1 + x2)
summary(lm_ref_b)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.24591    0.73770   1.689   0.1037
methoda     -0.32389    0.14025  -2.309   0.0295 *
methodc     -0.65705    0.13248  -4.960 4.14e-05 ***


In summary:
idea1:
a vs b:pvalue=0.0000001
a vs c:pvalue=0.0000002
b vs c:pvalue=0.8578386
idea2:
a vs b:pvalue=0.0683
a vs c:pvalue=0.2069
b vs c:pvalue<0.001
idea3:
a vs b:pvalue=0.0295
a vs c:pvalue=0.0978
b vs c:pvalue=4.14e-05

So the result of 3 ideas are different,and I don't know which one is correct.
Many thanks for your help.

My best
 
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.




	[[alternative HTML version deleted]]

