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]]


