[R-sig-ME] GLMM: tricking software to estimate normal error at lower level
bbolker at gmail.com
Mon Apr 18 03:55:37 CEST 2016
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
(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
[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
> quine$id <- 1:NROW(quine)
> write.dta(quine, file = "quine.dta12")
> ## Original model in MASS example(glm.nb)
> quine.pois1 <- glm(Days ~ Sex/(Age + Eth*Lrn), data = quine, family =
> ## 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)
> quine.nb1 <- glm.nb(Days ~ Sex*Age + Sex*Eth*Lrn, data = quine)
> quine$predpois <- predict(quine.pois2, type = "link")
> quine$prednb <- predict(quine.nb1, type = "link")
> plot(predpois ~ prednb, data = quine)
> with(quine, cor(predpois, prednb))
> # Stata uses nAGQ=14 equivalent:
> glmer1 <- glmer(Days ~ Sex*Age + Sex*Eth*Lrn + (1|id), data = quine,
> family = "poisson", nAGQ = 14)
> 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)))
> ## end of MRE
More information about the R-sig-mixed-models