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?