# Mixed effects negative binomial with robust standard errors (Huber-white) in R

I would like to fit a random effects model in R using the negative binomial distribution and reporting robust standard errors.

I was going to try using the sandwich package to compute the robust standard errors from the fitted model object:

lmtest::coeftest(my_me_model, vcov = sandwich::vcovHC(my_me_model, type = "HC0"))


lme4::glmer.nb() allows me to fit a mixed effects model however I am unable to calculate robust standard errors, it looks like the model returned by lme4::glmer.nb() is an s4 class.

lmtest::coeftest(my_me_model, vcov = sandwich::vcovHC(my_me_model, type = "HC0"))
Error in UseMethod("estfun") :
no applicable method for 'estfun' applied to an object of class "c('glmerMod', 'merMod')"


robustlmm package function rlmer() allows me to calculate robust standard errors "huberization of likelihood and DAS-Scale estimation" however I cannot see a way to use the negative binomial with this package.

ptmixed package allows me to fit a mixed effects negative binomial I but cannot see a way to compute robust standard errors. So the reverse issue I encountered with robustlmm.

I also came across glmTMB package which also allows me to fit a negative bionomial mixed effects model, but where I am also unable to use e.g. sandwich to compute robust standard errors.

How can I fit a mixed effects negative binomial regression and then compute robust standard errors (Huber-white)?

• Your request is exceedingly specialized and I am not sure it has ever been implemented in publicly available software. Standard errors are not very precise and/or meaningful for negative binomial regression at the best of times, even more so for mixed models, and even more so if you want a robust estimate. It there a reason why you think this is (a) a good approach and (b) should be readily available in publicly available software? Aug 21 '20 at 6:38
• Have Huber-White robust standard errors ever been defined for a non-normal count-based model? Aug 22 '20 at 0:44

Sandwich covariances are available for a wide range of standard maximum likelihood models in sandwich (including the output from glm(), glm.nb(), zeroinfl(), and hurdle()). However, for mixed-effects models this is less straightforward but there is the relatively recent work in merDeriv:
Wang T, Merkle EC (2018). "merDeriv: Derivative Computations for Linear Mixed Effects Models with Application to Robust Standard Errors." Journal of Statistical Software, Code Snippets, 87(1), 1–16. doi:10.18637/jss.v087.c01.
They include the necessary methods for plugging the output from lmer() and glmer() into the sandwich() function and also provide a vcov() method for these objects. However, vcovHC() is not applicable but sandwich() essentially corresponds to HC0 anyway.
As for objects fitted by glmer.nb(): The vcov() method in merDeriv works for these objects but sandwich() does not. I would recommend to get in touch with the merDeriv authors and ask them whether it would be possible to support glmer.nb() as well.