From: Michaël Coeurdassier <michael.coeurdassier_at_univ-fcomte.fr>

Date: Wed 15 Feb 2006 - 00:52:56 EST

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

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
*