Greetings R-ers, A colleague and I have been exploring the behaviour of glmmPQL in R and S-PLUS 6 and we appear to get different results using the same code and the same data set, which worries us. I have checked the behaviour in R 1.7.1 (MacOS 9.2) and R. 1.9.0 (Windows 2000) and the results are the same, but differ from S-PLUS 6 with the latest Mass and nlme libraries (Windows XP). Here is the R output:> fit <- glmmPQL(cbind(alive, setup-alive)~ inout, random=~1|family, >family=binomial, data=dat)Loading required package: nlme iteration 1> fit2 <- glmmPQL(cbind(alive, setup-alive)~ 1, random=~1|family, >family=binomial, data=dat)iteration 1 iteration 2> anova(fit)numDF denDF F-value p-value (Intercept) 1 17 4.711629 0.0444 inout 1 17 1.260877 0.2771> anova(fit2)numDF denDF F-value p-value (Intercept) 1 18 4.800814 0.0418> anova(fit, fit2)Model df AIC BIC logLik Test L.Ratio p-value fit 1 4 87.79085 94.12493 -39.89543 fit2 2 3 86.41927 91.16983 -40.20964 1 vs 2 0.628421 0.4279>Here is the S-PLUS output:> m1 <- glmmPQL(cbind(alive, setup - alive)~ inout, data=dat, >random=~1|family, family = binomial)iteration 1 iteration 2 iteration 3> anova(m1)numDF denDF F-value p-value (Intercept) 1 17 3.859334 0.0660 inout 1 17 1.971763 0.1783> m2 <- glmmPQL(cbind(alive, setup - alive) ~ 1, data = dat, random = >~ 1 | family, family = binomial)iteration 1 iteration 2 iteration 3 iteration 4> anova(m1, m2)Model df AIC BIC logLik Test L.Ratio p-value m1 1 4 87.00400 93.33808 -39.50200 m2 2 3 91.75447 96.50503 -42.87724 1 vs 2 6.75047 0.0094 Note that R and S-PLUS differ in the number of iterations. Also the logLikelihoods differ considerably too. Pinheiro and Bates argue that a likelihood-ratio test for fixed effects is not reliable ("anticonservative"), but I think that both packages should at least give the same answer!! I checked lmeControl() in R and S-PLUS and the settings for lme look the same on both platforms. Which output (if any) should we believe? Any insight would be greatly appreciated. Thanks in advance, Simon Blomberg. School of Botany and Zoology The Australian National University Canberra, Australia. [[alternative HTML version deleted]]
glmmPQL calls lme, and it is that which differs between S-PLUS and R. These are optimization problems with multiple local maxima, and like any complex statistical fitting problem you should not expect all programs to give the same answer. The short answer is to believe all of them, and to check the usefulnees of the answer for your criteria. On Thu, 19 Aug 2004, Simon wrote:> Greetings R-ers, > > A colleague and I have been exploring the behaviour of glmmPQL in R > and S-PLUS 6 and we appear to get different results using the same > code and the same data set, which worries us. I have checked the > behaviour in R 1.7.1 (MacOS 9.2) and R. 1.9.0 (Windows 2000) and the > results are the same, but differ from S-PLUS 6 with the latest Mass > and nlme libraries (Windows XP)....> Note that R and S-PLUS differ in the number of iterations. Also the > logLikelihoods differ considerably too. Pinheiro and Bates argue that > a likelihood-ratio test for fixed effects is not reliable > ("anticonservative"), but I think that both packages should at least > give the same answer!! I checked lmeControl() in R and S-PLUS and the > settings for lme look the same on both platforms. Which output (if > any) should we believe? Any insight would be greatly appreciated.-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595
Professor Ripley,> These are optimization problems with multiple local maxima, and > like any complex statistical fitting problem you should not expect all > programs to give the same answer.Examining the log-likelihoods from the anova statements, while fit1 and m1 are similar, the maximum for fit2 (due to R) is higher than the maximum for m2 (due to S-plus). R:> anova(fit, fit2)Model df AIC BIC logLik Test L.Ratio p-value fit 1 4 87.79085 94.12493 -39.89543 fit2 2 3 86.41927 91.16983 -40.20964 1 vs 2 0.628421 0.4279 S-PLUS:> anova(m1, m2)Model df AIC BIC logLik Test L.Ratio p-value m1 1 4 87.00400 93.33808 -39.50200 m2 2 3 91.75447 96.50503 -42.87724 1 vs 2 6.75047 0.0094 Can we interpret that in any way? For example, is it reasonable to infer that in this case the R results, suggesting no significant difference between the model, are preferable? Andrew Associate Professor, Forest Biometrics University of Idaho Moscow ID 83843