[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 T9
5 T7
6 T7
7 T9
8 T 12
9 T9
10 T 10
11 M9
12 M 10
13 M8
14 M8
15 M 10
16 M8
17 M8
18 M9
19 M 11
20 M7
21 F8
22 F8
23 F5
24 F9
25 F5
26 F7
27 F5
28 F7
29 F6
30 F3
*> 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 2916.3879
Trait 2 7.177027 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-M0.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
[EMAIL PROTECTED]: [EMAIL PROTECTED]
http://lbe.univ-fcomte.fr/
__
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
Re: [R] post-hoc comparisons following glmm
Dear Mr Graves,
this seems work on my data. Thank you very much for your help.
Sincerely
Michael Coeurdassier
Spencer Graves a écrit :
> The following appears to be an answer to your question, though
> I'd be pleased to receive critiques from others. Since your example
> is NOT self contained, I modified an example in the "glmmPQL" help file:
>
> (fit <- glmmPQL(y ~ factor(week)-1+trt, random = ~ 1 | ID,
> + family = binomial, data = bacteria))
> iteration 1
> iteration 2
> iteration 3
> iteration 4
> iteration 5
> iteration 6
> Linear mixed-effects model fit by maximum likelihood
> Data: bacteria
> Log-likelihood: -551.1184
> Fixed: y ~ factor(week) - 1 + trt
> factor(week)0 factor(week)2 factor(week)4 factor(week)6
> factor(week)11
> 3.3459650 3.5262521 1.9102037 1.7645881
> 1.7660845
>trtdrug trtdrug+
> -1.2527642 -0.7570441
>
> Random effects:
> Formula: ~1 | ID
> (Intercept) Residual
> StdDev:1.426534 0.7747477
>
> Variance function:
> Structure: fixed weights
> Formula: ~invwt
> Number of Observations: 220
> Number of Groups: 50
> > anova(fit)
> numDF denDF F-value p-value
> factor(week) 5 166 10.821682 <.0001
> trt 248 1.889473 0.1622
> > (denDF.week <- anova(fit)$denDF[1])
> [1] 166
> > (denDF.week <- anova(fit)$denDF[1])
> [1] 166
> > (par.week <- fixef(fit)[1:5])
> factor(week)0 factor(week)2 factor(week)4 factor(week)6
> factor(week)11
> 3.345965 3.526252 1.910204 1.764588
> 1.766085
> > (vc.week <- vcov(fit)[1:5, 1:5])
>factor(week)0 factor(week)2 factor(week)4 factor(week)6
> factor(week)0 0.3351649 0.1799365 0.1705898 0.1694884
> factor(week)2 0.1799365 0.3709887 0.1683038 0.1684096
> factor(week)4 0.1705898 0.1683038 0.2655072 0.1655673
> factor(week)6 0.1694884 0.1684096 0.1655673 0.2674647
> factor(week)11 0.1668450 0.1665177 0.1616748 0.1638169
>factor(week)11
> factor(week)0 0.1668450
> factor(week)2 0.1665177
> factor(week)4 0.1616748
> factor(week)6 0.1638169
> factor(week)11 0.2525962
> > CM <- array(0, dim=c(5*4/2, 5))
> > i1 <- 0
> > for(i in 1:4)for(j in (i+1):5){
> + i1 <- i1+1
> + CM[i1, c(i, j)] <- c(-1, 1)
> + }
> > CM
> [,1] [,2] [,3] [,4] [,5]
> [1,] -11000
> [2,] -10100
> [3,] -10010
> [4,] -10001
> [5,]0 -1100
> [6,]0 -1010
> [7,]0 -1001
> [8,]00 -110
> [9,]00 -101
> [10,]000 -11
> > library(multcomp)
> > csimint(par.week, df=denDF.week, covm=vc.week,cmatrix=CM)
>
> Simultaneous confidence intervals: user-defined contrasts
>
> 95 % confidence intervals
>
> Estimate 2.5 % 97.5 %
> [1,]0.180 -1.439 1.800
> [2,] -1.436 -2.838 -0.034
> [3,] -1.581 -2.995 -0.168
> [4,] -1.580 -2.967 -0.193
> [5,] -1.616 -3.123 -0.109
> [6,] -1.762 -3.273 -0.250
> [7,] -1.760 -3.244 -0.277
> [8,] -0.146 -1.382 1.091
> [9,] -0.144 -1.359 1.070
> [10,]0.001 -1.206 1.209
>
> > csimtest(par.week, df=denDF.week, covm=vc.week,cmatrix=CM)
>
> Simultaneous tests: user-defined contrasts
>
> Contrast matrix:
> [,1] [,2] [,3] [,4] [,5]
> [1,] -11000
> [2,] -10100
> [3,] -10010
> [4,] -10001
> [5,]0 -1100
> [6,]0 -1010
> [7,]0 -1001
> [8,]00 -110
> [9,]00 -101
> [10,]000 -11
>
> Adjusted P-Values
>
> p adj
> [1,] 0.011
> [2,] 0.013
> [3,] 0.014
> [4,] 0.015
> [5,] 0.020
> [6,] 0.024
> [7,] 0.985
> [8,] 0.985
> [9,] 0.985
> [10,] 0.997
> > sessionInfo()
> R version 2.2.1, 2005-12-20, i386-pc-mingw32
>
> attached base packages:
> [1] "methods" "stats" "graphics" "grDevices" "utils"
> "datasets"
> [7] "base"
>
> other attached packages:
> multcompmvtnorm MASSstatmod nlme
>"0.4-8""0.7-2" "7.2-24""1.2.4" "3.1-68.1"
>
> If this does NOT answer your question (or eve
[R] post-hoc comparisons following glmm
Dear R community, I performed a generalized linear mixed model using glmmPQL (MASS library) to analyse my data i.e : y is the response with a poisson distribution, t and Trait are the independent variables which are continuous and categorical (3 categories C, M and F) respectively, ind is the random variable. mydata<-glmmPQL(y~t+Trait,random=~1|ind,family=poisson,data=tab) Do you think it is OK? Trait is significant (p < 0.0001) and I would like to perform post-hoc comparisons to check where the difference among Trait categories but I did not find a solution in R help list or others. Thank you in advance for your help Michael __ 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
[R] multiple comparisons for lme using multcomp
0.027 0.135 0.110
control-Al2000.225 -1.5930.132 0.116 0.466 0.341
Al800-Al400 -0.217 -1.4750.152 0.145 0.466 0.341
Al600-Al400 -0.149 -1.0640.138 0.292 0.583 0.466
Al800-Al600 -0.068 -0.4490.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
-11
numDF denDF F-value p-value
1 112 2.538813 0.1371
> anova(lm1,L=c("treatmentcontrol"=1,"treatmentAl400"=-1))
F-test for linear combination(s)
treatmentAl400 treatmentcontrol
-11
numDF denDF F-value p-value
1 112 17.30181 0.0013
> anova(lm1,L=c("treatmentcontrol"=1,"treatmentAl600"=-1))
F-test for linear combination(s)
treatmentAl600 treatmentcontrol
-11
numDF denDF F-value p-value
1 112 25.72466 3e-04
> anova(lm1,L=c("treatmentcontrol"=1,"treatmentAl800"=-1))
F-test for linear combination(s)
treatmentAl800 treatmentcontrol
-11
numDF denDF F-value p-value
1 112 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
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
[EMAIL PROTECTED]: [EMAIL PROTECTED]
http://lbe.univ-fcomte.fr/
[[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
[R] multiple comparisons for lme using multcomp
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
1control a 1.8
2control a 2.3
3control a 1.3
4control a 0.8
5control a 2.0
6control b 1.1
7control b 2.2
8control b 1.3
9control 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
a1.546679 -0.1648540 -0.4895219 -0.6383375 -0.7066632
b1.546657 -0.1648540 -0.4895219 -0.6383375 -0.7066632
c1.546664 -0.1648540 -0.4895219 -0.6383375 -0.7066632
d1.546643 -0.1648540 -0.4895219 -0.6383375 -0.7066632
s1.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.060
>
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.4670.213 0.000 0.000 0.000
Al600-Al100 -2.185 -10.3890.207 0.000 0.000 0.000
Al400-Al100 -2.036 -9.8500.210 0.000 0.000 0.000
Al200-Al100 -1.712 -8.0510.215 0.000 0.000 0.000
control-Al100 -1.487 -7.2430.205 0.000 0.000 0.000
control-Al8000.767 -5.2820.143 0.000 0.000 0.000
control-Al6000.698 -5.0720.148 0.000 0.000 0.000
control-Al4000.550 -4.1600.155 0.000 0.001 0.001
Al800-Al200 -0.542 -3.4880.141 0.001 0.006 0.006
Al600-Al200 -0.473 -3.1910.140 0.002 0.014 0.012
Al400-Al200 -0.325 -2.2670.147 0.027 0.135 0.110
control-Al2000.225 -1.5930.132 0.116 0.466 0.341
Al800
[R] bootstrap function coefficients
Dear R community,
Please, can you help me with a problem concerning bootstrap. The data table
called «RMika», contained times (Tps) and corresponding concentration of a
chemical in a soil (SolA). I would like to get, by bootstraping, 10
estimations of the parameters C0 and k from the function: SolA =
C0*exp(-k*Tps).
# First, I fit the data and all is OK
> tabMika<-read.delim("RMika.txt")
> library(nls)
> attach(tabMika)
> Expon<-function(Tps,parm){
+ C0<-parm[1]
+ k<-parm[2]
+ }
> DegSA.nls<-nls(SolA~C0*exp(-k*Tps),start=c(C0=35, k=1),tabMika)
> summary(DegSA.nls)
Formula: SolA ~ C0 * exp(-k * Tps)
Parameters:
Estimate Std. Error t value Pr(>|t|)
C0 25.682104 1.092113 23.52 < 2e-16 ***
k 0.087356 0.007582 11.52 6.36e-13 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 2.584 on 32 degrees of freedom
Correlation of Parameter Estimates:
C0
k 0.7101
# Second, I try to use bootstrap to get the 10 estimates of k and C0
> library(bootstrap)
> theta<-function(tabMika){coef(eval(DegSA.nls$call))}
> bootSolA.nls<-bootstrap(tabMika,10,theta)
Warning message:
multi-argument returns are deprecated in: return(thetastar, func.thetastar,
jack.boot.val, jack.boot.se,
> bootSolA.nls
$thetastar
[,1][,2][,3][,4][,5][,6]
C0 25.68210358 25.68210358 25.68210358 25.68210358 25.68210358 25.68210358
k 0.08735615 0.08735615 0.08735615 0.08735615 0.08735615 0.08735615
[,7][,8][,9] [,10]
C0 25.68210358 25.68210358 25.68210358 25.68210358
k 0.08735615 0.08735615 0.08735615 0.08735615
# As you can notify, the 10 estimations have the same values for C0 and k.
Moreover, this correspond to the values of k and C0 determined by fitting
all the data without bootstrap!!!???
I cannot find what is wrong. Please, if you find a solution, thank you to
send it to me.
Sincerely
Michael Coeurdassier
Michaël COEURDASSIER, PhD
Department of Environmental Biology
UC 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
[EMAIL PROTECTED]: [EMAIL PROTECTED]
http://lbe.univ-fcomte.fr/
[[alternative HTML version deleted]]
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Basic questions on nls and bootstrap
Dear R community,
I have currently some problems with non linear regression analysis in R.
My data correspond to the degradation kinetic of a pollutant in two
different soil A and B, x data are time in day and y data are pollutant
concentration in soil.
In a first time, I want to fit the data for the soil A by using the Ct =
C0*exp(-k*Tpst) with Ct the concentration of pollutant at time t, C0 is the
initial concentration (i.e. t=0), Tps is the time in days.
The table containing data is called «tabMika»
Here, you can find my R program:
> tabMika<-read.delim("RMika.txt")
> tabMika
Tps SolA Solb
10 32.97 35.92
20 32.01 31.35
31 21.73 22.03
41 23.73 18.53
52 19.68 18.28
62 18.56 16.79
# and continue like that until 29 days.
> library(nls)
> attach(tabMika)
> Expon<-function(Tps,parm){
+ C0<-parm[1]
+ k<-parm[2]
+ }
> DegSA.nls<-nls(SolA~C0*exp(-k*Tps),start=c(C0=35, k=1),tabMika)
> summary(DegSA.nls)
Formula: SolA ~ C0 * exp(-k * Tps)
Parameters:
Estimate Std. Error t value Pr(>|t|)
C0 25.682104 1.092113 23.52 < 2e-16 ***
k 0.087356 0.007582 11.52 6.36e-13 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 2.584 on 32 degrees of freedom
Correlation of Parameter Estimates:
C0
k 0.7101
# until this step, I think (or may be hope) that it is OK.
#Then, I would like to use bootstrap function to obtain 10 estimations
(more of course but 10 is retained for this example) of the parameters C0
and t. I use this way:
> library(bootstrap)
> theta<-function(tabMika){coef(eval(DegSA.nls$call))}
> bootSolA.nls<-bootstrap(tabMika,10,theta)
Warning message:
multi-argument returns are deprecated in: return(thetastar, func.thetastar,
jack.boot.val, jack.boot.se,
> bootSolA.nls
$thetastar
[,1][,2][,3][,4][,5][,6]
C0 25.68210358 25.68210358 25.68210358 25.68210358 25.68210358 25.68210358
k 0.08735615 0.08735615 0.08735615 0.08735615 0.08735615 0.08735615
[,7][,8][,9] [,10]
C0 25.68210358 25.68210358 25.68210358 25.68210358
k 0.08735615 0.08735615 0.08735615 0.08735615
$func.thetastar
NULL
$jack.boot.val
NULL
$jack.boot.se
NULL
$call
bootstrap(x = tabMika, nboot = 10, theta = theta)
# as you can notify, the 10 estimations of C0 and k are the same value
(!!), so, my program is wrong but I cannot find where is the problem.
# In a second time, I would like to fit the same data with a new function
Ct = C0*Tpsb. I did as follow and obtained the following results:
> Pat<-function(Tps,parm){
+ a<-parm[1]
+ b<-parm[2]
+ }
> PatSA.nls<-nls(SolA~a*(Tps^b),start=c(a=35, b=-1),tabMika)
Error in numericDeriv(form[[3]], names(ind), env) : Missing value or an
Infinity produced when evaluating the model
> summary(PatSA.nls)
Error in summary(PatSA.nls) : Object "PatSA.nls" not found
# I cannot understand the mean of «Error in numericDeriv(form[[3]],
names(ind), env) : Missing value or an Infinity produced when evaluating
the model», please, can you help me.
Thank you very much in advance.
Michael C
Michaël COEURDASSIER, PhD
Department of Environmental Biology
UC 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
[EMAIL PROTECTED]: [EMAIL PROTECTED]
http://lbe.univ-fcomte.fr/
[[alternative HTML version deleted]]
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
