132 lines
No EOL
4.9 KiB
R
132 lines
No EOL
4.9 KiB
R
# =============================================================================
|
|
# 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")
|
|
} |