Many thanks
At 2013-03-14 00:23:43,"Richard M. Heiberger" <r...@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 <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.
