[R-sig-ME] lmer logistic nested random effects

[Joshua Wiley]

Hi Ben,

It seems as though you are correct that Stata uses more quadrature
points.  However, I do not believe that SAS does for more than a
single grouping factor.  Obviously this was not exhaustive, but I
tried when I wrote this page:

http://www.ats.ucla.edu/stat/sas/code/hdp.htm

and at least proc glimmix would not work with anything besides pseudo
likelihoods based on a linearization of the problem or the Laplace
approximation (quadrature with a single integration point).  I happen
to have also tried Stata and SAS with cross classified random effects
in a logistic model (this time with real data) and know that neither
of them use more integration points.

Another option you can try is to use the glmmADMB package which
implements a variety of models optimized in the backend with the AD
model builder.  Here is an example you can try and compare.  Note that
glmmadmb() was _substantially_ slower and took minutes to complete on
my machine with even this simple intercept only model.

hdp <- read.csv("http://www.ats.ucla.edu/stat/data/hdp.csv")
hdp$DID <- factor(hdp$DID)
hdp$HID <- factor(hdp$HID)
require(lme4)
m1 <- glmer(remission ~ 1 + (1 | DID) + (1 | HID), family = binomial,
data = hdp)
require(glmmADMB)
m2 <- glmmadmb(remission ~ 1 + (1 | DID) + (1 | HID), family =
"binomial", data = hdp)


which for me gives:

> summary(m1)
Generalized linear mixed model fit by maximum likelihood ['summary.mer']
 Family: binomial ( logit )
Formula: remission ~ 1 + (1 | DID) + (1 | HID)
   Data: hdp

      AIC       BIC    logLik  deviance
 7884.569  7905.721 -3939.284  7878.569

Random effects:
 Groups Name        Variance Std.Dev.
 DID    (Intercept) 3.4578   1.8595
 HID    (Intercept) 0.2015   0.4489
Number of obs: 8525, groups: DID, 407; HID, 35

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.3973     0.1269  -11.01   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> summary(m2)

Call:
glmmadmb(formula = remission ~ 1 + (1 | DID) + (1 | HID), data = hdp,
    family = "binomial")


Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)    -1.40       0.13   -10.7   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Number of observations: total=8525, DID=407, HID=35
Random effect variance(s):
Group=DID
            Variance StdDev
(Intercept)    3.456  1.859
Group=HID
            Variance StdDev
(Intercept)   0.2036 0.4513

Log-likelihood: -3939.28

Cheers,

Josh


On Mon, Aug 13, 2012 at 2:37 PM, Ben Pelzer <b.pelzer at maw.ru.nl> wrote:
> Dear list,
>
> I would like to fit a logistic model with lmer, for nested factors like
> classes nested in schools, and Y being 0 or 1. So I have dichotomous Y
> values of pupils in classes in schools. My null-model reads:
>
> lmer ( Y ~ (1|school/class), family=binomial(link=logit), nAGQ=7)
>
> The options nAGQ=7 causes the following error:
>
>    AGQ method requires a single grouping factor
>
> Omitting the nAGQ=7 works fine and produces Laplace approximation of the
> loglikelihood.
>
> Am I right in concluding that for hierarchical nested factors (like class
> nested in school) glmer can only use 1 quadrature point? Or am I overlooking
> something?
>
> Of course, more quadrature points would be nicer, and in e.g. xtmelogit in
> stata and glimmix in sas I believe more quad. points are possible (i.e., for
> hierarchically nested factors, not for crossed factors) Hence I was
> surprised that in lmer, it seems impossible. Could anyone confirm this,
> please? Kind regards,
>
> Ben.
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
Programmer Analyst II, Statistical Consulting Group
University of California, Los Angeles
https://joshuawiley.com/