[R] glmmPQL

[Christof Bigler]

Dear R users,

can anybody explain me why the function glmmPQL(.) behaves in different 
ways, depending on the number of measurements/individuals you use? To 
show you this, I generated two examples. The first one includes 20 
indivduals with each 100 repeated measurements (binary response), the 
second one includes 40 individuals. The 'individuals' differ only in 
different x values. I fitted logistic regression models with and without 
random intercepts.

The coefficients and p values between dummy.glm20 and dummy.glmm20 
differ. However, dummy.glm40 and dummy.glmm40 have the same coefficients 
and p values, respectively. Did the solution in the second example not 
converge (only one iteration step)?

Why does the AIC between e.g. dummy.glm20 and dummy.glmm20 differ so 
much?
And last question: how can dummy.glm20 and dummy.glmm20 be compared with 
an anova(.) function?

I use R 1.5.0 (Darwin version) on a Mac G4.

Thanks for any advice.

Christof

#######################################
#Example 1
 > y <- c(rep(0,40),sample(c(rep(0,20),rep(1,20)),40),rep(1,20))
 > dummy20 <- data.frame(ID=as.factor(c(rep(1:20,each=100))),
                     y=rep(y,20),
                     x=c(1:2000))
 > dummy.glm20 <- glm(y~x,data=dummy20,family=binomial)
 > summary(dummy.glm20)
...
Coefficients:
               Estimate Std. Error z value Pr(>|z|)
(Intercept) -5.330e-01  9.199e-02  -5.795 6.85e-09 ***
x            1.270e-04  7.915e-05   1.604    0.109
...
Null deviance: 2692.0  on 1999  degrees of freedom
Residual deviance: 2689.5  on 1998  degrees of freedom
AIC: 2693.5

 > dummy.glmm20 <- glmmPQL(y~x,random=~1 | ID,
+                       data=dummy20,family=binomial)
iteration 1
iteration 2
iteration 3
iteration 4
iteration 5

 > summary(dummy.glmm20)
Linear mixed-effects model fit by maximum likelihood
  Data: dummy20
        AIC      BIC    logLik
   10270.78 10293.19 -5131.392

Random effects:
  Formula: ~1 | ID
         (Intercept)  Residual
StdDev:    50.55603 0.7928198

Variance function:
  Structure: fixed weights
  Formula: ~invwt
Fixed effects: y ~ x
                 Value Std.Error   DF   t-value p-value
(Intercept) -88.62246 11.686506 1979 -7.583316  <.0001
x             0.08768  0.002913 1979 30.099686  <.0001
  Correlation:
   (Intr)
x -0.252

Standardized Within-Group Residuals:
        Min         Q1        Med         Q3        Max
-2.4592661 -0.4117634 -0.1389889  0.4223991  2.8757740

Number of Observations: 2000
Number of Groups: 20

#Example 2
 > dummy40 <- data.frame(ID=as.factor(c(rep(1:40,each=100))),
                     y=rep(y,40),
                     x=c(1:4000))
 > dummy.glm40 <- glm(y~x,data=dummy40,family=binomial)
 > summary(dummy.glm40)
...
Coefficients:
               Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.691e-01  6.478e-02  -7.241 4.45e-13 ***
x            3.173e-05  2.796e-05   1.135    0.256
...
Null deviance: 5384.1  on 3999  degrees of freedom
Residual deviance: 5382.8  on 3998  degrees of freedom
AIC: 5386.8

 > dummy.glmm40 <- glmmPQL(y~x,random=~1 | ID,
+                       data=dummy40,family=binomial)
iteration 1
 > summary(dummy.glmm40)
Linear mixed-effects model fit by maximum likelihood
  Data: dummy40
        AIC      BIC    logLik
   17068.15 17093.32 -8530.073

Random effects:
  Formula: ~1 | ID
         (Intercept)  Residual
StdDev: 0.003066688 0.9999229

Variance function:
  Structure: fixed weights
  Formula: ~invwt
Fixed effects: y ~ x
                  Value  Std.Error   DF   t-value p-value
(Intercept) -0.4690798 0.06479625 3959 -7.239305  <.0001
x            0.0000317 0.00002797 3959  1.134708  0.2566
  Correlation:
   (Intr)
x -0.867

Standardized Within-Group Residuals:
       Min        Q1       Med        Q3       Max
-0.842464 -0.820600 -0.798994  1.214569  1.263420

Number of Observations: 4000
Number of Groups: 40

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._