[R] multiple comparison of interaction of ANCOVA

Jinsong Zhao jszhao at yeah.net
Sun Dec 11 13:15:12 CET 2011


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



More information about the R-help mailing list