multiple comparison of interaction of ANCOVA
Inline below. -- Bert
On Sun, Dec 11, 2011 at 4:15 AM, Jinsong Zhao <jszhao at 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!
Graph the data sensibly to figure out what's going on. Statistical machinationsand anova tables with P values alone are not sufficient and can be opaque or misleading. If you do not know what "sensibly" is (or even if you do), consult: http://addictedtor.free.fr/graphiques/ Cheers, Bert
Regards, Jinsong
______________________________________________ R-help at 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.
Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm