# Adapted from Jörgens & Wassmer's supplied sims_fixed-variance.R.
# This teaching variant explicitly disables early efficacy for both methods.
# It generates normal group means directly (exact with known sigma = 1),
# uses identical data for both analyses, and needs only `rpact`.
library(rpact)

## ---- analyze-paired-trial
analyze_trial <- function(dataSet, designIN, designCD) {
  stageIN <- getStageResults(
    designIN,
    dataInput = dataSet,
    intersectionTest = "Dunnett",
    normalApproximation = TRUE
  )
  stageCD <- getStageResults(
    designCD,
    dataInput = dataSet,
    intersectionTest = "Dunnett",
    normalApproximation = TRUE
  )
  IN <- getClosedCombinationTestResults(stageIN)
  CD <- getClosedConditionalDunnettTestResults(stageCD)
  global <- which(rowSums(IN$indices) == 4)
  stopifnot(length(global) == 1L)
  list(IN = IN, CD = CD, global = global)
}
## ---- end-analysis

paired_simulation <- function(
  nsim = 1000L,
  keep = 1L,
  effects = c(0, 0, 0, 0.25),
  seed = 20261015L
) {
  stopifnot(
    length(effects) == 4L,
    all(is.finite(effects)),
    keep %in% 1:4,
    length(keep) == 1L,
    nsim >= 2,
    nsim == as.integer(nsim)
  )
  set.seed(seed)
  designIN <- getDesignInverseNormal(
    kMax = 2,
    alpha = 0.025,
    informationRates = c(0.5, 1),
    typeOfDesign = "noEarlyEfficacy"
  )
  designCD <- getDesignConditionalDunnett(
    alpha = 0.025,
    informationAtInterim = 0.5,
    secondStageConditioning = TRUE
  )
  means <- c(effects, 0) # rpact: control is last
  records <- replicate(nsim, {
    first <- rnorm(5, means, 1 / sqrt(100))
    selected <- order(first[1:4], decreasing = TRUE)[seq_len(keep)]
    retained <- c(selected, 5L)
    n2 <- rep(NA_real_, 5)
    n2[retained] <- ceiling(500 / (keep + 1))
    second <- rep(NA_real_, 5)
    second[retained] <- rnorm(
      keep + 1L,
      means[retained],
      1 / sqrt(n2[retained])
    )
    args <- list()
    for (j in 1:5) {
      args[[paste0("n", j)]] <- c(100, n2[j])
      args[[paste0("means", j)]] <- c(first[j], second[j])
      args[[paste0("stDevs", j)]] <- c(1, if (is.na(n2[j])) NA_real_ else 1)
    }
    dataSet <- do.call(getDataset, args)
    analysis <- analyze_trial(dataSet, designIN, designCD)
    IN <- analysis$IN
    CD <- analysis$CD
    global <- analysis$global
    c(
      IN = any(IN$rejected[, 2]),
      CD = any(CD$rejected[, 2]),
      globalIN = IN$rejectedIntersections[global, 2],
      globalCD = CD$rejectedIntersections[which(rowSums(CD$indices) == 4L), 2],
      trueIN = any(IN$rejected[effects > 0, 2]),
      trueCD = any(CD$rejected[effects > 0, 2])
    )
  })
  trials <- as.data.frame(t(records))
  delta <- as.numeric(trials$CD) - as.numeric(trials$IN)
  cdOnly <- trials$CD & !trials$IN
  list(
    summary = data.frame(
      nsim,
      keep,
      anyIN = mean(trials$IN),
      anyCD = mean(trials$CD),
      disjunctiveIN = mean(trials$trueIN),
      disjunctiveCD = mean(trials$trueCD),
      difference = mean(delta),
      mcse_difference = sd(delta) / sqrt(nsim),
      IN_only = mean(trials$IN & !trials$CD),
      CD_only = mean(cdOnly),
      CD_only_global_IN_fails = mean(cdOnly & !trials$globalIN),
      CD_only_global_IN_passes = mean(cdOnly & trials$globalIN)
    ),
    trials = trials
  )
}

# Partition each rejection bar on the denominator of ALL simulated trials.
# Black: this method rejects, the other passes its global test but rejects no arm.
rejection_components <- function(trials) {
  both <- trials$IN & trials$CD
  inOnly <- trials$IN & !trials$CD
  cdOnly <- trials$CD & !trials$IN
  data.frame(
    method = c("Inverse normal", "Conditional Dunnett"),
    concordant = rep(mean(both), 2),
    discordant_global_fails = c(
      mean(inOnly & !trials$globalCD),
      mean(cdOnly & !trials$globalIN)
    ),
    discordant_global_passes = c(
      mean(inOnly & trials$globalCD),
      mean(cdOnly & trials$globalIN)
    )
  )
}

plot_rejection_components <- function(components) {
  old <- par(
    mar = c(3.3, 4.2, 1.8, 0.5),
    family = "sans",
    las = 1,
    fg = "#173747",
    col.axis = "#173747",
    col.lab = "#173747"
  )
  on.exit(par(old))
  heights <- t(as.matrix(components[, c(
    "discordant_global_fails",
    "discordant_global_passes",
    "concordant"
  )]))
  totals <- colSums(heights)
  plot.new()
  plot.window(
    xlim = c(0.4, 2.6),
    ylim = c(0, max(totals) + 0.08),
    xaxs = "i",
    yaxs = "i"
  )
  ticks <- seq(0, max(totals) + 0.05, by = 0.1)
  abline(h = ticks, col = "#e4eaed")
  for (j in 1:2) {
    tops <- cumsum(heights[, j])
    rect(
      j - 0.30,
      c(0, head(tops, -1)),
      j + 0.30,
      tops,
      col = c("#bfc3c6", "#161616", c("#a9cfdf", "#b0d9a4")[j]),
      border = NA
    )
    text(j, totals[j] + 0.025, sprintf("%.1f%%", 100 * totals[j]), font = 2)
  }
  axis(1, at = 1:2, labels = components$method, tick = FALSE)
  axis(2, at = ticks, labels = paste0(round(ticks * 100), "%"))
  title(ylab = "Probability of any elementary rejection")
  box(bty = "l", col = "#8b9ba3")
  invisible(components)
}

# 10,000 iterations × six scenarios × four selection rules is intentionally
# not executed while rendering slides. Pass a count for a standalone run:
# Rscript R/paired-diagnostics.R 1000
if (sys.nframe() == 0L) {
  args <- commandArgs(trailingOnly = TRUE)
  nsim <- if (length(args)) as.integer(args[1]) else 1000L
  print(paired_simulation(nsim = nsim)$summary)
  print(sessionInfo())
}