thr3ads.net - R help - [R] multiple comparison of interaction of ANCOVA [Dec 2011]

If this information is useful, please help other people find it:
Share via:

Jinsong Zhao

2011-Dec-11 12:15 UTC

[R] multiple comparison of interaction of ANCOVA

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

Bert Gunter

2011-Dec-11 14:49 UTC

head link

[R] 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

Richard M. Heiberger

2011-Dec-11 17:48 UTC

head link

[R] multiple comparison of interaction of ANCOVA

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.


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.

Rich


On Sun, Dec 11, 2011 at 7:15 AM, Jinsong Zhao <jszhao@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
>
> ______________________________**________________
> R-help@r-project.org mailing list
>
https://stat.ethz.ch/mailman/**listinfo/r-help<https://stat.ethz.ch/mailman/listinfo/r-help>
> PLEASE do read the posting guide http://www.R-project.org/**
> posting-guide.html <http://www.r-project.org/posting-guide.html>
> and provide commented, minimal, self-contained, reproducible code.
>
	[[alternative HTML version deleted]]

Reasonably Related Threads

Search for more possibly parallel threads

R help - Dec 2011 - multiple comparison of interaction of ANCOVA

[R] multiple comparison of interaction of ANCOVA

[R] multiple comparison of interaction of ANCOVA

[R] multiple comparison of interaction of ANCOVA

Reasonably Related Threads