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)?
 A: 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.
