Skip to content

fitting glmm under lme4 and Gauss Hermite integration

3 messages · Jie Li, Ben Bolker

Original
Ben Bolker
#
On 10-10-30 01:07 PM, Jie Li wrote:
In principle, yes.
 1. You're essentially using the random effects here to account for
individual-level variance (i.e. a lognormal-Poisson model), rather than
accounting for grouping/correlation.
 2. This is a very small dataset. In particular, fitting k=6
fixed-effect parameters to n=12 data points means that you will be
far from reliable asymptotic rules of thumb (we would generally
prefer n/k >= 10 ...). I would strongly recommend bootstrapped
confidence intervals, or some other reasonably robust small-sample
procedure (permutation test, parametric bootstrap ...)
glmer uses the Laplace approximation (aGH with 1 quadrature point)
unless you set nAGQ>1.  I might recommend [Jiang, Jiming. 2007. Linear
and generalized linear mixed models and their applications. Springer].
<http://lme4.r-forge.r-project.org/book/Ch5.pdf> is Bates's draft of his
new book, but it only goes as far as defining the Laplace approximation.
 If you already have "a book" that talks about GH, does it give any more
information?

  I would suggest comparing the results of glmer() and glmm() for
different numbers of quadrature points; see if the two packages give
similar estimates, and see whether the estimates seem to be behaving
reasonably/converging to something sensible as you increase the number
of quadrature points. (For "extra credit" :-) or if this is really
important you could also try comparing the results with the glmmML package.)

  good luck,
    Ben Bolker
Ben Bolker
# 
PS it looks like the coefficients are quite stable under increasing
numbers of quadrature points ...


x=factor(rep(1:6,each=2))
y=c(14, 18, 10, 47, 13, 24, 10, 12,  1,  0,  6,  6)
adjust=c(5.505332,5.645447, 5.549076, 5.886104, 6.023448,
  6.177944, 6.077642, 6.186209, 5.030438, 5.117994, 5.170484, 5.493061)
eu=factor(1:12)
try=data.frame(x,y,eu,adjust)

library(lme4)
fit <- glmer(y ~ x + (1 | eu) + offset(log(adjust)),
              family = poisson, data =try)

nagqvec <- 1:20
fitlist1 <- lapply(nagqvec,
                   glmer,
                   formula=y ~ x + (1 | eu) + offset(log(adjust)),
                   family=poisson,
                   data=try,
                   start=NULL,
                   verbose=FALSE)

t(sapply(fitlist1,fixef))
On 10-10-30 01:07 PM, Jie Li wrote: