Dear R-help list,
I would like to perform multiple comparisons for lme. Can you report to me
if my way to is correct or not? Please, note that I am not nor a
statistician nor a mathematician, so, some understandings are sometimes
quite hard for me. According to the previous helps on the topic in R-help
list May 2003 (please, see Torsten Hothorn advices) and books such as
Venables & Ripley or Pinheiro & Bates. I need multcomp library. I followed
the different examples and get these results :
In this example, response is the mass of earthworms after one month of
exposure to different concentrations of pollutants (treatment). The
experimental design was three containers per treatment with five animals in
each container (or less if mortality occurred). Containers were referred as
box and considered as the random variable.
> tab<-read.delim("example1.txt")
> tab
treatment box response
1 control a 1.8
2 control a 2.3
3 control a 1.3
4 control a 0.8
5 control a 2.0
6 control b 1.1
7 control b 2.2
8 control b 1.3
9 control b 2.0
10 control b 1.3
11 control c 1.5
12 control c 1.4
13 control c 2.1
14 control c 1.7
15 control c 1.3
16 Al100 d 1.6
17 Al100 d 2.1
18 Al100 d 0.7
19 Al100 d 1.8
20 Al100 d 1.2
21 Al100 e 1.5
22 Al100 e 1.5
23 Al100 e 2.0
24 Al100 e 1.0
25 Al100 e 1.6
26 Al100 f 0.9
27 Al100 f 2.0
28 Al100 f 1.9
29 Al100 f 1.7
30 Al100 f 1.7
…
68 Al800 q 1.0
69 Al800 r 0.8
70 Al800 r 0.9
71 Al800 r 0.9
72 Al800 r 0.6
73 Al800 s 0.9
74 Al800 s 1.0
75 Al800 s 0.8
76 Al800 s 0.8
77 Al800 s 0.7
> attach(tab)
> library(nlme)
> lm1<-lme(response~treatment,random=~1|box)
> library(multcomp)
Loading required package: mvtnorm
> # first way to do (seem uncorrect)
> summary(csimtest(coef(lm1),vcov(lm1),cmatrix=contrMat(table(treatment),
type="Tukey"),df=59))
Error in csimtest(coef(lm1), vcov(lm1), cmatrix =
contrMat(table(treatment), : estpar not a vector
> #indeed
> coef(lm1)
(Intercept) treatmentAl200 treatmentAl400 treatmentAl600 treatmentAl800
a 1.546679 -0.1648540 -0.4895219 -0.6383375 -0.7066632
b 1.546657 -0.1648540 -0.4895219 -0.6383375 -0.7066632
c 1.546664 -0.1648540 -0.4895219 -0.6383375 -0.7066632
d 1.546643 -0.1648540 -0.4895219 -0.6383375 -0.7066632
…
s 1.546667 -0.1648540 -0.4895219 -0.6383375 -0.7066632
treatmentcontrol
a 0.06
b 0.06
c 0.06
d 0.06
…
s 0.06
> # second way to do could be to get a vector for lm1 coefficient removing
intercept(is it correct?)
> vect<-as.numeric(coef(lm1)[1,])
> vect
[1] 1.5466787 -0.1648540 -0.4895219 -0.6383375 -0.7066632 0.0600000
>
summary(csimtest(vect,vcov(lm1),cmatrix=contrMat(table(treatment),type="Tukey"),df=59))
Simultaneous tests: user-defined contrasts
user-defined contrasts for factor
Contrast matrix:
Al100 Al200 Al400 Al600 Al800 control
Al200-Al100 -1 1 0 0 0 0
Al400-Al100 -1 0 1 0 0 0
Al600-Al100 -1 0 0 1 0 0
Al800-Al100 -1 0 0 0 1 0
control-Al100 -1 0 0 0 0 1
Al400-Al200 0 -1 1 0 0 0
Al600-Al200 0 -1 0 1 0 0
Al800-Al200 0 -1 0 0 1 0
control-Al200 0 -1 0 0 0 1
Al600-Al400 0 0 -1 1 0 0
Al800-Al400 0 0 -1 0 1 0
control-Al400 0 0 -1 0 0 1
Al800-Al600 0 0 0 -1 1 0
control-Al600 0 0 0 -1 0 1
control-Al800 0 0 0 0 -1 1
Absolute Error Tolerance: 0.001
Coefficients:
Estimate t value Std.Err. p raw p Bonf p adj
Al800-Al100 -2.253 -10.467 0.213 0.000 0.000 0.000
Al600-Al100 -2.185 -10.389 0.207 0.000 0.000 0.000
Al400-Al100 -2.036 -9.850 0.210 0.000 0.000 0.000
Al200-Al100 -1.712 -8.051 0.215 0.000 0.000 0.000
control-Al100 -1.487 -7.243 0.205 0.000 0.000 0.000
control-Al800 0.767 -5.282 0.143 0.000 0.000 0.000
control-Al600 0.698 -5.072 0.148 0.000 0.000 0.000
control-Al400 0.550 -4.160 0.155 0.000 0.001 0.001
Al800-Al200 -0.542 -3.488 0.141 0.001 0.006 0.006
Al600-Al200 -0.473 -3.191 0.140 0.002 0.014 0.012
Al400-Al200 -0.325 -2.267 0.147 0.027 0.135 0.110
control-Al200 0.225 -1.593 0.132 0.116 0.466 0.341
Al800-Al400 -0.217 -1.475 0.152 0.145 0.466 0.341
Al600-Al400 -0.149 -1.064 0.138 0.292 0.583 0.466
Al800-Al600 -0.068 -0.449 0.145 0.655 0.655 0.655
> #a friend told me that it is possible to do multiple comparisons for lme
in a simplest way, i.e. :
> anova(lm1,L=c("treatmentcontrol"=1,"treatmentAl200"=-1))
F-test for linear combination(s)
treatmentAl200 treatmentcontrol
-1 1
numDF denDF F-value p-value
1 1 12 2.538813 0.1371
> anova(lm1,L=c("treatmentcontrol"=1,"treatmentAl400"=-1))
F-test for linear combination(s)
treatmentAl400 treatmentcontrol
-1 1
numDF denDF F-value p-value
1 1 12 17.30181 0.0013
> anova(lm1,L=c("treatmentcontrol"=1,"treatmentAl600"=-1))
F-test for linear combination(s)
treatmentAl600 treatmentcontrol
-1 1
numDF denDF F-value p-value
1 1 12 25.72466 3e-04
> anova(lm1,L=c("treatmentcontrol"=1,"treatmentAl800"=-1))
F-test for linear combination(s)
treatmentAl800 treatmentcontrol
-1 1
numDF denDF F-value p-value
1 1 12 27.9043 2e-04
> # however, p values are different that those obtained above. Is this way
OK or not?
Thank you in advance.
Sincerely
Michael Coeurdassier
[[alternative HTML version deleted]]
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
