`geex`

and sandwich
for robust covariance estimation`vignettes/v02_sandwich_comparison.Rmd`

`v02_sandwich_comparison.Rmd`

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:

```
##
## Attaching package: 'geex'
```

```
## The following object is masked from 'package:methods':
##
## show
```

```
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.