geex and
sandwich for robust covariance estimationThis 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.