[R-sig-ME] Function to get residuals & fitted values from Poisson GLMM with OLRE

[Paul Johnson]

Hi all,

I’ve written a function which some might find useful. It outputs “corrected” Pearson residuals from a glmer fit of the form

   fit <- glmer(y ~ . + (1 | obs), family = “poisson”)

where "(1 | obs)” fits an observation level random effect (OLRE), i.e. length(obs) == length(y).*

The function is called residfitted.olre and can be found here: 
https://github.com/pcdjohnson/miscR

The default behaviour of residuals(fit) and fitted(fit) is to treat the OLRE as part of the fitted values, but for most practical purposes I want the OLRE to be in the residuals and not in the fitted values. In particular, if I’m using a residuals-vs-fitted-values plot to assess the fit, the default residuals and fitted functions used by plot(fit) will generally produce nasty-looking plots showing a trend and severe heteroscedasticity even from well fitting models. The “corrected” residuals fix this problem.

There’s also example code below the function showing how this it allows interpretable residuals-vs-fitted plots to be produced from a Poisson GLMM with an OLRE (aka a Poisson-lognormal GLMM).

The same problem affects a binomial GLMM with an OLRE, but I haven’t (yet) made the function work for binomial GLMMs for the reason covered in this post:
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2014q1/021818.html
(basically I don’t know what the variance function is for a mixture of binomial + logit-normal — anyone know?).

Hope someone finds it useful,
Paul

*Refs:
https://peerj.com/articles/616/
http://www.ncbi.nlm.nih.gov/pubmed/11393830