18 Bootstrap
This vignette demonstrates how to use bootstrap methods for inference with the marginaleffects
package.
18.1 Parallelization
The bootstrap is an incredibly powerful and flexible method to quantify the uncertainty around estimates. Unfortunately, it can be quite computationally expensive. In some cases, parallelization can significantly speed up the process. This is most beneficial when bootstrap sample size (R
) is large (≥ 1000); model fitting is computationally expensive; you have multiple CPU cores available.
Under the hood, marginaleffects
uses the future
framework to enable parallel processing. This requires calling the plan()
function, and specifying a couple global options:
-
marginaleffects_parallel_inferences
: Set toTRUE
to enable parallel processing -
marginaleffects_parallel_packages
: Character vector of package names used to fit the model.
18.2 Simple example
All the objects produced by the marginaleffects
can be passed to the inferences()
function to compute bootstrap confidence intervals. For example, fit a simple model to the iris
dataset and compute average comparisons.
mod <- lm(Sepal.Length ~ Sepal.Width + Species, data = iris)
avg_comparisons(mod)
Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
Sepal.Width +1 0.804 0.106 7.56 <0.001 44.5 0.595 1.01
Species versicolor - setosa 1.459 0.112 13.01 <0.001 126.2 1.239 1.68
Species virginica - setosa 1.947 0.100 19.47 <0.001 277.9 1.751 2.14
Type: response
The default confidence intervals reported above were computed using the delta method. We can obtain bootstrap-based estimates by piping the command into inferences()
:
avg_comparisons(mod) |> inferences(method = "rsample")
Term Contrast Estimate 2.5 % 97.5 %
Sepal.Width +1 0.804 0.591 1.01
Species versicolor - setosa 1.459 1.266 1.65
Species virginica - setosa 1.947 1.764 2.14
Type: response
This also works with other marginaleffects
functions:
avg_predictions(mod, by = "Species") |> inferences(method = "rsample")
Species Estimate 2.5 % 97.5 %
setosa 5.01 4.91 5.11
versicolor 5.94 5.79 6.08
virginica 6.59 6.42 6.77
Type: response
hypotheses(mod, hypothesis = ratio ~ sequential) |> inferences(method = "rsample")
Hypothesis Estimate 2.5 % 97.5 %
(Sepal.Width) / ((Intercept)) 0.357 0.205 0.658
(Speciesversicolor) / (Sepal.Width) 1.815 1.519 2.260
(Speciesvirginica) / (Speciesversicolor) 1.335 1.189 1.503
The inferences()
function supports multiple bootstrap methods, and includes arguments to control some aspects of the procedure, including R
to determine the number of resamples. Additional (undocumented) arguments will be passed to the underlying bootstrap function, such as rsample::bootstraps()
. This allows users to supply arguments such as strata
or pool
. See ?rsample::bootstraps
for details.
18.3 Custom Estimation Procedures
Sometimes analysts need to use the bootstrap to compute uncertainy around more complex estimation procedures that go beyond the standard marginal effects workflow. For example, one may want to proceed in three steps:
- Fit a logistic regression model to the treatment in order to estimate propensity scores.
- Fit a linear model to the outcome using the estimated propensity scores as weights.
- Compute the average treatment effect using G-computation and weights with the
avg_comparisons()
function.
The bootstrap allows us to take into account the uncertainty in all stages of estimation. This can be useful in a variety of other scenarios, such as complex survey designs, or machine learning applications where analysts wish to bootstrap the entire ML pipeline.
To implement a custom estimation procedure, we define an estimator
argument in inferences()
function. estimator()
must
- Take a data frame as input
- Return an object of class
hypotheses
,predictions
,slopes
, orcomparisons
18.3.1 Example: Bootstrap with Propensity Score Weighting
library(marginaleffects)
lalonde <- get_dataset("lalonde")
# Define custom estimator function
# input: data frame
# output: marginaleffects object
estimator <- function(data) {
# Step 1: Estimate propensity scores
fit1 <- glm(treat ~ age + educ + race, family = binomial, data = data)
ps <- predict(fit1, type = "response")
# Step 2: Fit weighted outcome model
m <- lm(re78 ~ treat * (re75 + age + educ + race),
data = data, weight = ps
)
# Step 3: Compute average treatment effect by G-computation
avg_comparisons(m, variables = "treat", wts = ps, vcov = FALSE)
}
# Test the estimator on the full dataset
estimator(lalonde)
Estimate
1146
Term: treat
Type: response
Comparison: 1 - 0
# Bootstrap the entire procedure
cmp <- estimator(lalonde) |>
inferences(method = "rsample", estimator = estimator)
The benefit of this approach is that it abstracts away from the mechanics of bootstrapping. We do not need to interact with the rsample
or boot
packages directly. We do not need to post-process or reshape the draws to obtain the format of interest. Indeed, the output is a standard marginaleffects
object, that can be manipulated with all the standard tools, and comes with a nice pretty-print method:
cmp
Estimate 2.5 % 97.5 %
1146 -379 2623
Term: treat
Type: response
Comparison: 1 - 0
It also allows us to examine the bootstrap distribution easily, by extracting draws in a variety of convenient format using the get_draws()
helper function from marginaleffects
.