# Binary trial planning: standalone educational exercise, 18 October 2026.
# Base R only. All planning assumptions and participant data are synthetic.
# Run: Rscript --vanilla binary-planning.R OUTPUT_DIRECTORY
# In RStudio: source("binary-planning.R"); run_exercise("results")

power_binary <- function(p_control, p_treatment, n_control, n_treatment,
                         alpha = 0.05) {
  stopifnot(length(p_control) == 1L, length(p_treatment) == 1L,
    all(is.finite(c(p_control, p_treatment, n_control, n_treatment, alpha))),
    p_control > 0, p_control < 1, p_treatment > 0, p_treatment < 1,
    n_control > 0, n_treatment > 0,
    n_control == floor(n_control), n_treatment == floor(n_treatment),
    alpha > 0, alpha < 1)
  difference <- p_treatment - p_control
  pooled <- (n_treatment * p_treatment + n_control * p_control) /
    (n_treatment + n_control)
  v_null <- pooled * (1 - pooled) * (1 / n_control + 1 / n_treatment)
  v_alternative <- p_treatment * (1 - p_treatment) / n_treatment +
    p_control * (1 - p_control) / n_control
  critical <- stats::qnorm(1 - alpha / 2) * sqrt(v_null)
  stats::pnorm((difference - critical) / sqrt(v_alternative)) +
    stats::pnorm((-difference - critical) / sqrt(v_alternative))
}

plan_binary <- function(p_control, p_treatment, target = 0.90, alpha = 0.05,
                        allocation = c(1L, 1L), missing = 0) {
  # allocation is treatment:control; one integer block preserves this ratio.
  stopifnot(length(allocation) == 2L, all(is.finite(allocation)),
    all(allocation >= 1), all(allocation == floor(allocation)),
    length(target) == 1L, is.finite(target), target >= 0.5, target < 1,
    length(missing) == 1L, is.finite(missing), missing >= 0, missing <= 0.99,
    p_control != p_treatment)
  evaluate <- function(k) power_binary(p_control, p_treatment,
    allocation[2] * k, allocation[1] * k, alpha)
  maximum <- floor(1e6 / max(allocation))
  stopifnot(maximum >= 1, target > alpha)
  lower <- 0L
  upper <- 1L
  while (evaluate(upper) < target) {
    if (upper == maximum) stop("Target exceeds the one-million-per-arm limit")
    lower <- upper
    upper <- min(maximum, upper * 2L)
  }
  while (upper - lower > 1L) {
    middle <- floor((upper + lower) / 2)
    if (evaluate(middle) >= target) upper <- middle else lower <- middle
  }
  # Percentage form avoids cancellation at common whole-percentage boundaries.
  random_blocks <- ceiling(upper * 100 / (100 - missing * 100))
  if (random_blocks * max(allocation) > 1e6) stop("Randomized count exceeds limit")
  data.frame(p_control = p_control, p_treatment = p_treatment,
    rd = p_treatment - p_control, rr = p_treatment / p_control,
    alpha = alpha, target = target, missing = missing,
    n_control = allocation[2] * upper, n_treatment = allocation[1] * upper,
    random_control = allocation[2] * random_blocks,
    random_treatment = allocation[1] * random_blocks,
    total_randomized = sum(allocation) * random_blocks,
    power = evaluate(upper))
}

make_synthetic <- function() {
  # Fixed seed and generator; these are not SatraGO records or clinical evidence.
  set.seed(20261019, kind = "Mersenne-Twister", normal.kind = "Inversion",
    sample.kind = "Rejection")
  data.frame(id = sprintf("SYN-%03d", seq_len(400)),
    arm = rep(c("control", "treatment"), each = 200),
    response = c(stats::rbinom(200, 1, 0.25), stats::rbinom(200, 1, 0.40)))
}

compare_synthetic <- function() {
  # Independent synthetic normal estimates; neither is a SatraGO estimate.
  estimate <- c(0.12, 0.16)
  se <- c(0.07, 0.06)
  estimate <- c(estimate, estimate[2] - estimate[1])
  se <- c(se, sqrt(sum(se^2)))
  data.frame(comparison = c("synthetic_A", "synthetic_B", "B_minus_A"),
    estimate = estimate, se = se,
    lower95 = estimate - stats::qnorm(0.975) * se,
    upper95 = estimate + stats::qnorm(0.975) * se,
    p = 2 * stats::pnorm(-abs(estimate / se)))
}

run_exercise <- function(output = "results") {
  dir.create(output, recursive = TRUE, showWarnings = FALSE)
  scenarios <- expand.grid(p_control = c(0.15, 0.25, 0.35),
    rd = c(0.10, 0.15, 0.20))
  grid <- do.call(rbind, lapply(seq_len(nrow(scenarios)), function(i)
    plan_binary(scenarios$p_control[i], scenarios$p_control[i] + scenarios$rd[i],
      missing = 0.10)))
  grid <- grid[order(round(grid$rd, 10), grid$p_control), ]
  rownames(grid) <- NULL
  utils::write.csv(grid, file.path(output, "planning-scenarios.csv"), row.names = FALSE)
  allocation <- rbind(plan_binary(0.25, 0.40, missing = 0.10),
    plan_binary(0.25, 0.40, allocation = c(2L, 1L), missing = 0.10))
  utils::write.csv(allocation, file.path(output, "allocation-comparison.csv"), row.names = FALSE)
  synthetic <- make_synthetic()
  utils::write.csv(synthetic, file.path(output, "synthetic-participants.csv"), row.names = FALSE)
  demonstration <- compare_synthetic()
  utils::write.csv(demonstration, file.path(output, "synthetic-comparison.csv"), row.names = FALSE)
  satrago <- data.frame(study = c("SatraGO-1", "SatraGO-2"),
    response_treatment = c(49.0, 52.9), response_control = c(31.2, 23.4),
    adjusted_rd_pp = c(17.9, 28.9), lower9503_pp = c(-1.6, 11.5),
    upper9503_pp = c(37.3, 46.2), p = c(0.0715, 0.0011))
  utils::write.csv(satrago, file.path(output, "satrago-published-summary.csv"), row.names = FALSE)
  # Reported adjusted effects/95.03% CIs, not a reanalysis of trial participants.
  for (lang in c("es", "en")) {
    es <- lang == "es"
    grDevices::png(file.path(output, paste0("satrago-effects-", lang, ".png")),
      width = 1600, height = 750, res = 180)
    graphics::par(mar = c(5.2, 7, 4.5, 2), family = "sans")
    graphics::plot(satrago$adjusted_rd_pp, c(2, 1),
      xlim = c(-10, 55), ylim = c(0.5, 2.5), yaxt = "n", pch = 19,
      col = "#0b5563", cex = 1.5,
      xlab = if (es) "Diferencia ajustada de respuesta (puntos porcentuales)" else
        "Adjusted response difference (percentage points)", ylab = "",
      main = if (es) "SatraGO: estimaciones e incertidumbre" else
        "SatraGO: estimates and uncertainty")
    graphics::axis(2, at = c(2, 1), labels = satrago$study, las = 1)
    graphics::abline(v = 0, lty = 2, col = "#687582")
    graphics::arrows(satrago$lower9503_pp, c(2, 1), satrago$upper9503_pp,
      c(2, 1), angle = 90, code = 3, length = 0.06, lwd = 2, col = "#0b5563")
    graphics::mtext(if (es) "IC publicados del 95,03%; enfermedad activa, semana 24" else
      "Published 95.03% CIs; active disease, week 24", side = 3, line = 0.5, cex = 0.8)
    grDevices::dev.off()
    grDevices::png(file.path(output, paste0("planning-sensitivity-", lang, ".png")),
      width = 1600, height = 950, res = 180)
    graphics::par(mar = c(5, 5.5, 4.5, 1.5), family = "sans")
    graphics::plot(NA, xlim = c(15, 35), ylim = c(0, max(grid$total_randomized) * 1.12),
      xlab = if (es) "Respuesta asumida en control (%)" else "Assumed control response (%)",
      ylab = if (es) "Total que aleatorizar" else "Total to randomize",
      main = if (es) "El resultado cambia con los supuestos" else "Assumptions change the answer")
    colours <- c("#823c39", "#0b5563", "#46549b")
    for (j in seq_along(c(0.10, 0.15, 0.20))) {
      part <- grid[abs(grid$rd - c(0.10, 0.15, 0.20)[j]) < 1e-10, ]
      graphics::lines(100 * part$p_control, part$total_randomized,
        type = "b", pch = j + 14L, col = colours[j], lwd = 2)
    }
    graphics::legend("topleft", legend = paste0("RD +", c(10, 15, 20), " pp"),
      col = colours, pch = 15:17, lty = 1, bty = "n", cex = 0.9)
    graphics::mtext(if (es) "Escenarios sint\u00e9ticos: 90% potencia; alfa bilateral 0,05; 1:1; 10% ausentes" else
      "Synthetic scenarios: 90% power; two-sided alpha 0.05; 1:1; 10% missing", side = 3,
      line = 0.5, cex = 0.75)
    grDevices::dev.off()
  }
  treatment <- synthetic$response[synthetic$arm == "treatment"]
  control <- synthetic$response[synthetic$arm == "control"]
  analysis <- stats::prop.test(c(sum(treatment), sum(control)),
    c(length(treatment), length(control)), correct = FALSE)
  capture.output(analysis, file = file.path(output, "synthetic-analysis.txt"))
  capture.output(utils::sessionInfo(), file = file.path(output, "sessionInfo.txt"))
  print(grid[c("p_control", "p_treatment", "n_control", "total_randomized", "power")])
  print(allocation[c("n_control", "n_treatment", "total_randomized")])
  print(demonstration)
  invisible(grid)
}

if (sys.nframe() == 0L) {
  args <- commandArgs(trailingOnly = TRUE)
  run_exercise(if (length(args)) args[1] else "results")
}