Skip to content
Prev 82931 / 398506 Next

glmmADMB: Generalized Linear Mixed Models using AD Model Builder

On 12/15/05, Roel de Jong <dejongroel at gmail.com> wrote:
The current version of lmer takes method =  "PQL" (the default) or
"Laplace" or "AGQ" although AGQ is not available for vector-valued
random effects in that version so one must be content with "PQL" or
"Laplace"
Could you make the datasets, or the code that generates them,
available?  My code for such a simulation would be

genGLMM <- function(nobs, gsiz, fxd, Sigma, linkinv = binomial()$linkinv)
{
    ngrp <- nobs/gsiz
    ranef <- matrix(rnorm(ngrp * ncol(Sigma)), nr = ngrp) %*% chol(Sigma)
    pred2 <- rnorm(nobs)
    pred3 <- rnorm(nobs)
    mm <- model.matrix(~pred2 + pred3)
    rmm <- model.matrix(~pred2)
    grp <- gl(n = 1500/30, k = 30, len = 1500)
                                        # linear predictor
    lp <- as.vector(mm %*% fxd + rowSums(rmm * ranef[grp,]))
    resp <- as.integer(runif(nobs) < linkinv(lp))
    data.frame(resp = resp, pred2 = pred2, pred3 = pred3, grp = grp)
}

Running this function gives
Generalized linear mixed model fit using PQL
Formula: resp ~ pred2 + pred3 + (pred2 | grp)
   Data: sim1
 Family: binomial(logit link)
      AIC      BIC    logLik deviance
 1403.522 1440.714 -694.7609 1389.522
Random effects:
 Groups Name        Variance Std.Dev. Corr
 grp    (Intercept) 0.44672  0.66837
        pred2       0.55629  0.74585  0.070
# of obs: 1500, groups: grp, 50

Estimated scale (compare to 1)  0.9032712

Fixed effects:
             Estimate Std. Error z value  Pr(>|z|)
(Intercept) -1.081710   0.121640 -8.8927 < 2.2e-16
pred2        1.607273   0.141697 11.3430 < 2.2e-16
pred3        0.531071   0.072643  7.3107 2.657e-13
[1] 0.33 0.00 0.33 0.00 0.00
Generalized linear mixed model fit using Laplace
Formula: resp ~ pred2 + pred3 + (pred2 | grp)
   Data: sim1
 Family: binomial(logit link)
      AIC      BIC    logLik deviance
 1401.396 1438.588 -693.6979 1387.396
Random effects:
 Groups Name        Variance Std.Dev. Corr
 grp    (Intercept) 0.35248  0.59370
        pred2       0.46641  0.68294  0.077
# of obs: 1500, groups: grp, 50

Estimated scale (compare to 1)  0.9854841

Fixed effects:
             Estimate Std. Error z value  Pr(>|z|)
(Intercept) -1.119008   0.121640 -9.1993 < 2.2e-16
pred2        1.680916   0.141697 11.8627 < 2.2e-16
pred3        0.543548   0.072643  7.4825 7.293e-14
[1] 4.62 0.01 4.65 0.00 0.00

Fitting that model using glmmADMB gives
...
iteration output omitted
...

GLMM's in R powered by AD Model Builder:

  Family: binomial

Fixed effects:
  Log-likelihood: -602.035
  Formula: resp ~ pred2 + pred3
(Intercept)       pred2       pred3
   -1.11990     1.69030     0.54619

Random effects:
  Grouping factor: grp
  Formula: ~pred2
Structure: General positive-definite
               StdDev      Corr
(Intercept) 0.5890755
pred2       0.6712377 0.1023698

Number of Observations: 1500
Number of Groups: 50

The "Laplace" method in lmer and the default method in glmm.admb,
which according to the documentation is the Laplace approximation,
produce essentially the same model fit.  One difference is the
reported value of the log-likelihood, which we should cross-check, and
another difference is in the execution time
...
Iteration output omitted
...
[1]  0.23  0.02 21.44 19.45  0.24

Fitting this model takes about 4.7 seconds with the Laplace
approximation in lmer (and only 0.33 seconds for PQL, which is not
that far off) and about 20 seconds in glmm.admb
I agree.  I am particularly pleased that Otter Research allows access
to a Linux executable of their code (although I would, naturally,
prefer the code to be Open Source).