[R] multiple comparison of interaction of ANCOVA

[Jinsong Zhao]

On 2011-12-12 1:48, Richard M. Heiberger wrote:
> Thank you for you use of HH.
> I think the right graph for this data is the much simpler ancova function
> library(HH)
>
> ancova(y ~ year * Trt, data=mydata)
> where we see that the three treatments have totally different slopes.

Thank you very much for the advice. The plots seem to be that I hope to get.

> The WoodEnergy example doesn't apply here.  The WoodEnergy example
> illustrates
> a way of finding differences among treatments for a fixed value of the
> covariate when the
> slopes are similar.

Apart from the data set here, is there a way to do multiple comparison 
on the interaction of one way analysis of covariance?

BTW, when loading HH package, the TukeyHSD() from stats package can not 
be use. The following is the output:

 > mydata.aov <- aov(y~Trt, mydata)
 > TukeyHSD(mydata.aov)
   Tukey multiple comparisons of means
     95% family-wise confidence level

Fit: aov(formula = y ~ Trt, data = mydata)

$Trt
       diff       lwr        upr     p adj
B-A  5.004  3.275115  6.7328854 0.0000149
C-A -2.320 -4.048885 -0.5911146 0.0097872
C-B -7.324 -9.052885 -5.5951146 0.0000003

 > library(HH)
Loading required package: lattice
Loading required package: grid
Loading required package: multcomp
Loading required package: mvtnorm
Loading required package: survival
Loading required package: splines
Loading required package: leaps
Loading required package: RColorBrewer
Loading required package: latticeExtra
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

After detach(package:HH), TukeyHSD(mydata.aov) can output the correct 
results.

Thanks again.

Regards,
Jinsong



> Rich
> On Sun, Dec 11, 2011 at 7:15 AM, Jinsong Zhao <jszhao at yeah.net
> <mailto: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!
>
>     Regards,
>     Jinsong