[R] multiple comparisons for lme using multcomp

From: MichaŽl Coeurdassier <michael.coeurdassier_at_univ-fcomte.fr>
Date: Thu 10 Mar 2005 - 00:14:46 EST


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 Received on Thu Mar 10 00:19:29 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:30:42 EST