Reproducible research compendium and companion R package steCI for the methods
paper Cluster on the Subject, Not the Record: Confidence Intervals and Simultaneous
Bands for Additive-Hazards Sequential Trial Emulation (M. Ehsan Karim, 2026).
It derives, and evaluates the coverage of, confidence-interval procedures for the weighted-Aalen additive-hazards sequential-trial-emulation (STE) estimator of a causal marginal risk difference, answering — for the additive-hazards outcome model — the open problem stated by Limozin, Seaman & Su (2025): "inference procedures for STE based on additive hazard models warrant further research."
Two estimands. E1 — a constant hazard difference; E2 — the marginal risk-difference (MRD) curve over follow-up.
| You want… | Go to |
|---|---|
| Install and use the estimator + CI methods | steCI/ · Install · Quick start |
| The figures/tables and the code + data behind each | Where each result comes from |
| To re-run the whole simulation study | Reproduce the study |
| The applied illustration | applied/ · Reproduce the applied example |
| The exact-truth / sensitivity / robustness analyses | dev/ · the results map |
- Keogh et al. (2023, Stat. Med.) give the additive-hazards STE estimator of the marginal risk difference and recommend the nonparametric bootstrap for confidence intervals — without ever evaluating its coverage.
- Limozin et al. (2025, SMMR) compare sandwich / bootstrap / LEF-bootstrap / jackknife intervals, but only for a pooled-logistic marginal structural model, explicitly leaving the additive-hazards case to future work.
- Because each subject is reused across several stacked emulated trials, the default
row-level robust ("sandwich") standard error reported by
timereg::aalenis asymptotically anticonservative: it omits a positive within-subject cross-trial covariance, so its under-coverage does not vanish as the sample grows. This package derives the estimator's influence functions, proves the anticonservatism, and shows that clustering the robust variance on the subject (for the constant hazard difference) and a fast multiplier / LEF bootstrap with simultaneous bands (for the risk-difference curve) recover most of the deficit — leaving a small, quantified residual shortfall rather than fully restoring nominal coverage.
| Method | Estimand | Basis | Status |
|---|---|---|---|
| M1 model-based SE | E1 | Aalen model variance | anticonservative — avoid |
| M2 robust (row) SE | E1 | row-level sandwich (timereg default) |
inconsistent, anticonservative — avoid |
| M3 cluster-robust SE | E1 | sandwich clustered on subject | recommended (E1) — recovers most of the deficit |
| M4 nonparametric bootstrap | E1, E2 (MRD) | refit trials + weights + fit per resample | valid — reference standard |
| M5 jackknife | E1 | delete-one-subject | valid — optional |
| M6 multiplier / LEF | E2 (MRD) | perturb pre-computed influence functions; pointwise + sup-t band | recommended (E2) — recovers most of the deficit |
# install.packages("remotes")
remotes::install_github("ehsanx/steCI", subdir = "steCI")The installable package lives in the steCI/ subfolder; the rest of the
repository reproduces the simulation study and the applied example.
library(steCI)
set.seed(1)
## simulate from the additive-hazard MSM (Keogh et al. 2021)
d <- dat_sim(n = 1000, alpha.0 = 0.25, gamma.0 = -1, gamma.L = 1.5)
## point estimates: constant hazard difference (E1) and the marginal
## risk-difference curve (E2), by empirical standardisation of a weighted Aalen fit
est <- ste_estimate(d, estimand = "both")
## confidence intervals
ci_analytic(ste_prepare(d)) # M1 model / M2 row-robust / M3 cluster-robust analytic SEs
ci_lef(d, approach = "a1") # M6 multiplier/LEF: pointwise intervals + simultaneous band (E2)
ci_bootstrap(d, B = 200) # M4 nonparametric bootstrap (reference)
ci_jackknife(d) # M5 delete-one-subject jackknifeci_lef(..., approach = "a2") additionally propagates weight-estimation uncertainty
through the two logistic weight models (the effect is small with stabilised weights).
steCI/ the R package (estimator + all CI methods + dat_sim + tests + vignette)
R/ man/ tests/ vignettes/ DESCRIPTION NAMESPACE LICENSE
simulation/ the ADEMP simulation study (81-cell factorial design)
scenarios.R the factorial grid (n x event rate x confounding x prevalence)
run_sim.R one_replicate.R run one replicate / drive the design
aggregate.R metrics.R aggregate chunks -> ADEMP metrics (bias, coverage, SE ratio, ...)
figures.R regenerate every manuscript figure (writes to figures/)
report.R calibrate_dgm.R smoke_test.R
slurm/ UBC ARC (Sockeye) Slurm array-job scripts + submit_*.sh launchers
data/ the SMALL aggregated results the tables/figures are computed from
exact/results/ AUTHORITATIVE aggregation, scored against the exact closed-form truths
results/ companion aggregation (truth-free quantities such as SE ratios)
dev/ analysis, validation, sensitivity & robustness scripts backing specific
manuscript claims, plus their small result caches (*.rds) and the exact
closed-form truth cache (truth_cache_exact.rds). See the results map below.
applied/ Stanford heart-transplant illustration (public survival::jasa data)
figures/ the generated figures (fig1_e1_coverage.pdf, ...); regenerate via figures.R
The large raw simulation output (per-replicate objects, ~180 MB) and the local
simulation/output/ working directory are not tracked; they are regenerated by the
pipeline. The small aggregated metrics needed to reproduce every figure and table are
included, under data/.
Both are the same replicates, re-aggregated two ways:
data/exact/results/— scored against the exact closed-form marginal truths. This is authoritative for anything that depends on a truth (bias, MSE, coverage).simulation/figures.Rreadsdata/exact/results/metrics.rdsby default.data/results/— the companion aggregation; use it for truth-free quantities (e.g. empirical-vs-model SE ratios), which are identical across the two.
Each row: the reported object → the script that produces it → the input it reads.
(Paths are relative to the repo root; every script honours REPO_ROOT / RESULTS_DIR
environment overrides and defaults to the shipped locations.)
| Result | Script | Reads |
|---|---|---|
| Figures (E1 coverage, E2 coverage, MRD band, …) | simulation/figures.R |
data/exact/results/metrics.rds → writes figures/ |
| Main ADEMP metrics (bias, coverage, SE ratio by method × cell) | simulation/aggregate.R + simulation/metrics.R |
regenerates data/… from per-replicate chunks |
| Exact closed-form truths (E1 hazard diff., E2 MRD curve) | dev/rebuild_truth_multiseed.R, dev/apply_exact_truth.R |
→ dev/truth_cache_exact.rds |
max.clust sensitivity (M3 cluster SE under exact vs binned clustering) |
dev/maxclust_arc.R (runner) → dev/maxclust_collect.R (scoring) |
→ dev/maxclust_results.rds |
| Population-limit / finite-n bias of the Aalen fit | dev/plim_arc.R, dev/plim_closed_form.R → dev/plim_collect.R |
data/exact/results/, → dev/plim_results.rds |
| Jensen decomposition of the coverage shortfall | dev/jensen_decomposition.R |
data/results/, → dev/jensen_decomposition_results.rds |
| Weight-truncation sensitivity | dev/truncation_*.R |
→ dev/truncation_*.rds |
| Negative-control diagnostic | dev/negative_control.R |
→ dev/negative_control_results.rds |
| Robustness (misspecification, time-varying effect, 2nd confounder) | dev/robustness_*.R, dev/robust_arc.R → dev/robust_collect.R |
→ dev/robustness_*.rds |
| Variance-formula / M6 / Approach-2 validation | dev/validate_*.R, dev/verify_*.R |
(self-contained checks) |
| Applied Stanford heart-transplant analysis | applied/heart_analysis.R |
survival::jasa (public) → applied/heart_results.rds |
dev/ scripts are named by the claim they back (e.g. dev/failure_rates_and_nominal.R,
dev/negative_control.R); their *.rds caches hold the computed numbers.
# 1. Exact closed-form truths are cached in dev/truth_cache_exact.rds (already included;
# rebuild with dev/rebuild_truth_multiseed.R).
# 2. Run the 81-cell x 1000-replicate design (on a cluster; see simulation/slurm/).
# Locally, a quick smoke test of the pipeline end-to-end:
source("simulation/smoke_test.R")
# 3. Aggregate the per-replicate chunks, then regenerate the figures:
source("simulation/aggregate.R") # -> data/ (locally; the shipped data/ is the full run)
source("simulation/figures.R") # -> figures/On UBC ARC the full study runs as Slurm array jobs — see
simulation/slurm/submit_all.sh. The Slurm headers use
placeholder --account / --mail-user values; set your own before submitting.
source("applied/heart_analysis.R") # Stanford heart transplant; all six CI methodsPorted from Keogh et al. (2021, Biometrical Journal): an unmeasured frailty
U ~ N(0, 0.1^2); a time-varying confounder L_k ~ N(0.8 L_{k-1} - A_{k-1} + 0.1(k-1) + U, 1);
monotone treatment logit Pr(A_k = 1 | A_{k-1}=0, L_k) = gamma_0 + gamma_L L_k; and a
piecewise-constant additive hazard h = alpha_0 + alpha_A A + alpha_L L + alpha_U U
with alpha_A = -0.04, alpha_L = alpha_U = 0.015. The 81-cell grid varies n,
alpha_0 (event rate), gamma_L (confounding), and gamma_0 (prevalence).
Karim, M. E. Cluster on the Subject, Not the Record: Confidence Intervals and Simultaneous Bands for Additive-Hazards Sequential Trial Emulation. (Preprint, 2026.)
- Keogh RH, Gran JM, Seaman SR, Davies G, Vansteelandt S (2023). Causal inference in survival analysis using longitudinal observational data: Sequential trials and marginal structural models. Statistics in Medicine 42(13):2191–2225.
- Limozin JM, Seaman SR, Su L (2025). Inference procedures in sequential trial emulation with survival outcomes. Statistical Methods in Medical Research.
- Keogh RH, Seaman SR, Gran JM, Vansteelandt S (2021). Simulating longitudinal data from marginal structural models using the additive hazard model. Biometrical Journal 63(7):1526–1541.
- Lin DY, Ying Z (1994). Semiparametric analysis of the additive risk model. Biometrika 81(1):61–71.
MIT © 2026 M. Ehsan Karim. See LICENSE.