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