I would like to fit a mixed effects model using negative binomial and also compute robust standard errors (Huber-White).
I started out using lme4 package:
## modified from https://www.r-econometrics.com/methods/hcrobusterrors/ library(tidyverse) library(wooldridge) library(sandwich) library(lmtest) library(lme4) library(lmtest) # Load the sample saving <- wooldridge::saving %>% mutate(age_bin = cut(age, breaks = quantile(age))) # Only use positive values of saving, which are smaller than income saving <- saving %>% filter(sav > 0, inc < 20000, sav < inc) # make a mixed effects model mod.me.nb <- lme4::glmer.nb(sav ~ (1|age_bin) + scale(educ) + scale(inc), data = saving) summary(mod.me.nb) Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod'] Family: Negative Binomial(1.1842) ( log ) Formula: sav ~ (1 | age_bin) + scale(educ) + scale(inc) Data: saving AIC BIC logLik deviance df.resid 1259.4 1271.0 -624.7 1249.4 70 Scaled residuals: Min 1Q Median 3Q Max -1.0857 -0.6854 -0.2256 0.2799 3.4315 Random effects: Groups Name Variance Std.Dev. age_bin (Intercept) 4.406e-11 6.638e-06 Number of obs: 75, groups: age_bin, 4 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) 7.3375 0.1062 69.121 <2e-16 *** scale(educ) 0.1066 0.1181 0.903 0.3664 scale(inc) 0.2785 0.1178 2.363 0.0181 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation of Fixed Effects: (Intr) scl(d) scale(educ) 0.000 scale(inc) 0.000 -0.310 convergence code: 0 boundary (singular) fit: see ?isSingular
I then tried using the sandwich package to compute robust standard errors but this does not work:
## from https://www.r-econometrics.com/methods/hcrobusterrors/ > lmtest::coeftest(mod.me.nb, vcov = sandwich::vcovHC(mod.me.nb, type = "HC0")) Error in UseMethod("estfun") : no applicable method for 'estfun' applied to an object of class "c('glmerMod', 'merMod')" In addition: Warning message: In hatvalues.merMod(x) : the hat matrix may not make sense for GLMMs
A comment on here from earlier was that perhaps the returned object from glmer.nb is not what was expected by
sandwich::vcovHC(). mod.me.nb is an s4 object.
I then found the robustlmm package. Here there is a function
rlmer() with approach "huberization of likelihood and DAS-Scale estimation" however I cannot see a way to use the negative binomial with this package, instead it looks like it is based on the lmer function, but I'm not sure how to tell this function to use negative binomial.
I then found ptmixed. Here I have the inverse issue from robustlmm. I am able to fit a mixed effects negative binomial but cannot see a way to compute robust standard errors.
How can I fit a mixed effects negative binomial regression and then compute robust standard errors?