# Simulations for /guides/log-rank-vs-cox-regression/. R 4.3.3, survival.
# Seed 20261110, fixed before running (not chosen from results).
# Part 1: false-positive rate with no true difference, small groups. Exponential event times with
#   median 35 days, follow-up ends at day 60 (administrative censoring, about 30% censored).
# Part 2: crossing curves. Treated animals have a higher hazard before day 20 and a lower one after;
#   proportion of experiments with p < 0.05 for the log-rank test and an RMST(60) difference test.
library(survival)
set.seed(20261110)
rate <- log(2) / 35
tests_null <- function(n, reps) {
  hits <- matrix(0, reps, 4, dimnames = list(NULL, c("logrank", "cox_wald", "cox_lr", "cox_score")))
  for (i in seq_len(reps)) {
    t <- rexp(2 * n, rate); ev <- as.integer(t <= 60); t <- pmin(t, 60)
    g <- rep(0:1, each = n)
    if (sum(ev) < 2 || length(unique(g[ev == 1])) < 2) { hits[i, ] <- NA; next }
    s <- suppressWarnings(summary(coxph(Surv(t, ev) ~ g)))
    lr <- survdiff(Surv(t, ev) ~ g)
    hits[i, ] <- c(pchisq(lr$chisq, 1, lower.tail = FALSE), s$waldtest["pvalue"], s$logtest["pvalue"], s$sctest["pvalue"]) < 0.05
  }
  colMeans(hits, na.rm = TRUE)
}
rmst_p <- function(t, ev, g, tau = 60) {
  tab <- summary(survfit(Surv(t, ev) ~ g), rmean = tau)$table
  d <- tab[2, "rmean"] - tab[1, "rmean"]; se <- sqrt(sum(tab[, "se(rmean)"]^2))
  2 * pnorm(-abs(d / se))
}
cross_power <- function(n, reps) {
  hits <- matrix(0, reps, 2, dimnames = list(NULL, c("logrank", "rmst")))
  for (i in seq_len(reps)) {
    t0 <- rexp(n, rate)
    # treated: hazard 3x control before day 20, 0.25x control after day 20 (piecewise exponential)
    e1 <- rexp(n, 3 * rate); late <- e1 > 20
    e1[late] <- 20 + rexp(sum(late), 0.25 * rate)
    t <- c(t0, e1); ev <- as.integer(t <= 60); t <- pmin(t, 60); g <- rep(0:1, each = n)
    lr <- survdiff(Surv(t, ev) ~ g)
    hits[i, ] <- c(pchisq(lr$chisq, 1, lower.tail = FALSE), rmst_p(t, ev, g)) < 0.05
  }
  colMeans(hits)
}
reps <- 20000
cat("Part 1: false-positive rate (p < 0.05), no true difference,", reps, "runs per row\n")
for (n in c(6, 8, 10, 20)) {
  r <- tests_null(n, reps)
  cat(sprintf("n = %2d per group: log-rank %.4f  Cox Wald %.4f  Cox LR %.4f  Cox score %.4f\n", n, r[1], r[2], r[3], r[4]))
}
cat("Part 2: crossing curves (treated hazard 3x before day 20, 0.25x after), 5000 runs per row\n")
for (n in c(10, 20, 40)) {
  r <- cross_power(n, 5000)
  cat(sprintf("n = %2d per group: share with p < 0.05  log-rank %.4f  RMST(60) %.4f\n", n, r[1], r[2]))
}
# True RMST(60) difference under this scenario, by numerical integration
S0 <- function(x) exp(-rate * x)
S1 <- function(x) ifelse(x <= 20, exp(-3 * rate * x), exp(-3 * rate * 20) * exp(-0.25 * rate * (x - 20)))
cat(sprintf("True RMST(60): control %.3f, treated %.3f days\n", integrate(S0, 0, 60)$value, integrate(S1, 0, 60)$value))