Michaël Coeurdassier
2006-Mar-22 11:04 UTC
[R] post hoc comparison following with multcomp?
Dear R community,
I would like to check differences between treatment (Trait) following a
glm. From R help list, I tried in the following way using the multcom
library. Please, is it a correct manner to do post hoc comparison with
the data below.
Thank you in advance. Sincerely
*> tabp<-read.delim("ponteM1.txt")*
*> tabp*
Trait totponte
1 T 10
2 T 11
3 T 11
4 T 9
5 T 7
6 T 7
7 T 9
8 T 12
9 T 9
10 T 10
11 M 9
12 M 10
13 M 8
14 M 8
15 M 10
16 M 8
17 M 8
18 M 9
19 M 11
20 M 7
21 F 8
22 F 8
23 F 5
24 F 9
25 F 5
26 F 7
27 F 5
28 F 7
29 F 6
30 F 3
*> attach(tabp)*
*> glm1<-glm(totponte~Trait,family=poisson)* **
*> anova(glm1,test="F")*
Model: poisson, link: log
Response: totponte
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev F Pr(>F)
NULL 29 16.3879
Trait 2 7.1770 27 9.2109 3.5885 0.02764 *
* *
*> library(multcomp)*
*> (coefglm1<-coef(glm1))*
(Intercept) TraitM TraitT
1.8405496 0.3342021 0.4107422
*> (vc.trait<-vcov(glm1))*
(Intercept) TraitM TraitT
(Intercept) 0.01587301 -0.01587301 -0.01587301
TraitM -0.01587301 0.02723665 0.01587301
TraitT -0.01587301 0.01587301 0.02639933
*> (CM<-contrMat(table(tabp$Trait),type="Tukey"))*
F M T
M-F -1 1 0
T-F -1 0 1
T-M 0 -1 1
* *
*> csimint(coefglm1,df=27,covm=vc.trait,cmatrix=CM)*
Simultaneous confidence intervals: user-defined contrasts
95 % confidence intervals
Estimate 2.5 % 97.5 %
M-F -1.506 -2.177 -0.836
T-F -1.430 -2.097 -0.763
T-M 0.077 -0.286 0.439
---------------
Micha?l COEURDASSIER, PhD
Department of Environmental Biology
UsC INRA EA3184MRT
Institute for Environmental Sciences and Technology
University of Franche-Comte
Place Leclerc
25030 Besan?on cedex
FRANCE
Tel : +33 (0)381 665 741
Fax : +33 (0)381 665 797
E-m at il: michael.coeurdassier at univ-fcomte.fr
http://lbe.univ-fcomte.fr/
Your use of the multcomp package looks plausible to me (though I
certainly did not check every detail).
Do you have any specific concerns? [Also, have you reviewed
vignette("Rmc")?] And have you considered the plotting, e.g.:
plot(csimint(coefglm1,df=27,covm=vc.trait,cmatrix=CM))
hope this helps.
spencer graves
Micha?l Coeurdassier wrote:
> Dear R community,
>
> I would like to check differences between treatment (Trait) following a
> glm. From R help list, I tried in the following way using the multcom
> library. Please, is it a correct manner to do post hoc comparison with
> the data below.
>
> Thank you in advance. Sincerely
>
> *> tabp<-read.delim("ponteM1.txt")*
> *> tabp*
> Trait totponte
> 1 T 10
> 2 T 11
> 3 T 11
> 4 T 9
> 5 T 7
> 6 T 7
> 7 T 9
> 8 T 12
> 9 T 9
> 10 T 10
> 11 M 9
> 12 M 10
> 13 M 8
> 14 M 8
> 15 M 10
> 16 M 8
> 17 M 8
> 18 M 9
> 19 M 11
> 20 M 7
> 21 F 8
> 22 F 8
> 23 F 5
> 24 F 9
> 25 F 5
> 26 F 7
> 27 F 5
> 28 F 7
> 29 F 6
> 30 F 3
>
> *> attach(tabp)*
> *> glm1<-glm(totponte~Trait,family=poisson)* **
> *> anova(glm1,test="F")*
> Model: poisson, link: log
> Response: totponte
> Terms added sequentially (first to last)
>
> Df Deviance Resid. Df Resid. Dev F Pr(>F)
> NULL 29 16.3879
> Trait 2 7.1770 27 9.2109 3.5885 0.02764 *
> * *
> *> library(multcomp)*
> *> (coefglm1<-coef(glm1))*
> (Intercept) TraitM TraitT
> 1.8405496 0.3342021 0.4107422
>
> *> (vc.trait<-vcov(glm1))*
> (Intercept) TraitM TraitT
> (Intercept) 0.01587301 -0.01587301 -0.01587301
> TraitM -0.01587301 0.02723665 0.01587301
> TraitT -0.01587301 0.01587301 0.02639933
>
> *> (CM<-contrMat(table(tabp$Trait),type="Tukey"))*
> F M T
> M-F -1 1 0
> T-F -1 0 1
> T-M 0 -1 1
> * *
> *> csimint(coefglm1,df=27,covm=vc.trait,cmatrix=CM)*
> Simultaneous confidence intervals: user-defined contrasts
> 95 % confidence intervals
>
> Estimate 2.5 % 97.5 %
> M-F -1.506 -2.177 -0.836
> T-F -1.430 -2.097 -0.763
> T-M 0.077 -0.286 0.439
>
>
> ---------------
> Micha?l COEURDASSIER, PhD
> Department of Environmental Biology
> UsC INRA EA3184MRT
> Institute for Environmental Sciences and Technology
>
> University of Franche-Comte
> Place Leclerc
> 25030 Besan?on cedex
> FRANCE
> Tel : +33 (0)381 665 741
> Fax : +33 (0)381 665 797
>
> E-m at il: michael.coeurdassier at univ-fcomte.fr
> http://lbe.univ-fcomte.fr/
>
> ______________________________________________
> R-help at 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