From: Douglas Bates <dmbates_at_gmail.com>

Date: Sun 18 Dec 2005 - 07:51:54 EST

mm <- model.matrix(~pred2 + pred3)

rmm <- model.matrix(~pred2)

grp <- gl(n = 1500/30, k = 30, len = 1500)

1403.522 1440.714 -694.7609 1389.522

Random effects:

Groups Name Variance Std.Dev. Corr grp (Intercept) 0.44672 0.66837

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

1401.396 1438.588 -693.6979 1387.396

Random effects:

Groups Name Variance Std.Dev. Corr grp (Intercept) 0.35248 0.59370

*> system.time(fm2 <- lmer(resp ~ pred2 + pred3 + (pred2|grp), sim1, binomial, method = "Laplace"))
*

[1] 4.62 0.01 4.65 0.00 0.00

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Sun Dec 18 08:05:13 2005

Date: Sun 18 Dec 2005 - 07:51:54 EST

On 12/15/05, Roel de Jong <dejongroel@gmail.com> wrote:

*> Dear R-users,
**>
*

> because lme(r) & glmmpql, which are based on Penalized Quasi Likelihood,

*> are not very robust with Bernoulli responses,
*

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

genGLMM <- function(nobs, gsiz, fxd, Sigma, linkinv = binomial()$linkinv) {

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

*> (fm3 <- glmm.admb(resp ~ pred2 + pred3, ~ pred2, "grp", sim1, "binomial", "logit", "full"))
*

...

iteration output omitted

...

GLMM's in R powered by AD Model Builder:

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

> 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!
*

*>
**> 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.
**> >
**> > http://otter-rsch.com/admbre/examples/glmmadmb/glmmADMB.html
**> >
**> > Several people who used this package have requested
**> > additional features. We now have a new version ready.
**> > The major new feature is that glmmADMB allows Bernoulli responses
**> > with logistic and probit links. In addition there is
**> > a "ranef.glmm.admb()" function for getting the random effects.
**> >
**> > The download site is still:
**> >
**> > http://otter-rsch.com/admbre/examples/glmmadmb/glmmADMB.html
**> >
**> > 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
**> > about the software.
**> >
**> > 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@stat.math.ethz.ch mailing list
**> > https://stat.ethz.ch/mailman/listinfo/r-help
**> > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
**> >
**>
**> ______________________________________________
**> R-help@stat.math.ethz.ch mailing list
**> https://stat.ethz.ch/mailman/listinfo/r-help
**> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
**>
*

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Sun Dec 18 08:05:13 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:41:40 EST
*