## ----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")
}

## ----load-data-gene-----------------------------------------------------------
library(causaldef)
library(ggplot2)

data(gene_perturbation)
str(gene_perturbation)

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

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

## ----plot-gene, fig.height=4--------------------------------------------------
plot(deficiency_gene, type = "bar")

## ----nc-gene------------------------------------------------------------------
set.seed(123)
nc_test <- nc_diagnostic(
    spec_gene,
    method = "iptw",
    alpha = 0.05,
    n_boot = 100
)

print(nc_test)

## ----policy-gene--------------------------------------------------------------
bounds_gene <- policy_regret_bound(
    deficiency_gene,
    utility_range = c(0, 10)
)

print(bounds_gene)

## ----effect-gene--------------------------------------------------------------
effect_gene <- estimate_effect(
    deficiency_gene,
    target_method = "iptw"
)

print(effect_gene)

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

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

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

## ----estim-hct----------------------------------------------------------------
deficiency_hct <- estimate_deficiency(
    spec_hct,
    methods = c("unadjusted", "iptw"),
    n_boot = 50 # Use more for production
)

print(deficiency_hct)

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

## ----policy-hct---------------------------------------------------------------
# Utility = months of survival (horizon = 24)
bounds_hct <- policy_regret_bound(
    deficiency_hct,
    utility_range = c(0, 24)
)

print(bounds_hct)

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

## ----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.")

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