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 <laomen...@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]]
______________________________________________
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.
