[R-sig-ME] lmer logistic nested random effects
Joshua Wiley
jwiley.psych at gmail.com
Tue Aug 14 05:53:16 CEST 2012
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/
More information about the R-sig-mixed-models
mailing list