# 9Bootstrap & Simulation

marginaleffects offers an inferences() function to compute uncertainty estimates using the bootstrap and simulation-based inference.

WARNING: The inferences() function is experimental. It may be renamed, the user interface may change, or the functionality may migrate to arguments in other marginaleffects functions.

Consider a simple model:

mod <- lm(Sepal.Length ~ Petal.Width * Petal.Length + factor(Species), data = iris)

We will compute uncertainty estimates around the output of comparisons(), but note that the same approach works with the predictions() and slopes() functions as well.

## 9.1 Delta method

The default strategy to compute standard errors and confidence intervals is the delta method. This is what we obtain by calling:

avg_comparisons(mod, by = "Species", variables = "Petal.Width")
#>
#>         Term Contrast    Species Estimate Std. Error      z Pr(>|z|)   S  2.5 % 97.5 %
#>  Petal.Width mean(+1) setosa      -0.1103      0.285 -0.387    0.699 0.5 -0.669  0.449
#>  Petal.Width mean(+1) versicolor  -0.0201      0.160 -0.125    0.900 0.2 -0.334  0.293
#>  Petal.Width mean(+1) virginica    0.0216      0.169  0.128    0.898 0.2 -0.309  0.353
#>
#> Columns: term, contrast, Species, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted
#> Type:  response

Since this is the default method, we obtain the same results if we add the inferences() call in the chain:

avg_comparisons(mod, by = "Species", variables = "Petal.Width") |>
inferences(method = "delta")
#>
#>         Term Contrast    Species Estimate Std. Error      z Pr(>|z|)   S  2.5 % 97.5 %
#>  Petal.Width mean(+1) setosa      -0.1103      0.285 -0.387    0.699 0.5 -0.669  0.449
#>  Petal.Width mean(+1) versicolor  -0.0201      0.160 -0.125    0.900 0.2 -0.334  0.293
#>  Petal.Width mean(+1) virginica    0.0216      0.169  0.128    0.898 0.2 -0.309  0.353
#>
#> Columns: term, contrast, Species, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted
#> Type:  response

## 9.2 Bootstrap

marginaleffects supports three bootstrap frameworks in R: the well-established boot package, the newer rsample package, and the so-called “bayesian bootstrap” in fwb.

### 9.2.1boot

avg_comparisons(mod, by = "Species", variables = "Petal.Width") |>
inferences(method = "boot")
#>
#>         Term Contrast    Species Estimate Std. Error  2.5 % 97.5 %
#>  Petal.Width mean(+1) setosa      -0.1103      0.263 -0.624  0.412
#>  Petal.Width mean(+1) versicolor  -0.0201      0.161 -0.323  0.286
#>  Petal.Width mean(+1) virginica    0.0216      0.183 -0.332  0.380
#>
#> Columns: term, contrast, Species, estimate, predicted_lo, predicted_hi, predicted, std.error, conf.low, conf.high
#> Type:  response

All unknown arguments that we feed to inferences() are pushed forward to boot::boot():

est <- avg_comparisons(mod, by = "Species", variables = "Petal.Width") |>
inferences(method = "boot", sim = "balanced", R = 500, conf_type = "bca")
est
#>
#>         Term Contrast    Species Estimate Std. Error  2.5 % 97.5 %
#>  Petal.Width mean(+1) setosa      -0.1103      0.272 -0.633  0.435
#>  Petal.Width mean(+1) versicolor  -0.0201      0.159 -0.343  0.301
#>  Petal.Width mean(+1) virginica    0.0216      0.180 -0.369  0.354
#>
#> Columns: term, contrast, Species, estimate, predicted_lo, predicted_hi, predicted, std.error, conf.low, conf.high
#> Type:  response

We can extract the original boot object from an attribute:

attr(est, "inferences")
#>
#> BALANCED BOOTSTRAP
#>
#>
#> Call:
#> bootstrap_boot(model = model, FUN = FUN, newdata = ..1, vcov = ..2,
#>     variables = ..3, type = ..4, by = ..5, conf_level = ..6,
#>     cross = ..7, comparison = ..8, transform = ..9, wts = ..10,
#>     hypothesis = ..11, eps = ..12)
#>
#>
#> Bootstrap Statistics :
#>        original        bias    std. error
#> t1* -0.11025325  0.0039582054   0.2715939
#> t2* -0.02006005  0.0009499053   0.1591604
#> t3*  0.02158742 -0.0004392025   0.1803044

Or we can extract the individual draws with the posterior_draws() function:

posterior_draws(est) |> head()
#>   drawid        draw        term contrast    Species    estimate predicted_lo predicted_hi predicted std.error   conf.low conf.high
#> 1      1 -0.06570473 Petal.Width mean(+1)     setosa -0.11025325     5.013640     4.901389  4.957514 0.2715939 -0.6331331 0.4354420
#> 2      1  0.05048857 Petal.Width mean(+1) versicolor -0.02006005     6.330887     6.325011  6.327949 0.1591604 -0.3432209 0.3006181
#> 3      1  0.10414181 Petal.Width mean(+1)  virginica  0.02158742     6.997499     7.033528  7.015513 0.1803044 -0.3686909 0.3537212
#> 4      2 -0.07332859 Petal.Width mean(+1)     setosa -0.11025325     5.013640     4.901389  4.957514 0.2715939 -0.6331331 0.4354420
#> 5      2  0.02696760 Petal.Width mean(+1) versicolor -0.02006005     6.330887     6.325011  6.327949 0.1591604 -0.3432209 0.3006181
#> 6      2  0.07328022 Petal.Width mean(+1)  virginica  0.02158742     6.997499     7.033528  7.015513 0.1803044 -0.3686909 0.3537212

posterior_draws(est, shape = "DxP") |> dim()
#> [1] 500   3

### 9.2.2rsample

As before, we can pass arguments to rsample::bootstraps() through inferences(). For example, for stratified resampling:

est <- avg_comparisons(mod, by = "Species", variables = "Petal.Width") |>
inferences(method = "rsample", R = 100, strata = "Species")
est
#>
#>         Term Contrast    Species Estimate  2.5 % 97.5 %
#>  Petal.Width mean(+1) setosa      -0.1103 -0.593  0.359
#>  Petal.Width mean(+1) versicolor  -0.0201 -0.319  0.209
#>  Petal.Width mean(+1) virginica    0.0216 -0.297  0.304
#>
#> Columns: term, contrast, Species, estimate, predicted_lo, predicted_hi, predicted, conf.low, conf.high
#> Type:  response

attr(est, "inferences")
#> # Bootstrap sampling using stratification with apparent sample
#> # A tibble: 101 × 3
#>    splits           id           estimates
#>    <list>           <chr>        <list>
#>  1 <split [150/56]> Bootstrap001 <tibble [3 × 7]>
#>  2 <split [150/52]> Bootstrap002 <tibble [3 × 7]>
#>  3 <split [150/54]> Bootstrap003 <tibble [3 × 7]>
#>  4 <split [150/59]> Bootstrap004 <tibble [3 × 7]>
#>  5 <split [150/57]> Bootstrap005 <tibble [3 × 7]>
#>  6 <split [150/60]> Bootstrap006 <tibble [3 × 7]>
#>  7 <split [150/58]> Bootstrap007 <tibble [3 × 7]>
#>  8 <split [150/53]> Bootstrap008 <tibble [3 × 7]>
#>  9 <split [150/59]> Bootstrap009 <tibble [3 × 7]>
#> 10 <split [150/50]> Bootstrap010 <tibble [3 × 7]>
#> # ℹ 91 more rows

Or we can extract the individual draws with the posterior_draws() function:

posterior_draws(est) |> head()
#>   drawid        draw        term contrast    Species    estimate predicted_lo predicted_hi predicted   conf.low conf.high
#> 1      1 -0.14511321 Petal.Width mean(+1)     setosa -0.11025325     5.013640     4.901389  4.957514 -0.5926132 0.3586174
#> 2      1 -0.12316971 Petal.Width mean(+1) versicolor -0.02006005     6.330887     6.325011  6.327949 -0.3191469 0.2093705
#> 3      1 -0.11303711 Petal.Width mean(+1)  virginica  0.02158742     6.997499     7.033528  7.015513 -0.2965658 0.3035202
#> 4      2  0.10290137 Petal.Width mean(+1)     setosa -0.11025325     5.013640     4.901389  4.957514 -0.5926132 0.3586174
#> 5      2  0.07008285 Petal.Width mean(+1) versicolor -0.02006005     6.330887     6.325011  6.327949 -0.3191469 0.2093705
#> 6      2  0.05492862 Petal.Width mean(+1)  virginica  0.02158742     6.997499     7.033528  7.015513 -0.2965658 0.3035202

posterior_draws(est, shape = "PxD") |> dim()
#> [1]   3 100

### 9.2.3 Fractional Weighted Bootstrap (aka Bayesian Bootstrap)

The fwb package implements fractional weighted bootstrap (aka Bayesian bootstrap):

“fwb implements the fractional weighted bootstrap (FWB), also known as the Bayesian bootstrap, following the treatment by Xu et al. (2020). The FWB involves generating sets of weights from a uniform Dirichlet distribution to be used in estimating statistics of interest, which yields a posterior distribution that can be interpreted in the same way the traditional (resampling-based) bootstrap distribution can be.” -Noah Greifer

The inferences() function makes it easy to apply this inference strategy to marginaleffects objects:

avg_comparisons(mod) |> inferences(method = "fwb")
#>
#>          Term            Contrast Estimate Std. Error  2.5 % 97.5 %
#>  Petal.Length +1                    0.8929     0.0812  0.733  1.055
#>  Petal.Width  +1                   -0.0362     0.1591 -0.338  0.273
#>  Species      versicolor - setosa  -1.4629     0.3298 -2.118 -0.830
#>  Species      virginica - setosa   -1.9842     0.3875 -2.728 -1.209
#>
#> Columns: term, contrast, estimate, std.error, conf.low, conf.high
#> Type:  response

## 9.3 Simulation-based inference

This simulation-based strategy to compute confidence intervals was described in Krinsky & Robb (1986) and popularized by King, Tomz, Wittenberg (2000). We proceed in 3 steps:

1. Draw R sets of simulated coefficients from a multivariate normal distribution with mean equal to the original model’s estimated coefficients and variance equal to the model’s variance-covariance matrix (classical, “HC3”, or other).
2. Use the R sets of coefficients to compute R sets of estimands: predictions, comparisons, or slopes.
3. Take quantiles of the resulting distribution of estimands to obtain a confidence interval and the standard deviation of simulated estimates to estimate the standard error.

Here are a few examples:

avg_comparisons(mod, by = "Species", variables = "Petal.Width") |>
inferences(method = "simulation")
#>
#>         Term Contrast    Species Estimate Std. Error  2.5 % 97.5 %
#>  Petal.Width mean(+1) setosa      -0.1103      0.291 -0.693  0.457
#>  Petal.Width mean(+1) versicolor  -0.0201      0.154 -0.303  0.302
#>  Petal.Width mean(+1) virginica    0.0216      0.163 -0.291  0.356
#>
#> Columns: term, contrast, Species, estimate, std.error, conf.low, conf.high, predicted_lo, predicted_hi, predicted, tmp_idx
#> Type:  response

Since simulation based inference generates R estimates of the quantities of interest, we can treat them similarly to draws from the posterior distribution in bayesian models. For example, we can extract draws using the posterior_draws() function, and plot their distributions using packages likeggplot2 and ggdist:

avg_comparisons(mod, by = "Species", variables = "Petal.Width") |>
inferences(method = "simulation") |>
posterior_draws("rvar") |>
ggplot(aes(y = Species, xdist = rvar)) +

## 9.4 Multiple imputation and missing data

The same workflow and the same inferences() function can be used to estimate models with multiple imputation for missing data.