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
post-hoc comparisons following glmm
3 messages · Spencer Graves, Michaël Coeurdassier
3 days later
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 2 48 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,] -1 1 0 0 0
[2,] -1 0 1 0 0
[3,] -1 0 0 1 0
[4,] -1 0 0 0 1
[5,] 0 -1 1 0 0
[6,] 0 -1 0 1 0
[7,] 0 -1 0 0 1
[8,] 0 0 -1 1 0
[9,] 0 0 -1 0 1
[10,] 0 0 0 -1 1
> 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,] -1 1 0 0 0
[2,] -1 0 1 0 0
[3,] -1 0 0 1 0
[4,] -1 0 0 0 1
[5,] 0 -1 1 0 0
[6,] 0 -1 0 1 0
[7,] 0 -1 0 0 1
[8,] 0 0 -1 1 0
[9,] 0 0 -1 0 1
[10,] 0 0 0 -1 1
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:
multcomp mvtnorm MASS statmod 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 even if it does), PLEASE do
read the posting guide! "www.R-project.org/posting-guide.html". I'd
prefer not to have to guess whether you would think the example I chose
was relevant.
hope this helps,
spencer graves
Micha??l Coeurdassier wrote:
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 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
3 days later
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 2 48 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,] -1 1 0 0 0 [2,] -1 0 1 0 0 [3,] -1 0 0 1 0 [4,] -1 0 0 0 1 [5,] 0 -1 1 0 0 [6,] 0 -1 0 1 0 [7,] 0 -1 0 0 1 [8,] 0 0 -1 1 0 [9,] 0 0 -1 0 1 [10,] 0 0 0 -1 1
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,] -1 1 0 0 0
[2,] -1 0 1 0 0
[3,] -1 0 0 1 0
[4,] -1 0 0 0 1
[5,] 0 -1 1 0 0
[6,] 0 -1 0 1 0
[7,] 0 -1 0 0 1
[8,] 0 0 -1 1 0
[9,] 0 0 -1 0 1
[10,] 0 0 0 -1 1
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:
multcomp mvtnorm MASS statmod 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 even if it does),
PLEASE do read the posting guide!
"www.R-project.org/posting-guide.html". I'd prefer not to have to
guess whether you would think the example I chose was relevant.
hope this helps,
spencer graves
Micha??l Coeurdassier wrote:
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 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