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:

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:

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.