`geex`

`vignettes/v05_finite_sample_corrections.Rmd`

`v05_finite_sample_corrections.Rmd`

The empirical sandwich variance estimator is known to underestimate \(V(\theta)\) in small samples (Fay and Graubard 2001). Particularly in the context of GEE, many authors have proposed corrections that modify components of \(\hat{\Sigma}\) and/or by assuming \(\hat{\theta}\) follows a \(t\) (or \(F\)), as opposed to Normal, distribution with some estimated degrees of freedom. Many of the proposed corrections somehow modify a combination of the \(A_i\), \(A_m\), \(B_i\), or \(B_m\) matrices.

`geex`

provides an API that allows users to specify
functions that utilize these matrices to form corrections. A finite
sample correction function at a minimum takes the argument
`components`

, which is an object of class
`sandwich_components`

. For example,

```
correct_by_nothing <- function(components){
A <- grab_bread(components)
B <- grab_meat(components)
compute_sigma(A = A, B = B)
}
```

is a correctly formed function that does no corrections. Additional arguments may also be specified, as shown in the example.

`geex`

The `geex`

package includes the bias correction and
degrees of freedom corrections proposed by Fay
and Graubard (2001) in the `correct_by_fay_bias`

and
`correct_by_fay_df`

functions respectively. The following
demonstrates the construction and use of the bias correction. Fay and Graubard (2001) proposed the modified
variance estimator \(\hat{\Sigma}^{bc}(b) =
A_m^{-1} B_m^{bc}(b) \{A_m^{-1}\}^{\intercal}/m\), where:

\[\begin{equation} \label{eq:bc} B^{bc}_m(b) = \sum_{i = 1}^m H_i(b) B_i H_i(b)^{\intercal}, \end{equation}\]

\[\begin{equation} \label{eq:H} H_i(b) = \{1 - \min(b, \{A_i A^{-1}\}_{jj}) \}^{-1/2}, \end{equation}\]

and \(W_{jj}\) is the \((j, j)\) element of a matrix \(W\). When \(\{A_i A^{-1}\}_{jj}\) is close to 1, the adjustment to \(\hat{\Sigma}^{bc}(b)\) may be extreme, and the constant \(b\) is chosen by the analyst to limit over adjustments.

The bias corrected estimator \(\hat{\Sigma}^{bc}(b)\) can be implemented
in `geex`

by the following function:

```
bias_correction <- function(components, b){
A <- grab_bread(components)
A_i <- grab_bread_list(components)
B_i <- grab_meat_list(components)
Ainv <- solve(A)
H_i <- lapply(A_i, function(m){
diag( (1 - pmin(b, diag(m %*% Ainv) ) )^(-0.5) )
})
Bbc_i <- lapply(seq_along(B_i), function(i){
H_i[[i]] %*% B_i[[i]] %*% H_i[[i]]
})
Bbc <- apply(simplify2array(Bbc_i), 1:2, sum)
compute_sigma(A = A, B = Bbc)
}
```

The `compute_sigma`

function simply computes \(A^{-1} B \{A^{-1}\}^{\intercal}\). Note
that `geex`

computes \(A_m\)
and \(B_m\) as the sums of \(A_i\) and \(B_i\) rather than the means, hence the
appropriate function in the `apply`

call is `sum`

and not `mean`

. To use this bias correction, the
`m_estimate`

function accepts a named list of corrections to
perform. Each element of the list is also a list with two elements:
`correctFUN`

, the correction function; and
`correctFUN_control`

, a list of arguments passed to the
`correctFUN`

besides `A`

, `A_i`

,
`B`

, and `B_i`

.

Here we compare the `geex`

implementation of GEE with an
exchangeable correlation matrix to Fay’s `saws`

package.

The estimating functions are:

\[\begin{equation} \label{gee} \sum_{i= 1}^m \psi(\mathbf{Y}_i, \mathbf{X}_i, \beta) = \sum_{i = 1}^m \mathbf{D}_i^{\intercal} \mathbf{V}_i^{-1} (\mathbf{Y}_i - \mathbf{\mu}(\beta)) = 0 \end{equation}\]

where \(\mathbf{D}_i = \partial
\mathbf{\mu}/\partial \mathbf{\beta}\). The covariance matrix is
modeled by \(\mathbf{V}_i = \phi
\mathbf{A}_i^{0.5} \mathbf{R}(\alpha) \mathbf{A}_i^{0.5}\). The
matrix \(\mathbf{R}(\alpha)\) is the
“working” correlation matrix, which in this example is an exchangeable
matrix with off diagonal elements \(\alpha\). The matrix \(\mathbf{A}_i\) is a diagonal matrix with
elements containing the variance functions of \(\mu\). The equations in \(\eqref{gee}\) can be translated into an
`eeFUN`

as:

```
gee_eefun <- function(data, formula, family){
X <- model.matrix(object = formula, data = data)
Y <- model.response(model.frame(formula = formula, data = data))
n <- nrow(X)
function(theta, alpha, psi){
mu <- family$linkinv(X %*% theta)
Dt <- t(X) %*% diag(as.numeric(mu), nrow = n)
A <- diag(as.numeric(family$variance(mu)), nrow = n)
R <- matrix(alpha, nrow = n, ncol = n)
diag(R) <- 1
V <- psi * (sqrt(A) %*% R %*% sqrt(A))
Dt %*% solve(V) %*% (Y - mu)
}
}
```

This `eeFUN`

treats the correlation parameter \(\alpha\) and scale parameter \(\phi\) as fixed, though some estimation
algorithms use an iterative procedure that alternates between estimating
\(\beta\) and these parameters. By
customizing the root finding function, such an algorithm could be
implemented using `geex`

[see
`vignette("geex_root_solvers")`

for more information].

We use this example to compare covariance estimates obtained from the
`gee`

function, so root finding computations are turned off.
The `gee`

\(\beta\)
estimates are used instead. Estimates for \(\alpha\) and \(\phi\) are also extracted from the
`gee`

results in `m_estimate`

. This example shows
that an `eeFUN`

can accept additional arguments to be passed
to either the outer (data) function or the inner (theta) function.
Unlike previous examples, the independent units are the types of wool,
which is set in `m_estimate`

by the `units`

argument.

```
g <- gee::gee(breaks~tension, id=wool, data=warpbreaks, corstr="exchangeable")
guo <- saws::geeUOmega(g)
```

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

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

```
results <- m_estimate(
estFUN = gee_eefun, data = warpbreaks,
units = 'wool', roots = coef(g), compute_roots = FALSE,
outer_args = list(formula = breaks ~ tension,
family = gaussian()),
inner_args = list(alpha = g$working.correlation[1,2],
psi = g$scale),
corrections = list(
bias_correction_.1 = correction(bias_correction, b = .1),
bias_correction_.3 = correction(bias_correction, b = .3)))
```

In the `geex`

output, the item `corrections`

contains a list of the results of computing each item in the
`corrections_list`

. Comparing the `geex`

results
to the results of the `saws::geeUOmega`

function, the maximum
difference in the results for any of corrected estimated covariance
matrices is 1.1e-09.