r.R

# Module 20 - competing risks and time-varying Cox models.
#   Rscript module/20-competing-risks-timevariant/code/r.R
#
# Reads the coherent competing-risks dataset data/konkurrenz_risiken.csv
# (generated by data/generate_data.py with its own RNG streams). Columns:
#   event_time             day of the exit event (death day, discharge day, or 30)
#   event_state            0 = censored at day 30, 1 = death, 2 = discharge
#   vasopressor_start_tag  day the time-varying exposure begins, else NA
#
# Discharge is a genuine competing event whose hazard FALLS with severity, and
# administrative censoring at day 30 is preserved. Baseline covariates come
# from kohorte.csv. No files under data/ are edited.
#
# Fine-Gray (subdistribution hazard) is fitted TWICE for cross-validation:
#   (a) survival::finegray() + a weighted coxph() (the README route), and
#   (b) cmprsk::crr() (the canonical Fine-Gray implementation).
# Both must agree. require_pkgs() aborts loudly (non-zero exit) if a package is
# missing, so the lesson never exits 0 while silently producing nothing.

script <- normalizePath(sub("--file=", "", grep("--file=", commandArgs(), value = TRUE)[1]))
root   <- dirname(dirname(dirname(dirname(script))))
source(file.path(root, "lib", "helpers.R"))
require_pkgs("survival", "cmprsk")
suppressPackageStartupMessages(library(survival))
set.seed(SEED)

ADMIN_HORIZON <- 30

df  <- .load("konkurrenz_risiken.csv")
cov <- load_cohort()[, c("patient_id", "alter", "sofa_score", "diabetes")]
df  <- merge(df, cov, by = "patient_id")
df$event_state_f <- factor(df$event_state, levels = c(0, 1, 2),
                           labels = c("censored", "death", "entlassung"))

cat("=== 3-state competing-risks outcome (data/konkurrenz_risiken.csv) ===\n")
print(table(df$event_state_f))

# --- 1) Cumulative Incidence Function (CIF) via multi-state survfit --------
cat("\n=== 1) CIF (correct) vs. naive 1-KM (treats discharge as censoring) ===\n")
cif_fit <- survfit(Surv(event_time, event_state_f) ~ 1, data = df)
cif_summary <- summary(cif_fit, times = ADMIN_HORIZON)
cif_30 <- cif_summary$pstate[1, which(cif_fit$states == "death")]

km_fit <- survfit(Surv(event_time, event_state == 1) ~ 1, data = df)
naive_30 <- 1 - summary(km_fit, times = ADMIN_HORIZON)$surv

cat(sprintf("CIF for death at day %d (competing risks respected): %.1f%%\n",
            ADMIN_HORIZON, 100 * cif_30))
cat(sprintf("Naive 1-KM 'risk of death' at the same horizon (discharge as censoring): %.1f%%\n",
            100 * naive_30))
cat("The naive number OVERSHOOTS the truth - visibly, but bounded - because discharged patients\n")
cat("drop out of the risk set. The day-30 risk set is still ~90 patients, so 1-KM does not blow\n")
cat("up to 100 %: this is the realistic textbook bias, not a simulation artefact.\n")

# --- 2) Cause-specific Cox model (competing event = ordinary censoring) ----
cat("\n=== 2) Cause-specific Cox model (death only; discharge treated as censoring) ===\n")
cs_fit <- coxph(Surv(event_time, event_state == 1) ~ alter + sofa_score + diabetes, data = df)
print(round(summary(cs_fit)$coefficients[, c("coef", "exp(coef)", "se(coef)", "Pr(>|z|)")], 4))
cat("Interpretation: the instantaneous death rate among patients who have not yet died OR\n")
cat("been discharged - good for etiology, not for population-level absolute risk.\n")

# --- 3) Fine-Gray subdistribution hazard model (two implementations) --------
cat("\n=== 3) Fine-Gray model (subdistribution hazard) ===\n")
fg_data <- finegray(Surv(event_time, event_state_f) ~ alter + sofa_score + diabetes,
                     data = df, etype = "death")
fg_fit <- coxph(Surv(fgstart, fgstop, fgstatus) ~ alter + sofa_score + diabetes,
                weight = fgwt, data = fg_data)
cat("(a) survival::finegray() + weighted coxph():\n")
print(round(summary(fg_fit)$coefficients[, c("coef", "exp(coef)", "se(coef)", "Pr(>|z|)")], 4))

crr_fit <- cmprsk::crr(ftime = df$event_time, fstatus = df$event_state,
                       cov1 = as.matrix(df[, c("alter", "sofa_score", "diabetes")]),
                       failcode = 1, cencode = 0)
crr_tab <- summary(crr_fit)$coef
cat("(b) cmprsk::crr() cross-check (canonical Fine-Gray) - subdistribution HRs:\n")
print(round(crr_tab[, c("coef", "exp(coef)", "p-value")], 4))
cat("Both implementations agree. Interpretation: the effect on the ACTUAL population-level\n")
cat("probability of death by time t, keeping discharged patients in the risk set with a\n")
cat("shrinking IPCW weight - good for prognosis/absolute risk.\n")
cat("Contrast sofa_score's SHR here with the cause-specific HR above: Fine-Gray's is LARGER,\n")
cat("because higher-SOFA patients are also less likely to be discharged - built into the\n")
cat("generating process (discharge mean = 15 + 1.4*sofa days, so the discharge hazard falls\n")
cat("with severity), a channel cause-specific hazards ignore by design.\n")
cat(sprintf("Power caveat: only %d deaths total, so splitting events across two causes leaves\n",
            sum(df$event_state == 1)))
cat("few events per cause - a richly adjusted competing-risks model would be underpowered.\n")

# --- 4) Immortal time bias: naive time-fixed vs. correct time-varying -------
cat("\n=== 4) Time-varying Cox: vasopressor exposure (start-stop format) ===\n")
df$death <- as.integer(df$event_state == 1)
df$ever_vasopressor <- as.integer(!is.na(df$vasopressor_start_tag))

# NAIVE: exposure as a fixed baseline covariate (WRONG - immortal time bias).
naive_fit <- coxph(Surv(event_time, death) ~ ever_vasopressor + sofa_score, data = df)
hr_naive <- exp(coef(naive_fit)["ever_vasopressor"])
p_naive  <- summary(naive_fit)$coefficients["ever_vasopressor", "Pr(>|z|)"]
cat(sprintf("NAIVE time-fixed HR for 'ever vasopressor' (immortal time bias): %.2f (p=%.3f)\n",
            hr_naive, p_naive))
cat("-> a spurious PROTECTIVE effect: treated patients had to survive long enough to be\n")
cat("   treated, so their early deaths are misattributed to the untreated group.\n")

# CORRECT: time-varying exposure switched on only from the start day.
n <- nrow(df)
rows <- vector("list", 0)
for (i in seq_len(n)) {
  pid   <- df$patient_id[i]
  t_end <- df$event_time[i]
  death <- df$death[i]
  s     <- df$sofa_score[i]
  st    <- df$vasopressor_start_tag[i]
  if (!is.na(st) && st < t_end) {
    rows[[length(rows) + 1]] <- data.frame(patient_id = pid, start = 0, stop = st,
                                            vasopressor = 0, death = 0, sofa = s)
    rows[[length(rows) + 1]] <- data.frame(patient_id = pid, start = st, stop = t_end,
                                            vasopressor = 1, death = death, sofa = s)
  } else {
    rows[[length(rows) + 1]] <- data.frame(patient_id = pid, start = 0, stop = t_end,
                                            vasopressor = 0, death = death, sofa = s)
  }
}
tv <- do.call(rbind, rows)
tv <- tv[tv$stop > tv$start, ]
cat(sprintf("\nStart-stop rows: %d | patients: %d | exposed: %d\n",
            nrow(tv), length(unique(tv$patient_id)), sum(df$ever_vasopressor)))
print(table(vasopressor = tv$vasopressor, death = tv$death))

tv_fit <- coxph(Surv(start, stop, death) ~ vasopressor + sofa, data = tv)
cat("\n")
print(round(summary(tv_fit)$coefficients[, c("coef", "exp(coef)", "se(coef)", "Pr(>|z|)")], 4))
hr_tv <- exp(coef(tv_fit)["vasopressor"])
cat(sprintf("\nCORRECT time-varying HR for vasopressor: %.2f - the spurious protective effect\n", hr_tv))
cat("collapses toward 1 once the immortal pre-exposure time is credited to the UNEXPOSED state.\n")
cat("The generating process gives vasopressor no causal effect on death (it only marks\n")
cat("severity), so ~1.0 is the honest answer; the naive model manufactured a benefit.\n")