16  R packages for estimation of the causal (in)direct effects

We’ll now turn to working through a few examples of estimating the natural and interventional direct and indirect effects. Note that we will be using the medoutcon R package, which supports multiple mediators and a single binary intermediate confounder, but, if your data scenario includes multiple mediators and multiple intermediate confounders, you should consider using the HDmediation R package instead.

As our running example, we’ll use a simple data set from an observational study of the relationship between BMI and kids’ behavior, freely distributed with the mma R package on CRAN. First, let’s load the packages we’ll be using and set a seed; then, load this data set and take a quick look.

library(tidyverse)
library(sl3)
library(medoutcon)
library(medshift)
library(mma)
set.seed(429153)
# load and examine data
data(weight_behavior)
dim(weight_behavior)

# drop missing values
weight_behavior <- weight_behavior %>%
  drop_na() %>%
  as_tibble()
weight_behavior

The documentation for the data set describes it as a “database obtained from the Louisiana State University Health Sciences Center, New Orleans, by Dr. Richard Scribner. He explored the relationship between BMI and kids’ behavior through a survey at children, teachers and parents in Grenada in 2014. This data set includes 691 observations and 15 variables.” Note that the data set contained several observations with missing values, which we removed above to simplify the demonstration of our analytic methods. In practice, we recommend instead using appropriate corrections (e.g., imputation, inverse weighting) to fully take advantage of the observed data.

Following the motivation of the original study, we focus on the causal effects of participating in a sports team (sports) on the BMI of children (bmi), taking into consideration several mediators (snack, exercises, overweigh); all other measured covariates are taken to be potential baseline confounders.

16.1 medoutcon: Natural and interventional (in)direct effects

The data on a single observational unit can be represented \(O = (W, A, M, Y)\), with the data pooled across all participants denoted \(O_1, \ldots, O_n\), for a of \(n\) i.i.d. observations of \(O\). Recall the DAG from an earlier chapter, which represents the data-generating process:

16.1.1 Natural (in)direct effects

To start, we will consider estimation of the natural direct and indirect effects, which, we recall, are defined as follows

\[ \E[Y_{1,M_1} - Y_{0,M_0}] = \underbrace{\E[Y_{\color{red}{1},\color{blue}{M_1}} - Y_{\color{red}{1},\color{blue}{M_0}}]}_{\text{natural indirect effect}} + \underbrace{\E[Y_{\color{blue}{1},\color{red}{M_0}} - Y_{\color{blue}{0},\color{red}{M_0}}]}_{\text{natural direct effect}}. \]

  • Our medoutcon R package (1,2), which accompanies (3), implements one-step and TML estimators of both the natural and interventional (in)direct effects.
  • Both types of estimators are capable of accommodating flexible modeling strategies (e.g., ensemble machine learning) for the initial estimation of nuisance parameters.
  • The medoutcon R package uses cross-validation in initial estimation: this results in cross-validated (or “cross-fitted”) one-step and TML estimators (4–6), which exhibit greater robustness than their non-sample-splitting analogs.
  • To this end, medoutcon integrates with the sl3 R package (7), which is extensively documented in this book chapter (8,9).

16.1.2 Interlude: sl3 for nuisance parameter estimation

  • To fully take advantage of the one-step and TML estimators, we’d like to rely on flexible, data adaptive strategies for nuisance parameter estimation.

  • Doing so minimizes opportunities for model misspecification to compromise our analytic conclusions.

  • Choosing among the diversity of available machine learning algorithms can be challenging, so we recommend using the Super Learner algorithm for ensemble machine learning (10), which is implemented in the sl3 R package (7).

  • Below, we demonstrate the construction of an ensemble learner based on a limited library of algorithms, including n intercept model, a main terms GLM, Lasso (\(\ell_1\)-penalized) regression, and random forest (ranger).

    # instantiate learners
mean_lrnr <- Lrnr_mean$new()
fglm_lrnr <- Lrnr_glm_fast$new()
lasso_lrnr <- Lrnr_glmnet$new(alpha = 1, nfolds = 3)
mars_lrnr <- Lrnr_earth$new()
rf_lrnr <- Lrnr_ranger$new(num.trees = 200)

# create learner library and instantiate super learner ensemble
lrnr_lib <- Stack$new(mean_lrnr, fglm_lrnr, lasso_lrnr, mars_lrnr, rf_lrnr)
sl_lrnr <- Lrnr_sl$new(learners = lrnr_lib)

  • Of course, there are many alternatives for learning algorithms to be included in such a modeling library. Feel free to explore!

16.1.3 Efficient estimation of the natural (in)direct effects

  • Estimation of the natural direct and indirect effects requires estimation of a few nuisance parameters. Recall that these are

    • \(g(a\mid w)\), which denotes \(\P(A=a \mid W=w)\)
    • \(h(a\mid m, w)\), which denotes \(\P(A=a \mid M=m, W=w)\)
    • \(b(a, m, w)\), which denotes \(\E(Y \mid A=a, M=m, W=w)\)

  • While we recommend the use of Super Learning, we opt to instead estimate all nuisance parameters with Lasso regression below (to save computational time).

  • Now, let’s use the medoutcon() function to estimate the natural direct effect:

    # compute one-step estimate of the natural direct effect
nde_onestep <- medoutcon(
  W = weight_behavior[, c("age", "sex", "race", "tvhours")],
  A = (as.numeric(weight_behavior$sports) - 1),
  Z = NULL,
  M = weight_behavior[, c("snack", "exercises", "overweigh")],
  Y = weight_behavior$bmi,
  g_learners = lasso_lrnr,
  h_learners = lasso_lrnr,
  b_learners = lasso_lrnr,
  effect = "direct",
  estimator = "onestep",
  estimator_args = list(cv_folds = 5)
)
summary(nde_onestep)

  • We can similarly call medoutcon() to estimate the natural indirect effect:

# compute one-step estimate of the natural indirect effect
nie_onestep <- medoutcon(
  W = weight_behavior[, c("age", "sex", "race", "tvhours")],
  A = (as.numeric(weight_behavior$sports) - 1),
  Z = NULL,
  M = weight_behavior[, c("snack", "exercises", "overweigh")],
  Y = weight_behavior$bmi,
  g_learners = lasso_lrnr,
  h_learners = lasso_lrnr,
  b_learners = lasso_lrnr,
  effect = "indirect",
  estimator = "onestep",
  estimator_args = list(cv_folds = 5)
)
summary(nie_onestep)
  • From the above, we can conclude that the effect of participation on a sports team on BMI is primarily mediated by the variables snack, exercises, and overweigh, as the natural indirect effect is several times larger than the natural direct effect.
  • Note that we could have instead used the TML estimators, which have improved finite-sample performance, instead of the one-step estimators. Doing this is as simple as setting the estimator = "tmle" in the relevant argument.

16.1.4 Interventional (in)direct effects

Since our knowledge of the system under study is incomplete, we might worry that one (or more) of the measured variables are not mediators, but, in fact, intermediate confounders affected by treatment. While the natural (in)direct effects are not identified in this setting, their interventional (in)direct counterparts are, as we saw in an earlier section. Recall that both types of effects are defined by static interventions on the treatment. The interventional effects are distinguished by their use of a stochastic intervention on the mediator to aid in their identification.

Recall that the interventional (in)direct effects are defined via the decomposition:

\[ \E[Y_{1,G_1} - Y_{0,G_0}] = \underbrace{\E[Y_{\color{red}{1},\color{blue}{G_1}} - Y_{\color{red}{1},\color{blue}{G_0}}]}_{\text{interventional indirect effect}} + \underbrace{\E[Y_{\color{blue}{1},\color{red}{G_0}} - Y_{\color{blue}{0},\color{red}{G_0}}]}_{\text{interventional direct effect}} \]

  • In our data example, we’ll consider the eating of snacks as a potential intermediate confounder, since one might reasonably hypothesize that participation on a sports team might subsequently affect snacking, which then could affect mediators like the amount of exercises and overweight status.
  • The interventional direct and indirect effects may also be easily estimated with the medoutcon R package (1,2).
  • Just as for the natural (in)direct effects, medoutcon implements cross-validated one-step and TML estimators of the interventional effects.

16.1.5 Efficient estimation of the interventional (in)direct effects

  • Estimation of these effects is more complex, so a few additional nuisance parameters arise when expressing the (more general) EIF for these effects:

    • \(q(z \mid a, w)\), the conditional density of the intermediate confounders, conditional only on treatment and baseline covariates;
    • \(r(z \mid a, m, w)\), the conditional density of the intermediate confounders, conditional on mediators, treatment, and baseline covariates.

  • To estimate the interventional effects, we only need to set the argument Z of medoutcon to a value other than NULL.

  • Note that the implementation in medoutcon is currently limited to settings with only binary intermediate confounders, i.e., \(Z \in \{0, 1\}\).

  • Let’s use medoutcon() to estimate the interventional direct effect:

    # compute one-step estimate of the interventional direct effect
interv_de_onestep <- medoutcon(
  W = weight_behavior[, c("age", "sex", "race", "tvhours")],
  A = (as.numeric(weight_behavior$sports) - 1),
  Z = (as.numeric(weight_behavior$snack) - 1),
  M = weight_behavior[, c("exercises", "overweigh")],
  Y = weight_behavior$bmi,
  g_learners = lasso_lrnr,
  h_learners = lasso_lrnr,
  b_learners = lasso_lrnr,
  effect = "direct",
  estimator = "onestep",
  estimator_args = list(cv_folds = 5)
)
summary(interv_de_onestep)

  • We can similarly estimate the interventional indirect effect:

    # compute one-step estimate of the interventional indirect effect
interv_ie_onestep <- medoutcon(
  W = weight_behavior[, c("age", "sex", "race", "tvhours")],
  A = (as.numeric(weight_behavior$sports) - 1),
  Z = (as.numeric(weight_behavior$snack) - 1),
  M = weight_behavior[, c("exercises", "overweigh")],
  Y = weight_behavior$bmi,
  g_learners = lasso_lrnr,
  h_learners = lasso_lrnr,
  b_learners = lasso_lrnr,
  effect = "indirect",
  estimator = "onestep",
  estimator_args = list(cv_folds = 5)
)
summary(interv_ie_onestep)

  • From the above, we can conclude that the effect of participation on a sports team on BMI is largely through the interventional indirect effect (i.e., through the pathways involving the mediating variables) rather than via its direct effect.

  • Just as before, we could have instead used the TML estimators, instead of the one-step estimators. Doing this is as simple as setting the estimator = "tmle" in the relevant argument.