library(marginaleffects)
mod <- lm(Sepal.Length ~ Petal.Width * Petal.Length + factor(Species), data = iris)9 Bootstrap & 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:
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: responseSince 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: response9.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.1 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.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: responseAll 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: responseWe 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.1803044Or 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.2 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.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 rowsOr 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 1009.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: response9.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:
- Draw
Rsets 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
Rsets of coefficients to computeRsets 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.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: responseSince 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()
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.