Bootstrap and Simulation-Based Inference (Experimental) • marginaleffects

marginaleffects offers an experimental 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:

library(marginaleffects)

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.

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|)  2.5 % 97.5 %
#>  Petal.Width mean(+1) setosa      -0.1103      0.285 -0.387    0.699 -0.669  0.449
#>  Petal.Width mean(+1) versicolor  -0.0201      0.160 -0.125    0.900 -0.334  0.293
#>  Petal.Width mean(+1) virginica    0.0216      0.169  0.128    0.898 -0.309  0.353
#> 
#> Columns: term, contrast, Species, estimate, std.error, statistic, p.value, conf.low, conf.high, predicted, predicted_hi, predicted_lo

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|)  2.5 % 97.5 %
#>  Petal.Width mean(+1) setosa      -0.1103      0.285 -0.387    0.699 -0.669  0.449
#>  Petal.Width mean(+1) versicolor  -0.0201      0.160 -0.125    0.900 -0.334  0.293
#>  Petal.Width mean(+1) virginica    0.0216      0.169  0.128    0.898 -0.309  0.353
#> 
#> Columns: term, contrast, Species, estimate, std.error, statistic, p.value, conf.low, conf.high, predicted, predicted_hi, predicted_lo

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.

`boot`

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.636  0.418
#>  Petal.Width mean(+1) versicolor  -0.0201      0.162 -0.358  0.307
#>  Petal.Width mean(+1) virginica    0.0216      0.186 -0.350  0.372
#> 
#> Columns: term, contrast, Species, estimate, predicted, predicted_hi, predicted_lo, std.error, conf.low, conf.high

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.273 -0.722  0.416
#>  Petal.Width mean(+1) versicolor  -0.0201      0.166 -0.338  0.305
#>  Petal.Width mean(+1) virginica    0.0216      0.188 -0.334  0.426
#> 
#> Columns: term, contrast, Species, estimate, predicted, predicted_hi, predicted_lo, std.error, conf.low, conf.high

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, 
#>     comparison = ..7, transform = ..8, wts = ..9, hypothesis = ..10, 
#>     eps = ..11)
#> 
#> 
#> Bootstrap Statistics :
#>        original       bias    std. error
#> t1* -0.11025325 0.0016255817   0.2731068
#> t2* -0.02006005 0.0010577310   0.1658978
#> t3*  0.02158742 0.0007955211   0.1883568

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

posterior_draws(est) |> head()
#>   drawid        draw        term contrast    Species    estimate predicted predicted_hi predicted_lo std.error   conf.low conf.high
#> 1      1 -0.17362042 Petal.Width mean(+1)     setosa -0.11025325  4.957514     4.901389     5.013640 0.2731068 -0.7221460 0.4159604
#> 2      1  0.21983670 Petal.Width mean(+1) versicolor -0.02006005  6.327949     6.325011     6.330887 0.1658978 -0.3379984 0.3052733
#> 3      1  0.40151883 Petal.Width mean(+1)  virginica  0.02158742  7.015513     7.033528     6.997499 0.1883568 -0.3340813 0.4256430
#> 4      2  0.06297702 Petal.Width mean(+1)     setosa -0.11025325  4.957514     4.901389     5.013640 0.2731068 -0.7221460 0.4159604
#> 5      2  0.11957405 Petal.Width mean(+1) versicolor -0.02006005  6.327949     6.325011     6.330887 0.1658978 -0.3379984 0.3052733
#> 6      2  0.14570820 Petal.Width mean(+1)  virginica  0.02158742  7.015513     7.033528     6.997499 0.1883568 -0.3340813 0.4256430

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

`rsample`

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.703  0.382
#>  Petal.Width mean(+1) versicolor  -0.0201 -0.253  0.302
#>  Petal.Width mean(+1) virginica    0.0216 -0.288  0.385
#> 
#> Columns: term, contrast, Species, estimate, predicted, predicted_hi, predicted_lo, conf.low, conf.high

attr(est, "inferences")
#> # Bootstrap sampling using stratification with apparent sample 
#> # A tibble: 101 × 3
#>    splits           id           estimates       
#>    <list>           <chr>        <list>          
#>  1 <split [150/64]> Bootstrap001 <tibble [3 × 7]>
#>  2 <split [150/53]> Bootstrap002 <tibble [3 × 7]>
#>  3 <split [150/57]> Bootstrap003 <tibble [3 × 7]>
#>  4 <split [150/55]> Bootstrap004 <tibble [3 × 7]>
#>  5 <split [150/54]> Bootstrap005 <tibble [3 × 7]>
#>  6 <split [150/54]> Bootstrap006 <tibble [3 × 7]>
#>  7 <split [150/51]> Bootstrap007 <tibble [3 × 7]>
#>  8 <split [150/57]> Bootstrap008 <tibble [3 × 7]>
#>  9 <split [150/60]> 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 predicted_hi predicted_lo   conf.low conf.high
#> 1      1 -0.3102611 Petal.Width mean(+1)     setosa -0.11025325  4.957514     4.901389     5.013640 -0.7032414 0.3821446
#> 2      1  0.1065051 Petal.Width mean(+1) versicolor -0.02006005  6.327949     6.325011     6.330887 -0.2530083 0.3024675
#> 3      1  0.2989504 Petal.Width mean(+1)  virginica  0.02158742  7.015513     7.033528     6.997499 -0.2882844 0.3848567
#> 4      2 -0.1564778 Petal.Width mean(+1)     setosa -0.11025325  4.957514     4.901389     5.013640 -0.7032414 0.3821446
#> 5      2 -0.0142775 Petal.Width mean(+1) versicolor -0.02006005  6.327949     6.325011     6.330887 -0.2530083 0.3024675
#> 6      2  0.0513847 Petal.Width mean(+1)  virginica  0.02158742  7.015513     7.033528     6.997499 -0.2882844 0.3848567

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

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.Width  +1                   -0.0362      0.159 -0.339  0.274
#>  Petal.Length +1                    0.8929      0.079  0.742  1.048
#>  Species      versicolor - setosa  -1.4629      0.316 -2.100 -0.863
#>  Species      virginica - setosa   -1.9842      0.372 -2.724 -1.281
#> 
#> Columns: term, contrast, estimate, std.error, conf.low, conf.high

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:

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).
Use the R sets of coefficients to compute R sets of estimands: predictions, comparisons, or slopes.
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:

library(ggplot2)
library(ggdist)

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.1155      0.282 -0.685  0.447
#>  Petal.Width mean(+1) versicolor  -0.0174      0.163 -0.345  0.296
#>  Petal.Width mean(+1) virginica    0.0310      0.172 -0.311  0.351
#> 
#> Columns: term, contrast, Species, estimate, std.error, conf.low, conf.high, predicted, predicted_hi, predicted_lo, tmp_idx

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)) +
  stat_slabinterval()

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.