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 <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
______________________________________________
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.
