Inverse probability weighting with causatr

Inverse probability weighting (IPW) estimates causal effects by reweighting the observed data so that the distribution of confounders is balanced across treatment groups. causatr ships a self-contained density-ratio engine that fits the conditional treatment density via the user’s propensity_model_fn (default stats::glm) and builds per-intervention Hájek weights for an intercept-only marginal structural model (MSM) refit inside contrast().

This vignette demonstrates IPW in causatr using the NHEFS dataset from Hernán & Robins (2025), covering every supported combination of treatment type, outcome type, contrast scale, and inference method for time-fixed treatments.

For non-static interventions (shift, scale, dynamic, IPSI), see vignette("interventions").

Note: IPW and matching support effect modification via A:modifier terms in confounders. IPW expands the per-intervention MSM to include modifier main effects (Y ~ 1 + modifier); matching expands to a saturated MSM (Y ~ A + modifier + A:modifier). Use by = "modifier" in contrast() to obtain stratum-specific estimates.

Setup

The examples in this vignette use the NHEFS (National Health and Nutrition Examination Survey — Epidemiologic Follow-up Study) dataset, the running example in Hernán & Robins (2025). The assumed causal structure is:

  • Treatment (\(A\)): qsmk — quit smoking between 1971 and 1982 (binary: 0/1).
  • Outcome (\(Y\)): wt82_71 — weight change in kg over the same period (continuous). A derived binary version, gained_weight (\(\mathbf{1}\{Y > 0\}\)), is used in the binary outcome sections.
  • Confounders (\(L\)): sex, age, race, education, smoking intensity (cigarettes/day), years smoked, exercise, physical activity, and baseline weight. These jointly affect both the probability of quitting and the subsequent weight change.

Structural causal model (assumed):

\[ \begin{aligned} L &\sim F_L \\ A &\sim \text{Bernoulli}\!\bigl(\text{logit}^{-1}(\alpha_0 + \alpha_L^\top L)\bigr) \\ Y &= \beta_0 + \beta_A A + \beta_L^\top L + \beta_{AL}^\top (A \cdot L) + \varepsilon, \quad \varepsilon \sim N(0, \sigma^2) \end{aligned} \]

The confounders \(L\) are common causes of \(A\) and \(Y\); adjusting for \(L\) blocks all backdoor paths, identifying \(\mathbb{E}[Y^a]\) via the Horvitz–Thompson reweighting estimator. IPW models the treatment mechanism \(P(A \mid L)\) directly and reweights outcomes to remove confounding, rather than modelling \(E[Y \mid A, L]\) as g-computation does.

Code 
library(causatr)
library(tinytable)
library(tinyplot)

data("nhefs")

nhefs_complete <- nhefs[!is.na(nhefs$wt82_71) & !is.na(nhefs$education), ]

nhefs_complete$gained_weight <- as.integer(nhefs_complete$wt82_71 > 0)

nhefs_complete$sex <- factor(
  nhefs_complete$sex,
  levels = 0:1,
  labels = c("Male", "Female")
)

Binary treatment, continuous outcome

This replicates the IPW analysis from Chapter 12 of Hernán & Robins: the average causal effect of quitting smoking (qsmk) on weight change (wt82_71). By default, confounders specifies the same covariate set for the propensity model. Use confounders_treatment to specify a different set for the treatment model when needed (e.g. for AIPW double robustness).

ATE with sandwich SE

Code 
fit_ipw <- causat(
  nhefs_complete,
  outcome = "wt82_71",
  treatment = "qsmk",
  confounders = ~ sex + age + I(age^2) + race + factor(education) +
    smokeintensity + I(smokeintensity^2) + smokeyrs + I(smokeyrs^2) +
    factor(exercise) + factor(active) + wt71 + I(wt71^2),
  estimator = "ipw",
  estimand = "ATE",
  model_fn = stats::glm,
  propensity_model_fn = stats::glm
)

res_ate_sw <- contrast(
  fit_ipw,
  interventions = list(quit = static(1), continue = static(0)),
  reference = "continue",
  type = "difference",
  ci_method = "sandwich"
)
res_ate_sw

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: difference  ·  CI method: sandwich  ·  N: 1566

Intervention means
intervention estimate se ci_lower ci_upper
quit 5.288 0.447 4.412 6.163
continue 1.788 0.217 1.363 2.214
Contrasts
comparison estimate se ci_lower ci_upper
quit vs continue 3.499 0.49 2.538 4.46

The book reports an IPW estimate of ATE \(\approx\) 3.4 kg (95% CI: 2.4, 4.5).

ATE with bootstrap SE

Code 
res_ate_bs <- contrast(
  fit_ipw,
  interventions = list(quit = static(1), continue = static(0)),
  reference = "continue",
  type = "difference",
  ci_method = "bootstrap",
  n_boot = 50L
)
res_ate_bs

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: difference  ·  CI method: bootstrap  ·  N: 1566

Intervention means
intervention estimate se ci_lower ci_upper
quit 5.288 0.525 4.258 6.317
continue 1.788 0.227 1.343 2.234
Contrasts
comparison estimate se ci_lower ci_upper
quit vs continue 3.499 0.564 2.393 4.606

ATT estimand

For IPW, the estimand is fixed at fitting time because it determines the weights. To estimate the ATT, refit with estimand = "ATT".

Code 
fit_ipw_att <- causat(
  nhefs_complete,
  outcome = "wt82_71",
  treatment = "qsmk",
  confounders = ~ sex + age + I(age^2) + race + factor(education) +
    smokeintensity + I(smokeintensity^2) + smokeyrs + I(smokeyrs^2) +
    factor(exercise) + factor(active) + wt71 + I(wt71^2),
  estimator = "ipw",
  estimand = "ATT",
  model_fn = stats::glm,
  propensity_model_fn = stats::glm
)

res_att <- contrast(
  fit_ipw_att,
  interventions = list(quit = static(1), continue = static(0)),
  reference = "continue",
  type = "difference",
  ci_method = "sandwich"
)
res_att

Estimator: ipw  ·  Estimand: ATT  ·  Contrast: difference  ·  CI method: sandwich  ·  N: 403

Intervention means
intervention estimate se ci_lower ci_upper
quit 4.525 0.435 3.672 5.378
continue 1.222 0.296 0.642 1.803
Contrasts
comparison estimate se ci_lower ci_upper
quit vs continue 3.303 0.498 2.326 4.28

ATC estimand

Code 
fit_ipw_atc <- causat(
  nhefs_complete,
  outcome = "wt82_71",
  treatment = "qsmk",
  confounders = ~ sex + age + I(age^2) + race + factor(education) +
    smokeintensity + I(smokeintensity^2) + smokeyrs + I(smokeyrs^2) +
    factor(exercise) + factor(active) + wt71 + I(wt71^2),
  estimator = "ipw",
  estimand = "ATC",
  model_fn = stats::glm,
  propensity_model_fn = stats::glm
)

res_atc <- contrast(
  fit_ipw_atc,
  interventions = list(quit = static(1), continue = static(0)),
  reference = "continue",
  type = "difference",
  ci_method = "sandwich"
)
res_atc

Estimator: ipw  ·  Estimand: ATC  ·  Contrast: difference  ·  CI method: sandwich  ·  N: 1163

Intervention means
intervention estimate se ci_lower ci_upper
quit 5.554 0.489 4.596 6.511
continue 1.984 0.218 1.557 2.412
Contrasts
comparison estimate se ci_lower ci_upper
quit vs continue 3.569 0.528 2.534 4.604

ATT with bootstrap SE

Code 
res_att_bs <- contrast(
  fit_ipw_att,
  interventions = list(quit = static(1), continue = static(0)),
  reference = "continue",
  type = "difference",
  ci_method = "bootstrap",
  n_boot = 50L
)
res_att_bs

Estimator: ipw  ·  Estimand: ATT  ·  Contrast: difference  ·  CI method: bootstrap  ·  N: 403

Intervention means
intervention estimate se ci_lower ci_upper
quit 4.525 0.425 3.693 5.357
continue 1.222 0.312 0.61 1.835
Contrasts
comparison estimate se ci_lower ci_upper
quit vs continue 3.303 0.537 2.249 4.356

Extracting results programmatically

Code 
coef(res_ate_sw)
#>     quit continue 
#> 5.287648 1.788473
confint(res_ate_sw)
#>             lower    upper
#> quit     4.412464 6.162831
#> continue 1.362560 2.214387
vcov(res_ate_sw)
#>                 quit    continue
#> quit     0.199389310 0.003155553
#> continue 0.003155553 0.047222294
tidy(res_ate_sw)
#>               term estimate std.error     type conf.low conf.high
#> 1 quit vs continue 3.499174 0.4902045 contrast 2.538391  4.459958
glance(res_ate_sw)
#>   estimator estimand contrast_type ci_method    n n_interventions
#> 1       ipw      ATE    difference  sandwich 1566               2

Binary treatment, binary outcome

Using the gained_weight indicator as a binary outcome. The IPW engine refits an intercept-only weighted MSM per intervention; contrast() computes marginal risks and derives contrasts via the unified influence-function variance engine (variance_if()), which for IPW uses compute_ipw_if_self_contained_one() to build the per-individual IF on the stacked propensity + MSM M-estimation system.

Risk difference (sandwich)

Code 
fit_ipw_bin <- causat(
  nhefs_complete,
  outcome = "gained_weight",
  treatment = "qsmk",
  confounders = ~ sex + age + I(age^2) + race + factor(education) +
    smokeintensity + I(smokeintensity^2) + smokeyrs + I(smokeyrs^2) +
    factor(exercise) + factor(active) + wt71 + I(wt71^2),
  estimator = "ipw",
  model_fn = stats::glm,
  propensity_model_fn = stats::glm
)

res_rd <- contrast(
  fit_ipw_bin,
  interventions = list(quit = static(1), continue = static(0)),
  reference = "continue",
  type = "difference",
  ci_method = "sandwich"
)
res_rd

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: difference  ·  CI method: sandwich  ·  N: 1566

Intervention means
intervention estimate se ci_lower ci_upper
quit 0.776 0.021 0.735 0.818
continue 0.639 0.014 0.611 0.666
Contrasts
comparison estimate se ci_lower ci_upper
quit vs continue 0.138 0.025 0.089 0.187

Risk difference (bootstrap)

Code 
res_rd_bs <- contrast(
  fit_ipw_bin,
  interventions = list(quit = static(1), continue = static(0)),
  reference = "continue",
  type = "difference",
  ci_method = "bootstrap",
  n_boot = 50L
)
res_rd_bs

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: difference  ·  CI method: bootstrap  ·  N: 1566

Intervention means
intervention estimate se ci_lower ci_upper
quit 0.776 0.022 0.733 0.82
continue 0.639 0.012 0.615 0.662
Contrasts
comparison estimate se ci_lower ci_upper
quit vs continue 0.138 0.026 0.088 0.188

Risk ratio

Code 
res_rr <- contrast(
  fit_ipw_bin,
  interventions = list(quit = static(1), continue = static(0)),
  reference = "continue",
  type = "ratio",
  ci_method = "sandwich"
)
res_rr

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: ratio  ·  CI method: sandwich  ·  N: 1566

Intervention means
intervention estimate se ci_lower ci_upper
quit 0.776 0.021 0.735 0.818
continue 0.639 0.014 0.611 0.666
Contrasts
comparison estimate se ci_lower ci_upper
quit vs continue 1.216 0.042 1.137 1.301

Odds ratio

Code 
res_or <- contrast(
  fit_ipw_bin,
  interventions = list(quit = static(1), continue = static(0)),
  reference = "continue",
  type = "or",
  ci_method = "sandwich"
)
res_or

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: or  ·  CI method: sandwich  ·  N: 1566

Intervention means
intervention estimate se ci_lower ci_upper
quit 0.776 0.021 0.735 0.818
continue 0.639 0.014 0.611 0.666
Contrasts
comparison estimate se ci_lower ci_upper
quit vs continue 1.967 0.264 1.511 2.559

Risk ratio (bootstrap)

Code 
res_rr_bs <- contrast(
  fit_ipw_bin,
  interventions = list(quit = static(1), continue = static(0)),
  reference = "continue",
  type = "ratio",
  ci_method = "bootstrap",
  n_boot = 50L
)
res_rr_bs

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: ratio  ·  CI method: bootstrap  ·  N: 1566

Intervention means
intervention estimate se ci_lower ci_upper
quit 0.776 0.023 0.731 0.822
continue 0.639 0.015 0.609 0.668
Contrasts
comparison estimate se ci_lower ci_upper
quit vs continue 1.216 0.044 1.133 1.305

Parallel bootstrap

Bootstrap replicates can be parallelised by passing parallel and ncpus straight through to boot::boot():

Code 
res_par <- contrast(
  fit_ipw,
  interventions = list(quit = static(1), continue = static(0)),
  reference = "continue",
  ci_method = "bootstrap",
  n_boot = 200L,
  parallel = "multicore",  # or "snow" on Windows
  ncpus = 4L
)

These default to getOption("boot.parallel", "no") and getOption("boot.ncpus", 1L), so a session-wide options(boot.parallel = "multicore", boot.ncpus = 4L) flips parallelism on for every contrast() call without per-call plumbing.

Dynamic intervention (binary treatment)

A dynamic() intervention assigns each individual to a treatment value based on their covariates. Under IPW, dynamic rules on binary treatment use Horwitz–Thompson indicator weights: \(w_i = \mathbb{1}\{A_i = d(L_i)\} / f(d(L_i) \mid L_i)\).

Code 
rule <- function(data, trt) {
  ifelse(data$smokeintensity > 20, 1L, 0L)
}

res_dyn <- contrast(
  fit_ipw,
  interventions = list(
    targeted = dynamic(rule),
    all_quit = static(1)
  ),
  reference = "all_quit",
  type = "difference"
)
res_dyn

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: difference  ·  CI method: sandwich  ·  N: 1566

Intervention means
intervention estimate se ci_lower ci_upper
targeted 2.803 0.355 2.107 3.499
all_quit 5.288 0.447 4.412 6.163
Contrasts
comparison estimate se ci_lower ci_upper
targeted vs all_quit -2.485 0.385 -3.239 -1.73

Heavy smokers (> 20 cigs/day) are assigned to quit, others are left untreated. The contrast with static(1) (all quit) tells us how much is lost by targeting only heavy smokers.

Code 
tt(tidy(res_dyn), digits = 3)
term estimate std.error type conf.low conf.high
targeted vs all_quit -2.48 0.385 contrast -3.24 -1.73

IPSI — incremental propensity score intervention

ipsi(delta) multiplies each individual’s odds of treatment by delta (Kennedy 2019). It is a smooth stochastic intervention that avoids positivity violations — no one is forced to a treatment value they would never take.

Code 
res_ipsi <- contrast(
  fit_ipw,
  interventions = list(
    double_odds = ipsi(2),
    observed = NULL
  ),
  reference = "observed",
  type = "difference"
)
res_ipsi

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: difference  ·  CI method: sandwich  ·  N: 1566

Intervention means
intervention estimate se ci_lower ci_upper
double_odds 3.125 0.218 2.697 3.552
observed 2.638 0.199 2.248 3.028
Contrasts
comparison estimate se ci_lower ci_upper
double_odds vs observed 0.486 0.066 0.357 0.615

Varying delta traces out a dose–response curve:

Code 
deltas <- c(0.25, 0.5, 1, 2, 4)
ipsi_results <- lapply(deltas, function(d) {
  r <- contrast(
    fit_ipw,
    interventions = list(ipsi = ipsi(d), obs = NULL),
    reference = "obs",
    type = "difference"
  )
  data.frame(
    delta = d,
    estimate = r$contrasts$estimate[1],
    ci_lower = r$contrasts$ci_lower[1],
    ci_upper = r$contrasts$ci_upper[1]
  )
})
ipsi_df <- do.call(rbind, ipsi_results)

tinyplot(
  estimate ~ delta,
  data = ipsi_df,
  type = "pointrange",
  ymin = ipsi_df$ci_lower,
  ymax = ipsi_df$ci_upper,
  xlab = expression(delta ~ "(odds multiplier)"),
  ylab = "E[Y(delta)] - E[Y(obs)]",
  main = "IPSI dose-response"
)
abline(h = 0, lty = 2, col = "grey40")

IPSI dose-response curve showing the effect of varying the odds multiplier delta.

At delta = 1 the intervention equals the natural course, so the contrast is zero. Values below 1 reduce the odds of treatment; values above 1 increase them.

Categorical treatment (k > 2 levels)

The density-ratio engine auto-selects a multinomial propensity model (via nnet::multinom) when the treatment is a factor or character column. The per-intervention MSM is intercept-only, so contrast() computes any pairwise difference against a user-chosen reference level.

DGM. A single confounder \(L\) influences both the three-level categorical treatment \(A \in \{A, B, C\}\) and the outcome \(Y\). Treatment assignment follows a softmax (multinomial logistic) model with category A as the reference level. The outcome is linear in the treatment indicators and \(L\).

\[ \begin{aligned} L &\sim N(0, 1) \\ P(A = k \mid L) &= \frac{\exp(\gamma_k L)}{\sum_{j} \exp(\gamma_j L)}, \quad \gamma_A = 0,\; \gamma_B = 0.5,\; \gamma_C = -0.3 \\ Y &= 2 + 1.5 \cdot \mathbf{1}\{A = B\} - 0.8 \cdot \mathbf{1}\{A = C\} + 0.5\, L + \varepsilon, \quad \varepsilon \sim N(0, 1) \end{aligned} \]

True contrasts (vs. reference A): \(\mathbb{E}[Y^B] - \mathbb{E}[Y^A] = 1.5\) and \(\mathbb{E}[Y^C] - \mathbb{E}[Y^A] = -0.8\).

Code 
set.seed(42)
n <- 2000
L <- rnorm(n)
probs <- exp(cbind(0, 0.5 * L, -0.3 * L))
probs <- probs / rowSums(probs)
A <- factor(
  apply(probs, 1, function(p) sample(c("A", "B", "C"), 1, prob = p)),
  levels = c("A", "B", "C")
)
Y <- 2 + (A == "B") * 1.5 + (A == "C") * -0.8 + 0.5 * L + rnorm(n)
cat_df <- data.frame(Y = Y, A = A, L = L)

fit_ipw_cat <- causat(
  cat_df,
  outcome = "Y",
  treatment = "A",
  confounders = ~ L,
  estimator = "ipw",
  model_fn = stats::glm,
  propensity_model_fn = nnet::multinom
)

res_ipw_cat <- contrast(
  fit_ipw_cat,
  interventions = list(
    arm_a = static("A"),
    arm_b = static("B"),
    arm_c = static("C")
  ),
  reference = "arm_a",
  type = "difference"
)
res_ipw_cat

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: difference  ·  CI method: sandwich  ·  N: 2000

Intervention means
intervention estimate se ci_lower ci_upper
arm_a 1.99 0.041 1.91 2.07
arm_b 3.451 0.046 3.361 3.54
arm_c 1.141 0.042 1.058 1.223
Contrasts
comparison estimate se ci_lower ci_upper
arm_b vs arm_a 1.461 0.06 1.344 1.577
arm_c vs arm_a -0.85 0.057 -0.961 -0.738

Continuous treatment (shift / scale MTPs)

For continuous treatments the engine fits a Gaussian treatment density model and uses density-ratio weights for modified treatment policies. static() interventions on continuous treatments are rejected (nobody is observed exactly at the target value, so the HT indicator is zero almost surely). Use shift() or scale_by() instead:

DGM. A single confounder \(L\) drives both a continuous (Gaussian) treatment \(A\) and a continuous outcome \(Y\). The treatment density is \(A \mid L \sim N(1 + 0.5\, L,\; 1)\) and the outcome is linear in \(A\) and \(L\) with additive noise.

\[ \begin{aligned} L &\sim N(0, 1) \\ A \mid L &\sim N(1 + 0.5\, L,\; 1) \\ Y &= 1 + 2\, A + L + \varepsilon, \quad \varepsilon \sim N(0, 1) \end{aligned} \]

Because \(\partial Y / \partial A = 2\), a unit additive shift \(A \mapsto A + 1\) has true causal effect \(\mathbb{E}[Y(A + 1)] - \mathbb{E}[Y(A)] = 2\).

Code 
set.seed(5)
n <- 2000
L <- rnorm(n)
A <- 1 + 0.5 * L + rnorm(n)
Y <- 1 + 2 * A + L + rnorm(n)
cont_df <- data.frame(Y = Y, A = A, L = L)

fit_ipw_cont <- causat(
  cont_df,
  outcome = "Y",
  treatment = "A",
  confounders = ~ L,
  estimator = "ipw",
  model_fn = stats::glm,
  propensity_model_fn = stats::glm
)

res_shift <- contrast(
  fit_ipw_cont,
  interventions = list(up1 = shift(1), observed = NULL),
  reference = "observed",
  type = "difference"
)
res_shift

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: difference  ·  CI method: sandwich  ·  N: 2000

Intervention means
intervention estimate se ci_lower ci_upper
up1 5.05 0.132 4.79 5.309
observed 3.083 0.068 2.949 3.217
Contrasts
comparison estimate se ci_lower ci_upper
up1 vs observed 1.967 0.112 1.748 2.186
Code 

res_scale <- contrast(
  fit_ipw_cont,
  interventions = list(doubled = scale_by(2), observed = NULL),
  reference = "observed",
  type = "difference"
)
res_scale

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: difference  ·  CI method: sandwich  ·  N: 2000

Intervention means
intervention estimate se ci_lower ci_upper
doubled 4.697 0.87 2.993 6.402
observed 3.083 0.068 2.949 3.217
Contrasts
comparison estimate se ci_lower ci_upper
doubled vs observed 1.615 0.863 -0.077 3.306

Structural truth: E[Y(A+1)] - E[Y(A)] = 2 (the coefficient on A).

See vignette("interventions") for a full tour of shift, scale, dynamic, and IPSI interventions with side-by-side g-comp and IPW examples.

Count treatment (Poisson / negative binomial)

Integer-valued treatments — dose levels, event counts, visit counts — default to a Gaussian treatment density, which works when the support is dense enough that the Normal approximation is harmless (age in years, cigarettes/day). For sparse count data (0–5 doses), the Gaussian pdf evaluates between integers that carry zero probability, producing badly calibrated density ratios.

Pass propensity_family = "poisson" or propensity_family = "negbin" to opt into a Poisson or negative binomial treatment density model:

DGM. A confounder \(L\) affects both a Poisson-distributed count treatment \(A\) and a continuous outcome \(Y\). The treatment is generated from a log-linear Poisson model and the outcome is linear in \(A\) and \(L\).

\[ \begin{aligned} L &\sim N(0, 1) \\ A \mid L &\sim \text{Poisson}\!\bigl(\exp(0.5\, L)\bigr) \\ Y &= 2 + 1.5\, A + L + \varepsilon, \quad \varepsilon \sim N(0, 1) \end{aligned} \]

Because \(\partial Y / \partial A = 1.5\), a unit shift \(A \mapsto A + 1\) has true causal effect \(\mathbb{E}[Y(A + 1)] - \mathbb{E}[Y(A)] = 1.5\).

Code 
set.seed(6)
n <- 3000
L <- rnorm(n)
A <- rpois(n, exp(0.5 * L))
Y <- 2 + 1.5 * A + L + rnorm(n)
count_df <- data.frame(Y = Y, A = A, L = L)

fit_pois <- causat(
  count_df,
  outcome = "Y",
  treatment = "A",
  confounders = ~ L,
  estimator = "ipw",
  propensity_family = "poisson",
  model_fn = stats::glm,
  propensity_model_fn = stats::glm
)

res_pois <- contrast(
  fit_pois,
  interventions = list(plus1 = shift(1), observed = NULL),
  reference = "observed",
  type = "difference"
)
res_pois

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: difference  ·  CI method: sandwich  ·  N: 3000

Intervention means
intervention estimate se ci_lower ci_upper
plus1 5.267 0.069 5.133 5.402
observed 3.692 0.049 3.595 3.788
Contrasts
comparison estimate se ci_lower ci_upper
plus1 vs observed 1.576 0.05 1.478 1.674

Structural truth: E[Y(A+1)] - E[Y(A)] = 1.5.

For negative binomial (which nests Poisson when theta is large), pass propensity_family = "negbin". MASS::glm.nb is auto-selected as the propensity fitter:

Code 
fit_nb <- causat(
  count_df,
  outcome = "Y",
  treatment = "A",
  confounders = ~ L,
  estimator = "ipw",
  propensity_family = "negbin",
  model_fn = stats::glm,
  propensity_model_fn = MASS::glm.nb
)
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached
#> Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
#> control$trace > : iteration limit reached

res_nb <- contrast(
  fit_nb,
  interventions = list(plus1 = shift(1), observed = NULL),
  reference = "observed",
  type = "difference"
)
res_nb

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: difference  ·  CI method: sandwich  ·  N: 3000

Intervention means
intervention estimate se ci_lower ci_upper
plus1 5.261 0.068 5.127 5.395
observed 3.692 0.049 3.595 3.788
Contrasts
comparison estimate se ci_lower ci_upper
plus1 vs observed 1.57 0.05 1.472 1.667

scale_by() also works on count treatments, but with a stricter constraint: the density ratio evaluates at A / factor, so A / factor must be a non-negative integer for every observed value. In practice, this means scale_by() on counts is only feasible when all observed counts are multiples of the factor.

ImportantInteger constraints on count interventions

Only shift() with an integer delta and scale_by() where A / factor is integer for every observation are allowed. static(), dynamic(), threshold(), and ipsi() are rejected because the count pmf is zero at non-integer arguments and the HT indicator / IPSI closed form require binary structure. Use estimator = "gcomp" for those intervention types on count data.

GAM propensity model

Pass propensity_model_fn = mgcv::gam to use a GAM for the treatment density model while keeping the default stats::glm for the outcome MSM:

Code 
fit_gam_ps <- causat(
  nhefs_complete,
  outcome = "wt82_71",
  treatment = "qsmk",
  confounders = ~ sex + s(age) + race + factor(education) +
    s(smokeintensity) + s(smokeyrs) + factor(exercise) +
    factor(active) + s(wt71),
  estimator = "ipw",
  model_fn = stats::glm,
  propensity_model_fn = mgcv::gam
)

res_gam_ps <- contrast(
  fit_gam_ps,
  interventions = list(quit = static(1), continue = static(0)),
  reference = "continue",
  type = "difference"
)
res_gam_ps

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: difference  ·  CI method: sandwich  ·  N: 1566

Intervention means
intervention estimate se ci_lower ci_upper
quit 5.308 0.44 4.446 6.17
continue 1.811 0.216 1.387 2.234
Contrasts
comparison estimate se ci_lower ci_upper
quit vs continue 3.497 0.482 2.553 4.441

External (survey) weights

Pre-computed survey weights are passed via weights =; causatr validates them via check_weights() and multiplies them element-wise with the estimated density-ratio weights. The combined weight drives both the per-intervention MSM fit and the influence-function sandwich variance.

DGM. A simple confounded binary treatment with a continuous outcome. Random survey weights (independent of the data-generating process) are attached to illustrate external weight handling.

\[ \begin{aligned} L &\sim N(0, 1) \\ A \mid L &\sim \text{Bernoulli}\!\bigl(\text{logit}^{-1}(0.5\, L)\bigr) \\ Y &= 2 + 3\, A + 1.5\, L + \varepsilon, \quad \varepsilon \sim N(0, 1) \\ w_i &\stackrel{\text{iid}}{\sim} \text{Uniform}(0.5,\; 2.0) \end{aligned} \]

True ATE: \(\mathbb{E}[Y^1] - \mathbb{E}[Y^0] = 3\).

Code 
set.seed(6)
n <- 2000
L <- rnorm(n)
A <- rbinom(n, 1, plogis(0.5 * L))
Y <- 2 + 3 * A + 1.5 * L + rnorm(n)
w <- runif(n, 0.5, 2.0)
sv_df <- data.frame(Y = Y, A = A, L = L)

fit_ipw_sv <- causat(
  sv_df,
  outcome = "Y",
  treatment = "A",
  confounders = ~ L,
  estimator = "ipw",
  weights = w,
  model_fn = stats::glm,
  propensity_model_fn = stats::glm
)

res_ipw_sv <- contrast(
  fit_ipw_sv,
  interventions = list(a1 = static(1), a0 = static(0)),
  reference = "a0",
  type = "difference"
)
res_ipw_sv

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: difference  ·  CI method: sandwich  ·  N: 2000

Intervention means
intervention estimate se ci_lower ci_upper
a1 5.079 0.051 4.979 5.18
a0 1.928 0.051 1.829 2.028
Contrasts
comparison estimate se ci_lower ci_upper
a1 vs a0 3.151 0.049 3.054 3.248

Direct survey::svydesign integration

causat(weights = ...) also accepts a survey::svydesign object directly — stats::weights(design) is used to extract the sampling weights and the design’s first-stage PSU is auto-adopted as the contrast-time cluster (see below). Explicit cluster = overrides the adoption; single-PSU designs (svydesign(ids = ~1, ...)) propagate only the weights.

Code 
library(survey)
des <- svydesign(ids = ~psu, weights = ~pw, data = d)
fit <- causat(d, outcome = "Y", treatment = "A",
              confounders = ~ L, estimator = "ipw",
              weights = des,
              model_fn = stats::glm, propensity_model_fn = stats::glm)

Cluster-robust sandwich

For design clusters (site, household, PSU, repeated measures) the IPW sandwich should account for within-cluster comovement of the per-individual IFs. causat(cluster = "col") preserves the column through prepare_data() and stashes it on the fit; contrast() defaults to that cluster and runs the sum-within-cluster-then-square IF aggregation. Numerically equivalent to sandwich::vcovCL(type = "HC0") on the final MSM-fit → predict-then-average step, up to the usual \(G/(G-1)\) small-cluster factor.

Code 
fit_cl <- causat(d, outcome = "Y", treatment = "A",
                 confounders = ~ L, estimator = "ipw",
                 cluster = "site",
                 model_fn = stats::glm, propensity_model_fn = stats::glm)
contrast(fit_cl, list(a1 = static(1), a0 = static(0)))

Comparing estimands

Code 
est_df <- data.frame(
  estimand = c("ATE", "ATT", "ATC"),
  estimate = c(
    res_ate_sw$contrasts$estimate[1],
    res_att$contrasts$estimate[1],
    res_atc$contrasts$estimate[1]
  ),
  ci_lower = c(
    res_ate_sw$contrasts$ci_lower[1],
    res_att$contrasts$ci_lower[1],
    res_atc$contrasts$ci_lower[1]
  ),
  ci_upper = c(
    res_ate_sw$contrasts$ci_upper[1],
    res_att$contrasts$ci_upper[1],
    res_atc$contrasts$ci_upper[1]
  )
)

tinyplot(
  estimate ~ estimand,
  data = est_df,
  type = "pointrange",
  ymin = est_df$ci_lower,
  ymax = est_df$ci_upper,
  xlab = "Estimand",
  ylab = "Effect on weight change (kg)",
  main = "IPW estimates by estimand"
)
abline(h = 0, lty = 2, col = "grey40")

Point estimates and confidence intervals for ATE, ATT, and ATC estimated by IPW.

Effect modification with by

When confounders includes an A:modifier interaction term, IPW expands the per-intervention MSM from Y ~ 1 to Y ~ 1 + modifier. The modifier main effects let predict() return stratum-specific counterfactual means; the treatment itself is still absorbed by the density-ratio weights. Use by = "modifier" in contrast() to recover heterogeneous treatment effects.

Here we add a qsmk:sex interaction to examine how the effect of quitting smoking on weight change differs by sex. The propensity model automatically strips the interaction term (it only enters the MSM, not the treatment model).

Code 
fit_ipw_hte <- causat(
  nhefs_complete,
  outcome = "wt82_71",
  treatment = "qsmk",
  confounders = ~ sex + age + I(age^2) + race + factor(education) +
    smokeintensity + I(smokeintensity^2) + smokeyrs + I(smokeyrs^2) +
    factor(exercise) + factor(active) + wt71 + I(wt71^2) +
    qsmk:sex,
  estimator = "ipw",
  model_fn = stats::glm,
  propensity_model_fn = stats::glm
)

res_ipw_by_sex <- contrast(
  fit_ipw_hte,
  interventions = list(quit = static(1), continue = static(0)),
  reference = "continue",
  type = "difference",
  ci_method = "sandwich",
  by = "sex"
)
res_ipw_by_sex

Estimator: ipw  ·  Estimand: ATE  ·  Contrast: difference  ·  CI method: sandwich  ·  N: 1566

Intervention means
intervention estimate se ci_lower ci_upper by n_by
quit 5.367 0.562 4.266 6.468 Male 762
continue 1.798 0.302 1.206 2.389 Male 762
quit 5.21 0.717 3.805 6.615 Female 804
continue 1.78 0.32 1.153 2.406 Female 804
Contrasts
comparison estimate se ci_lower ci_upper by n_by
quit vs continue 3.569 0.633 2.328 4.81 Male 762
quit vs continue 3.43 0.777 1.907 4.953 Female 804
Note

The A:modifier term does not enter the propensity model — it would create a post-treatment variable problem. build_ps_formula() strips interaction terms involving the treatment automatically. The modifier enters only the per-intervention MSM as a main effect, enabling stratum-specific counterfactual mean estimation under the density-ratio weights.

Tidy and glance

causatr results work with the broom ecosystem via tidy() and glance():

Code 
tidy(res_ate_sw)
#>               term estimate std.error     type conf.low conf.high
#> 1 quit vs continue 3.499174 0.4902045 contrast 2.538391  4.459958
glance(res_ate_sw)
#>   estimator estimand contrast_type ci_method    n n_interventions
#> 1       ipw      ATE    difference  sandwich 1566               2

Forest plot

The plot() method produces a forest plot using the forrest package.

Code 
plot(res_ate_sw)

Forest plot of the IPW-estimated ATE of quitting smoking on weight change.

Diagnostics

After fitting an IPW model, use diagnose() to assess positivity and covariate balance. The balance diagnostics use cobalt to compute standardised mean differences (SMD) before and after weighting.

Code 
diag <- diagnose(fit_ipw)
#> Note: `s.d.denom` not specified; assuming "pooled".
diag
#> <causatr_diag>
#>  Estimator:ipw
#>  Treatment: binary
#>  Estimand:  ATE
#> 
#> Positivity:
#>       statistic        value
#>          <char>        <num>
#>             min 2.667709e-07
#>             q25 1.739280e-01
#>          median 2.356540e-01
#>             q75 3.232473e-01
#>             max 8.723380e-01
#>   n_below_lower 6.000000e+00
#>   n_above_upper 0.000000e+00
#>    n_violations 6.000000e+00
#>  pct_violations 3.800000e-01
#> 
#> Covariate balance:
#> Balance Measures
#>                         Type Diff.Un     M.Threshold.Un V.Ratio.Un
#> sex_Female            Binary -0.1603 Not Balanced, >0.1          .
#> age                  Contin.  0.2820 Not Balanced, >0.1     1.0731
#> I(age^2)             Contin.  0.2816 Not Balanced, >0.1     1.1632
#> race                  Binary -0.1771 Not Balanced, >0.1          .
#> factor(education)_0   Binary  0.0439     Balanced, <0.1          .
#> factor(education)_1   Binary -0.1018 Not Balanced, >0.1          .
#> factor(education)_2   Binary  0.0797     Balanced, <0.1          .
#> factor(education)_3   Binary  0.0234     Balanced, <0.1          .
#> factor(education)_4   Binary -0.0129     Balanced, <0.1          .
#> factor(education)_5   Binary -0.0615     Balanced, <0.1          .
#> factor(education)_6   Binary  0.0496     Balanced, <0.1          .
#> factor(education)_7   Binary  0.0025     Balanced, <0.1          .
#> factor(education)_8   Binary  0.0397     Balanced, <0.1          .
#> factor(education)_9   Binary  0.0173     Balanced, <0.1          .
#> factor(education)_10  Binary -0.0826     Balanced, <0.1          .
#> factor(education)_11  Binary -0.1044 Not Balanced, >0.1          .
#> factor(education)_12  Binary -0.0472     Balanced, <0.1          .
#> factor(education)_13  Binary  0.0089     Balanced, <0.1          .
#> factor(education)_15  Binary -0.0665     Balanced, <0.1          .
#> factor(education)_16  Binary  0.1354 Not Balanced, >0.1          .
#> factor(education)_17  Binary  0.0890     Balanced, <0.1          .
#> smokeintensity       Contin. -0.2167 Not Balanced, >0.1     1.1679
#> I(smokeintensity^2)  Contin. -0.1289 Not Balanced, >0.1     1.1519
#> smokeyrs             Contin.  0.1589 Not Balanced, >0.1     1.1846
#> I(smokeyrs^2)        Contin.  0.1789 Not Balanced, >0.1     1.3279
#> factor(exercise)_0    Binary -0.1237 Not Balanced, >0.1          .
#> factor(exercise)_1    Binary  0.0398     Balanced, <0.1          .
#> factor(exercise)_2    Binary  0.0568     Balanced, <0.1          .
#> factor(active)_0      Binary -0.0718     Balanced, <0.1          .
#> factor(active)_1      Binary  0.0268     Balanced, <0.1          .
#> factor(active)_2      Binary  0.0740     Balanced, <0.1          .
#> wt71                 Contin.  0.1332 Not Balanced, >0.1     1.0606
#> I(wt71^2)            Contin.  0.1272 Not Balanced, >0.1     1.0876
#> 
#> Balance tally for mean differences
#>                    count
#> Balanced, <0.1        19
#> Not Balanced, >0.1    14
#> 
#> Variable with the greatest mean difference
#>  Variable Diff.Un     M.Threshold.Un
#>       age   0.282 Not Balanced, >0.1
#> 
#> Sample sizes
#>     Control Treated
#> All    1163     403
#> 
#> Weight distribution:
#>    group     n     mean        sd      min       max       ess
#>   <char> <int>    <num>     <num>    <num>     <num>     <num>
#>  treated   403 3.867604 2.0214552 1.146345 18.317409  316.6996
#>  control  1163 1.346257 0.2502549 1.000000  3.517515 1124.1873
#>  overall  1566 1.995110 1.5204893 1.000000 18.317409  990.8646

Love plot

A Love plot visualises covariate balance. Covariates with absolute SMD below the threshold (dashed line at 0.1) are considered well-balanced.

Code 
plot(diag)
#> Note: `s.d.denom` not specified; assuming "pooled".

Love plot showing covariate balance before and after IPW weighting.

Weight distribution

Extreme weights can inflate variance. The weight summary shows the effective sample size (ESS) — a large drop from the nominal sample size suggests influential weights.

Code 
diag$weights
#>      group     n     mean        sd      min       max       ess
#>     <char> <int>    <num>     <num>    <num>     <num>     <num>
#> 1: treated   403 3.867604 2.0214552 1.146345 18.317409  316.6996
#> 2: control  1163 1.346257 0.2502549 1.000000  3.517515 1124.1873
#> 3: overall  1566 1.995110 1.5204893 1.000000 18.317409  990.8646

Weight truncation

When propensity models yield extreme density ratios — due to near-positivity violations, heavy confounding, or the cumulative product over many components or time periods — the resulting IPW estimator can be unstable with inflated variance. The trim argument on contrast() winsorizes density-ratio weights at a given quantile (Cole & Hernán 2008; Spreafico et al. 2025). This is a per-density-ratio operation applied before any multivariate or longitudinal product, following the approach used internally by lmtp (default 0.999).

Code 
# Without truncation
res_raw <- contrast(
  fit_ipw, interventions = list(a1 = static(1), a0 = static(0)),
  ci_method = "sandwich"
)

# With truncation at the 99.9th percentile of the density ratio
res_trim <- contrast(
  fit_ipw, interventions = list(a1 = static(1), a0 = static(0)),
  ci_method = "sandwich", trim = 0.999
)
Truncation Estimate SE
None (trim = 1) -3.5 0.49
trim = 0.999 -3.5 0.486

Truncation clips individual density ratios at their trim-th quantile (upper-tail only). Key design choices:

  • Per-component, pre-product: for multivariate IPW, each of the \(K\) per-component ratios is truncated before the \(\prod_k\) product. For longitudinal IPW, per-period ratios are truncated before the cumulative \(\prod_t\) product.
  • Not applied to sampling weights (transport) or IPCW weights (censoring) — those define the target population, not causal reweighting.
  • Sandwich SE reflects the truncated weights: the threshold is held fixed at the point estimate, so numDeriv perturbations evaluate the IF on the truncated weight surface.
  • Bootstrap recomputes the truncation threshold on each replicate, capturing the full variance including the effect of truncation.

The default trim = 1 applies no truncation and reproduces the standard IPW estimator.

Multivariate (joint) treatment

estimator = "ipw" accepts multiple treatment columns treatment = c("A1", "A2", ...). The engine factorises the joint conditional density via the chain rule

\[f(A_1, \ldots, A_K \mid L) = \prod_{k=1}^{K} f_k(A_k \mid A_{1..k-1}, L)\]

and fits one univariate density model per component (sequential factorisation). The joint density-ratio weight is the product of \(K\) per-component ratios, with both numerator and denominator of the \(k\)-th factor conditioning on observed upstream treatment values — only the \(k\)-th argument (\(A_k\) vs \(d_k^{-1}(A_k)\)) changes:

\[w_k(\alpha_k) = \frac{f_k\bigl(d_k^{-1}(A_k) \mid A_{1..k-1}^{\mathrm{obs}}, L; \alpha_k\bigr) \cdot |\mathrm{Jac}|}{f_k\bigl(A_k \mid A_{1..k-1}^{\mathrm{obs}}, L; \alpha_k\bigr)}.\]

This is the sequential MTP estimand of Díaz, Williams, Hoffman, Schenck (2023), the standard in the causal-inference MTP literature. Sandwich variance routes through a stacked propensity bread that is block-diagonal across the \(K\) models (each component is fit independently), so the propensity correction sums \(K\) single-model corrections.

Multivariate gcomp implements a different estimand. Multivariate gcomp computes \(E[Y^d] = E\bigl\{E[Y \mid A_1 = d_1(A_1^{\mathrm{obs}}), A_2 = d_2(A_2^{\mathrm{obs}}), L]\bigr\}\) under deterministic joint transformation: each individual’s counterfactual is \((d_1(A_1^{\mathrm{obs}}), d_2(A_2^{\mathrm{obs}}))\), applied simultaneously per individual. For static-only interventions this coincides with the sequential MTP estimand; for non-static interventions on non-final components the two can differ by the upstream \(\to\) downstream cross-dependence. Users who need gcomp-IPW triangulation on non-static multivariate policies should be aware of this estimand divergence.

The intervention is supplied as a named list, one entry per treatment column, each entry a causatr_intervention:

Code 
set.seed(42)
n <- 3000
L  <- rnorm(n)
A1 <- rbinom(n, 1, plogis(0.3 * L))
A2 <- rbinom(n, 1, plogis(-0.3 * L))
Y  <- 2 + 1.5 * A1 + 1.0 * A2 - 0.5 * L + rnorm(n)
df_mv <- data.frame(L = L, A1 = A1, A2 = A2, Y = Y)

fit_mv <- causat(
  df_mv, "Y", c("A1", "A2"), ~ L,
  estimator = "ipw",
  model_fn = stats::glm,
  propensity_model_fn = stats::glm
)
res_mv <- contrast(
  fit_mv,
  interventions = list(
    both    = list(A1 = static(1), A2 = static(1)),
    a1_only = list(A1 = static(1), A2 = static(0)),
    a2_only = list(A1 = static(0), A2 = static(1)),
    neither = list(A1 = static(0), A2 = static(0))
  ),
  reference = "neither"
)
tt(tidy(res_mv))
term estimate std.error type conf.low conf.high
both vs neither 2.5763591 0.05232820 contrast 2.473798 2.678920
a1_only vs neither 1.5722688 0.05622049 contrast 1.462079 1.682459
a2_only vs neither 0.9732751 0.05569496 contrast 0.864115 1.082435

Truth on this DGP: \(E[Y(1,1)] = 4.5\), \(E[Y(0,0)] = 2.0\), \(\mathrm{ATE}(\text{both vs neither}) = 2.5\). Each component intervention uses any of the regular constructors (static(), shift(), scale_by(), dynamic()); IPSI in any component is rejected because Kennedy’s closed form has no joint analogue. Component types can mix freely: binary × continuous, binary × categorical (via nnet::multinom), binary × Poisson / negative binomial (via propensity_family = c("", "poisson") / c("", "negbin")), continuous × continuous, and up to \(K\)-component combinations. Effect modification A_k:modifier is supported — per-component propensity formulas strip all treatment-touching terms and the MSM expands to \(Y \sim 1 + \mathrm{modifier\_main\_effects}\). Multivariate matching remains rejected because MatchIt is binary-only.

Stabilized weights

Pass stabilize = "marginal" to swap the per-component numerator for a separately-fit model \(g_k(A_k \mid A_{1..k-1})\) that drops \(L\) from the conditioning (Robins, Hernán, Brumback 2000, extended to multivariate). The per-component weight becomes

\[ w_k^{s}(\alpha_k) = \frac{g_k\bigl(d_k^{-1}(A_k) \mid A_{1..k-1}\bigr) \cdot |\mathrm{Jac}|}{f_k(A_k \mid A_{1..k-1}, L; \alpha_k)}. \]

Code 
fit_mv_stab <- causat(
  df_mv, "Y", c("A1", "A2"), ~ L,
  estimator = "ipw", stabilize = "marginal",
  model_fn = stats::glm, propensity_model_fn = stats::glm
)
res_mv_stab <- contrast(
  fit_mv_stab,
  interventions = list(
    both    = list(A1 = static(1), A2 = static(1)),
    neither = list(A1 = static(0), A2 = static(0))
  ),
  reference = "neither"
)
tt(tidy(res_mv_stab))
term estimate std.error type conf.low conf.high
both vs neither 2.576359 0.0523282 contrast 2.473798 2.67892

Whether stabilization reduces or inflates the sandwich SE depends on whether the marginal density of \(A_k\) has lower variance than the full-\(L\) conditional — it’s not universally beneficial, just a common variance-reduction lever when \(L\) strongly predicts \(A_k\). Numerator-model uncertainty is held fixed in the sandwich (nuisance-fixed convention, same as \(\sigma\) / \(\theta\)); bootstrap refits both models and captures full variance.

Mixed-type example

Mixed-type treatment (binary + continuous) just works:

Code 
set.seed(7)
n <- 3000
L <- rnorm(n)
A1 <- rbinom(n, 1, 0.5)
A2 <- 5 + 0.5 * L + rnorm(n)
Y  <- 1 + 2 * A1 + 0.3 * A2 - 0.5 * L + rnorm(n)
df_mix <- data.frame(L = L, A1 = A1, A2 = A2, Y = Y)

fit_mix <- causat(
  df_mix, "Y", c("A1", "A2"), ~ L,
  estimator = "ipw",
  model_fn = stats::glm,
  propensity_model_fn = stats::glm
)
res_mix <- contrast(
  fit_mix,
  interventions = list(
    treat   = list(A1 = static(1), A2 = shift(-2)),
    control = list(A1 = static(0), A2 = shift(-2))
  ),
  reference = "control"
)
tt(tidy(res_mix))
term estimate std.error type conf.low conf.high
treat vs control 1.810289 0.2155841 contrast 1.387752 2.232826

The A2 shift is applied uniformly under both arms, so the contrast recovers the marginal A1 effect (truth = 2).

Multivariate + truncation

With \(K\) components, the product weight \(\prod_k w_k\) can grow rapidly. trim truncates each per-component ratio before the product:

Code 
res_mv_trim <- contrast(
  fit_mv,
  interventions = list(
    both    = list(A1 = static(1), A2 = static(1)),
    neither = list(A1 = static(0), A2 = static(0))
  ),
  reference = "neither",
  trim = 0.99
)
tt(tidy(res_mv_trim))
term estimate std.error type conf.low conf.high
both vs neither 2.570902 0.05261264 contrast 2.467784 2.674021

Summary of covered combinations

Legend. ✅ covered and truth-pinned in tests · 🟡 smoke test only · ⛔ rejected with an informative error.

Treatment Outcome Intervention Estimand Contrast Inference Weights Status
Binary Continuous Static ATE Difference Sandwich none ✅
Binary Continuous Static ATE Difference Bootstrap none ✅
Binary Continuous Static ATT Difference Sandwich none ✅
Binary Continuous Static ATC Difference Sandwich none ✅
Binary Continuous Static ATT Difference Bootstrap none ✅
Binary Continuous Static ATE Difference Sandwich survey ✅
Binary Binary Static ATE Difference Sandwich none ✅
Binary Binary Static ATE Difference Bootstrap none ✅
Binary Binary Static ATE Ratio Sandwich none ✅
Binary Binary Static ATE OR Sandwich none ✅
Binary Binary Static ATE Ratio Bootstrap none ✅
Binary Continuous Dynamic ATE Difference Sandwich none ✅
Binary Continuous IPSI(\(\delta\)) ATE Difference Sandwich none ✅
Continuous Continuous Shift ATE Difference Sandwich none ✅
Continuous Continuous Scale ATE Difference Sandwich none ✅
Categorical (k>2) Continuous Static ATE Difference Sandwich none ✅
Categorical (k>2) Continuous Dynamic ATE Difference Sandwich none 🟡
Count (Poisson) Continuous Shift ATE Difference Sandwich none ✅
Count (NB) Continuous Shift ATE Difference Sandwich none ✅
Any — Threshold — — — — ⛔ (use gcomp)
Continuous — Static — — — — ⛔ (use shift/scale)
Continuous — Dynamic — — — — ⛔ (use shift/scale)
Binary Gaussian Static (with A:modifier) ATE Difference Sandwich / Bootstrap none ✅
Multi (bin × bin) Gaussian Static ATE Difference Sandwich / Bootstrap none ✅
Multi (bin × cont) Gaussian Static + Shift ATE Difference Sandwich none ✅
Multi (cont × cont) Gaussian Shift + Shift ATE Difference Sandwich none ✅ (seq-MTP)
Multi (bin × bin) Binary Static ATE Diff / Ratio / OR Sandwich none ✅
Multi (3+ binary) Gaussian Static ATE Difference Sandwich none ✅
Multi (bin × bin) Gaussian Static (with A_k:modifier) by-stratified Difference Sandwich none ✅
Multi (bin × categorical) Gaussian Static ATE Difference Sandwich none ✅
Multi (bin × Poisson) Gaussian Static + Shift ATE Difference Sandwich none ✅
Multi (bin × negbin) Gaussian Static + Shift ATE Difference Sandwich none ✅
Multi (any) — Any + stabilize = "marginal" ATE — Sandwich / Bootstrap none ✅
Multi + IPSI component — — — — — — ⛔
Longitudinal Continuous Static / Shift / Scale / Dynamic ATE Difference Sandwich / Bootstrap none ✅ (see vignette("longitudinal"))
Longitudinal Continuous Static + stabilize=“marginal” ATE Difference Sandwich none ✅
Longitudinal Any IPSI — — — — ⛔ (use point IPW)
Multivariate longitudinal Continuous Static / Shift ATE Difference Sandwich / Bootstrap none / marginal ✅ (see vignette("longitudinal"))
Multivariate longitudinal Gaussian Static + EM (A_k:modifier) by-stratified Difference Sandwich none / marginal ✅ vs ICE
Multivariate longitudinal Any IPSI component — — — — ⛔ (use ICE)

See FEATURE_COVERAGE_MATRIX.md for the authoritative coverage status of every method × treatment × outcome × intervention × variance combination.

References

Hernán MA, Robins JM (2025). Causal Inference: What If. Chapman & Hall/CRC. Chapter 12: IP weighting and marginal structural models.

Kennedy EH (2019). Nonparametric causal effects based on incremental propensity score interventions. Journal of the American Statistical Association 114:645–656.

Díaz I, Williams N, Hoffman KL, Schenck EJ (2023). Non-parametric causal effects based on longitudinal modified treatment policies. Journal of the American Statistical Association 118:846–857.