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