Re: [R] post-hoc comparisons following glmm

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Sat 11 Feb 2006 - 15:05:31 EST

          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 Sat Feb 11 15:08:26 2006

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:42:27 EST