# Kaplan-Meier con survival: mismo CSV y ajuste que Survival Lab.
# Descarga synthetic-data.csv en la misma carpeta que este script.
# event: 1 = event, 0 = right censoring; time is in months.
library(survival)
d <- read.csv("synthetic-data.csv")
stopifnot(all(c("id", "group", "time", "event", "censor_reason") %in% names(d)),
          nrow(d) == 300L, !anyNA(d), !anyDuplicated(d$id),
          setequal(unique(d$group), c("A", "B")),
          all(table(d$group) == 150L), is.numeric(d$time),
          all(is.finite(d$time) & d$time >= 0 & d$time <= 24),
          all(d$event %in% c(0L, 1L)),
          all(d$censor_reason %in% c("event", "dropout", "administrative")),
          all((d$event == 1L) == (d$censor_reason == "event")),
          all(d$time[d$censor_reason == "administrative"] == 24),
          all(d$time[d$censor_reason == "dropout"] < 24))
d$group <- factor(d$group, levels = c("A", "B"))

fit <- survfit(Surv(time, event) ~ group, data = d,
               conf.type = "log-log", conf.int = 0.95)
at12 <- summary(fit, times = 12, extend = FALSE)
results <- data.frame(group = sub("^group=", "", as.character(at12$strata)),
                      S12_pct = 100 * at12$surv,
                      lower95_pct = 100 * at12$lower,
                      upper95_pct = 100 * at12$upper,
                      n_risk = at12$n.risk)
cat("\nResultados a 12 meses (%)\n")
print(results, digits = 4, row.names = FALSE)
risk <- summary(fit, times = c(0, 6, 12, 18, 24), extend = FALSE)
cat("\nPersonas en riesgo justo antes de cada mes\n")
risk_table <- matrix(risk$n.risk, nrow = 2, byrow = TRUE,
                     dimnames = list(c("A", "B"), c("0 m", "6 m", "12 m", "18 m", "24 m")))
print(risk_table)

plot(fit, col = c("#147D82", "#163B53"), lwd = 2,
     conf.int = TRUE, mark.time = TRUE,
     xlab = "Tiempo (meses)", ylab = "Probabilidad de permanecer sin evento",
     main = "Kaplan-Meier e IC puntuales del 95%", xlim = c(0, 24), ylim = c(0, 1))
legend("bottomleft", legend = c("Grupo A", "Grupo B"),
       col = c("#147D82", "#163B53"), lwd = 2, bty = "n")
cat("\n", R.version.string, "; survival ",
    as.character(packageVersion("survival")), "\n", sep = "")