# =============================================================================
# Retrospective Sample Size Evaluation — Event-Free Survival (2 Groups)
# Freedman Method | works directly from hazard ratio and observed events
# Appropriate for: unequal group sizes, single censored-at-event designs
# =============================================================================

# --- 0. Install & Load Packages ----------------------------------------------

# Uncomment to install if needed:
# install.packages("gsDesign")

library(gsDesign)


# --- 1. Observed Data --------------------------------------------------------

# n           : number of participants per group
# events      : number of observed events per group
# person_time : total person-years of follow-up per group
#               = sum of each individual's follow-up duration
#               naturally accounts for censoring and variable follow-up
# max_follow_up : maximum follow-up time (years)

n_A           <- 323;   events_A <- 80;   person_time_A <- 2089
n_B           <- 56;    events_B <- 32;   person_time_B <- 294
max_follow_up <- 9      # years

# Derived quantities
hazard_A <- events_A / person_time_A   # events per person-year
hazard_B <- events_B / person_time_B
hr       <- hazard_B / hazard_A        # observed hazard ratio (B vs A)

overall_event_rate <- (events_A + events_B) / (n_A + n_B)
pi_A     <- n_A / (n_A + n_B)
pi_B     <- n_B / (n_A + n_B)

cat("=== Observed Data Summary ===\n")
cat(sprintf("Group A : N=%-4d  Events=%-3d  Person-years=%6.0f  Hazard=%.4f\n",
            n_A, events_A, person_time_A, hazard_A))
cat(sprintf("Group B : N=%-4d  Events=%-3d  Person-years=%6.0f  Hazard=%.4f\n",
            n_B, events_B, person_time_B, hazard_B))
cat(sprintf("Hazard ratio (B vs A)  : %.3f\n", hr))
cat(sprintf("Allocation (A / B)     : %.1f%% / %.1f%%\n", pi_A*100, pi_B*100))
cat(sprintf("Overall event rate     : %.3f\n", overall_event_rate))


# --- 2. Required Events — Freedman Method ------------------------------------
# Freedman (1982): works directly from the hazard ratio.
# More stable than Schoenfeld when group sizes are unequal, because it does
# not depend on a variance term V that collapses under unequal allocation.
#
# Formula: E = (z_alpha/2 + z_beta)^2 * (1 + hr)^2 / (hr - 1)^2
#
# Note: hr must not equal 1 (no difference). If hr < 1, invert it so
# the formula always uses the ratio > 1.

required_events_freedman <- function(alpha = 0.05, power = 0.80, hr) {
  if (hr == 1) stop("Hazard ratio must not equal 1 (no detectable difference).")
  hr      <- ifelse(hr < 1, 1/hr, hr)   # ensure hr > 1
  z_alpha <- qnorm(1 - alpha / 2)
  z_beta  <- qnorm(power)
  E       <- (z_alpha + z_beta)^2 * (1 + hr)^2 / (hr - 1)^2
  return(ceiling(E))
}

E_required <- required_events_freedman(alpha = 0.05, power = 0.80, hr = hr)
N_required <- ceiling(E_required / overall_event_rate)

cat("\n=== Freedman Method (alpha=0.05, power=0.80) ===\n")
cat("Required total events :", E_required, "\n")
cat("Observed total events :", events_A + events_B, "\n")
cat("Required total N      :", N_required, "\n")
cat("Observed total N      :", n_A + n_B, "\n")
cat("Difference (req-obs)  :", N_required - (n_A + n_B), "\n")


# --- 3. Validation via gsDesign::nSurv ---------------------------------------

cat("\n=== gsDesign Validation ===\n")
tryCatch({
  gs <- nSurv(
    lambdaC = hazard_A,
    hr      = hr,
    sided   = 2,
    alpha   = 0.05,
    beta    = 0.20,
    T       = max_follow_up,
    minfup  = 0
  )
  print(gs)
}, error = function(e) {
  cat("gsDesign error:", conditionMessage(e), "\n")
})


# --- 4. Sensitivity Analysis Across Power Levels -----------------------------

power_levels <- c(0.70, 0.75, 0.80, 0.85, 0.90)

sensitivity <- data.frame(
  power      = power_levels,
  events_req = sapply(power_levels, function(pw)
    required_events_freedman(0.05, pw, hr))
)
sensitivity$n_req    <- ceiling(sensitivity$events_req / overall_event_rate)
sensitivity$adequate <- ifelse(sensitivity$n_req <= (n_A + n_B), "Yes", "No")

cat("\n=== Sensitivity Analysis by Power Level ===\n")
print(sensitivity)


# --- 5. Interpretation -------------------------------------------------------

cat("\n=== Interpretation ===\n")
cat(sprintf("Hazard ratio: %.3f — Group B events occur %.1fx faster than Group A\n",
            hr, hr))

if (E_required <= (events_A + events_B)) {
  cat("Event count : ADEQUATE — observed events meet Freedman requirement.\n")
} else {
  cat(sprintf("Event count : INSUFFICIENT — %d observed, %d required.\n",
              events_A + events_B, E_required))
}

if (N_required <= (n_A + n_B)) {
  cat("Sample size : ADEQUATELY POWERED at 80%.\n")
} else {
  cat(sprintf(paste0("Sample size : UNDERPOWERED at 80%% — ",
                     "N=%d observed, N=%d required (shortfall: %d).\n"),
              n_A + n_B, N_required, N_required - (n_A + n_B)))
  cat("Type II error risk is elevated — interpret null results with caution.\n")
}