gapclosing: An R package to estimate gap closing estimands

knitr::opts_chunk$set(dpi = 300, fig.width = 6.5, fig.height = 3, out.width = "650px")
t0 <- Sys.time()

Gaps across social categories like race, class, and gender are important to understand. We would like to know whether there is anything we can do to close these gaps. What if we intervened to reduce incarceration or increase access to education? Would those interventions close gaps across categories of race, class, or gender?

These types of questions are at the core of a growing literature in epidemiology addresses these questions with techniques for causal decomposition analysis (@vanderweele2014, @jackson2018, @jackson2020). The paper which accompanies this software package introduces these methods to a broader set of social scientists.

Lundberg, Ian. Forthcoming. "The gap-closing estimand: A causal approach to study interventions that close disparities across social categories." Sociological Methods and Research. Preprint available at

This package provides software to support inquiry into gap-closing estimands. A guiding principle is to partition the tasks the human user must do from the tasks that the package can automate.

The human user must carry out these tasks:

Given those inputs, the package will automatically:

The structure of the vignette is as follows.

  1. Data structure
  2. Quick overview: Functionality
  3. Detailed procedures: Customize your call
    • Define the intervention
    • Make causal assumptions for identification
    • Specify a treatment model and/or an outcome model
  4. Estimate by machine learning: Two examples
  5. Estimate with stochastic assignment rules
  6. Conclusion

1. Data structure

In a data frame data, we have a gap-defining category such as race, gender, or class. We have a binary treatment variable that could have counterfactually been different for any individual. We want to know the degree to which an intervention to change the treatment would close gaps across the categories.

These boxes will present an example.

Example. Suppose we have the following data.

  • $X$ (category): Category of interest, taking values {A, B, C}
  • $T$ (treatment): Binary treatment variable, taking values 0 and 1
  • $L$ (confounder): A continuous confounding variable, Uniform(-1,1)
  • $Y$ (outcome): A continuous outcome variable, conditionally normal
simulated_data <- generate_simulated_data(n = 1000)

2. Quick overview: Functionality

The gapclosing() function estimates gaps across categories and the degree to which they would close under the specified counterfactual_assignments of the treatment.

estimate <- gapclosing(
  data = simulated_data,
  counterfactual_assignments = 1,
  outcome_formula = formula(outcome ~ confounder + category*treatment),
  treatment_formula = formula(treatment ~ confounder + category),
  category_name = "category",
  se = TRUE,
  # You can process the bootstrap in parallel with as many cores as available
  parallel_cores = 2

By default, this function will do the following:

In this example, the plot(estimate) function produces the following visualization. The factual outcomes are unequal across categories, but the counterfactual outcomes are roughly equal. In this simulated setting, the intervention almost entirely closes the gaps across the categories.

plots <- plot(estimate, return_plots = TRUE)
print(plots[[1]] +
        ggtitle("First result of a call to plot()"))

The disparityplot() function lets us zoom in on the factual and counterfactual disparity between two categories, of interest. In this case, we see that the intervention lifts outcomes in category A to be more comparable to category B. A disparityplot is a ggplot2 object and can be customized by passing additional layers.

disparityplot(estimate, category_A = "A", category_B = "B") +
  ggtitle("A disparityplot()")

The summary function will print estimates, standard errors, and confidence intervals for all of these results.

old <- options()
options(width = 300)

3. Detailed procedures: Customize your call

This section provides a more detailed overview of the use of gapclosing(). It is structured from the perspective of the three key tasks for the researcher: define the intervention, make causal assumptions for identification, and specify a treatment and/or outcome model. Along the way, this section introduces many of the possible arguments in a gapclosing() call.

Step 1. Define the intervention

To answer a gap-closing question, we first need to define what that intervention would be. To what treatment value would units be counterfactually assigned? There are several options.

Example. We are interested in the disparity across populations defined by category that would persist under counterfactual assignment to set treatment to the value 1. counterfactual_assignments = 1

Step 2. Make causal assumptions for identification

The package does help you with this step. Gap-closing estimands involve unobserved potential outcomes. Because they are unobserved, the data cannot tell us which variables are needed for estimation. Instead, that is a conceptual choice to be carried out with tools like Directed Acyclic Graphs (DAGs). See the accompanying @lundberg2021 paper for more on identification.

Example. Assume that the set of variables ${X,L}$ is a sufficient conditioning set to identify the gap-closing estimand. Formally, this requires us to assume that within each stratum of $X$ and $L$ the expected value of the potential outcome $Y(1)$ is the same as the expected value among units who factually have $T = 1$ within those strata. $$\mathbb{E}(Y(1)\mid X, L) = \mathbb{E}(Y\mid X, L, T = 1)$$ DAGs are a good way to reason about this assumption: in this example, conditioning (depicted by boxes) on category and confounder is sufficient to identify the causal effect of treatment (blue edge in the DAG), because doing so blocks all backdoor paths between the treatment and the outcome. Notably, the gap-closing estimand makes no claims about the causal effect of category since the counterfactual is defined over treatment only.

for_DAG <- data.frame(label = c("Category","Confounder","Treatment","Outcome","Unobserved"),
                      x = c(1,2,3,4,1),
                      y = c(1,1,1,1,0))
for_DAG %>%
  ggplot(aes(x = x, y = y, label = label)) +
  # Edges coming out of category
  annotate(geom = "segment",
         x = 1, y = 1, xend = 1.55, yend = 1,
         arrow = arrow(length = unit(.1,"in"))) +
  annotate(geom = "curve",
           x = 1, y = 1, xend = c(2.9,3.9), yend = c(1.2,1.2),
           arrow = arrow(length = unit(.1,"in")),
           curvature = -.3) +
  # Edges coming out of confounder
  annotate(geom = "segment",
           x = 2, y = 1, xend = 2.65, yend = 1,
           arrow = arrow(length = unit(.1,"in"))) +
  annotate(geom = "curve",
           x = 2, y = 1, xend = 3.75, yend = .82,
           arrow = arrow(length = unit(.1,"in")),
           curvature = .25) +
  # Edges coming out of treatment
  annotate(geom = "segment",
           x = 3, y = 1, xend = 3.65, yend = 1,
           color = "blue", size = 1.05,
           arrow = arrow(length = unit(.13,"in"))) +
  # Edges coming out of unobserved
  annotate(geom = "segment",
           x = 1, y = 0, xend = c(1,1.9), yend = .8,
           arrow = arrow(length = unit(.1,"in"))) +
  annotate(geom = "curve",
           x = 1, y = 0, xend = c(3.9), yend = .8,
           arrow = arrow(length = unit(.1,"in")),
           curvature = .25) +
  geom_label(color = "white") +
  geom_label(data = for_DAG %>% filter(label %in% c("Category","Confounder"))) +
  geom_text(data = for_DAG %>% filter(label %in% c("Unobserved","Treatment","Outcome"))) +
  ylim(c(-.3,2)) +
  xlim(c(0,5)) +

Once we select a sufficient conditioning set, those predictors will appear in both the treatment and/or the outcome model used for estimation.

Step 3. Specify a treatment model and/or an outcome model

We need to estimate one or both of (1) the probability of treatment given confounders and (2) the conditional mean of the outcome given treatment and confounders. We do that by providing one or both of the following.

treatment_formula = formula(outcome ~ confounder + category)

outcome_formula = formula(outcome ~ confounder + category*treatment)

Whether treatment_formula or outcome_formula is left NULL will determine the estimation procedure.

Choose an estimation algorithm

The treatment_formula and outcome_formula are handed to treatment_algorithm and outcome_algorithm, which can take the following values.

treatment_algorithm = "glm"
outcome_algorithm = "lm"

If the data are a sample from a population selected with unequal probabilities, you can also use the weight_name option to pass estimation functions the name of the sampling weight (a variable in data proportional to the inverse probability of sample inclusion). If omitted, a simple random sample is assumed.

Why doubly robust? A side note

Doubly-robust estimation yields advantages that can be conceptualized in two ways.

Although double robustness has strong mathematical properties, in any given application with a finite sample it is possible that treatment or outcome modeling could outperform doubly-robust estimation. Therefore, the package supports all three approaches.

Why sample splitting? Another side note

Taking the bias-correction view of double robustness above, it is clear that sample splitting affords a further opportunity for improvement: if you learn an outcome model and estimate its average bias on the same sample, you might get a poor estimate of the bias. For this reason, one should consider using one sample (which I call data_learn) to learn the prediction functions and another sample (which I call data_estimate) to estimate the bias and aggregate to an estimate of the estimand.

In particular, the option sample_split = "cross_fit" allows the user to specify that estimation should proceed by a cross-fitting procedure which is analogous to cross-validation.

  1. Split the sample into folds $f = 1,\dots,$n_folds (default here is n_folds = 2)
  2. Use all folds except $f'$ to estimate the treatment and outcome models
  3. Aggregate to an estimate using the predictions in $f'$
  4. Average the estimate that results from (2) and (3) repeated n_folds times with each fold playing the role of $f'$ in turn

This is the procedure that @chernozhukov2018 argue is critical to double machine learning for causal estimation, although this type of sample splitting is not new [@bickel1982].

If sample_split = "cross_fit", the default is to conduct 2-fold cross-fitting, but this can be changed with the n_folds argument. The user can also specify their own vector folds of fold assignments of length nrow(data), if there is something about the particular setting that would make a manual fold assignment preferable.

Example. (this is the default and can be left implicit)
sample_split = "one_sample"

Produce standard errors

The package supports bootstrapped standard error estimation. The procedure with bootstrap_method = "simple" (the default) is valid when the data are a simple random sample from the target population. In this case, each bootstrap iteration conducts estimation on a resampled dataset selected with replacement with equal probabilities. The standard error is calculated as the standard deviation of the estimate across bootstrap samples, and confidence intervals are calculated by a normal approximation.

Example. (line 2 is the default and can be left implicit)
se = T bootstrap_samples = 1000

In some settings, the sample size may be small and categories or treatments of interest may be rare. In these cases, it is possible for one or more simple bootstrap samples to contain zero cases in some (treatment $\times$ category) cell of interest. To avoid this problem, bootstrap_method = "stratified" conducts bootstrap resampling within blocks defined by (treatment $\times$ category). This procedure is valid if you assume that the data are selected at random from the population within these strata, so that across repeated samples from the true population the proportion in each stratum would remain the same.

Many samples are not simple random samples. In complex sample settings, users should implement their own standard error procedures to accurately capture sampling variation related to how their data were collected. The way the data were collected could motivate a resampling strategy to mimic the sources of variation in that sampling process, which the user can implement manually by calling gapclosing to calculate a point estimate on each resampled dataset with se = FALSE.

4. Estimate by machine learning: Two examples

Suppose you want to relax parametric functional form assumptions by plugging in a machine learning estimator. As the user, this simply involves changing the arguments to the gapclosing() function.

Estimation with Generalized Additive Models (GAMs) via mgcv

Perhaps you are concerned about linearity assumptions: the continuous confounder, for instance, might actually have a nonlinear association with the outcome. We know the truth is linear in this simulated example, but in practice you would never know. You can address this concern by estimating with a GAM, using the s() operator from mgcv for smooth terms (see @mgcv).

estimate_gam <- gapclosing(
  data = simulated_data,
  counterfactual_assignments = 1,
  outcome_formula = formula(outcome ~ s(confounder) + category*treatment),
  treatment_formula = formula(treatment ~ s(confounder) + category),
  category_name = "category",
  treatment_algorithm = "gam",
  outcome_algorithm = "gam",
  sample_split = "cross_fit"
  # Note: Standard errors with `se = TRUE` are supported.
  # They are omitted here only to speed vignette build time.

Estimation with random forests via ranger

Perhaps you are concerned that the true treatment probability and expected outcome functions have many interactions among the predictors. You can set treatment_algorithm and outcome_algorithm to "ranger" to estimate via the ranger function in the ranger package [@ranger].

One aspect of the way gapclosing() operationalizes ranger() is unique out of all the estimation algorithm options. When you choose a random forest, it is because you believe there are many important interactions. Some of the most important interactions may be between the treatment and the other predictors. Therefore, outcome_algorithm = ranger enforces those interactions by estimating the outcome model separately for treated and control units. For this reason, when outcome_algorithm = ranger there is no need to include the treatment variable explicitly in the outcome_formula.

estimate_ranger <- gapclosing(
  data = simulated_data,
  counterfactual_assignments = 1,
  outcome_formula = formula(outcome ~ confounder + category),
  treatment_formula = formula(treatment ~ confounder + category),
  category_name = "category",
  treatment_algorithm = "ranger",
  outcome_algorithm = "ranger",
  sample_split = "cross_fit"
  # Note: Standard errors with `se = TRUE` are supported.
  # They are omitted here only to speed vignette build time.

Estimates from these three algorithms are roughly the same

In this simulation, the GLM models are correctly specified and there are no nonlinearities or interactions for the machine learning approaches to learn. In this case, the sample size is large enough that those approaches correctly learn the linear functional form, and all three estimation strategies yield similar estimates.

glm_plot <- plot(estimate, return_plots = T)[[1]]
gam_plot <- plot(estimate_gam, return_plots = T)[[1]]
ranger_plot <- plot(estimate_ranger, return_plots = T)[[1]]
glm_plot +
    ggtitle("GLM estimate") +
    ylim(c(-1.5,1.5)) +
  theme(legend.position = "bottom")
gam_plot +
    ggtitle("GAM estimate") +
    ylim(c(-1.5,1.5)) +
  theme(legend.position = "bottom")
ranger_plot +
  ggtitle("Random forest estimate") +
  ylim(c(-1.5,1.5)) +
  theme(legend.position = "bottom")

Note that confidence intervals for GAM and random forest can also be generated with SE = TRUE, which is turned off here only to speed vignette build time.

A warning about complex estimators

The assumptions of a parametric model are always doubtful, leading to a common question of whether one should always use a more flexible machine learning approach like ranger. In a very large sample, a flexible learner would likely be the correct choice. In the sample sizes of social science settings, the amount of data may sometimes be insufficient for these algorithms to discover a complex functional form. When the parametric assumptions are approximately true, the parametric estimators may have better performance in small sample sizes. What counts as "small" and "large" is difficult to say outside of any specific setting.

5. Estimate with stochastic assignment rules

The examples above focused on estimation for a fixed treatment assignment: assign to treatment 1 with probability 1. But we might also want to know about the gap-closing estimand if we assigned people to treatment stochastically with some probability between 0 and 1. The counterfactual_assignments argument can handle this possibility.

For example, consider the gap-closing estimand if assigned to treatment 1 with each probability .75.

estimate_stochastic <- gapclosing(
  data = simulated_data,
  counterfactual_assignments = .75,
  outcome_formula = formula(outcome ~ confounder + category*treatment),
  treatment_formula = formula(treatment ~ confounder + category),
  category_name = "category"

The disparity between categories A and B under that stochastic intervention (0.75 probability of treatment = 1) is estimated to be r round((estimate_stochastic$counterfactual_disparities %>% filter(category == "B - A"))$estimate,2), whereas under the previous deterministic intervention to assign treatment to the value 1 the disparity would be r round((estimate$counterfactual_disparities %>% filter(category == "B - A"))$estimate,2). This illustrates an important point: the gap-closing estimand can be different depending on the counterfactual assignment rule, as the figure below shows for counterfactuals in which treatment is assigned with probabilities ranging from 0 to 1.

counterfactual_assignments_values <- seq(0,1,.25)
many_stochastic_estimates <- rep(NA, length(counterfactual_assignments_values))
for (i in 1:length(counterfactual_assignments_values)) {
  estimate_case <- gapclosing(
    data = simulated_data,
    counterfactual_assignments = counterfactual_assignments_values[i],
    outcome_formula = formula(outcome ~ confounder + category*treatment),
    treatment_formula = formula(treatment ~ confounder + category),
    category_name = "category"
  estimate_we_want <- estimate_case$change_disparities %>%
    filter(change_type == "proportional") %>%
    filter(category == "B - A")
  many_stochastic_estimates[i] <- estimate_we_want$estimate
data.frame(counterfactual_assignments = counterfactual_assignments_values,
           estimate = many_stochastic_estimates) %>%
  ggplot(aes(x = counterfactual_assignments,
             y = estimate)) +
  geom_hline(yintercept = c(0,1),
             color = "gray", linetype = "dashed") +
  geom_line(color = "gray") +
  geom_point(aes(color = factor(counterfactual_assignments == .75))) +
  xlab("Counterfactual Treatment Assignment Probability") +
  scale_y_continuous(name = "Percent of B - A Gap Closed",
                     labels = function(x) paste0(round(100*x),"%")) +
  annotate(geom = "text",
           x = 1, y = c(0,1),
           label = c("No change to disparity",
                     "Disparity completely eliminated"),
           color = "darkgray", size = 3, hjust = 1, vjust = 1.5) +
  annotate(geom = "text",
           x = .1, y = .3, 
           label = "Disparities eliminated under various\nstochastic\ninterventions",
           color = "gray", size = 3, hjust = 0, vjust = 1) +
  annotate(geom = "text",
           x = .75, y = .3,
           label = "Gap-Closing Estimand\nCalculated Above",
           color = "blue", size = 3, vjust = -.5) +
  scale_color_manual(values = c("gray","blue")) +
  theme_bw() +
  theme(panel.grid = element_blank(),
        legend.position = "none")

Your stochastic assignments can also be different for different people. For example, suppose we assign those in Category A to treatment 1 with probability .5, those in Category B to treatment with probability .4, and those in Category C to treatment with probability .3. In this case, counterfactual_assignments will be set to a vector of length nrow(data).

our_assignments <- case_when(simulated_data$category == "A" ~ .5,
                             simulated_data$category == "B" ~ .4, 
                             simulated_data$category == "C" ~ .3)
estimate_stochastic <- gapclosing(
  data = simulated_data,
  counterfactual_assignments = our_assignments,
  outcome_formula = formula(outcome ~ confounder + category*treatment),
  treatment_formula = formula(treatment ~ confounder + category),
  category_name = "category"

That intervention would close the B - A gap by r paste0(round(100 * (estimate_stochastic$change_disparities %>% filter(category == "B - A" & change_type == "proportional"))$estimate),"%").

6. Conclusion

The gapclosing package is designed to support inquiry into gap closing estimands, thus promoting new understanding about interventions that can close gaps across social categories. The goal of the package is to automate technical tasks (sample splitting, aggregation to doubly robust estimates, visualization), thus freeing the researcher to devote more attention to scientific tasks like defining the intervention and making causal assumptions.

If you use this package and find a bug, it would be most helpful if you would create an issue on GitHub. Suggestions for additional features are also welcome.

This vignette was compiled on r Sys.time().


Try the gapclosing package in your browser

Any scripts or data that you put into this service are public.

gapclosing documentation built on Oct. 11, 2021, 9:07 a.m.