A comparison of geex and sandwich for robust covariance estimation

Bradley Saul

2022-07-24

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

## [1] -0.36869683 -0.02037916

coef(mglm) # from the GLM function

## (Intercept)          X1 
## -0.36869683 -0.02037916

# Compare variance estimates
vcov(estimates)

##               [,1]          [,2]
## [1,]  0.0028345579 -0.0007476536
## [2,] -0.0007476536  0.0003870030

sandwich::sandwich(mglm)

##               (Intercept)            X1
## (Intercept)  0.0028345579 -0.0007476536
## X1          -0.0007476536  0.0003870030

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.