...
On Mon, Dec 12, 2011 at 7:00 AM, Jinsong Zhao <jsz...@yeah.net> wrote:
> Apart from the data set here, is there a way to do multiple comparison on
> the interaction of one way analysis of covariance?
>
Interaction in ancova usually means the slopes are different.
...
BTW, when loading HH package, the TukeyHSD() from stats package can not be
use. The following is the output:
> library(HH)
Warning messages:
1: package 'HH' was built under R version 2.13.2
2: package 'mvtnorm' was built under R version 2.13.2
3: package 'leaps' was built under R version 2.13.2
4: package 'RColorBrewer' was built under R version 2.13.2
5: package 'latticeExtra' was built under R version 2.13.2
> TukeyHSD(mydata.aov)
height
All those warnings say your version of R is very out of date. R-2.14.0 is
current.
R-2.14.1 is scheduled to be released at the end of December.
...
After detach(package:HH), TukeyHSD(mydata.aov) can output the correct
results.
Thank you for alterting me to the conflict.
The problem arises because the returned object from TukeyHSD has class
[1] "multicomp" "TukeyHSD"
and the print.multicomp method in HH picks it up. The stats package itself
has no "multicomp" methods. The workaround to get the printed TukeyHSD
output when HH is loaded is to use the long statements
stats:::print.TukeyHSD(TukeyHSD(mydata.aov))
plot.TukeyHSD(TukeyHSD(mydata.aov))
I will detect that conflict in the next release of HH and prevent it from
happening.
A better solution is to use glht in the multcomp package instead of
TukeyHSD.
Then the statement is
confint(glht(mydata.aov, linfct=mcp(Trt="Tukey")))
glht is better because it does not have the limitations of TukeyHSD.
TukeyHSD is not applicable in your example because TukeyHSD doesn't work
with covariates.
The documentation ?TukeyHSD says "non-factor terms these will be dropped
with a warning".
With glht and the HH package loaded, you can get a better plot.
glht.mmc(mydata.aov, linfct=mcp(Trt="Tukey"))
plot(glht.mmc(mydata.aov, linfct=mcp(Trt="Tukey")))
With the MMC plot ?HH:::MMC, "Each contrast is plotted at a height which is
the weighted
average of the means being compared"
Thanks again.
Regards,
Jinsong
Rich
>> On Sun, Dec 11, 2011 at 7:15 AM, Jinsong Zhao <jsz...@yeah.net
>> <mailto:jsz...@yeah.net>> wrote:
>>
>> Hi there,
>>
>> The following data is obtained from a long-term experiments.
>>
>> > mydata <- read.table(textConnection("
>> + y year Trt
>> + 9.37 1993 A
>> + 8.21 1995 A
>> + 8.11 1999 A
>> + 7.22 2007 A
>> + 7.81 2010 A
>> + 10.85 1993 B
>> + 12.83 1995 B
>> + 13.21 1999 B
>> + 13.70 2007 B
>> + 15.15 2010 B
>> + 5.69 1993 C
>> + 5.76 1995 C
>> + 6.39 1999 C
>> + 5.73 2007 C
>> + 5.55 2010 C"), header = TRUE)
>> > closeAllConnections()
>>
>> The experiments is designed without replication, thus I have to use
>> ANCOVA or linear mixed effect model to analyze the data. In the
>> model, variable year is coded as a continuous variable, and Trt is
>> factor variable.
>>
>> > mydata.aov <- aov(y~Trt*year, mydata)
>> > anova(mydata.aov)
>> Analysis of Variance Table
>>
>> Response: y
>> Df Sum Sq Mean Sq F value Pr(>F)
>> Trt 2 140.106 70.053 197.9581 3.639e-08 ***
>> year 1 0.610 0.610 1.7246 0.221600
>> Trt:year 2 8.804 4.402 12.4387 0.002567 **
>> Residuals 9 3.185 0.354
>> ---
>> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>>
>> As you have seen, the interaction effect is significant. I hope to
>> use TukeyHSD() or glht() to do multiple comparison on Trt:year.
>> However, for variable year is not a factor, they all give error
>> messages.
>>
>> I try to follow the demo("MMC.WoodEnergy") in HH package, as follwoing:
>>
>> > library(HH)
>> > mca.1993 <- mcalinfct(mydata.aov, "Trt")
>> > non.zero <- mca.1993[,5:6] != 0
>> > mca.1993[,5:6][non.zero] <- 1993 * sign(mca.1993[,5:6][non.zero])
>> > summary(glht(mydata.aov, linfct=mca.1993))
>>
>> Simultaneous Tests for General Linear Hypotheses
>>
>> Fit: aov(formula = y ~ Trt * year, data = mydata)
>>
>> Linear Hypotheses:
>> Estimate Std. Error t value Pr(>|t|)
>> B - A == 0 2.8779 0.5801 4.961 0.00215 **
>> C - A == 0 -2.8845 0.5801 -4.972 0.00191 **
>> C - B == 0 -5.7624 0.5801 -9.933 < 0.001 ***
>> ---
>> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>> (Adjusted p values reported -- single-step method)
>>
>> It can give comparison between levels of Trt within one year, e.g.,
>> 1993.
>>
>> Is it possible to do multiple comparison for the following pairs:
>>
>> A.1995 - A.1993
>> B.1995 - A.1993
>> C.1995 - A.1993
>>
>> Any suggestions or comments will be really appreciated. Thanks in
>> advance!
>>
>> Regards,
>> Jinsong
>>
>
[[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.
