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)
)

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(), or ipsi().
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: Non-static interventions (shift(), scale_by(), threshold(), dynamic(), ipsi()) are only supported for estimator = “gcomp”. For IPW and matching, the weights/matched sets were estimated under the observed treatment distribution; applying a different regime requires re-calling causat() with updated data. ipsi() additionally requires a fitted propensity model and currently aborts for all estimators.

Survival contrasts are not implemented. contrast() on a causatr_fit returned by causat_survival() aborts with an informative error. The pooled logistic fit itself works.

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). Passed directly to boot::boot(). Ignored when ci_method = “sandwich”.

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

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	binary: `static()` / `dynamic()` / `ipsi()`; continuous: `shift()` / `scale_by()`
`“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 (2011). 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"
)