Introduction

Although regression models are frequently used in empirical research to study relationships among variables, often the quantity of substantive interest is not one of the coefficients of the model, but rather a quantity derived from the coefficients, such as predicted values or average marginal effects. The usual method for estimating the uncertainty of the derived quantities is an approximation known as the “delta method”. The delta method involves two approximations: 1) that the variance of the derived quantity can be represented as a first-order Taylor series, and 2) that the resulting estimate is normally distributed¹.In many cases, especially with nonlinear models, these approximation can fail badly. clarify implements an alternative to the delta method—simulation-based inference—which involves simulating the sampling distributions of the derived quantities.

The methodology clarify relies on is described in King, Tomz, and Wittenberg (2000) and Rainey (2023). Similar functionality exists in the CLARIFY package in Stata [Tomz, Wittenberg, and King (2003)]² and used to be available in the Zelig R package (Imai, King, and Lau 2008), though there are differences in these implementations. clarify provides additional flexibility by allowing the user to request any derived quantity, in addition to providing shortcuts for common quantities, including predictions at representative values, average marginal effects, and average dose-response functions (described below). clarify relies on and can be seen as a companion to the marginaleffects package, which offers similar functionality but uses the delta method for calculating uncertainty.

Using clarify

There are four steps to using clarify:

1. Fit the model to the data using modeling functions in supported packages

2. Use sim() to take draws from the sampling distribution of the estimated model coefficients

3. Use sim_apply() or its wrappers sim_setx(), sim_ame(), and sim_adrf() to compute derived quantities using each simulated set of coefficients

4. Use summary() and plot() to summarize and visualize the distribution of the derived quantities and perform inference on them

In the sections below, we’ll describe how to implement these steps in detail. First, we’ll load clarify using library().

library(clarify)

For a running example, we’ll use the lalonde dataset in the MatchIt package, which contains data on 614 participants enrolled in a job training program or sampled from a survey. The treatment variable is treat and the outcome is re78, and all other variables are confounders. Although the original intent was to estimate the effect of treat on re78, we’ll use it more generally to demonstrate all of clarify’s capabilities. In addition, we’ll use a transformation of the outcome variable to demonstrate applications to nonlinear models.

data("lalonde", package = "MatchIt")

lalonde$re78_0 <- ifelse(lalonde$re78 == 0, 1, 0)

head(lalonde)
#>      treat age educ   race married nodegree re74 re75       re78 re78_0
#> NSW1     1  37   11  black       1        1    0    0  9930.0460      0
#> NSW2     1  22    9 hispan       0        1    0    0  3595.8940      0
#> NSW3     1  30   12  black       0        0    0    0 24909.4500      0
#> NSW4     1  27   11  black       0        1    0    0  7506.1460      0
#> NSW5     1  33    8  black       0        1    0    0   289.7899      0
#> NSW6     1  22    9  black       0        1    0    0  4056.4940      0

1. Fitting the model

The first step is to fit the model. clarify can operate on a large set of models (those supported by marginaleffects), including generalized linear models, multinomial models, multivariate models, and instrumental variable models, many of which are available in other R packages. Even if clarify does not offer direct support for a given model, there are ways to use its functionality regardless (explained in more detail below).

Because we are computing derived quantities, it is not critical to parameterize the model in such a way that the coefficients are interpretable. Below, we’ll fit a probit regression model for the outcome given the treatment and confounders. Coefficients in probit regression do not have a straightforward interpretation, but that’s okay; our quantities of interest can be expressed as derived quantities–functions of the model parameters, such as predictions, counterfactual predictions, and averages and contrasts of them.

fit <- glm(re78_0 ~ treat + age + educ + race + married +
             nodegree + re74 + re75, data = lalonde,
           family = binomial("probit"))

2. Drawing from the coefficient distribution

After fitting the model, we will use sim() to draw coefficients from their sampling distribution. The sampling distribution is assumed to be multivariate normal or multivariate t with appropriate degrees of freedom, with a mean vector equal to the coefficient vector and a covariance matrix equal to the asymptotic covariance matrix extracted from the model. The arguments to sim() are listed below:

sim(fit = , n = , vcov = , coefs = , dist = )

• fit – the fitted model object, the output of the call to the fitting function (e.g., glm())

• n – the number of simulated values to draw; by default, 1000. More values will yield more replicable and precise results at the cost of speed.

• vcov – either the covariance matrix of the estimated coefficients, a function used to extract it from the model (e.g., sandwich::vcovHC() for the robust covariance matrix), or a string or formula giving a code for extracting the covariance matrix (see marginaleffects::get_vcov() for details). If left unspecified, the default covariance matrix will be extracted from the model.

• coefs – either a vector of coefficients to be sampled or a function to extract them from the fitted model. If left unspecified, the default coefficients will be extracted from the model. Typically this does not need to be specified.

• dist – the name of the distribution from which to draw the sampled coefficients. Can be "normal" for a normal distribution or t(#) for a t-distribution, where # represents the degrees of freedom. If left unspecified, sim() will decide on which distribution makes sense given the characteristics of the model (the decision is made by insight::get_df() with type = "wald"). Typically this does not need to be specified.

If your model is not supported by clarify, you can omit the fit argument and just specify the vcov and coefs argument, which will draw the coefficients from the distribution named in dist ("normal" by default).

sim() uses a random number generator to draw the sampled coefficients from the sampling distribution, so a seed should be set using set.seed() to ensure results are replicable across sessions.

The output of the call to sim() is a clarify_sim object, which contains the sampled coefficients, the original model fit object if supplied, and the coefficients and covariance matrix used to sample.

set.seed(1234)

# Drawing simulated coefficients using an HC2 robust
# covariance matrix
s <- sim(fit, vcov = "HC2")

s
#> A `clarify_sim` object
#>  - 10 coefficients, 1000 simulated values
#>  - sampled distribution: multivariate normal
#>  - original fitting function call:
#> 
#> glm(formula = re78_0 ~ treat + age + educ + race + married + 
#>     nodegree + re74 + re75, family = binomial("probit"), data = lalonde)

3. Computing derived quantities

After sampling the coefficients, we will compute derived quantities on each set of sampled coefficient and store the result, which represents a “posterior” distribution of the derived quantity, as well as on the original coefficients, which are used as the final estimates. The core functionality is provided by sim_apply(), which accepts a clarify_sim object from sim() and a function to compute and return one or more derived quantities, then applies that function to each set of simulated coefficients. The arguments to sim_apply() are below:

sim_apply(sim = , FUN = , verbose = , cl = , ...)

• sim – a clarify_sim object; the output of a call to sim().

• FUN – a function that takes in either a model fit object or a vector of coefficients and returns one or more derived quantities. The first argument should be named fit to take in a model fit object or coefs to take in coefficients.

• verbose – whether to display a progress bar.

• cl – an argument that controls parallel processing, which can be the number of cores to use or a cluster object resulting from parallel::makeCluster().

• ... - further arguments to FUN.

The FUN argument can be specified in one of two ways: either as a function that takes in a model fit object or a function that takes in a vector of coefficients. The latter will always work but the former only works for supported models. When the function takes in a model fit object, sim_apply() will first insert each set of sampled coefficients into the model fit object and then supply the modified model to FUN.

For example, let’s say our derived quantity of interest is the predicted probability of the outcome for participant PSID1. We would specified our FUN function as follows:

sim_fun1 <- function(fit) {
  predict(fit, newdata = lalonde["PSID1",], type = "response")
}

The fit object supplied to this function will be one in which the coefficients have been set to their values in a draw from their sampling distribution as generated by sim(). We then supply the function to sim_apply() to simulate the sampling distribution of the predicted value of interest:

est1 <- sim_apply(s, FUN = sim_fun1, verbose = FALSE)

est1
#> A `clarify_est` object (from `sim_apply()`)
#>  - 1000 simulated values
#>  - 1 quantity estimated:                 
#>  PSID1 0.02377322

The resulting clarify_est object contains the simulated estimates in matrix form as well as the estimate computed on the original coefficients. We’ll examine the sampling distribution shortly, but first we’ll demonstrate computing a derived quantity from the coefficients directly.

The race variable is a factor, and the black category is used as the reference level, so it’s not immediately clear whether there is a difference between the coefficients racehispan and racewhite, which represent the non-reference categories hispan and white. To compare these two directly, we can use sim_apply() to compute a derived quantity that corresponds to the difference between them.

sim_fun2 <- function(coefs) {
  hispan <- unname(coefs["racehispan"])
  white <- unname(coefs["racewhite"])
  
  c("w - h" = white - hispan)
}

est2 <- sim_apply(s, FUN = sim_fun2, verbose = FALSE)

est2
#> A `clarify_est` object (from `sim_apply()`)
#>  - 1000 simulated values
#>  - 1 quantity estimated:                
#>  w - h 0.1089068

The function supplied to FUN can be arbitrarily complicated and return as many derived quantities as you want, though the slower each run of FUN is, the longer it will take to simulate the derived quantities. Using parallel processing by supplying an argument to cl can sometimes dramatically speed up evaluation.

There are several functions in clarify that serve as convenience wrappers for sim_apply() to automate some common derived quantities of interest. These include

• sim_setx() – computing predicted values and first differences at representative or user-specified values of the predictors

• sim_ame() – computing average marginal means, contrasts of average marginal means, and average marginal effects

• sim_adrf() – computing average dose-response functions and average marginal effects functions

These are described in their own sections below. In addition, there are functions that have methods for clarify_est objects, including cbind() for combining two clarify_est objects together and transform() for computing quantities that are derived from the already-computed derived quantities. These are also described in their own sections below.

4. Summarize and visualize the simulated distribution

To examine the uncertainty around and perform inference on our estimated quantities, we can use plot() and summary() on the clarify_est object.

plot() displays a density plot of the resulting estimates across the simulations, with markers identifying the point estimate (computed using the original model coefficients as recommended by Rainey (2023) ) and, optionally, uncertainty bounds (which function like confidence or credible interval bounds). The arguments to plot() are below:

plot(x = , parm = , ci = , level = , method = , reference =)

• x – the clarify_est object (the output of a call to sim_apply()).

• parm – the names or indices of the quantities to be plotted if more than one was estimated in sim_apply(); if unspecified, all will be plotted.

• ci – whether to display lines at the uncertainty bounds. The default is TRUE to display them.

• level – if ci is TRUE, the desired two-sided confidence level. The default is .95 so that that the bounds are at the .025 and .975 quantiles when method (see below) is "quantile".

• method – if ci is TRUE, the method used to compute the bounds. Allowable methods include a Normal approximation ("wald") or using the quantiles of the resulting distribution ("quantile"). The Normal approximation involves multiplying the standard deviation of the estimates (i.e., which functions like the standard error of the sampling distribution) by the critical Z-statistic computed using (1-level)/2 to create a symmetric margin of error around the point estimate. The default is "quantile" to instead use quantile-based bounds, which are more appropriate when the distribution is non-Normal. However, quantile-based bounds may require more simulations to stabilize.

• reference – whether to display a normal density over the plot for each estimate, The default is FALSE to omit them.

Below, we plot the first estimate we computed above, the predicted probability for participant PSID1:

plot(est1, reference = TRUE, ci = FALSE)

From the plot, we can see that the distribution of simulated values is non-Normal, asymmetrical, and not centered around the estimate, with no values falling below 0 because the outcome is a predicted probability. Given its non-Normality, the quantile-based bounds are clearly more appropriate, as the bounds computed from the Normal approximation would be outside the bounds of the estimate. The plot itself is a ggplot object that can be modified using ggplot2 syntax.

We can use summary() to display the value of the point estimate, the uncertainty bounds, and other statistics that describe the distribution of estimates. The arguments to summary() are below:

summary(object = , parm = , level = , method = , null = )

• object – the clarify_est object (the output of a call to sim_apply()).

• parm – the names or indices of the quantities to be displayed if more than one was estimated in sim_apply(); if unspecified, all will be displayed.

• level – the desired two-sided confidence level. The default is .95 so that that the bounds are at the .025 and .975 quantiles when method (see below) is "quantile".

• method – the method used to compute the uncertainty bounds. Allowable methods include a Normal approximation ("wald") or using the quantiles of the resulting distribution ("quantile"). See plot() above.

• null – an optional argument specifying the null value in a hypothesis test for the estimates. If specified, a p-value will be computed using either a standard Z-test (if method is "quantile") or an inversion of the uncertainty interval. The default is not to display any p-values.

We’ll use summary() with the default arguments on our first clarify_est object to view the point estimate and quantile-based uncertainty bounds.

summary(est1)
#>       Estimate   2.5 %  97.5 %
#> PSID1  0.02377 0.00339 0.10539

Our second estimated quantity, the difference between two regression coefficients, is closer to Normally distributed, as the plot below demonstrates (and would be expected theoretically), so we’ll use the Normal approximation to test the hypothesis that difference differs from 0.

plot(est2, reference = TRUE, ci = FALSE)

summary(est2, method = "wald", null = 0)
#>       Estimate  2.5 % 97.5 % Std. Error Z value P-value
#> w - h    0.109 -0.307  0.525      0.212    0.51    0.61

The uncertainty intervals and p-values in the summary() output are computed using the Normal approximation because we set method = "wald", and the p-value for the test that our estimate is equal to 0 is returned because we set null = 0. The computed bounds are very close to the quantile-based bound (displayed in the plot, and can be requested in summary() by setting method = "quantile") because the simulated sampling distribution is close to Normal. Note that the Normal approximation should be used only when the simulated sampling distribution is both close to Normally distributed and centered around the estimate (i.e., when the mean of the simulated values [red vertical line] coincides with the estimate computed on the original coefficients [black vertical line]).

sim_apply() wrappers: sim_setx(), sim_ame(), sim_adrf()

sim_apply() can be used to compute the simulated sampling distribution for any arbitrary derived quantity of interest, but there are some quantities that are common in applied research and may otherwise be somewhat challenging to program on their own, so clarify provides shortcut functions to make computing these quantities simple. These functions include sim_setx(), sim_ame(), and sim_adrf(). Each of these can only be used when regression models compatible with clarify are supplied to the original call to sim().

Like sim_apply(), each of these functions is named sim_*(), which signifies that they are to be used on an object produced by sim() (i.e., a clarify_sim object). (Multiple calls to these functions can be applied to the same clarify_sim object and combined; see the cbind() section below.) These functions are described below.

sim_setx(sim = , x = , x1 = , outcome = , type = , verbose = , cl = )

est3 <- sim_setx(s,
                 x = list(treat = 0:1,
                          re75 = c(0, 20000),
                          race = "black"),
                 verbose = FALSE)

attr(est3, "setx")
#>                         treat      age     educ  race married nodegree     re74
#> treat = 0, re75 = 0         0 27.36319 10.26873 black       0        1 4557.547
#> treat = 1, re75 = 0         1 27.36319 10.26873 black       0        1 4557.547
#> treat = 0, re75 = 20000     0 27.36319 10.26873 black       0        1 4557.547
#> treat = 1, re75 = 20000     1 27.36319 10.26873 black       0        1 4557.547
#>                          re75
#> treat = 0, re75 = 0         0
#> treat = 1, re75 = 0         0
#> treat = 0, re75 = 20000 20000
#> treat = 1, re75 = 20000 20000

plot(est3, var = "re75", ci = FALSE)

est4 <- sim_setx(s,
                 x = list(treat = 0:1,
                          re75 = seq(0, 20000, by = 2000),
                          race = "black"),
                 verbose = FALSE)

plot(est4)

est5 <- sim_setx(s,
                 x = list(treat = 0, re75 = 0),
                 x1 = list(treat = 1, re75 = 0),
                 verbose = FALSE)

summary(est5)
#>           Estimate   2.5 %  97.5 %
#> treat = 0   0.2316  0.1658  0.3117
#> treat = 1   0.1667  0.0930  0.2728
#> FD         -0.0649 -0.1346  0.0209

sim_ame(sim = , var = , subset = , contrast = , outcome = ,
        type = , eps = , verbose = , cl = )

est6 <- sim_ame(s,
                var = "treat", subset = treat == 1,
                contrast = "rr", verbose = FALSE)

summary(est6, null = c(`RR` = 1))
#>         Estimate 2.5 % 97.5 % P-value
#> E[Y(0)]    0.319 0.250  0.400       .
#> E[Y(1)]    0.244 0.189  0.306       .
#> RR         0.765 0.559  1.082    0.13

est7 <- sim_ame(s, var = "age", verbose = FALSE)

summary(est7)
#>           Estimate   2.5 %  97.5 %
#> dY/d(age)  0.00618 0.00257 0.00965

sim_adrf(sim = , var = , subset = , contrast = , at = ,
         n = , outcome = , type = , eps = , verbose = ,
         cl = )

est8 <- sim_adrf(s, var = "age", contrast = "adrf",
                 at = seq(18, 50, by = 2),
                 verbose = FALSE)

plot(est8)

summary(est8, parm = 1:4)
#>          Estimate 2.5 % 97.5 %
#> E[Y(18)]    0.178 0.143  0.227
#> E[Y(20)]    0.189 0.155  0.233
#> E[Y(22)]    0.200 0.169  0.240
#> E[Y(24)]    0.211 0.182  0.249

est9 <- sim_adrf(s, var = "age", contrast = "amef",
                 at = seq(18, 50, by = 2),
                 verbose = FALSE)

plot(est9)

Transforming and combining estimates

Often, our quantities of interest are not just the outputs of the functions above, but comparisons between them. For example, to test for moderation of a treatment effect, we may want to compare AMEs in multiple treatment groups. Or, it might be that we are interested in an effect described using a different effect measure than the one originally produced; for example, we may decide we want the risk difference AME after computing the risk ratio AME. The functions transform() and cbind() allow users to transform quantities in a single clarify_est object and combine two clarify_est objects. These are essential for computing quantities that themselves are derived from the derived quantities computed by the sim_*() functions.

lalonde <- transform(lalonde, re78_0 = ifelse(re78 == 0, 1, 0))

est6 <- transform(est6, RD = `E[Y(1)]` - `E[Y(0)]`)

summary(est6, null = c(`RR` = 1, `RD` = 0))
#>         Estimate   2.5 %  97.5 % P-value
#> E[Y(0)]   0.3195  0.2501  0.4000       .
#> E[Y(1)]   0.2443  0.1892  0.3062       .
#> RR        0.7646  0.5592  1.0822    0.13
#> RD       -0.0752 -0.1692  0.0216    0.13

# AME of treat with race = "black"
est10b <- sim_ame(s, var = "treat", subset = race == "black",
                  contrast = "diff", verbose = FALSE)
summary(est10b)
#>         Estimate   2.5 %  97.5 %
#> E[Y(0)]   0.3370  0.2562  0.4279
#> E[Y(1)]   0.2596  0.2004  0.3258
#> Diff     -0.0774 -0.1743  0.0215

# AME of treat with race = "hispan"
est10h <- sim_ame(s, var = "treat", subset = race == "hispan",
                  contrast = "diff", verbose = FALSE)
summary(est10h)
#>         Estimate   2.5 %  97.5 %
#> E[Y(0)]   0.1745  0.1052  0.2875
#> E[Y(1)]   0.1225  0.0602  0.2376
#> Diff     -0.0520 -0.1161  0.0169

names(est10b) <- paste(names(est10b), "b", sep = "_")
names(est10h) <- paste(names(est10h), "h", sep = "_")

est10 <- cbind(est10b, est10h)
summary(est10)
#>           Estimate   2.5 %  97.5 %
#> E[Y(0)]_b   0.3370  0.2562  0.4279
#> E[Y(1)]_b   0.2596  0.2004  0.3258
#> Diff_b     -0.0774 -0.1743  0.0215
#> E[Y(0)]_h   0.1745  0.1052  0.2875
#> E[Y(1)]_h   0.1225  0.0602  0.2376
#> Diff_h     -0.0520 -0.1161  0.0169

est10 <- transform(est10, `Dh - Db` = Diff_h - Diff_b)
summary(est10, parm = "Dh - Db")
#>         Estimate    2.5 %   97.5 %
#> Dh - Db  0.02540 -0.00406  0.07855

Using clarify with multiply imputed data

Simulation-based inference in multiply imputed data is relatively straightforward. Simulated coefficients are drawn from the model estimated in each imputed dataset separately, and then the simulated coefficients are pooled into a single set of simulated coefficients. In Bayesian terms, this would be considered “mixing draws” and is the recommended approach for Bayesian analysis with multiply imputed data (Zhou and Reiter 2010).

Using clarify with multiply imputed data is simple. Rather than using sim(), we use the function misim(). misim() functions just like sim() except that it takes in a list of model fits (i.e, containing a model fit to each imputed dataset) or an object containing such a list (e.g., a mira object from mice::with() or a mimira object from MatchThem::with()). misim() simulates coefficient distributions within each imputed dataset and then appends them together to a form a single combined set of coefficient draws.

sim_apply() and its wrappers accept the output of misim() and compute the desired quantity using each set of coefficients. When these function rely on using a dataset (e.g., sim_ame(), which averages predicted outcomes across all units in the dataset used to fit the model), they automatically know to associate a given coefficient draw with the imputed dataset that was used to fit the model that produced that draw. In user-written functions supplied to the FUN argument of sim_apply(), it is important to extract the dataset from the model fit. This is demonstrated below.

The final estimates of the quantity of interest is computed as the mean of the estimates computed in each imputed dataset (i.e., using the original coefficients, not the simulated ones), which is the same quantity that would be computed using standard pooling rules. This is not always valid for noncollapsible estimates, like ratios, and so care should be taken to ensure the mean of the resulting estimates has a valid interpretation.

The arguments to misim() are as follows:

misim(fitlist = , n = , vcov = , coefs = , dist = )

• fitslist – a list of model fits or an accepted object containing them (e.g., mira object from mice::with())

• n – the number of simulations to run for each imputed dataset. The default is 1000, but fewer can be used because the total number of simulated quantities will be m*n, where m is the number of imputed datasets.

• vcov, coefs, dist – the same as with sim(), except that a list of such arguments can be supplied to be applied to each imputed dataset.

Below we illustrate using misim() and sim_apply() with multiply imputed data. We’ll use the africa dataset from the Amelia package.

library(Amelia)
data("africa", package = "Amelia")

# Multiple imputation
a.out <- amelia(x = africa, m = 10, cs = "country",
                ts = "year", logs = "gdp_pc", p2s = 0)

# Fit model to each dataset
model.list <- with(a.out, lm(gdp_pc ~ infl * trade))

# Simulate coefficients
si <- misim(model.list, n = 200)

si
#> A `clarify_misim` object
#>  - 4 coefficients, 10 imputations with 200 simulated values each
#>  - sampled distributions: multivariate t(116)

The function we’ll be applying to each imputed dataset will be one that computes the average marginal effect of infl. (We’ll run the same analysis afterward using sim_ame().)

sim_fun <- function(fit) {
  #Extract the original dataset using get_predictors()
  X <- insight::get_predictors(fit)
  
  p0 <- predict(fit, newdata = X)
  
  #Perturb infl slightly
  p1 <- predict(fit, newdata = transform(X, infl = infl + 1e-5))
  
  return(c(AME = mean((p1 - p0) / 1e-5)))
}

est_mi <- sim_apply(si, FUN = sim_fun, verbose = FALSE)

summary(est_mi)
#>     Estimate 2.5 % 97.5 %
#> AME    -5.75 -8.91  -2.28

Note that sim_apply() “knows” which imputation produced each set of simulated coefficients, so using insight::get_predictors() on the fit supplied to sim_fun() will use the right dataset. Care should be taken when analyses restrict each imputed dataset in a different way (e.g. when matching with a caliper in each one), as the resulting imputations may not refer to a specific target population and mixing the draws may be invalid.

Below, we can use sim_ame():

est_mi2 <- sim_ame(si, var = "infl", verbose = FALSE)

We get the same results, as expected.