A comparison of `geex` and sandwich for robust covariance estimation
In geex: An API for M-Estimation

This examples uses the vaccinesim dataset from the inferference package to compare the estimated covariance matrix obtained from geex and sandwich. An example $\psi$ function written in R.

This function computes the score functions for a GLM.

eefun <- function(data, model){
  X <- model.matrix(model, data = data)
  Y <- model.response(model.frame(model, data = data))
  function(theta){
    lp  <- X %*% theta
    rho <- plogis(lp)

    score_eqns <- apply(X, 2, function(x) sum((Y - rho) * x))
    score_eqns
  }
}

Compare sandwich variance estimators to sandwich treating individuals as units:

library(geex)
library(inferference)
mglm    <- glm(A ~ X1, data = vaccinesim, family = binomial)
estimates <- m_estimate(
  estFUN = eefun,
  data = vaccinesim,
  root_control = setup_root_control(start = c(-.35, 0)),
  outer_args = list(model = mglm))

# Compare point estimates
coef(estimates) # from GEEX
coef(mglm) # from the GLM function

# Compare variance estimates
vcov(estimates)
sandwich::sandwich(mglm)

Pretty darn good! Note that the geex method is much slower than sandwich (especially using method = 'Richardson' for numDeriv), but this is because sandwich uses the closed form of the score equations, while geex compute them numerically. However, geex's real utility comes when you have more complicated estimating equations. Also, the analyst has the ability to code faster $\psi$ functions by optimizing their code or using Rccp, for example.