Diagnostics workflow

library(CausalState)
library(SuperLearner)

Overview

CausalState produces a lot of diagnostic output. This vignette walks through a stable, objective workflow for using it. The core idea is that several checks and tuning decisions should happen before the final estimation, and a smaller set of post-estimation checks confirm the fit behaved sensibly.

The recommended flow:

  1. Density ratios — fit the weights first; use weight_diagnostics() and policy_change_summary to judge whether the intervention is supported and to pick trim.
  2. Q calibration — fit sdr(), itmle(), or qreg() under the natural-course policy; use branch_cal_summary() to tune the SuperLearner libraries for the g- and Q-branches. For itmle() this is also where you tune the targeting library sl_tmle.
  3. Final estimation — once weights and Q-branches look reasonable, run sdr() / itmle() under the actual policy of interest.
  4. Post-estimation sanity checksrecursion_diag, target_cal (iTMLE only), and ic_df. These describe what the estimator did with a final fit; they are informative but not part of the tuning loop.

We use sim_bin() throughout for illustration. sim_cont() (continuous treatment) and sim_multi() (any mix of binary + continuous treatments) follow the same shape.

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

1. Density ratios: is the intervention supportable?

Fit the density ratios before touching any outcome model. This step answers two independent questions:

sl_lib <- c("SL.mean", "SL.glm")

policy_bin <- function(D_block, t, a_names) {
  out <- D_block[, ..a_names, drop = FALSE]
  out[[a_names[1]]] <- pmax(
    D_block[[a_names[1]]], as.integer(D_block[["L1"]] > 1.0)
  )
  out
}

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         = 3L,
  v               = 3L,
  seed            = 1L,
  id              = "id",
  time            = "time",
  policy_spec_fun = policy_bin
)

weight_diagnostics()

weight_diagnostics(wr)

The most important columns:

Secondary columns, useful for tuning and plausibility checks:

policy_change_summary

wr$policy_change_summary

Shows the fraction of subjects whose treatment actually changed under the policy at each time-step, plus the mean/quantiles of the observed-vs-shifted difference among those who changed. A policy that changes 100% of subjects by large amounts will have extreme weights by construction — the diagnostics above are then just confirming a design choice rather than a modelling problem.

Picking a trim

weight_diagnostics() accepts trim (default 1, i.e. no trimming) so you can see how a candidate trim reshapes the weight distribution before committing:

weight_diagnostics(wr, trim = 0.99)   # cap Rt_t at the 99th percentile
weight_diagnostics(wr, trim = 0.95)

sdr() and itmle() take the same trim argument. Trimming reduces variance at the cost of a small (usually negligible) bias in the weighted correction — the tradeoff is worth it when a small number of extreme weights are dominating the influence curve.

2. Q calibration: are the outcome models trustworthy?

Once weights look reasonable, fit the estimator under the natural course — i.e. with a policy that returns the observed treatment unchanged. This diagnostic run has two properties that make it uniquely useful for tuning:

Under the actual intervention policy, calibration is unmeasurable in the usual sense — the counterfactual outcome is never observed. That is why you should tune here, under the natural course, and then lock the SL library choice before the intervention run.

policy_nat <- function(D_block, t, a_names) {
  D_block[, ..a_names, drop = FALSE]
}

wr_nat <- 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 = 3L, v = 3L, seed = 1L,
  id = "id", time = "time",
  policy_spec_fun = policy_nat
)

res_nat <- sdr(
  df = df, weight_object = wr_nat, 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,
  k = 1L, inner_v = 3L, parallel = FALSE, seed = 1L,
  policy_spec_fun = policy_nat
)

branch_cal_summary()

bc <- branch_cal_summary(res_nat)
print(bc)

Focus on:

A slope of 0.7-1.3 and tgt_vl ≈ pred_vl per branch is a reasonable target. When a branch is systematically off, revise the corresponding sl_* library — usually by adding smoother learners (SL.glm, SL.glmnet) if calibration slope is < 1, or richer learners (SL.xgboost, spline learners) if slope is > 1 and the target is poorly matched.

Alternative: pure Q-view with qreg()

sdr() computes both g- and Q-branches, and both are reflected in branch_cal. If you want to isolate the Q-side without the g-branch noise (or without needing weights), qreg() runs the pure Q-recursion under the natural course:

res_qreg <- qreg(
  df = df, 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,
  k = 1L, inner_v = 3L, parallel = FALSE, seed = 1L,
  policy_spec_fun = policy_nat
)
branch_cal_summary(res_qreg)

Tuning the targeting library (itmle() only)

itmle() adds a fluctuation/targeting step on top of the Q-mixture. The targeting model has its own SuperLearner library (sl_tmle / tgt_lib). A natural-course itmle() run lets you tune this library too:

tgt_lib <- c("SL.tmle_empty", "SL.tmle_intercept", "SL.tmle_glm")

res_itmle_nat <- itmle(
  df = df, weight_object = wr_nat, 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 = tgt_lib,
  k = 1L, inner_v = 3L, v_target_itmle = 3L, v_sl_inner_itmle = 3L,
  parallel = FALSE, seed = 1L, policy_spec_fun = policy_nat
)
res_itmle_nat$diagnostics$target_sl

target_sl shows which targeting wrappers were selected per (fold, t). If a wrapper never gets picked, drop it from tgt_lib. If the fit looks unstable, prefer covariate-adaptive wrappers (SL.tmle_glm, SL.tmle_glmnet_*) over intercept-only ones.

3. Final estimation

Only after the weight and calibration diagnostics are satisfactory should you run the estimator under the actual policy of interest. Keep the SL library choice locked from the tuning step.

res <- 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,
  k = 1L, inner_v = 3L, parallel = FALSE, seed = 1L,
  policy_spec_fun = policy_bin, trim = 0.99
)
res

4. Post-estimation sanity checks

The remaining diagnostics describe what the estimator did with the final fit. They cannot easily be used to tune before running the estimator because they depend on the joint behaviour of the weights, the Q-models, and the specific intervention. Treat them as sanity checks rather than a feedback loop.

diagnostics$recursion_diag

One row per fold × time-step, ordered from t = tmax down to 1. The columns are grouped by concept.

Sample counts and fit-status flags

Predicted branch probabilities and means (natural vs shifted, training side)

Differences between _nat and _shf columns reflect how much the intervention shifts the branch predictions. Very small differences on a policy that meaningfully changes treatment suggest the models are insensitive to the treatment covariates — worth investigating.

Full mixture Q on training data

EIF update / targeting magnitudes

Model residuals

Validation-side sanity

diagnostics$target_cal (iTMLE only)

One row per (fold, outer iteration, inner iteration) of the Luedtke targeting loop. Tracks:

Useful for spotting a targeting loop that is not converging (EIF magnitude not shrinking, or coefficients oscillating). If this happens under the actual policy, it usually points to weight-driven instability in the fluctuation step.

diagnostics$sl_summary

Per-fold, per-time, per-component SuperLearner weight table. Consistently zero-weight learners can be pruned from the library on the next run — this is a legitimate feedback loop back into the tuning step.

ic_df

Per-subject influence curve values, used by contrast() to build risk-difference / risk-ratio / odds-ratio contrasts with valid standard errors. Not typically inspected directly; consumed by:

# res_nat and res are the natural-course and shifted fits from above.
contrast(res, res_nat)

Workflow summary

  1. density_ratio()weight_diagnostics() + policy_change_summary → pick trim, revise policy or sl_g if weights or shift look unreasonable.
  2. NAT-run sdr() (or itmle(), or qreg()) → branch_cal_summary() → revise sl_remain / sl_death / sl_recursive / sl_y (and for iTMLE, sl_tmle) until each branch calibrates.
  3. Final sdr() / itmle() under the intervention policy.
  4. Inspect recursion_diag (and target_cal for iTMLE) for sanity; use sl_summary to prune dead learners for the next run; use ic_df via contrast() for downstream comparisons.