Getting started with CausalState

Overview

CausalState estimates the causal effect of a longitudinal modified treatment policy (MTP) on an outcome in care-episode data where patients can transition out of an active state (e.g. ICU discharge or death). The distinctive feature is that the MTP can shift the transition dynamics themselves — not only the terminal outcome. The MTP framework is due to Díaz & van der Laan (2012) and Haneuse & Rotnitzky (2013).

Two main estimators are provided, both from Luedtke et al. (2017/2018) and both sequentially doubly robust (SDR): consistent whenever, at each time point $t$, either the treatment model $g_t$ or the outcome model $Q_t$ is correctly specified ($2^K$-robust, Definition 2 of Luedtke et al.):

sdr() — Sequential Doubly Robust estimator. Applies an EIF-based pseudo-outcome update within the backward Q-regression (Díaz et al. 2021).
itmle() — Infinite-dimensional TMLE (iTMLE), Algorithm 4 of Luedtke et al. Shares the same backward Q-regression as SDR; differs only in the update step, which uses an infinite-dimensional TMLE fluctuation.

qreg() is a pure Q-recursion plug-in (no update step, no DR guarantees) included as a sensitivity check when density-ratio weights are extreme.

The workflow is always:

density_ratio()  →  sdr() / itmle() / qreg()

density_ratio() must be run first; all downstream estimators inherit its fold structure and weights.

Data structure

Data must be in long format: one row per subject per time point, covering only periods while the subject is in the active state. The last row for each subject carries their transition outcome.

Column	Type	Meaning
`id`	integer / character	Subject identifier
`time`	integer 1, 2, …	Time index
`alive`	0 / 1	1 = alive at end of period
`in_state`	0 / 1	1 = still in active state (e.g. ICU)
treatment	numeric	One or more treatment variables (`a_names`)
outcome	numeric	Outcome variable (`y`)
covariates	numeric	Baseline and time-varying covariates

A subject exits when in_state == 0 on their last row (alive = 1 for discharge, alive = 0 for death). Subjects still in the active state at tmax have in_state == 1 on their last row.

Simulated example data

The DGP below generates a simple ICU panel with binary treatment and a binary outcome. Copy it directly or adapt it for your own simulation studies.

sim_panel <- function(n = 500L, tmax = 4L, seed = 42L) {
  set.seed(seed)
  rows <- vector("list", n)
  for (i in seq_len(n)) {
    age <- round(rnorm(1, 65, 10))
    sex <- rbinom(1, 1, 0.5)
    L1  <- rnorm(1, 0, 1)
    L2  <- rbinom(1, 1, 0.4)
    pat <- list()
    for (t in seq_len(tmax)) {
      A     <- rbinom(1, 1, plogis(0.3 * L1 - 0.4 + 0.2 * sex))
      p_die <- plogis(-4.0 + 0.3 * L1 - 0.1 * age / 10)
      p_dc  <- plogis(-2.5 + 0.5 * A  - 0.2 * L2)
      u     <- runif(1)
      if      (u < p_die)        { alive <- 0L; in_state <- 0L }
      else if (u < p_die + p_dc) { alive <- 1L; in_state <- 0L }
      else                       { alive <- 1L; in_state <- 1L }
      py <- if (!alive)      plogis(-3.0 + 0.1 * L1)
            else if (!in_state) plogis(1.5 + 0.2 * L1 - 0.1 * age / 10 + 0.3 * A)
            else             plogis(-0.5 + 0.4 * A - 0.2 * L1 + 0.1 * L2)
      Y <- rbinom(1, 1, py)
      pat[[length(pat) + 1L]] <- data.frame(
        id = i, time = t, age = age, sex = sex,
        alive = alive, in_state = in_state,
        L1 = L1, L2 = L2, A = A, Y = Y
      )
      if (in_state == 0L) break
      if (t < tmax) {
        L1 <- L1 + rnorm(1, -0.1 * A, 0.3)
        L2 <- rbinom(1, 1, plogis(0.5 * L2 + 0.3 * A - 0.5))
      }
    }
    rows[[i]] <- do.call(rbind, pat)
  }
  do.call(rbind, rows)
}

df <- sim_panel(n = 2000L, tmax = 5L, seed = 1L)
head(df)

Step 1 — Define a policy

A policy is a function (D_block, t, a_names) that returns a data.table with the shifted treatment values for all subjects at time t. Here we apply a soft upward shift: increase treatment probability by 0.3, capped at 1.

policy_up <- function(D_block, t, a_names) {
  out <- D_block[, ..a_names, drop = FALSE]
  out[[a_names[1]]] <- pmin(D_block[[a_names[1]]] + 0.3, 1)
  out
}

Step 2 — Density ratios

density_ratio() fits per-time-point treatment models and returns the instantaneous ratio $r_t = d\tilde{P}(A_t \mid H_t) / dP(A_t \mid H_t)$ comparing the MTP to the natural course.

# For real analyses replace with a richer library:
# e.g. c("SL.glm", "SL.earth", "SL.xgboost", "SL.dbarts")
sl_lib <- c("SL.mean", "SL.glm")

wr <- density_ratio(
  df              = df,
  a_names         = "A",
  tmax            = 5L,
  baseline        = c("age", "sex"),
  tv_names        = c("L1", "L2"),
  sl_g            = sl_lib,
  k               = 1L,
  inner_v         = 5L,
  v               = 5L,
  seed            = 1L,
  id              = "id",
  time            = "time",
  policy_spec_fun = policy_up
)

wr$weights_dt contains one row per subject per time with Rt_t (instantaneous ratio) and global_fold (cross-fitting fold). Inspect wr$sl_summary to check SuperLearner ensemble weights.

Bypass for identical natural and shifted treatment. At any time point where the policy produces no shift for any subject (natural treatment = shifted treatment in every row), density_ratio() skips model fitting entirely and fills $r_t = 1$ for all subjects at that time. This covers policies with a finite intervention window: time points outside the window are filled with 1 automatically.

Weights, intervention windows, and trim

Intervention window shorter than tmax. A common pattern is an MTP that intervenes only during times 1 to $x < \texttt{tmax}$ (e.g. treatment is feasible only in the first few ICU days). density_ratio() is called with this tmax, fits models only where the shift is non-trivial, and fills $r_t = 1$ outside the intervention window. The downstream estimators sdr() and itmle() are still called with the full tmax:

The Q-models covering time points beyond $x$ are estimated under the natural course (weights = 1), but the risk set (who is still in the ICU) continues to evolve and is modelled correctly. For SDR this means the backward recursion propagates through those time points using natural-course predictions; for iTMLE the targeting step uses unit weights there.
Modelling beyond the intervention window is therefore still useful: it correctly accounts for transitions that occur after the policy has ended.

Trimming. Trim is applied globally across the entire weights_dt object at the point it is consumed by sdr() or itmle(). The trim quantile is computed over all time points present in the weight object, not just those within the intervention window or the current tmax. This means:

The trim threshold is identical regardless of whether you call the estimator with tmax = 3 or tmax = 10, as long as they share the same weight object — estimates across horizons are directly comparable.
The same trim value must be passed consistently to all estimators that share a weight object.

Using `dr_sl = TRUE` (Wu-Benkeser metalearner)

By default (dr_sl = FALSE) treatment models are binary classifiers and the density ratio is recovered as $\hat{r} = \hat{p} / (1 - \hat{p})$. Setting dr_sl = TRUE switches to the Wu-Benkeser (2024) metalearner, which minimises a density-ratio loss directly on the simplex.

This requires custom SuperLearner wrappers that return density ratios directly, not class probabilities. The package does not export such wrappers because the right design choices (learner type, regularisation, bandwidth) are analysis-specific. The density_ratio() help page documents the interface a custom wrapper must satisfy.

Step 3 — SDR estimator

res_sdr <- sdr(
  df              = df,
  weight_object   = wr,
  tmax            = 5L,
  id              = "id",
  time            = "time",
  alive           = "alive",
  in_state        = "in_state",
  y               = "Y",
  baseline        = c("age", "sex"),
  tv_names        = c("L1", "L2"),
  a_names         = "A",
  sl_remain       = sl_lib,
  sl_death        = sl_lib,
  sl_recursive    = sl_lib,
  sl_y            = sl_lib,
  outcome_family  = "binomial",
  k               = 1L,
  inner_v         = 5L,
  seed            = 1L,
  policy_spec_fun = policy_up
)

cat(sprintf(
  "SDR  psi = %.3f  (natural = %.3f)  RD = %.3f  SE = %.3f  95%% CI [%.3f, %.3f]\n",
  res_sdr$psi, res_sdr$psi_nat, res_sdr$rd, res_sdr$se,
  res_sdr$psi - 1.96 * res_sdr$se,
  res_sdr$psi + 1.96 * res_sdr$se
))

The return list contains psi (point estimate under MTP), psi_nat / psi_shf (plug-in estimates under natural course and MTP), rd (risk difference psi_shf − psi_nat), se (from the efficient influence curve), ic (per-subject influence curve values), sl_summary, and fold_diag.

Step 4 — Infinite-dimensional TMLE (iTMLE)

iTMLE applies an infinite-dimensional TMLE fluctuation as the update step. It requires sl_tmle: a SuperLearner library of targeting wrappers from the sl_itmle family that handle the logit offset passed as column ._sl_offset in the design matrix. These are the only SL wrappers the package exports, because they must accommodate the offset-as-column structure that standard SuperLearner wrappers do not handle (see ?sl_itmle).

For a quick start, SL.tgt.intercept (standard one-parameter TMLE update) and SL.tgt.glm (GLM fluctuation) are sufficient. The full sl_tmle default vector adds penalised regression and gradient boosting.

res_itmle <- itmle(
  df               = df,
  weight_object    = wr,
  tmax             = 5L,
  id               = "id",
  time             = "time",
  alive            = "alive",
  in_state         = "in_state",
  y                = "Y",
  baseline         = c("age", "sex"),
  tv_names         = c("L1", "L2"),
  a_names          = "A",
  sl_remain        = sl_lib,
  sl_death         = sl_lib,
  sl_recursive     = sl_lib,
  sl_y             = sl_lib,
  sl_tmle          = c("SL.tgt.intercept", "SL.tgt.glm"),
  outcome_family   = "binomial",
  k                = 1L,
  inner_v          = 5L,
  v_target_itmle   = 5L,
  v_sl_inner_itmle = 5L,
  seed             = 1L,
  policy_spec_fun  = policy_up
)

cat(sprintf(
  "iTMLE  psi = %.3f  SE = %.3f  95%% CI [%.3f, %.3f]  targeting gap = %.2e\n",
  res_itmle$psi, res_itmle$se,
  res_itmle$ci[1], res_itmle$ci[2],
  res_itmle$targeting_gap
))

targeting_gap (mean EIF after targeting) should be near zero; a large value indicates the targeting step did not converge.

Step 5 — Q-recursion sensitivity check (`qreg`)

qreg() runs the same backward Q-regression without any update step. The point estimate $\hat\Psi_Q = n^{-1}\sum_i \hat{Q}_1(H_{i1})$ depends only on the Q-models and is independent of density-ratio weights.

It is not sequentially doubly robust and carries first-order bias when Q-models are misspecified. Use it as a sensitivity check: if qreg and sdr agree, the DR update and density-ratio weighting are not decisive. A large discrepancy suggests extreme weights or Q-model misspecification and warrants closer inspection of ESS diagnostics.

Passing weight_object is optional but enables an EIF-based SE alongside the always-computed naive SE (which treats the fitted Q as fixed and underestimates uncertainty).

res_qreg <- qreg(
  df              = df,
  weight_object   = wr,
  tmax            = 5L,
  id              = "id",
  time            = "time",
  alive           = "alive",
  in_state        = "in_state",
  y               = "Y",
  baseline        = c("age", "sex"),
  tv_names        = c("L1", "L2"),
  a_names         = "A",
  sl_remain       = sl_lib,
  sl_death        = sl_lib,
  sl_recursive    = sl_lib,
  sl_y            = sl_lib,
  outcome_family  = "binomial",
  k               = 1L,
  inner_v         = 5L,
  seed            = 1L,
  policy_spec_fun = policy_up
)

cat(sprintf(
  "qreg  estimate = %.3f  SE (naive) = %.3f  SE (EIF) = %.3f\n",
  res_qreg$estimate, res_qreg$se_naive, res_qreg$se_eif
))

Diagnostics

All estimators return diagnostic objects alongside the point estimate.

SuperLearner weights — `sl_summary`

A data.table of ensemble weights per fold, time point, and model component (g_remain, g_death_exit, Q_exit, Q_rem). Check that no single learner dominates unexpectedly and that weights are stable across folds.

head(res_sdr$sl_summary)

Per-fold summaries — `fold_diag`

Counts and mean Q predictions per fold and time point. Useful for spotting time points with very few at-risk subjects, which can destabilise model fits.

res_sdr$fold_diag

Branch calibration — `diagnostics$branch_cal` (SDR / qreg)

Compares in-fold training vs. held-out validation performance (AUC, Brier score, MSE) for each model component at each time point. A large training–validation gap signals overfitting.

head(res_sdr$diagnostics$branch_cal)

Targeting diagnostics — `target_diag`, `targeting_gap` (iTMLE only)

target_diag tracks the mean EIF and the fluctuation parameter $\epsilon$ across targeting iterations. targeting_gap summarises convergence — values near zero confirm the targeting step has solved the efficient score equation.

res_itmle$targeting_gap
head(res_itmle$target_diag)

Density-ratio ESS

wr$weights_dt provides Rt_t (instantaneous ratio) and cum_ratio (cumulative product) per subject per time. Monitor the effective sample size (ESS) of the cumulative weights — ESS collapse across time is the main practical failure mode of density-ratio-based estimators and will inflate standard errors even when the point estimate looks stable.

Parallelism

All estimators default to parallel = FALSE and parallel_t = FALSE and run on a single core out of the box.

Two independent levels of parallelism are available, both implemented via parallel::mclapply() (process forking):

Fold-level (parallel = TRUE, fold_workers): each cross-fitting fold runs in a separate forked process.
Within-fold regression (reg_workers): per-time-point regression tasks within each fold are distributed across worker processes.
Time-point level in density_ratio() (parallel_t = TRUE, t_workers): time points are parallelised independently of fold workers.

Because forking is used, this does not work on Windows. On Linux and macOS, set BLAS, xgboost, and dbarts thread counts to 1 before enabling process-level parallelism to avoid oversubscription:

# Example: 5-fold outer parallelism on Linux/macOS
res_sdr_par <- sdr(
  ...,
  parallel     = TRUE,
  fold_workers = 5L,
  reg_workers  = 1L      # keep inner regressions single-threaded
)

References

Bang H, Robins JM (2005). Doubly Robust Estimation in Missing Data and Causal Inference Models. Biometrics 61(4):962–973.

Díaz I, van der Laan MJ (2012). Population Intervention Causal Effects Based on Stochastic Interventions. Biometrics 68(2):541–549.

Haneuse S, Rotnitzky A (2013). Estimation of the Effect of Interventions that Modify the Received Treatment. Statistics in Medicine 32(30):5260–5277.

Luedtke AR, Sofrygin O, van der Laan MJ, Carone M (2017/2018). Sequential Double Robustness in Right-Censored Longitudinal Models. arXiv:1705.02459.

Rotnitzky A, Robins J, Babino L (2017). On the Multiply Robust Estimation of the Mean of the G-Functional. arXiv:1705.08582.

Díaz I, Williams N, Hoffman KL, Schenck EJ (2021). Nonparametric Causal Effects Based on Longitudinal Modified Treatment Policies. JASA 118(542):846–857.

Williams NT, Díaz I (2023). lmtp: An R package for estimating the causal effects of modified treatment policies. Observational Studies.

Wu C, Benkeser D (2024). Nonparametric Efficient Estimation of Marginal Structural Models using Targeted Machine Learning. arXiv:2408.10847.