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


