Fit a mixed effects model using negative binomial and also compute robust standard errors

[dougfir]

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?

[dougfir]

Adding to the list of relevant packages that I have found: CRAN - Package glmmTMB

Looks like I can fit a mixed effects negbin there but I still cannot see a way to get robust standard errors. Relevant Function vcovHC from package sandwich does not work on glmmTMB object · Issue #443 · glmmTMB/glmmTMB · GitHub

At present this just isn't available. In order to make this work we'd need to figure out a convenient way implement the pieces needed for robust estimation, following this document: it looks like we could get the "HC0" and "HC1" estimators described there relatively easily, but for the "HC2-HC4" estimators we would need the leverages/"hat values"/diagonals of the hat matrix, which are not easily available from glmmTMB objects at the moment (see #391).