[R-sig-ME] GLMM: tricking software to estimate normal error at lower level

[Ben Bolker]

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