gsDesignNB

Negative Binomial / Recurrent Event Designs

Guide users through sample size calculation, group sequential design, and simulation for clinical trials with negative binomial (recurrent event) outcomes using the gsDesignNB R package. Use this skill when the user asks about: negative binomial sample size, recurrent event trials, overdispersed counts, event gaps, rate ratios, Wald test for count data, seasonal event rates, blinded or unblinded sample size re-estimation, group sequential designs for negative binomial endpoints, or the Zhu-Lakkis method.

Sample Size and Simulation for Negative Binomial Outcomes with gsDesignNB

API reference

Full function docs: references/llms.txt (built from local man pages)
Workflow patterns: references/code_patterns.md

Key functions

Sample size calculation

sample_size_nbinom() - Sample size/power for negative binomial outcomes (Zhu-Lakkis Method 3)

Group sequential design

gsNBCalendar() - Group sequential design with calendar-time analysis schedule
compute_info_at_time() - Statistical information at a given calendar time
toInteger() - Round sample sizes to integers preserving allocation ratio
check_gs_bound() - Check if group sequential bounds are crossed
summarize_gs_sim() - Summarize operating characteristics from simulations

Simulation

nb_sim() - Simulate recurrent events (Gamma-Poisson mixture)
nb_sim_seasonal() - Simulate recurrent events with seasonal variation
sim_gs_nbinom() - Simulate multiple group sequential trials

Data cutting and analysis timing

cut_data_by_date() - Cut simulated data at a calendar date
get_analysis_date() - Find date when target event count is reached
get_cut_date() - Find earliest date satisfying multiple analysis criteria
cut_date_for_completers() - Find date when target completers are reached
cut_completers() - Cut data for completers analysis

Statistical testing and estimation

mutze_test() - Wald test for treatment rate ratio (NB or Poisson)
estimate_nb_mom() - Method of moments estimation for NB parameters
calculate_blinded_info() - Blinded information and dispersion estimation

Sample size re-estimation

blinded_ssr() - Blinded SSR using Friede & Schmidli method
unblinded_ssr() - Unblinded SSR using observed group rates

Workflow patterns

For detailed code templates, read references/code_patterns.md.

Topics covered: - Fixed sample size calculation with piecewise accrual, dropout, event gaps - Power calculation from a fixed design - Non-inferiority and super-superiority designs (rr0 parameter) - Group sequential design with calendar-time analysis schedule - Simulation of recurrent events and group sequential trials - Seasonal event rate simulation - Data cutting at interim/final analyses - Wald test (mutze_test) for treatment rate ratio - Blinded and unblinded information estimation at interim - Blinded and unblinded sample size re-estimation - Completers-based interim analysis - Verification of theoretical vs. simulated operating characteristics

Important design considerations

Dispersion parameter k: Controls overdispersion; k = 0 reduces to Poisson. Larger k means more overdispersion. Can be scalar (common) or length-2 vector (group-specific).
Event gaps: After each event, patients are “off risk” for event_gap time units. This reduces effective exposure: lambda_eff = lambda / (1 + lambda * gap). Specified in the same time units as rates.
Variance inflation factor Q: When follow-up varies across patients, Q = E[t^2] / E[t]^2 inflates the variance. sample_size_nbinom() handles this automatically.
Rate parameterization: Rates lambda1 (control) and lambda2 (experimental) are events per unit time. The treatment effect is the rate ratio RR = lambda2/lambda1.
Wald test (Mütze et al.): mutze_test() fits a negative binomial GLM with offset for log exposure. Falls back to Poisson when the NB dispersion estimate is very large (> poisson_threshold).
Calendar-time analysis schedule: gsNBCalendar() takes analysis_times as calendar months. Information at each analysis depends on enrollment pattern, dropout, and follow-up.
Spending time vs. information fraction: For group sequential designs, usTime/lsTime control alpha spending and may differ from the information fraction. This allows calendar-based or event-based spending schedules.
Blinded information: calculate_blinded_info() uses the blinded (pooled) rate and dispersion to estimate information. Can produce extreme values when the NB MLE is unstable — bound dispersion or use planning values as a fallback.
SSR: Blinded SSR (Friede & Schmidli) maintains the blind; unblinded SSR uses observed group rates for more accurate re-estimation but requires unblinding.
check_gs_bound() info_scale: Use “blinded” (default) or “unblinded” to select which information drives bound updates at analysis time.

Code Patterns

Code Patterns for gsDesignNB

Fixed sample size calculation
Power calculation
Event gaps and piecewise accrual
Non-inferiority and super-superiority
Group sequential design
Simulation of recurrent events
Group sequential simulation
Data cutting and analysis timing
Wald test for treatment rate ratio
Blinded and unblinded information
Sample size re-estimation
Completers-based interim analysis
Seasonal event simulation
Verification of operating characteristics

Fixed sample size calculation

sample_size_nbinom() implements the Zhu-Lakkis Method 3 for negative binomial outcomes, accounting for variable accrual, dropout, maximum follow-up, and event gaps.

Basic calculation

library(gsDesignNB)

design <- sample_size_nbinom(
  lambda1 = 0.5,          # Control event rate (events/time unit)
  lambda2 = 0.3,          # Experimental event rate
  dispersion = 0.5,       # Overdispersion parameter k
  power = 0.9,            # Target power
  alpha = 0.025,          # One-sided significance level
  accrual_rate = 10,      # Enrollment rate (patients/time unit)
  accrual_duration = 12,  # Total enrollment duration
  trial_duration = 24     # Total trial duration
)
design
summary(design)

Piecewise accrual ramp-up

design <- sample_size_nbinom(
  lambda1 = 0.5, lambda2 = 0.3, dispersion = 0.5,
  power = 0.9, alpha = 0.025,
  accrual_rate = c(5, 10),        # Rate per period
  accrual_duration = c(3, 3),     # Duration per period (total: 6 months)
  trial_duration = 12
)

With dropout and maximum follow-up

design <- sample_size_nbinom(
  lambda1 = 0.5, lambda2 = 0.3, dispersion = 0.5,
  power = 0.9, alpha = 0.025,
  accrual_rate = 10, accrual_duration = 12,
  trial_duration = 24,
  dropout_rate = 0.05,       # Annual dropout rate (scalar: same for both groups)
  max_followup = 12          # Cap follow-up at 12 months
)

Group-specific parameters

design <- sample_size_nbinom(
  lambda1 = 0.5, lambda2 = 0.3,
  dispersion = c(0.5, 0.3),       # Different dispersion per group
  power = 0.9, alpha = 0.025,
  accrual_rate = 10, accrual_duration = 12, trial_duration = 24,
  dropout_rate = c(0.10, 0.05),   # Different dropout per group
  max_followup = c(12, 12)        # Can differ per group
)

Unequal allocation

design <- sample_size_nbinom(
  lambda1 = 0.5, lambda2 = 0.3, dispersion = 0.5,
  power = 0.9, alpha = 0.025, ratio = 2,   # 2:1 experimental:control
  accrual_rate = 10, accrual_duration = 12, trial_duration = 24
)
# design$n1 = control, design$n2 = experimental (n2 = ratio * n1)

Key output fields

design$n_total        # Total sample size
design$n1             # Control group N
design$n2             # Experimental group N
design$power          # Achieved power
design$total_events   # Expected total events
design$events_n1      # Expected control events
design$events_n2      # Expected experimental events
design$exposure       # Average exposure per group
design$variance       # Variance of log rate ratio

Power calculation

Set power = NULL to compute power for a given sample size:

power_result <- sample_size_nbinom(
  lambda1 = 0.5, lambda2 = 0.4,
  dispersion = 0.5,
  power = NULL,                     # Request power calculation
  alpha = 0.025,
  accrual_rate = design$accrual_rate,
  accrual_duration = design$inputs$accrual_duration,
  trial_duration = design$inputs$trial_duration,
  dropout_rate = design$inputs$dropout_rate,
  max_followup = design$inputs$max_followup
)
power_result$power

Event gaps and piecewise accrual

Event gaps model time “off risk” after each event (e.g., recovery period).

design <- sample_size_nbinom(
  lambda1 = 2.0, lambda2 = 1.0,    # Higher rates (events/year)
  dispersion = 0.5,
  power = 0.9, alpha = 0.025,
  accrual_rate = 10, accrual_duration = 12, trial_duration = 24,
  event_gap = 30 / 365.25          # 30-day gap after each event (in years)
)

# Event gaps reduce effective exposure:
# lambda_eff = lambda / (1 + lambda * gap)
# Exposure_at_risk < Exposure_total
design$exposure_at_risk_n1  # Exposure at risk (excludes gap time)

Non-inferiority and super-superiority

The rr0 parameter specifies the null hypothesis rate ratio.

Non-inferiority (rr0 > 1)

# H0: RR >= 1.1, H1: RR < 1.1
ni_design <- sample_size_nbinom(
  lambda1 = 0.1, lambda2 = 0.09,   # Experimental slightly better
  rr0 = 1.1,                        # Non-inferiority margin
  dispersion = 0.5,
  power = 0.9, alpha = 0.025,
  accrual_rate = 100, accrual_duration = 1, trial_duration = 2
)

Super-superiority (rr0 < 1)

# H0: RR >= 0.5, H1: RR < 0.5 (must show > 50% reduction)
ss_design <- sample_size_nbinom(
  lambda1 = 0.1, lambda2 = 0.02,
  rr0 = 0.5,
  dispersion = 0.5,
  power = 0.9, alpha = 0.025,
  accrual_rate = 100, accrual_duration = 1, trial_duration = 2,
  max_followup = 1
)

Group sequential design

gsNBCalendar() creates a group sequential design for negative binomial outcomes with calendar-time analysis schedule, compatible with gsDesign spending functions.

Basic group sequential design

# Step 1: fixed design
nb_ss <- sample_size_nbinom(
  lambda1 = 1.5/12, lambda2 = 1.0/12, dispersion = 0.5,
  power = 0.9, alpha = 0.025,
  accrual_rate = 1, accrual_duration = 12, trial_duration = 24,
  max_followup = 12, event_gap = 20/30.4375
)

# Step 2: group sequential design
gs_nb <- gsNBCalendar(
  x = nb_ss,
  k = 3,                              # 3 analyses
  test.type = 4,                       # 2-sided, non-binding futility
  alpha = 0.025,
  sfu = gsDesign::sfLDOF,             # Upper spending: Lan-DeMets O'Brien-Fleming
  sfl = gsDesign::sfHSD, sflpar = -8, # Lower spending: Hwang-Shih-DeCani
  analysis_times = c(10, 18, 24)       # Calendar months
)

# Step 3: integer sample sizes
gs_nb <- toInteger(gs_nb)

# Display bounds
gsDesign::gsBoundSummary(gs_nb, deltaname = "RR", logdelta = TRUE)

With custom spending times

gs_nb <- gsNBCalendar(
  x = nb_ss,
  k = 3, test.type = 1,               # One-sided (efficacy only)
  sfu = gsDesign::sfLinear,
  sfupar = c(0.5, 0.5),               # Linear spending parameters
  usTime = c(0.1, 0.18, 1),           # Spending time (may differ from info fraction)
  analysis_times = c(10, 18, 24)
)

Key output fields

gs_nb$n.I           # Information at each analysis
gs_nb$n_total       # Total N at each analysis (vector)
gs_nb$n1            # Control N at each analysis
gs_nb$n2            # Experimental N at each analysis
gs_nb$events        # Expected total events at each analysis
gs_nb$exposure      # Expected exposure at each analysis
gs_nb$T             # Calendar times of analyses
gs_nb$upper$bound   # Upper (efficacy) Z bounds
gs_nb$lower$bound   # Lower (futility) Z bounds

Simulation of recurrent events

nb_sim() simulates a single trial using a Gamma-Poisson mixture model.

Basic simulation

sim_data <- nb_sim(
  enroll_rate = data.frame(rate = 20, duration = 12),
  fail_rate = data.frame(
    treatment = c("Control", "Experimental"),
    rate = c(0.5, 0.3),
    dispersion = c(0.5, 0.5)
  ),
  dropout_rate = data.frame(
    treatment = c("Control", "Experimental"),
    rate = c(0.05, 0.05),
    duration = c(100, 100)
  ),
  max_followup = 12,
  n = 200
)

# Output: one row per event/censoring per patient
head(sim_data)
# id, treatment, enroll_time, tte, calendar_time, event

With enrollment ramp-up and event gaps

sim_data <- nb_sim(
  enroll_rate = data.frame(
    rate = c(5, 10, 20),
    duration = c(2, 2, 8)
  ),
  fail_rate = data.frame(
    treatment = c("Control", "Experimental"),
    rate = c(1.5/12, 1.0/12),
    dispersion = c(0.5, 0.5)
  ),
  dropout_rate = data.frame(
    treatment = c("Control", "Experimental"),
    rate = c(0.001, 0.001),
    duration = c(100, 100)
  ),
  max_followup = 12,
  event_gap = 5/365.25     # 5-day gap after each event
)

Group sequential simulation

sim_gs_nbinom() runs many group sequential trials and returns per-analysis results.

Basic usage

sim_results <- sim_gs_nbinom(
  n_sims = 1000,
  enroll_rate = data.frame(rate = 20, duration = 12),
  fail_rate = data.frame(
    treatment = c("Control", "Experimental"),
    rate = c(1.5/12, 1.0/12),
    dispersion = c(0.5, 0.5)
  ),
  dropout_rate = data.frame(
    treatment = c("Control", "Experimental"),
    rate = c(0.001, 0.001),
    duration = c(100, 100)
  ),
  max_followup = 12,
  event_gap = 20/30.4375,
  analysis_times = c(10, 18, 24),
  n_target = gs_nb$n_total[3],
  design = gs_nb
)

# Check bounds
sim_checked <- check_gs_bound(sim_results, design = gs_nb, info_scale = "blinded")

# Summarize operating characteristics
oc <- summarize_gs_sim(sim_checked)
oc$power                 # Overall power
oc$analysis_summary      # Per-analysis cumulative crossing probabilities

With custom cutting criteria

sim_results <- sim_gs_nbinom(
  n_sims = 500,
  enroll_rate = data.frame(rate = 20, duration = 12),
  fail_rate = fail_rate_df,
  dropout_rate = dropout_df,
  max_followup = 12,
  event_gap = 5/365.25,
  n_target = 240,
  design = gs_nb,
  cuts = list(
    # Analysis 1: calendar time 10 OR 80 events, whichever later
    list(planned_calendar = 10, target_events = 80),
    # Analysis 2: calendar time 18 OR 160 events
    list(planned_calendar = 18, target_events = 160),
    # Analysis 3: calendar time 24 (final)
    list(planned_calendar = 24)
  )
)

Data cutting and analysis timing

Cut at a calendar date

# Aggregate to one row per subject at cut date
cut_summary <- cut_data_by_date(sim_data, cut_date = 18, event_gap = 5/365.25)
# Output: id, treatment, enroll_time, tte, tte_total, events

Find date when target events reached

event_date <- get_analysis_date(sim_data, planned_events = 150, event_gap = 5/365.25)

Find date satisfying multiple criteria

cut_date <- get_cut_date(
  data = sim_data,
  planned_calendar = 18,       # At least 18 months
  target_events = 150,          # At least 150 events
  target_completers = 100,      # At least 100 completers
  event_gap = 5/365.25
)

Find date for completers

comp_date <- cut_date_for_completers(sim_data, target_completers = 80)
cut_comp <- cut_completers(sim_data, cut_date = comp_date, event_gap = 5/365.25)

Wald test for treatment rate ratio

mutze_test() performs a Wald test based on negative binomial (or Poisson) GLM.

# Cut data first
cut_summary <- cut_data_by_date(sim_data, cut_date = 24, event_gap = 5/365.25)

# Negative binomial Wald test
result <- mutze_test(cut_summary)
result$z              # Z statistic
result$p_value        # P-value
result$rate_ratio     # Estimated rate ratio (exp(estimate))
result$conf_int       # Confidence interval for rate ratio
result$dispersion     # Estimated NB dispersion (1/theta)
result$group_summary  # Events, exposure by group

# Poisson fallback
result_poisson <- mutze_test(cut_summary, method = "poisson")

# Two-sided test
result_2sided <- mutze_test(cut_summary, sided = 2, conf_level = 0.95)

Extracting information from the test

# Statistical information = 1 / SE^2
info <- 1 / result$se^2

# Log rate ratio estimate
log_rr <- result$estimate

Blinded and unblinded information

Blinded information at interim

blinded_res <- calculate_blinded_info(
  data = cut_summary,
  ratio = 1,                         # 1:1 randomization
  lambda1_planning = 0.5,            # Planned control rate
  lambda2_planning = 0.3             # Planned experimental rate
)
blinded_res$blinded_info             # Estimated information
blinded_res$dispersion_blinded       # Estimated dispersion
blinded_res$lambda_blinded           # Blinded overall rate

Method of moments estimation

# Blinded (pooled) estimate
mom_blinded <- estimate_nb_mom(cut_summary)
mom_blinded$lambda      # Overall rate
mom_blinded$dispersion  # Dispersion k

# Unblinded (by group) estimate
mom_unblinded <- estimate_nb_mom(cut_summary, group = "treatment")
mom_unblinded$lambda    # Named vector: Control, Experimental
mom_unblinded$dispersion  # Common dispersion

Information at a planning time

info <- compute_info_at_time(
  analysis_time = 18,
  accrual_rate = 20,
  accrual_duration = 12,
  lambda1 = 0.5, lambda2 = 0.3,
  dispersion = 0.5,
  ratio = 1,
  dropout_rate = 0.05,
  event_gap = 5/365.25,
  max_followup = 12
)

Sample size re-estimation

Blinded SSR (Friede & Schmidli)

blinded_ssr_res <- blinded_ssr(
  data = cut_summary,
  ratio = 1,
  lambda1_planning = 0.5,
  lambda2_planning = 0.3,
  power = 0.9,
  alpha = 0.025,
  accrual_rate = 20,
  accrual_duration = 12,
  trial_duration = 24,
  max_followup = 12,
  event_gap = 5/365.25
)
blinded_ssr_res$n_total_blinded    # Re-estimated total N
blinded_ssr_res$info_fraction      # Current information fraction

Unblinded SSR

unblinded_ssr_res <- unblinded_ssr(
  data = cut_summary,
  ratio = 1,
  lambda1_planning = 0.5,
  lambda2_planning = 0.3,
  power = 0.9,
  alpha = 0.025,
  accrual_rate = 20,
  accrual_duration = 12,
  trial_duration = 24,
  max_followup = 12,
  event_gap = 5/365.25
)
unblinded_ssr_res$n_total_unblinded
unblinded_ssr_res$lambda1_unblinded  # Observed control rate
unblinded_ssr_res$lambda2_unblinded  # Observed experimental rate

Updating group sequential bounds after SSR

# Observed information fraction
info_frac <- blinded_ssr_res$info_fraction

# Update bounds with observed information
gs_update <- gsDesign::gsDesign(
  k = 2,
  test.type = 1,
  alpha = 0.025,
  sfu = gsDesign::sfLDOF,
  timing = info_frac         # Observed information fraction
)
# Use gs_update$upper$bound for updated efficacy bounds

Completers-based interim analysis

An interim analysis triggered by the number of patients completing full follow-up.

target_completers <- 80

# Find date when 80 patients have completed max_followup
comp_date <- cut_date_for_completers(sim_data, target_completers)

# Cut data at that date
interim_data <- cut_completers(sim_data, cut_date = comp_date, event_gap = 5/365.25)

# Test
result <- mutze_test(interim_data)
z_interim <- result$z
info_interim <- 1 / result$se^2

# Final analysis
final_date <- max(sim_data$calendar_time)
final_data <- cut_data_by_date(sim_data, cut_date = final_date, event_gap = 5/365.25)
result_final <- mutze_test(final_data)
info_final <- 1 / result_final$se^2

# Adaptive bounds based on observed information fraction
gs_adaptive <- gsDesign::gsDesign(
  k = 2, test.type = 1, alpha = 0.025,
  sfu = gsDesign::sfLDOF,
  timing = info_interim / info_final
)

Seasonal event simulation

nb_sim_seasonal() simulates events with season-dependent rates.

fail_rate <- data.frame(
  treatment = rep(c("Control", "Experimental"), each = 4),
  season = rep(c("Winter", "Spring", "Summer", "Fall"), 2),
  rate = c(2.0, 0.5, 0.2, 1.5,       # Control by season
           1.4, 0.35, 0.14, 1.05),    # Experimental (30% reduction)
  dispersion = 0.5
)

sim_seasonal <- nb_sim_seasonal(
  enroll_rate = data.frame(rate = 40, duration = 0.5),
  fail_rate = fail_rate,
  dropout_rate = data.frame(
    treatment = c("Control", "Experimental"),
    rate = c(0.05, 0.05),
    duration = c(100, 100)
  ),
  max_followup = 1,
  randomization_start_date = as.Date("2024-01-01"),
  n = 200
)

# Output: one row per season-interval per patient
# Columns: id, treatment, season, enroll_time, start, end, event,
#           calendar_start, calendar_end

# Cut and analyze
cut_seasonal <- cut_data_by_date(sim_seasonal, cut_date = 1.25, event_gap = 7/365.25)

# Fit NB GLM with seasonal adjustment
fit <- MASS::glm.nb(
  events ~ treatment + season + offset(log(tte)),
  data = cut_seasonal[cut_seasonal$tte > 0, ]
)
rr <- exp(coef(fit)["treatmentExperimental"])

Verification of operating characteristics

Compare theoretical (formula-based) predictions to simulation results.

# Theoretical predictions from design
design <- sample_size_nbinom(
  lambda1 = 0.4, lambda2 = 0.3, dispersion = 0.5,
  power = 0.9, alpha = 0.025,
  accrual_rate = c(1, 2), accrual_duration = c(6, 6),
  trial_duration = 24, dropout_rate = 0.1/12,
  max_followup = 12, event_gap = 20/30.42
)

# Simulate many trials
sim_results <- sim_gs_nbinom(
  n_sims = 3600,   # SE of power estimate ~ 0.005
  enroll_rate = data.frame(
    rate = design$accrual_rate * design$n_total,
    duration = design$inputs$accrual_duration
  ),
  fail_rate = data.frame(
    treatment = c("Control", "Experimental"),
    rate = c(0.4, 0.3),
    dispersion = c(0.5, 0.5)
  ),
  dropout_rate = data.frame(
    treatment = c("Control", "Experimental"),
    rate = c(0.1/12, 0.1/12),
    duration = c(100, 100)
  ),
  max_followup = 12,
  event_gap = 20/30.42,
  analysis_times = design$inputs$trial_duration,
  n_target = design$n_total,
  design = design
)

# Compare theoretical vs. simulated
sim_final <- sim_results[sim_results$analysis == 1, ]
comparison <- data.frame(
  Metric = c("Events (Control)", "Events (Experimental)",
             "Exposure at Risk (Control)", "Power"),
  Theoretical = c(design$events_n1, design$events_n2,
                   design$exposure_at_risk_n1 * design$n1, design$power),
  Simulated = c(mean(sim_final$events_ctrl), mean(sim_final$events_exp),
                 mean(sim_final$exposure_at_risk_ctrl), mean(sim_final$z_stat > qnorm(1 - 0.025)))
)

Sample Size and Simulation for Negative Binomial Outcomes with gsDesignNB

API reference

Key functions

Sample size calculation

Group sequential design

Simulation

Data cutting and analysis timing

Statistical testing and estimation

Sample size re-estimation

Workflow patterns

Important design considerations

Code Patterns

Code Patterns for gsDesignNB

Table of Contents

Fixed sample size calculation

Basic calculation

Piecewise accrual ramp-up

With dropout and maximum follow-up

Group-specific parameters

Unequal allocation

Key output fields

Power calculation

Event gaps and piecewise accrual

Non-inferiority and super-superiority

Non-inferiority (rr0 > 1)

Super-superiority (rr0 < 1)

Group sequential design

Basic group sequential design

With custom spending times

Key output fields

Simulation of recurrent events

Basic simulation

With enrollment ramp-up and event gaps

Group sequential simulation

Basic usage

With custom cutting criteria

Data cutting and analysis timing

Cut at a calendar date

Find date when target events reached

Find date satisfying multiple criteria

Find date for completers

Wald test for treatment rate ratio

Extracting information from the test

Blinded and unblinded information

Blinded information at interim

Method of moments estimation

Information at a planning time

Sample size re-estimation

Blinded SSR (Friede & Schmidli)

Unblinded SSR

Updating group sequential bounds after SSR

Completers-based interim analysis

Seasonal event simulation

Verification of operating characteristics