contrast – causatr

Compute causal contrasts from a fitted model

Description

Standardises outcome model predictions under each named intervention and reports pairwise causal contrasts (differences, ratios, or odds ratios) with uncertainty estimates.

For all three causal estimators, contrast() computes E[Y^a] by setting each individual’s treatment to the intervened value and averaging the fitted outcome model predictions over the target population. The target population is controlled by estimand (or subset for subgroup effects). The estimators differ only in how the outcome model was fitted; variance is computed by the unified influence-function engine variance_if() (R/variance_if.R), which handles all four cases (g-comp, IPW, matching, ICE) via one entry point. See vignettes/variance-theory.qmd for the derivation and the R/variance_if.R roxygen header for the architecture.

“gcomp”: standard glm/gam on the full data; one-model IF correction via prepare_model_if() + apply_model_correction().
“ipw”: density-ratio weighted MSM refit per intervention; IF Channel 2 is the stacked-sandwich MSM-plus-propensity correction assembled by compute_ipw_if_self_contained_one(), with the cross-derivative $A_{\beta\alpha}$ computed numerically by numDeriv::jacobian().
“matching”: glm() on the matched sample with match weights; IF is computed on the matched sample and aggregated cluster-robustly on matched$subclass via vcov_from_if(cluster = …).

Usage

contrast(
  fit,
  interventions,
  type = c("difference", "ratio", "or"),
  estimand = NULL,
  subset = NULL,
  reference = NULL,
  ci_method = c("sandwich", "bootstrap"),
  n_boot = 500L,
  conf_level = 0.95,
  by = NULL,
  parallel = getOption("boot.parallel", "no"),
  ncpus = getOption("boot.ncpus", 1L),
  cluster = NULL,
  treatment_values = NULL,
  trim = 1
)

Arguments

fit A causatr_fit object returned by causat().

interventions

A named list of interventions. Each element must be one of:

A causatr_intervention object created by static(), shift(), dynamic(), scale_by(), threshold(), ipsi(), or stochastic().
NULL, meaning the natural course (observed treatment values are used as-is). The natural course is the standard reference for modified treatment policies on continuous treatments (e.g. shift(-10) vs NULL; Diaz et al. 2023) and for longitudinal dynamic regimes (Hernan & Robins Ch. 21). For binary treatments, the natural-course marginal mean equals E[Y] under the observed treatment mechanism (by consistency); use static(1) vs static(0) for the conventional ATE. Supported for all estimators.
A named list of causatr_intervention objects, one per treatment variable, for multivariate (joint) treatments. Each sub-list must name every treatment variable specified in causat().

Note: intervention support is estimator-specific. gcomp supports static() / shift() / scale_by() / threshold() / dynamic() / stochastic(). IPW supports static() / shift() / scale_by() / dynamic(), plus ipsi() on binary treatment (point and univariate longitudinal) and stochastic() when a density is supplied. matching supports static() only. ipsi() requires a fitted propensity model and a binary treatment; multivariate / stabilized longitudinal IPSI and longitudinal AIPW IPSI are rejected with classed errors.

type Character. The contrast scale: “difference” (default), “ratio”, or “or” (odds ratio). All pairwise contrasts are reported.

estimand Character or NULL. The target estimand: “ATE”, “ATT”, or “ATC”. For estimator = “gcomp”, overrides the estimand stored in fit (allowing one fit to produce multiple estimands). For estimator = “ipw” or “matching”, must match the estimand used at fitting time – changing it aborts with an informative message. If NULL, defaults to fit$estimand. Mutually exclusive with subset.

subset A quoted expression (quote(…)) defining a subgroup to average over instead of an estimand. Evaluated in the context of the fitted data. Works for any treatment type. For example, subset = quote(age > 50) or subset = quote(A == 1) (the latter is equivalent to estimand = “ATT” for binary treatments). Mutually exclusive with estimand.

reference Character or NULL. Name of the reference intervention (the denominator/subtracted value for pairwise contrasts). Default: the first intervention in the list. Only relevant when type = “difference” or “ratio” and there are more than two interventions. For categorical treatments, use this to specify the reference level.

ci_method Character. The variance/CI approach: “sandwich” (default) or “bootstrap”.

n_boot Integer. Number of bootstrap replications when ci_method = “bootstrap”. Default 500.

conf_level Numeric. Confidence level for intervals. Default 0.95.

by Character or NULL. Name of a variable to stratify estimates by (effect modification). If provided, E[Y^a] is computed within each level of by.

parallel Character. Parallelisation backend for bootstrap (ci_method = “bootstrap” only). “no” (default) runs sequentially; “multicore” uses forked processes via parallel::mclapply() (Unix only); “snow” uses socket clusters via parallel::parLapply() (cross-platform); “future” dispatches replicates through future.apply::future_lapply() so any active future::plan() is honoured (local multisession, remote cluster, future.batchtools, etc.). The “no” / “multicore” / “snow” paths go through boot::boot(); “future” bypasses it to let future own the scheduling, which means the ncpus argument is ignored for “future” (the plan’s worker count governs). Ignored when ci_method = “sandwich”.

ncpus Integer. Number of CPU cores for parallel bootstrap. Default getOption(“boot.ncpus”, 1L). Passed directly to boot::boot().

cluster Character or NULL. Name of a column in fit$data identifying cluster membership (e.g. site, household, PSU id). When non-NULL, the sandwich variance is aggregated cluster-robustly: per-individual influence functions are summed within each cluster before squaring (Liang & Zeger 1986), and the resulting $(1/n^2)\sum_c (\sum_{i \in c} \mathrm{IF}_i)^2$ estimator is numerically equivalent to the stacked-EE cluster-robust sandwich and to sandwich::vcovCL() applied on the final predict-then-average step. Works for estimator = “gcomp”, “ipw”, and “ice”; aborts under “matching” (matching already clusters on its own subclass, and layering a user cluster on top would double-count). For ICE with a longitudinal fit, the cluster vector is read from the first-time-point rows. Cluster is unused when ci_method = “bootstrap” for point fits; ICE bootstrap already resamples entire individual trajectories.

treatment_values Numeric vector of length 2 or NULL. Only for estimator = “snm”. When supplied, contrast() returns the population-averaged blip effect (1/n)_i [(a_1, L_i; ) - (a_0, L_i; )] at treatment values c(a0, a1), with delta-method SE. When NULL (the default), returns the raw blip parameter table. Rejected for non-SNM estimators.

trim Numeric scalar in (0, 1]. Density-ratio weight truncation quantile. Per-component density ratios above the trim-th quantile of the fitted weight distribution are winsorized (capped) before entering the Hajek estimator. Default 1 (no truncation). Recommended values: 0.999 (lmtp default) or 0.99 (Spreafico et al. 2025). Only affects IPW and AIPW estimators; ignored for gcomp, matching, and SNM. Reference: Cole & Hernan (2008).

Details

Estimands and standardisation

For each intervention a, contrast() evaluates the intervention function on the target rows to obtain the intervened treatment vector $a(L_i)$, then computes:

$E[Y^a] = (1/|S|) \sum_{i \in S} \hat{E}[Y | A = a(L_i), L_i]$

where S is the set of rows determined by the estimand:

“ATE”: all rows (full population)
“ATT”: rows where A == 1 (observed treated)
“ATC”: rows where A == 0 (observed controls)
subset: rows satisfying the quoted expression

For “gcomp”, one fit supports multiple estimands because the outcome model is the same – only the rows averaged over change. For “ipw” and “matching”, the estimand is baked into the weights/matching and cannot be changed after fitting.

Estimand applicability

Treatment type	ATE	ATT/ATC	subset
Binary point	Yes	Yes	Yes
Continuous point	Yes	No (abort)	Yes
Categorical point	Yes	No (abort)	Yes
Multivariate point	Yes	No (abort)	Yes
Longitudinal	Yes	No (abort)	Yes

Treatment types and intervention support

Method	Treatment types	Intervention types
`“gcomp”`	binary, categorical, continuous, multivariate	all
`“ipw”`	binary, continuous	all deterministic
`“matching”`	binary	`static()` only

Variance estimation

“sandwich”: The unified influence-function engine variance_if(). Per-individual IF = Channel 1 (direct covariate-sampling term) + Channel 2 (one correction per nuisance model). Aggregated via vcov_from_if(), with cluster-robust aggregation for matching. Mathematically equivalent to the stacked M-estimation sandwich; see vignettes/variance-theory.qmd Sections 3-4.
“bootstrap”: Resamples individuals (entire pipeline refitted n_boot times). Respects cluster structure for longitudinal data.
The delta method is applied internally for ratio / odds-ratio contrasts on top of the “sandwich” vcov; it is not a separate ci_method option.

Value

A causatr_result object with slots:

estimates: data.table with one row per intervention: intervention, estimate, se, ci_lower, ci_upper.
contrasts: data.table with one row per pairwise comparison: comparison, estimate, se, ci_lower, ci_upper.
type: Contrast scale.
estimand: “ATE”, “ATT”, “ATC”, or “subset”.
ci_method: Inference approach.
reference: Name of the reference intervention.
interventions: The intervention list.
n: Number of individuals averaged over.
estimator: Causal estimator from the fit.
vcov: Full variance-covariance matrix for all E[Y^a].
call: The original call.

References

Hernan MA, Robins JM (2025). Causal Inference: What If. Chapman & Hall/CRC. Chapters 12-13.

Imai K, King G, Stuart EA (2008). Misunderstandings between experimentalists and observationalists about causal inference. Journal of the Royal Statistical Society Series A 171:481-502.

Zivich PN, Ross RK, Shook-Sa BE, Cole SR, Edwards JK (2024). Empirical sandwich variance estimator for iterated conditional expectation g-computation. Statistics in Medicine 43:5562-5572.

Examples

library("causatr")

data("nhefs", package = "causatr")
fit <- causat(nhefs, outcome = "wt82_71", treatment = "qsmk",
              confounders = ~ sex + age + wt71)

# ATE: mean difference with sandwich SE
result <- contrast(fit,
  interventions = list(quit = static(1), continue = static(0)),
  type = "difference"
)

# ATT: override estimand in contrast() (gcomp only)
result_att <- contrast(fit,
  interventions = list(quit = static(1), continue = static(0)),
  estimand = "ATT"
)

# Subgroup effect: age > 50
result_sub <- contrast(fit,
  interventions = list(quit = static(1), continue = static(0)),
  subset = quote(age > 50)
)

# Continuous treatment: shift with NULL reference
fit_cont <- causat(nhefs, outcome = "wt82_71",
                   treatment = "smokeintensity",
                   confounders = ~ sex + age + wt71)
contrast(fit_cont,
  interventions = list(reduce10 = shift(-10), observed = NULL),
  type = "difference"
)

# Categorical treatment: three arms, reference = "radio"
contrast(fit_cat,
  interventions = list(chemo = static("A"), radio = static("B"),
                       combo = static("C")),
  type = "difference",
  reference = "radio"
)

# Stochastic intervention: random treatment assignment
sampler <- function(data, trt) {
  rbinom(nrow(data), 1, plogis(0.3 * data$age))
}
contrast(fit,
  interventions = list(
    stoch = stochastic(sampler, n_mc = 100L),
    always = static(1)
  ),
  type = "difference"
)