GLMM: tricking software to estimate normal error at lower level

Paul Johnson <pauljohn32 at ...> writes:

Thanks for the help last week on GLMM estimators.  I'm gradually
figuring this out.

    Some *very* quick responses to a good question (if I wait this
will just get buried in the e-mail pile).

1. The general strategy you're describing (lognormal-Poisson error
via observation-level random effect) has been discussed on this list
before, I think: it's in the GLMM FAQ, the new version of which is at

https://rawgit.com/bbolker/mixedmodels-misc/master/glmmFAQ.html

(search for "overdispersion")

2.  When you predict from the model with the observation-level
random effects, you should leave out those terms: use re.form=~0
or re.form=NULL if the observation-level RE is the only random
effect in your model, otherwise use an appropriate re.form that
includes only the higher-level groups.  That should make the predictions
make a lot more sense.

3. I'm not sure your convergence problems are actually that bad
(i.e. false positive); do the results from different optimizers
agree reasonably well?  (I would compare "bobyqa" and "nloptwrap")

4. It doesn't make sense to specify a family argument with glmer.nb

   cheers
     Ben Bolker

  [I'm deleting a huge amount below because I'm posting via gmane, which
doesn't like too much quoted stuff -- argh.)

Today I'm writing notes comparing NB with Poisson models. Recall, NB
is Poisson with a log(gamma) error variable added to the linear
predictor.  I've often wondered how different that would be from
adding a normal error, but did not have way to estimate.

In Rabe-Hesketh and Skrondal Multilevel / Longitudinal Modeling Using
Stata, 3ed, they show a way to estimate a normally distributed random
error inserted into a Poisson GLM in a one level model by "tricking"
an estimator for a multilevel model (p. 707). I've checked, can do
with xtpoisson or meglm in Stata 14).  Code is below.

I have some trouble with glmer trying to do same. Estimates are close
to the "seems to work" answer from Stata, but I get glmer convergence
warning and different numbers, grossly different predicted values. I
understand glmer was not intended for this purpose, but wouldn't it be
fun if we could make it work?

Here's the basic plan in Stata that RHS demonstrate. Declare an id
variable with values 1, ..., N. This is non-replicated data, there
will be just one observation per group.

generate id = _n
xtset id
xtpoisson y x1 x2 x3, normal

The Stata xtpoisson defaults the additive random error as log(gamma),
but you can ask for normal. I've confirmed this runs, and it also
works in the newer Stata meglm fitting function (meglm is only for
normal errors).  Then you can compare side by side xtpoisson with
log(gamma), xtpoisson with normal, and meglm with family(poisson) and
normal mixed effects.

I fiddled around with glmer trying to do same, and the results are not
completely different, but I hit convergence errors. Since I don't
understand the fitting process, as I confessed last week, I'm a little
blocked here.

My minimal running example R code:

## Paul Johnson <pauljohn <at> ku.edu>
## 20160417

library(MASS)
quine$id <- 1:NROW(quine)

library(foreign)
write.dta(quine, file = "quine.dta12")
## Original model in MASS example(glm.nb)
quine.pois1 <- glm(Days ~ Sex/(Age + Eth*Lrn), data = quine, family =
poisson(log))
summary(quine.pois1)
## Don't know how to write Stata formula to compare to that, so
## rewrite without "/"
quine.pois2 <- glm(Days ~ Sex*Age + Sex*Eth*Lrn, family =
poisson(link=log), data = quine)
summary(quine.pois2)

quine.nb1 <- glm.nb(Days ~ Sex*Age + Sex*Eth*Lrn, data = quine)
summary(quine.nb1)
quine$predpois <- predict(quine.pois2, type = "link")
quine$prednb <- predict(quine.nb1, type = "link")

plot(predpois ~ prednb, data = quine)
with(quine, cor(predpois, prednb))

library(lme4)
# Stata uses  nAGQ=14 equivalent:
glmer1 <- glmer(Days ~ Sex*Age + Sex*Eth*Lrn + (1|id), data = quine,
                 family = "poisson", nAGQ = 14)
summary(glmer1)
quine$predglmer1 <- predict(glmer1)
plot(predglmer1 ~ prednb, data = quine)

glmer2 <- glmer.nb(Days ~ Sex*Age + Sex*Eth*Lrn + (1|id), data =
quine, family = "poisson",
                 control = glmerControl(optCtrl = list(maxfun = 10000)))

glmer3 <- glmer(Days ~ Sex*Age + Sex*Eth*Lrn + (1|id), data = quine,
family = "poisson",
                      control = glmerControl(optimizer = "Nelder_Mead",
optCtrl = list(maxfun = 10000)))
summary(glmer3)

## end of MRE