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)?

  • $\begingroup$ 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? $\endgroup$ Aug 21, 2020 at 6:38
  • $\begingroup$ Have Huber-White robust standard errors ever been defined for a non-normal count-based model? $\endgroup$ Aug 22, 2020 at 0:44

1 Answer 1


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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.