Re: [R] post-hoc comparisons following glmm

From: MichaŽl Coeurdassier <michael.coeurdassier_at_univ-fcomte.fr>
Date: Wed 15 Feb 2006 - 00:52:56 EST

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@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-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 Wed Feb 15 01:00:22 2006

This archive was generated by hypermail 2.1.8 : Wed 15 Feb 2006 - 14:35:08 EST