R/correct_by_funs.R

In geex: An API for M-Estimation

#------------------------------------------------------------------------------#
# correct_by_** description:
# Functions that correct the output of compute_matrices *by* some manipulation.
# Also includes necessary (non-exported) internal functions.
#------------------------------------------------------------------------------#

#------------------------------------------------------------------------------#
# Make corrections to basic sandwich estimators
#
# @param mats list of matrices resulting from \code{\link{compute_matrices}}
# @param corrections a list with one sublist for each correction to make. The
# sublist must contain at least \code{correctFUN} and optionally additional
# arguments in \code{correctFUN_control}.
# @return a list with corrected values. The list retains the names of the input
#  \code{corrections}
# @export
#------------------------------------------------------------------------------#
make_corrections <- function(components, corrections){
  lapply(corrections, function(correction){
     correct_by(.components = components, .correct_control = correction)
  })
}

###############################################################################
# correct_by functions provided by geex ####
###############################################################################

#------------------------------------------------------------------------------#
#' Correct sandwich components
#'
#' Modifies the matrices in a \code{\linkS4class{sandwich_components}} object
#' using the function and options in a \code{\linkS4class{correct_control}} object.
#' The function \code{\link{correction}} is a utility for creating
#' \code{\linkS4class{correct_control}} objects.
#'
#' See the finite sample corrections vignette for further examples.
#'
#' @param .components an object of class \code{\linkS4class{sandwich_components}}
#' @param .correct_control an object of class \code{\linkS4class{correct_control}}
#'
#' @return the result of \code{.FUN} in \code{.correct_control}.
#' @seealso \code{\link{fay_bias_correction}} and \code{\link{fay_df_correction}}
#' for corrections provided by \code{geex}
#' @examples
#' myee <- function(data){
#'    function(theta){
#'     c(data$Y1 - theta[1],
#'      (data$Y1 - theta[1])^2 - theta[2])
#'    }
#'  }
#' mybasis <- create_basis(
#'    estFUN = myee,
#'    data   = geexex)
#' mats <- estimate_sandwich_matrices(mybasis, .theta = c(5.04, 10.04))
#' correct_by(mats,
#'    .correct_control =  correction(fay_bias_correction, b = .75))
#' @export
#------------------------------------------------------------------------------#

correct_by <- function(.components, .correct_control){
  f <- grab_FUN(.correct_control)
  opts <- grab_options(.correct_control)
  do.call(f, args = append(list(components = .components), opts))
}

#------------------------------------------------------------------------------#
#' Creates a correct_control object
#'
#' @param FUN a correction to perform. \code{components} must be the
#' first argument
#' @param ... additional arguments passed to \code{FUN}
#' @return a \code{\linkS4class{correct_control}} object
#' @export
#' @examples
#' correction(FUN = fay_bias_correction, b = 0.75)
#------------------------------------------------------------------------------#

correction <- function(FUN, ...){
  dots <- list(...)
  methods::new(Class="correct_control",
      .FUN = FUN,
      .options   = dots)
}

#------------------------------------------------------------------------------#
#' Correct sandwich variance estimator byFay's bias correction
#'
#' Computes the bias corrected sandwich covariance matrix described in Fay and
#' Graubard (2001). See \code{vignette("05_finite_sample_corrections", package = "geex")}
#' for further information.
#'
#' @param components an object of class \code{\linkS4class{sandwich_components}}
#' @param b a numeric value < 1. Defaults to 0.75 as in Fay.
#' @return a corrected covariance matrix
#' @references Fay, M. P., & Graubard, B. I. (2001). Small-Sample adjustments for
#' Wald-type tests using sandwich estimators. Biometrics, 57(4), 1198-1206
#' @export
#' @examples
#' # This example demonstrates usage of the corrections, not a meaningful application
#' myee <- function(data){
#'  function(theta){
#'    c(data$Y1 - theta[1],
#'    (data$Y1 - theta[1])^2 - theta[2])
#'   }
#' }
#'
#' results <- m_estimate(
#'    estFUN = myee,
#'    data = geexex,
#'    root_control = setup_root_control(start = c(1,1)),
#'    corrections  = list(
#'      bias_correction_.1 = correction(fay_bias_correction, b = .1),
#'      bias_correction_.3 = correction(fay_bias_correction, b = .3))
#'    )
#'
#' get_corrections(results)
#------------------------------------------------------------------------------#

fay_bias_correction <- function(components, b = 0.75){
  fay_bias_correction_partial(components, b = b) ->  corrected_matrices

  compute_sigma(A = grab_bread(components), B = corrected_matrices$Bbc)
}

#------------------------------------------------------------------------------#
#' Correct sandwich variance inference by Fay's degrees of freedom correction
#'
#' Computes the degrees of freedom correction described in Fay and
#' Graubard (2001). See \code{vignette("05_finite_sample_corrections", package = "geex")}
#' for further information.
#'
#' @inheritParams fay_bias_correction
#' @param L a k x p matrix where p is the dimension of theta
#' @param version either 1 or 2, corresponding to hat(d) or tilde(d), respectively
#' @references Fay, M. P., & Graubard, B. I. (2001). Small-Sample adjustments for
#' Wald-type tests using sandwich estimators. Biometrics, 57(4), 1198-1206
#' @return a scalar corresponding to the estimated degrees of freedom
#' @export
#' @examples
#'
#' # This example demonstrates usage of the corrections, not a meaningful application
#' myee <- function(data){
#'  function(theta){
#'    c(data$Y1 - theta[1],
#'    (data$Y1 - theta[1])^2 - theta[2])
#'   }
#' }
#'
#' results <- m_estimate(
#'    estFUN = myee,
#'    data = geexex,
#'    root_control = setup_root_control(start = c(1,1)),
#'    corrections  = list(
#'      df_correction1 = correction(fay_df_correction,
#'                         b = .75, L = c(0, 1), version = 1 ),
#'      df_correction2 = correction(fay_df_correction,
#'                         b = .75, L = c(0, 1), version = 2 ))
#'    )
#'
#' get_corrections(results)
#------------------------------------------------------------------------------#

fay_df_correction <- function(components, b = .75, L, version){

  ## Prepare necessary matrices ##
  bias_mats <- fay_bias_correction_partial(components, b = b)
  A   <- grab_bread(components)
  A_i <- grab_bread_list(components)
  df_prep   <- df_correction_prep(L = L, A = A, A_i = A_i, H_i = bias_mats$H_i)

  ## Compute DF corrections ##
  if(version == 1){
    out <- df_correction_1(df_prep$A_d, df_prep$C)
  } else if(version == 2){
    out <- df_correction_2(A = A, A_i = A_i, C = df_prep$C, L = L, Bbc = bias_mats$Bbc)
  }
  out
}

###############################################################################
# functions internal to provided correct_by ####
###############################################################################


#------------------------------------------------------------------------------#
# Compute the matrices necessary for Fay's bias correction
#
# @inheritParams fay_bias_correction
# @return a list with
# \itemize{
# \item H_i - the H_i matrix
# \item Bbc - the 'bias corrected' meat matrix
# }
# @export
#------------------------------------------------------------------------------#

fay_bias_correction_partial <- 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(Ai){
    diag( (1 - pmin(b, diag(Ai %*% Ainv) ))^(-0.5) )
  })

  stopifnot(length(B_i) == length(H_i))

  Bbc_i <- lapply(seq_along(B_i), function(i){
    H_i[[i]] %*% B_i[[i]] %*% H_i[[i]]
  })
  Bbc <- compute_sum_of_list(Bbc_i) #apply(simplify2array(Bbc_i), 1:2, sum)

  list(H_i = H_i, Bbc = Bbc)
}

#------------------------------------------------------------------------------#
# Preparations for Fay's degrees of freedom corrections
#
# @param H_i a list of bias adjusted matrices
# @inheritParams fay_bias_correction
# @inheritParams fay_df_correction
# @export
#------------------------------------------------------------------------------#
df_correction_prep <- function(L, A, A_i, H_i){
  p <- ncol(A)
  m <- length(A_i)

  Ainv <- solve(A)
  II   <- diag(1, p*m)
  AA   <- do.call(rbind, A_i)
  calI <- do.call(cbind, args = lapply(1:m, function(i) diag(1, p) ))
  G    <- II - (AA %*% solve(A) %*% calI)

  M_i  <- lapply(H_i, function(mat){
    mat %*% Ainv %*% L %*% t(L) %*% t(Ainv) %*% mat
  })
  M    <- Matrix::bdiag(M_i)

  C    <- t(G) %*% M %*% G

  A_diag  <- Matrix::bdiag(A_i)

  list(A_d = A_diag, C = C)
}

#------------------------------------------------------------------------------#
# Estimate Fay's degrees of freedom corrections 1
#
# @param A_d see \code{\link{df_correction_prep}}
# @param C see \code{\link{df_correction_prep}}
#
# @export
#------------------------------------------------------------------------------#
df_correction_1 <- function(A_d, C){
  estimate_df(A = A_d, C = C)
}

#------------------------------------------------------------------------------#
# Estimate Fay's degrees of freedom correction 2
#
# @inheritParams fay_bias_correction
# @inheritParams fay_df_correction
# @param C same as the C matrix in the Fay (2001) notation Section 2.3
# @param Bbc the bias corrected "B" matrix
#
# @export
#------------------------------------------------------------------------------#
df_correction_2 <- function(A, A_i, C, L, Bbc){
  w_i  <- lapply(seq_along(A_i), function(i) {
    # exclude the ith element
    Oi <- apply(simplify2array(A_i[-i]), 1:2, sum)
    t(L) %*% (solve(Oi) - solve(A) ) %*% L
  })
  wbar <- sum(unlist(w_i))

  Abc_i <- lapply(w_i, function(w){
    as.numeric(w/wbar) * Bbc
  })
  Abc  <- Matrix::bdiag(Abc_i)

  estimate_df(A = Abc, C = C)
}

#------------------------------------------------------------------------------#
# Estimate Fay's degrees of freedom corrections
#
# @inheritParams df_correction_2
# @inheritParams fay_df_correction
# @export
#------------------------------------------------------------------------------#

estimate_df <- function(A, C){
  AC <- A %*% C
  (sum(Matrix::diag(AC)))^2 / sum(Matrix::diag(AC %*% AC))
}