[R] post-hoc comparisons following glmm

[Michaël Coeurdassier]

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
>