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

[Jarrod Hadfield]

Hi,

Regarding 2) you also need to specify type="terms" in the glmer predict 
function in orer to make them comparable to the previous link function 
predictions.

Cheers,

Jarrod


On 18/04/2016 02:55, Ben Bolker wrote:
>
>
> 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
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>


-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.