Finite sample correction API in geex

Introduction

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.

Corrections included with 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.

Bias correction example

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.

Comparision to saws package

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 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)

library(geex)
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.