# Permutation test of the RMST(60) difference for small groups, for /guides/log-rank-vs-cox-regression/.
# R 4.3.3. Seed 20261112, fixed before running. Same null and crossing scenarios as simulation.R.
# The permutation p-value is (1 + number of permutations with |diff| >= observed) / (1 + B), B = 999.
# Also: the worked example's permutation p (seed 1, B = 9999).
rate <- log(2) / 35
rmst_km <- function(t, ev, tau = 60) {
  o <- order(t); t <- t[o]; ev <- ev[o]
  ut <- unique(t[ev == 1 & t < tau])
  s <- 1; prev <- 0; area <- 0
  for (u in ut) {
    area <- area + s * (u - prev)
    n <- sum(t >= u); d <- sum(t == u & ev == 1)
    s <- s * (1 - d / n); prev <- u
  }
  area + s * (tau - prev)
}
perm_p <- function(t, ev, g, B = 999) {
  obs <- rmst_km(t[g == 1], ev[g == 1]) - rmst_km(t[g == 0], ev[g == 0])
  cnt <- 0
  for (b in seq_len(B)) {
    gp <- sample(g)
    dd <- rmst_km(t[gp == 1], ev[gp == 1]) - rmst_km(t[gp == 0], ev[gp == 0])
    if (abs(dd) >= abs(obs) - 1e-12) cnt <- cnt + 1
  }
  (1 + cnt) / (1 + B)
}
set.seed(20261112)
reps <- 2000
for (n in c(10, 20)) {
  hit0 <- logical(reps); hit1 <- logical(reps)
  for (i in seq_len(reps)) {
    g <- rep(0:1, each = n)
    t <- rexp(2 * n, rate); ev <- as.integer(t <= 60); t <- pmin(t, 60)
    hit0[i] <- perm_p(t, ev, g) < 0.05
    t0 <- rexp(n, rate); 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)
    hit1[i] <- perm_p(t, ev, g) < 0.05
  }
  cat(sprintf("n = %2d per group, %d runs, B = 999: permutation RMST test p < 0.05 with no true difference %.4f; crossing curves %.4f\n",
              n, reps, mean(hit0), mean(hit1)))
}
d <- read.csv("survival.csv"); g <- as.integer(d$group == "treated")
set.seed(1)
cat(sprintf("Worked example: RMST difference %.3f days, permutation p = %.4f (B = 9999)\n",
            rmst_km(d$day[g == 1], d$event[g == 1]) - rmst_km(d$day[g == 0], d$event[g == 0]),
            perm_p(d$day, d$event, g, B = 9999)))