R/est_sub.R

In nhejazi/medshift: Causal mediation analysis for stochastic interventions

#' Substitution estimator
#'
#' @param data A \code{data.table} containing the observed data, with columns
#'  in the order specified by the NPSEM (Y, Z, A, W), with column names set
#'  appropriately based on the original input data. Such a structure is merely
#'  a convenience utility to passing data around to the various core estimation
#'  routines and is automatically generated by \code{\link{medshift}}.
#' @param delta A \code{numeric} value indicating the degree of shift in the
#'  intervention to be used in defining the causal quantity of interest. In the
#'  case of binary interventions, this takes the form of an incremental
#'  propensity score shift, acting as a multiplier of the odds with which a
#'  given observational unit receives the intervention (EH Kennedy, 2018, JASA;
#'  <doi:10.1080/01621459.2017.1422737>).
#' @param g_learners A \code{\link[sl3]{Stack}} (or other learner class that
#'   inherits from \code{\link[sl3]{Lrnr_base}}), containing a single or set of
#'   instantiated learners from \pkg{sl3}, to be used in fitting the propensity
#'   score, i.e., g = P(A | W).
#' @param m_learners A \code{\link[sl3]{Stack}} (or other learner class that
#'   inherits from \code{\link[sl3]{Lrnr_base}}), containing a single or set of
#'   instantiated learners from \pkg{sl3}, to be used in fitting the outcome
#'   regression, i.e., m(A, Z, W).
#' @param w_names A \code{character} vector of the names of the columns that
#'  correspond to baseline covariates (W). The input for this argument is
#'  automatically generated by a call to the wrapper function \code{medshift}.
#' @param z_names A \code{character} vector of the names of the columns that
#'  correspond to mediators (Z). The input for this argument is automatically
#'  generated by a call to the wrapper function \code{medshift}.
#' @param ... Other arguments currently ignored.
est_substitution <- function(data,
                             delta,
                             g_learners,
                             m_learners,
                             w_names,
                             z_names,
                             ...) {
  # estimate propensity score
  g_out <- fit_g_mech(
    data = data, delta = delta,
    learners = g_learners, w_names = w_names
  )

  # fit regression for incremental propensity score intervention
  m_out <- fit_m_mech(
    data = data, learners = m_learners,
    z_names = z_names, w_names = w_names
  )

  # compute Dzw component of EIF using convenience function
  Dzw_groupwise <- compute_Dzw(g_output = g_out, m_output = m_out)

  # compute estimator
  estim_sub <- mean(Dzw_groupwise$dzw_cntrl) + mean(Dzw_groupwise$dzw_treat)

  # output
  estim_sub_out <- list(theta = estim_sub, type = "substitution")
  return(estim_sub_out)
}

###############################################################################

#' Construct joint Z-W component of efficient influence function
#'
#' @param g_output Object containing results from fitting the propensity score
#'  regression, as produced by a call to \code{\link{fit_g_mech}}.
#' @param m_output Object containing results from fitting the outcome
#'  regression, as produced by a call to \code{\link{fit_m_mech}}.
#'
#' @keywords internal
compute_Dzw <- function(g_output, m_output) {
  # get g components from output for that nuisance parameter
  g_shifted_A1 <- g_output$g_est$g_pred_shifted_A1
  g_shifted_A0 <- g_output$g_est$g_pred_shifted_A0

  # get m components from output for that nuisance parameter
  m_pred_A1 <- m_output$m_est$m_pred_A1
  m_pred_A0 <- m_output$m_est$m_pred_A0

  # compute component Dzw from nuisance parameters
  Dzw_A1 <- g_shifted_A1 * m_pred_A1
  Dzw_A0 <- g_shifted_A0 * m_pred_A0

  # output as simple list
  return(list(
    dzw_cntrl = Dzw_A0,
    dzw_treat = Dzw_A1
  ))
}