Analyzing Binary Outcomes in a Trial With Small Sample Sizes

mathematical statistics
sap
trial design
R
A pilot Monte Carlo simulation study to evaluate and compare the performance of three statistical methods for analyzing binary outcomes in a randomized clinical trial with small sample sizes and stratification.
Author

Yupng L.

Published

January 18, 2026

We will be conducting a pilot Monte Carlo simulation study to evaluate and compare the performance of three statistical methods for analyzing binary outcomes in a randomized clinical trial with small sample sizes and stratification (maybe seen in small biotech trials???).

The Statistical Objective

The Risk Difference (RD) will be estimated between an Investigational Product (IP/TRT) and Placebo (PLB), adjusting for a stratification factor (LEVEL1 vs. LEVEL2).

We will be comparing three specific approaches:

  1. Miettinen-Nurminen (MN): Stratified Risk Difference using CMH weights (via metalite.ae in R).

  2. Mantel-Haenszel (MH): Stratified Risk Difference and Sato’s confidence intervals (via risks and stats).

  3. Fisher’s Exact Test: Unstratified analysis (ignoring the stratum).

Simulation Design

Seven scenarios (NULLS, SCEN1–SCEN6) are defined varying the following parameters:

  • Sample Size: Fixed at a very small number (perarm = 16, total N = 32).

  • Stratum Prevalence: Varying from 10% to 40%.

  • Effect Sizes:

    • Treatment (TRT): Base effect (effectip) + an additional benefit if the subject is allocated to a specific stratum (effectplus).

    • Placebo (PLB): Fixed effect (effectplb).

The Workflow

Data Generation (simulate_once)

For each iteration, we generate a clinical trial dataset:

  • Allocation: Split the 16 subjects per arm into different strata.

  • Outcome Generation: Generate binary responses (resp) using runif based on the defined probabilities.

    • Note: Assume the Placebo rate is constant across strata, but the Treatment rate increases by effectplus in one specific stratum.

  • True Estimand: Calculate the true_rd (True Risk Difference) as a weighted average of the risk differences in the two strata.

Statistical Analysis

Inside the loop, run the three tests on the generated data:

  1. Fisher’s Exact Test:

    • Calculates a p-value ignoring the strata.

    • Metric: Rejection rate (p < 0.05).

  2. Mantel-Haenszel (MH):

    • Calculates the stratified RD point estimate and CI.

    • Metrics: CI Width, Coverage Probability (does CI contain true_rd?), and Rejection rate (mantelhaen.test p < 0.05).

  3. Miettinen-Nurminen (MN; not exact the naïve methods proposed in 1985):

    • Calculates stratified RD using rate_compare with weight = “cmh” (proposed by Kaifeng Lu).

    • Metrics: CI Width, Coverage Probability, and Rejection rate.

Execution (run_scenario_parallel)

  • Use Parallel Computing (parallel package) to run >=1000 simulations (n_iter) for each scenario.

  • Correctly handle Random Number Generation (RNG) in parallel by using clusterSetRNGStream and passing a unique seed i to every worker.

Aggregation

  • Finally, take the column means of results. This gives us the average performance (e.g., “Power” or “Type I Error” depending on the scenario, and average CI width) for each method under each scenario.

Measuring

  1. MH_cover / MN_cover: How often the 95% CI captures the true difference. (Is the method valid?)

  2. MH_width / MN_width: How precise are the estimates? (Narrower is better, provided coverage is maintained).

  3. MH_rej_h0 / MN_rej_h0 / Fish_rej_h0:

    • If True RD = 0 (Scenario NULLS): This measures Type I Error (False Positive Rate).

    • If True RD \(\neq\) 0 (Scenarios SCEN1-6): This measures Power.

R Program

rm(list=ls())

library(parallel)
library(dplyr)
library(metalite.ae)
library(risks)

n_cores <- max(1, detectCores() - 1)  # leave one core free

scenarios <- list(
  NULLS = list(scen=1, perarm=16, stratprevalence=0.2, effectip=0.45,effectplus=0.0,effectplb=0.45),
  SCEN1 = list(scen=2, perarm=16, stratprevalence=0.1, effectip=0.45,effectplus=0.1,effectplb=0.1),
  SCEN2 = list(scen=3, perarm=16, stratprevalence=0.2, effectip=0.45,effectplus=0.1,effectplb=0.1),
  SCEN3 = list(scen=4, perarm=16, stratprevalence=0.3, effectip=0.45,effectplus=0.1,effectplb=0.1),
  SCEN4 = list(scen=5, perarm=16, stratprevalence=0.4, effectip=0.45,effectplus=0.1,effectplb=0.1),
  SCEN5 = list(scen=6, perarm=16, stratprevalence=0.3, effectip=0.4,effectplus=0.2,effectplb=0.15),
  SCEN6 = list(scen=7, perarm=16, stratprevalence=0.3, effectip=0.35,effectplus=0.1,effectplb=0.15)
)

simulate_once <- function(params,
                          seed){
  
  if (!is.null(seed)) set.seed(seed)
  
  scen <- params$scen
  perarm <- params$perarm
  stratprevalence <- params$stratprevalence
  effectip <- params$effectip
  effectplus <- params$effectplus
  effectplb <- params$effectplb
  
  # sample size per arm
  samparm <- perarm
  
  # prevalence of stratification
  stratum <- stratprevalence
  
  # --- Simulation Parameters ---
  # n[1]: stratum 1
  # n[2]: stratum 2
  # IP and PBO sizes per stratum
  n1 <- c(samparm-ceiling(samparm*stratum), ceiling(samparm*stratum)); 
  n2 <- c(samparm-ceiling(samparm*stratum), ceiling(samparm*stratum)); 
  
  # effect size for IP and Placebo
  # subjects with biology use history have better effect + ~x%
  effectplb <- effectplb
  effectip <- effectip
  effectplus <- effectplus
  p1 <- c(effectip, effectip+effectplus); 
  p2 <- rep(effectplb,2)
  true_rd <- p1[1]*(1-stratum)+p1[2]*stratum-effectplb
  
  # Storage for widths
  widths <- data.frame(MH_cover = numeric(1), MH_width = numeric(1), MH_rej_h0 = numeric(1), 
                       MN_cover = numeric(1), MN_width = numeric(1), MN_rej_h0 = numeric(1),
                       Fish_rej_h0 = numeric(1))
  
  
  if (exists("simdata")) {rm(list="simdata")}
  
  simdata <- data.frame(
    treatment  = c(rep("TRT", n1[1]),  
                   rep("TRT", n1[2]), 
                   rep("PLB", n2[1]),
                   rep("PLB", n2[2])),
    
    strat      = c(rep("LEVEL 1", n1[1]),  
                   rep("LEVEL 2", n1[2]), 
                   rep("LEVEL 1", n2[1]),
                   rep("LEVEL 2", n2[2])),
    
    probuni = runif(samparm*2)
    
  )
  
  simdata <- simdata %>% mutate(resp=case_when(
    treatment=="TRT" & strat == "LEVEL 1" & probuni <= p1[1] ~ 1,
    treatment=="TRT" & strat == "LEVEL 2" & probuni <= p1[2] ~ 1,
    treatment=="PLB" & probuni <= p2[1] ~ 1,
    TRUE ~ 0
  )
  )
  
  # prepare frequency table for fisher
  # Create a contingency table
  table_data <- table(simdata$treatment, simdata$resp)
  
  # Perform Fisher's Exact Test
  test_result <- fisher.test(table_data)
  
  fshp <- test_result$p.value<0.05
  
  # extend to datasets for calling R package --- MH
  
  mh <- risks::rr_rd_mantel_haenszel(
    data=simdata,
    exposure=treatment,
    outcome=resp,
    confounders=strat,
    estimand = c("rd"),
    conf.level = 0.95
  )
  
  mhcover <- mh$conf.low<=true_rd & true_rd<=mh$conf.high
  
  cmhtest = mantelhaen.test(x = simdata$treatment, y = simdata$resp, z = simdata$strat)
  
  mhrej <- cmhtest$p.value<0.05
  
  # extend to datasets for calling R package --- MN with CMH weights 
  
  mn <- rate_compare(
    data=simdata,
    resp ~ factor(treatment, levels = c("PLB","TRT")),
    strata = strat,
    delta = 0,
    weight = "cmh",
    test = "two.sided",
    alpha = 0.05
  )
  
  mncover <- mn$lower<=true_rd & true_rd<=mn$upper  
  
  mnrej <- mn$p<0.05
  
  # MN
  widths$MN_width[1] <- mn$upper - mn$lower
  widths$MN_cover[1] <- mncover
  widths$MN_rej_h0[1] <- mnrej
  
  # MH Sato
  widths$MH_width[1] <- mh$conf.high - mh$conf.low
  widths$MH_cover[1] <- mhcover
  widths$MH_rej_h0[1] <- mhrej
  
  widths$Fish_rej_h0 <- fshp
  
  widths$TrueRD <- true_rd
  
  return(widths)
}

run_scenario_parallel <- function(params, n_iter, n_cores) {
  cl <- makeCluster(n_cores)
  on.exit(stopCluster(cl))
  
  # Export the function to worker nodes
  clusterExport(cl, varlist = c("simulate_once"), envir = environment())
  
  # Load dplyr in each worker environment
  clusterEvalQ(cl, {
    library(dplyr)
    library(risks)
    library(metalite.ae)
  })  
  
  # Set different RNG streams for each worker
  clusterSetRNGStream(cl, 123)  # 123 is the seed to initialize RNG, ensuring variation for each worker
  
  # Run simulations in parallel (each worker gets a different stream)
  res <- parLapply(
    cl,
    X = seq_len(n_iter),
    fun = function(i) simulate_once(params, seed = i)  # passing a unique seed per iteration
  )
  
  # Combine results into a matrix
  do.call(rbind, res)
}

n_iter <- 1000

results <- lapply(
  scenarios,
  run_scenario_parallel,
  n_iter = n_iter,
  n_cores = n_cores
)

calculate_col_means <- function(df, id) {
  col_means <- colMeans(df, na.rm = TRUE)  # Calculate column means
  col_means <- as.data.frame(t(col_means))  # Transpose to match the row format
  col_means$ID <- id  # Add the ID as a column
  return(col_means)
}

# Apply the function to each dataframe in the list and bind the results together
all <- lapply(names(results), function(id) {
  calculate_col_means(results[[id]], id)  # Pass the dataframe and its ID
})

# Combine all results into one final data frame
final_results <- do.call(rbind, all)

print(final_results)

Reference

Miettinen, C.W. (2008). Nurminen Method for Confidence Intervals of the Difference in Binomial Proportions from Stratified 2 x 2 Samples.

Appendix: An Update Regarding PROC FREQ on SAS VIYA

LTS 2023.10 (NOV 2023)

The new COMMONRISKDIFF(CL=MN) option in the TABLES statement provides Miettinen-Nurminen confidence limits for the common risk difference, which is the overall difference in binomial proportions for stratified 2 × 2 tables.

The COMMONRISKDIFF(CL=MNMH) option provides confidence limits that are computed by the Miettinen-Nurminen method with Mantel-Haenszel weights (Lu 2008). The PLOTS=RISKDIFFPLOT(COMMON=MN) option enables you to display Miettinen-Nurminen confidence limits for the common risk difference in risk difference plots.