Roel de Jong dejongroel at gmail.com
Mon Dec 19 00:16:13 CET 2005

```Dear professor Bates,

thank you for your reaction. To make sure that no errors occur in the
data generation process I used the elegant function you so neatly
provided to generate a couple of datasets under the model specification
specified earlier. Running lmer with a Laplace approximation to the
high-dimensional integral in the likelihood gives me a warning an then
this show-stopper:

Warning: IRLS iterations for PQL did not converge
Error in objective(.par, ...) : Unable to invert singular factor of
downdated X'X

Fitting the dataset with glmmADMB gives no apparent problems and
reasonable estimates. I attached the particular dataset to the email.

The difference in computation time can be attributed to the fact that
glmmadmb uses a generic technique called automatic differentiation with
the Laplace approximation. The same technique can be employed to fit
much more complex nonlinear models, but I'm sure Hans & Dave can tell

Best regards,
Roel de Jong

Douglas Bates wrote:
> On 12/15/05, Roel de Jong <dejongroel at gmail.com> wrote:
>
>>Dear R-users,
>>
>>because lme(r) & glmmpql, which are based on Penalized Quasi Likelihood,
>>are not very robust with Bernoulli responses,
>
>
> 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"
>
>
>>I wanted to test glmmADMB.  I run the following simulation study:
>
>
>>500 samples are drawn with the model specification:
>>y = (intercept*f1+pred2*f2+pred3*f3)+(intercept*ri+pred2*rs)
>>     where pred2 and pred3 are predictors distributed N(0,1)
>>     f1..f3 are fixed effects, f1=-1, f2=1.5, f3=0.5
>>     ri is random intercept with associated variance var_ri=0.2
>>     rs is random slope with associated variance var_rs=0.4
>>     the covariance between ri and rs "covr"=0.1
>>
>>1500 units/dataset, class size=30
>
>
> Could you make the datasets, or the code that generates them,
> available?  My code for such a simulation would be
>
> {
>     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
>
>>nobs <- 1500
>>gsiz <- 30
>>fxd <- c(-1, 1.5, 0.5)
>>Sigma <- matrix(c(0.2, 0.1, 0.1, 0.4), nc = 2)
>>set.seed(123454321)
>>sim1 <- genGLMM(nobs, gsiz, fxd, Sigma)
>>(fm1 <- lmer(resp ~ pred2 + pred3 + (pred2|grp), sim1, binomial))
>
> Generalized linear mixed model fit using PQL
> Formula: resp ~ pred2 + pred3 + (pred2 | grp)
>    Data: sim1
>       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
>
>>system.time(fm1 <- lmer(resp ~ pred2 + pred3 + (pred2|grp), sim1, binomial))
>
> [1] 0.33 0.00 0.33 0.00 0.00
>
>>(fm2 <- lmer(resp ~ pred2 + pred3 + (pred2|grp), sim1, binomial, method = "Laplace"))
>
> Generalized linear mixed model fit using Laplace
> Formula: resp ~ pred2 + pred3 + (pred2 | grp)
>    Data: sim1
>       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
>
>>system.time(fm2 <- lmer(resp ~ pred2 + pred3 + (pred2|grp), sim1, binomial, method = "Laplace"))
>
> [1] 4.62 0.01 4.65 0.00 0.00
>
> Fitting that model using glmmADMB gives
>
>>(fm3 <- glmm.admb(resp ~ pred2 + pred3, ~ pred2, "grp", sim1, "binomial", "logit", "full"))
>
> ...
> iteration output omitted
> ...
>
>
>   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
>
>
>>system.time(fm3 <- glmm.admb(resp ~ pred2 + pred3, ~ pred2, "grp", sim1, "binomial", "logit", "full"))
>
> ...
> 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
>
>
>
>
>>convergence:
>>~~~~~~~~~~~~
>>No crashes.
>>5/500 Datasets had on exit a gradient of the log-likelihood > 0.001
>>though. Removing the datasets with questionable convergence doesn't seem
>>to effect the simulation analysis.
>>
>>bias:
>>~~~~~~
>>f1=-1.00531376
>>f2= 1.49891060
>>f3= 0.50211520
>>ri= 0.20075947
>>covr=0.09886267
>>rs= 0.38948382
>>
>>Only the random slope "rs" is somewhat low, but i don't think it is of
>>significance
>>
>>coverage alpha=.95: (using asymmetric confidence intervals)
>>~~~~~~~~~~~~~~~~~~~~~~~~
>>f1=0.950
>>f2=0.950
>>f3=0.966
>>ri=0.974
>>covr=0.970
>>rs=0.970
>>
>>While some coverages are somewhat high, confidence intervals based on
>>asymptotic theory will not have exactly the nominal coverage level, but
>>with simulations (parametric bootstrap) that can be corrected for.
>>
>>I can highly recommend this excellent package to anyone fitting these
>>kinds of models, and want to thank Hans Skaug & Dave Fournier for their
>>hard work!
>
>
> 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).
>
>
>>Roel de Jong.
>>
>>
>>Hans Julius Skaug wrote:
>>
>>>Dear R-users,
>>>
>>>Half a year ago we put out the R package "glmmADMB" for fitting
>>>overdispersed count data.
>>>
>>>
>>>Several people who used this package have requested
>>>The major new feature is that glmmADMB allows Bernoulli responses
>>>a "ranef.glmm.admb()" function for getting the random effects.
>>>
>>>
>>>
>>>The package is based on the software ADMB-RE, but the full
>>>unrestricted R-package is made freely available by Otter Research Ltd
>>>and does not require ADMB-RE to run. Versions for Linux and Windows
>>>are available.
>>>
>>>We are still happy to get feedback for users, and to get suggestions
>>>for improvement.
>>>
>>>We have set up a forum at http://www.otter-rsch.ca/phpbb/ for discussions
>>>
>>>Regards,
>>>
>>>Hans
>>>
>>>_____________________________
>>>Hans Julius Skaug
>>>
>>>Department of Mathematics
>>>University of Bergen
>>>Johannes Brunsgate 12
>>>5008 Bergen
>>>Norway
>>>ph. (+47) 55 58 48 61
>>>
>>>______________________________________________
>>>R-help at stat.math.ethz.ch mailing list
>>>https://stat.ethz.ch/mailman/listinfo/r-help
>>>
>>
>>______________________________________________
>>R-help at stat.math.ethz.ch mailing list
>>https://stat.ethz.ch/mailman/listinfo/r-help