--- title: "Complete Workflow: From Data to Decision" author: "Deniz Akdemir" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Complete Workflow: From Data to Decision} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} eval_surv_effect <- getRversion() >= "4.0.0" && requireNamespace("survival", quietly = TRUE) knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5 ) if (!exists("deparse1", envir = baseenv())) { deparse1 <- function(expr, collapse = " ", width.cutoff = 500L, ...) { paste(deparse(expr, width.cutoff, ...), collapse = collapse) } } fmt_num <- function(x, digits = 3) { formatC(x, digits = digits, format = "f") } ``` ## Overview This vignette demonstrates the complete **causaldef** workflow, from data specification through to policy decision-making. We show how deficiency theory translates abstract statistical concepts into actionable clinical insights. The workflow consists of four stages: 1. **Specify** → Define the causal problem 2. **Estimate** → Compute deficiency for different adjustment strategies 3. **Diagnose** → Validate assumptions using negative controls and sensitivity analysis 4. **Decide** → Compute policy regret bounds and make informed decisions --- ## Part 1: Gene Perturbation Study (Continuous Outcome) We begin with the `gene_perturbation` dataset, which simulates a CRISPR knockout experiment. This illustrates the core workflow for continuous outcomes. ### 1.1 Data Description ```{r load-data-gene} library(causaldef) library(ggplot2) data(gene_perturbation) str(gene_perturbation) ``` **Variables:** - `knockout_status`: Treatment (Control vs. Knockout) - `target_expression`: Primary outcome (gene expression level) - `housekeeping_gene`: Negative control outcome (shouldn't be affected by knockout) - `batch`, `library_size`: Technical confounders **Causal Structure:** ``` [Batch, Library Size] | v [Knockout] -----> [Target Expression] \ \---X--> [Housekeeping Gene] (no causal effect) ``` The housekeeping gene is affected by the same technical variations but NOT by the knockout, making it an ideal negative control. ### 1.2 Step 1: Specification ```{r spec-gene} spec_gene <- causal_spec( data = gene_perturbation, treatment = "knockout_status", outcome = "target_expression", covariates = c("batch", "library_size"), negative_control = "housekeeping_gene", estimand = "ATE", outcome_type = "continuous" ) print(spec_gene) ``` ### 1.3 Step 2: Deficiency Estimation We compare three adjustment strategies: 1. **Unadjusted**: Ignores technical confounders 2. **IPTW**: Reweights samples to balance batch and library size 3. **AIPW**: Augmented IPTW (doubly robust) ```{r estim-gene} deficiency_gene <- estimate_deficiency( spec_gene, methods = c("unadjusted", "iptw", "aipw"), n_boot = 100 # Use more for production (e.g., 1000) ) print(deficiency_gene) ``` **Interpretation:** In this run, the unadjusted PS-TV proxy is about `r fmt_num(deficiency_gene$estimates["unadjusted"])`, while IPTW and AIPW reduce it to about `r fmt_num(deficiency_gene$estimates["iptw"])`. That is a substantial reduction in the observational-to-interventional gap, but it should still be interpreted through the downstream regret bounds rather than as literal proof of randomized equivalence. ```{r plot-gene, fig.height=4} plot(deficiency_gene, type = "bar") ``` ### 1.4 Step 3: Diagnose with Negative Control The negative control diagnostic tests whether our adjustment removes ALL confounding, not just the measured confounders. ```{r nc-gene} set.seed(123) nc_test <- nc_diagnostic( spec_gene, method = "iptw", alpha = 0.05, n_boot = 100 ) print(nc_test) ``` **Decision Logic:** | Result | Interpretation | Action | |--------|----------------|--------| | `falsified = FALSE` | The screen does not find strong evidence against the current adjustment, but this is not proof that all confounding is gone | Proceed with explicit uncertainty qualification | | `falsified = TRUE` | Residual confounding is detected by the screen | Add covariates, reconsider the design, or report the limitation | ### 1.5 Step 4: Policy Decision Suppose we're deciding whether to pursue this gene target for drug development. The utility is measured on a scale where: - 0 = no promise (no effect on expression) - 10 = maximum promise (strong effect) ```{r policy-gene} bounds_gene <- policy_regret_bound( deficiency_gene, utility_range = c(0, 10) ) print(bounds_gene) ``` **Regret Bounds:** `policy_regret_bound()` reports: - **Transfer penalty** \(M\cdot\delta\): additive worst-case regret inflation term, and - **Minimax safety floor** \((M/2)\cdot\delta\): irreducible worst-case regret when \(\delta>0\). **Decision Rule:** - If `transfer_penalty` is small relative to the practical stakes of the decision → the observational evidence may be usable with caveats - If `transfer_penalty` is large relative to those stakes → seek more evidence or narrow the decision scope ### 1.6 Effect Estimation Finally, we estimate the causal effect using the best-performing method: ```{r effect-gene} effect_gene <- estimate_effect( deficiency_gene, target_method = "iptw" ) print(effect_gene) ``` **Complete Report:** ```{r report-gene, results='asis'} conclusion_gene <- if (nc_test$falsified) { "Residual confounding was detected by the negative control screen. Treat the effect estimate as exploratory." } else if (bounds_gene$transfer_penalty < 0.5) { "Adjustment reduced the proxy gap substantially and the negative control did not falsify the analysis, but residual uncertainty remains." } else { "Adjustment improved balance, but the remaining proxy gap is still material on the chosen utility scale." } cat(sprintf( " ## Gene Perturbation Analysis Report **Treatment Effect (IPTW-adjusted):** %.2f log2 expression units **Deficiency (δ):** %.3f **Negative Control Screen:** %s (p = %.3f; delta_NC = %.3f) **Transfer Penalty:** %.3f on [0, 10] scale **Minimax Safety Floor:** %.3f on [0, 10] scale **Conclusion:** %s ", effect_gene$estimate, deficiency_gene$estimates["iptw"], ifelse(nc_test$falsified, "Falsified", "Not falsified"), nc_test$p_value, nc_test$delta_nc, bounds_gene$transfer_penalty, bounds_gene$minimax_floor, conclusion_gene )) ``` --- ## Part 2: Hematopoietic Cell Transplantation (Survival Outcome) Next, we analyze the `hct_outcomes` dataset, which mimics a retrospective registry study comparing conditioning regimens in HCT. ### 2.1 Data Description ```{r load-data-hct} data(hct_outcomes) str(hct_outcomes) # Summarize key variables summary(hct_outcomes[, c("age", "kps", "time_to_event")]) table(hct_outcomes$conditioning_intensity, hct_outcomes$event_status) ``` **Clinical Context:** - **Myeloablative conditioning**: High-intensity chemotherapy (younger, healthier patients) - **Reduced-intensity conditioning**: Lower dose (older, sicker patients) - **Outcome in this vignette**: 24-month restricted mean survival using death as the binary event The key challenge is **confounding by indication**: doctors assign treatment based on patient status, making naive comparisons biased. ### 2.2 Step 1: Survival Specification The deficiency, sensitivity, and regret-bound calculations in this section run on all supported runtimes. The final RMST effect-estimation chunk requires a compatible `survival` runtime, which in the current support matrix means `R >= 4.0`. ```{r spec-hct} # Create binary event indicator for death hct_outcomes$event_death <- as.integer(as.character(hct_outcomes$event_status) == "Death") spec_hct <- causal_spec_survival( data = hct_outcomes, treatment = "conditioning_intensity", time = "time_to_event", event = "event_death", covariates = c("age", "disease_status", "kps", "donor_type"), estimand = "RMST", horizon = 24 # 24-month restricted mean survival time ) print(spec_hct) ``` For a full competing-risks analysis of relapse versus death, prefer `causal_spec_competing()` together with `estimate_deficiency_competing()`. Here we keep the workflow aligned with the current reproducible death-endpoint survival interface. ### 2.3 Step 2: Deficiency Estimation ```{r estim-hct} deficiency_hct <- estimate_deficiency( spec_hct, methods = c("unadjusted", "iptw"), n_boot = 50 # Use more for production ) print(deficiency_hct) ``` **Clinical Interpretation:** The deficiency tells us how much our observational evidence differs from what an RCT would provide. In this run, the unadjusted proxy is about `r fmt_num(deficiency_hct$estimates["unadjusted"])`, while IPTW reduces it to about `r fmt_num(deficiency_hct$estimates["iptw"])`. That is a meaningful improvement, but the remaining gap is still large enough to matter on a 24-month clinical utility scale. ### 2.4 Step 3: Confounding Frontier Beyond point estimates, we can map a *sensitivity analysis* showing how deficiency varies with hypothetical unmeasured confounding: ```{r frontier-hct, fig.height=5} frontier <- confounding_frontier( spec_hct, alpha_range = c(-2, 2), # Confounding path: U → Treatment gamma_range = c(-2, 2), # Confounding path: U → Outcome grid_size = 30 ) print(frontier) plot(frontier) ``` **Reading the Frontier Map:** - **Center** (α = 0 or γ = 0): No unmeasured confounding → δ = 0 - **Corners**: Strong confounding on both paths → high δ - **Observed covariates** (dots): Benchmark strengths of measured confounders If an unmeasured confounder would need extreme strength (beyond observed benchmarks) to substantially increase δ, conclusions are robust. ### 2.5 Step 4: Policy Regret and RMST Effect ```{r policy-hct} # Utility = months of survival (horizon = 24) bounds_hct <- policy_regret_bound( deficiency_hct, utility_range = c(0, 24) ) print(bounds_hct) ``` **Clinical Regret Bounds:** In the current example, the IPTW proxy implies a transfer penalty of about `r fmt_num(bounds_hct$transfer_penalty, 2)` months and a minimax safety floor of about `r fmt_num(bounds_hct$minimax_floor, 2)` months on the 0--24 month utility scale. That means the residual observational-to-interventional gap is materially improved relative to the unadjusted analysis, but still not negligible for treatment decisions. ```{r effect-hct, eval = eval_surv_effect} # Estimate RMST difference effect_hct <- estimate_effect( deficiency_hct, target_method = "iptw", contrast = c("Myeloablative", "Reduced") ) print(effect_hct) ``` ```{r effect-hct-note, results='asis', eval = !eval_surv_effect} cat("Effect-estimation chunks are skipped on runtimes without the required survival support. The deficiency, frontier, and regret-bound calculations above still provide the main observational-versus-interventional diagnostics for this example.") ``` ### 2.6 Complete Decision Framework ```{r decision-hct, results='asis', eval = eval_surv_effect} delta_iptw <- deficiency_hct$estimates["iptw"] transfer_penalty <- bounds_hct$transfer_penalty minimax_floor <- bounds_hct$minimax_floor rmst_diff <- effect_hct$estimate # Decision logic if (transfer_penalty < 1) { evidence_quality <- "LOW residual decision risk" } else if (transfer_penalty < 2) { evidence_quality <- "MODERATE residual decision risk" } else { evidence_quality <- "HIGH residual decision risk" } # Benefit-to-risk ratio if (!is.na(rmst_diff) && !is.na(transfer_penalty) && transfer_penalty > 0) { benefit_to_risk <- abs(rmst_diff) / transfer_penalty recommendation <- ifelse(benefit_to_risk > 2, "The estimated effect exceeds the transfer penalty, but the recommendation should still be qualified by residual confounding risk.", "The estimated effect is not clearly separated from the transfer penalty; avoid strong recommendations from observational evidence alone." ) } else { benefit_to_risk <- NA recommendation <- "Unable to calculate benefit-to-risk ratio" } cat(sprintf( " ## HCT Treatment Decision Report **RMST Difference (IPTW):** %.2f months (%s favored) **Deficiency:** %.3f **Residual Decision Risk:** %s **Transfer Penalty:** %.2f months **Minimax Safety Floor:** %.2f months **Benefit-to-Risk Ratio:** %.1f:1 **Recommendation:** %s ### Clinical Translation The observational evidence suggests %s conditioning provides approximately %.1f months difference in restricted mean survival within the first 24 months. However, the transfer penalty is %.1f months and the minimax safety floor is %.1f months on the 0--24 month utility scale. This is a decision-relevant amount of residual uncertainty, so clinicians should weigh it against individual patient factors and avoid over-interpreting the observational comparison. ", abs(rmst_diff), ifelse(rmst_diff > 0, "Myeloablative", "Reduced"), delta_iptw, evidence_quality, transfer_penalty, minimax_floor, benefit_to_risk, recommendation, ifelse(rmst_diff > 0, "myeloablative", "reduced-intensity"), abs(rmst_diff), transfer_penalty, minimax_floor )) ``` --- ## Part 3: Comparative Analysis Across Studies ### 3.1 When Is Observational Evidence Sufficient? | Study | IPTW proxy | Additional diagnostic | Transfer penalty | Interpretation | |-------|------------|---------------------|------------------|----------------| | Gene Perturbation | `r fmt_num(deficiency_gene$estimates["iptw"])` | NC not falsified; `r fmt_num(nc_test$delta_nc)` observable residual-association proxy | `r fmt_num(bounds_gene$transfer_penalty, 2)` on [0, 10] | Residual uncertainty appears modest but nonzero | | HCT Death-Endpoint RMST | `r fmt_num(deficiency_hct$estimates["iptw"])` | Confounding frontier recommended; no negative control in this example | `r fmt_num(bounds_hct$transfer_penalty, 2)` months | Residual uncertainty remains material on the 24-month scale | ### 3.2 General Workflow Summary ``` ┌─────────────────────────────────────────────────────────────────┐ │ SPECIFY: causal_spec() / causal_spec_survival() │ │ ↓ Define treatment, outcome, covariates, NC │ ├─────────────────────────────────────────────────────────────────┤ │ ESTIMATE: estimate_deficiency() │ │ ↓ Compare unadjusted, IPTW, AIPW, TMLE, etc. │ │ ↓ Select method with lowest δ │ ├─────────────────────────────────────────────────────────────────┤ │ DIAGNOSE: nc_diagnostic() + confounding_frontier() │ │ ↓ Test whether assumptions are falsified │ │ ↓ Map sensitivity to unmeasured confounding │ ├─────────────────────────────────────────────────────────────────┤ │ DECIDE: policy_regret_bound() + estimate_effect() │ │ ↓ Compute transfer penalty / minimax floor │ │ ↓ Report effect with uncertainty qualification │ └─────────────────────────────────────────────────────────────────┘ ``` --- ## References 1. Akdemir, D. (2026). Constraints on Causal Inference as Experiment Comparison: A Framework for Identification, Transportability, and Policy Learning. DOI: 10.5281/zenodo.21877511 2. Le Cam, L., & Yang, G. L. (2000). Asymptotics in Statistics: Some Basic Concepts. Springer. 3. VanderWeele, T. J., & Ding, P. (2017). Sensitivity Analysis in Observational Research: Introducing the E-value. Annals of Internal Medicine.