Count-Outcome Equivalence: Parallel Design with 3 Arms and 3 Endpoints
In SimTOST: Sample Size Estimation for Bio-Equivalence Trials Through Simulation

knitr::opts_chunk$set(echo = TRUE, collapse = TRUE, comment = "#>",
                      warning = FALSE, message = FALSE)
library(SimTOST)

Overview

This example illustrates how to plan a parallel count-outcome equivalence study with three treatment arms and three count endpoints. The arms are a test product (TEST) and two reference products (EU_REF and US_REF). We use three clinically interpretable, illustrative endpoint names: Exacerbations, Hospitalizations, and RescueEvents (rescue-medication events). These labels represent event counts collected over the same follow-up period. The general distributional assumptions used by all SimTOST outcomes are described in the companion vignette methodological_assumptions.Rmd.

The count module supports scalar or vector-valued endpoint inputs. A joint simulation can require equivalence for k of the m endpoints and can apply Bonferroni, Mielke's weak k-out-of-m (adjust = "k"), Mielke's strong k-out-of-m (adjust = "t"), or Šidák adjustment to the endpoint-wise one-sided significance level.

The joint count kernel can also simulate correlated endpoints through the cor_mat argument. This dependence is generated with a Gaussian-copula construction [@nelsen2006].

Study assumptions

The assumed event rates are:

rate_list <- list(
  TEST = c(Exacerbations = 0.19,
           Hospitalizations = 0.13,
           RescueEvents = 0.08),
  EU_REF = c(Exacerbations = 0.20,
             Hospitalizations = 0.13,
             RescueEvents = 0.08),
  US_REF = c(Exacerbations = 0.21,
             Hospitalizations = 0.12,
             RescueEvents = 0.08)
)

comparisons <- list(
  EU_comparison = c("TEST", "EU_REF"),
  US_comparison = c("TEST", "US_REF")
)

lower_margin <- 0.80
upper_margin <- 1.25
exposure <- 5

The rate-ratio equivalence interval is 0.80 to 1.25. We use a Poisson model for the primary example and a one-sided significance level of 0.05 for each TOST component. The Poisson and negative-binomial planning framework follows published count-outcome equivalence methods [@chang_sample_2017; @zhu_sample_2017].

Power for a fixed sample size

Before estimating sample size, power can be examined for each endpoint and comparison. Here we use 100 participants per treatment arm and a small number of simulations to keep the vignette fast.

run_fixed_power <- function(comparison_name, endpoint_name) {
  arms <- comparisons[[comparison_name]]
  result <- simPower(
    n = 100,
    distribution = "pois",
    rate_list = setNames(list(
      setNames(rate_list[[arms[1]]][[endpoint_name]], endpoint_name),
      setNames(rate_list[[arms[2]]][[endpoint_name]], endpoint_name)
    ), arms),
    list_comparator = setNames(list(arms), comparison_name),
    list_lequi.tol = setNames(list(lower_margin), comparison_name),
    list_uequi.tol = setNames(list(upper_margin), comparison_name),
    exposure = exposure,
    dtype = "parallel",
    nsim = 1000,
    seed = 1234
  )
  data.frame(
    comparison = comparison_name, endpoint = endpoint_name,
    power = result$power, power_LCI = result$power_LCI,
    power_UCI = result$power_UCI
  )
}

# Each object is one transparent comparison-endpoint calculation.
fixed_power_EU_Exacerbations <- run_fixed_power("EU_comparison", "Exacerbations")
fixed_power_EU_Hospitalizations <- run_fixed_power("EU_comparison", "Hospitalizations")
fixed_power_EU_RescueEvents <- run_fixed_power("EU_comparison", "RescueEvents")
fixed_power_US_Exacerbations <- run_fixed_power("US_comparison", "Exacerbations")
fixed_power_US_Hospitalizations <- run_fixed_power("US_comparison", "Hospitalizations")
fixed_power_US_RescueEvents <- run_fixed_power("US_comparison", "RescueEvents")

fixed_power <- rbind(
  fixed_power_EU_Exacerbations,
  fixed_power_EU_Hospitalizations,
  fixed_power_EU_RescueEvents,
  fixed_power_US_Exacerbations,
  fixed_power_US_Hospitalizations,
  fixed_power_US_RescueEvents
)

fixed_power

Sample-size estimation for all comparisons and endpoints

We now estimate the smallest sample size per arm that reaches 80% power for each individual comparison and endpoint.

run_sample_size <- function(comparison_name, endpoint_name) {
  arms <- comparisons[[comparison_name]]
  result <- sampleSize(
    power = 0.80,
    distribution = "pois",
    rate_list = setNames(list(
      setNames(rate_list[[arms[1]]][[endpoint_name]], endpoint_name),
      setNames(rate_list[[arms[2]]][[endpoint_name]], endpoint_name)
    ), arms),
    list_comparator = setNames(list(arms), comparison_name),
    list_lequi.tol = setNames(list(lower_margin), comparison_name),
    list_uequi.tol = setNames(list(upper_margin), comparison_name),
    exposure = exposure,
    dtype = "parallel",
    nsim = 1000,
    seed = 1234,
    lower = 10,
    upper = 2000
  )
  data.frame(
    comparison = comparison_name, endpoint = endpoint_name,
    n_per_arm = result$n_per_arm, n_total_for_pair = result$n_total,
    achieved_power = result$power
  )
}

# Again, keep the six calculations as named objects so each result can be
# inspected independently before combining them into one table.
sample_size_EU_Exacerbations <- run_sample_size("EU_comparison", "Exacerbations")
sample_size_EU_Hospitalizations <- run_sample_size("EU_comparison", "Hospitalizations")
sample_size_EU_RescueEvents <- run_sample_size("EU_comparison", "RescueEvents")
sample_size_US_Exacerbations <- run_sample_size("US_comparison", "Exacerbations")
sample_size_US_Hospitalizations <- run_sample_size("US_comparison", "Hospitalizations")
sample_size_US_RescueEvents <- run_sample_size("US_comparison", "RescueEvents")

sample_size_results <- rbind(
  sample_size_EU_Exacerbations,
  sample_size_EU_Hospitalizations,
  sample_size_EU_RescueEvents,
  sample_size_US_Exacerbations,
  sample_size_US_Hospitalizations,
  sample_size_US_RescueEvents
)

sample_size_results

For a simple conservative planning rule, select the maximum required sample size per arm across all six comparison-endpoint combinations:

required_per_arm <- max(sample_size_results$n_per_arm)
required_total <- 3 * required_per_arm

c(required_per_arm = required_per_arm, required_total = required_total)

The total is multiplied by three because the study has three treatment arms. This rule ensures that each individual comparison and endpoint has at least the target simulated power under its own assumptions. It is not a substitute for a multiplicity-adjusted joint power calculation.

Joint sample-size estimation with `k = 3`

The joint sample-size function receives all three arms, both comparison families, and all three endpoints at once. Here, k = 3 requires all three endpoints to demonstrate equivalence for every comparison. Because k = m = 3, all three endpoints must pass; no endpoint can be selected or omitted. An endpoint-wise adjustment is needed when the rule allows the study to pass by choosing only some endpoints, such as k = 2 of m = 3, because there are then several possible successful endpoint subsets. With k = m, there is only one acceptable outcome: all three endpoints pass. The joint success rule therefore already defines the required criterion, and no endpoint-wise Bonferroni adjustment is required. The returned sample size is based on one joint simulated success criterion rather than the maximum of separate searches.

endpoint_corr <- matrix(c(
  1.0, 0.40, 0.25,
  0.40, 1.0, 0.35,
  0.25, 0.35, 1.0
), nrow = 3, byrow = TRUE)

joint_result <- sampleSize(
  power = 0.80,
  distribution = "pois",
  rate_list = rate_list,
  list_comparator = comparisons,
  list_lequi.tol = list(
    EU_comparison = rep(lower_margin, 3),
    US_comparison = rep(lower_margin, 3)
  ),
  list_uequi.tol = list(
    EU_comparison = rep(upper_margin, 3),
    US_comparison = rep(upper_margin, 3)
  ),
  exposure = rep(exposure, 3),
  cor_mat = endpoint_corr,
  dtype = "parallel",
  nsim = 500,
  seed = 1234,
  lower = 10,
  upper = 3000,
  k = 3,
  adjust = "none"
)

joint_sample_size <- data.frame(
  n_per_arm = joint_result$n_per_arm,
  n_total = joint_result$n_total,
  achieved_power = joint_result$power,
  k = joint_result$k,
  adjustment = joint_result$adjust
)

joint_sample_size

For a three-arm allocation, the reported total is three times the selected number per arm. Both comparison families share the simulated test-arm counts, and the endpoint correlation is generated through a Gaussian-copula count model [@nelsen2006]. This is therefore a joint count simulation rather than a maximum of separate comparison-specific searches.

Why joint planning is important

The separate calculation above targets 80% power for each comparison-endpoint combination. That does not mean that the complete trial has 80% probability of passing all six requirements. For example, if two requirements each have 80% power and are independent, their joint success probability is only $0.80 \times 0.80 = 0.64$.

The following comparison evaluates the separate result under the actual joint criterion. The joint calculations require all three endpoints for both comparison families (k = 3). Since k = m, no endpoint-wise adjustment is needed here; Bonferroni could be added only as a conservative sensitivity analysis. The independent-endpoint scenario uses an identity correlation matrix; the correlated scenario uses the matrix specified above.

joint_power_at_separate <- simPower(
  n = required_per_arm,
  distribution = "pois",
  rate_list = rate_list,
  list_comparator = comparisons,
  list_lequi.tol = list(
    EU_comparison = rep(lower_margin, 3),
    US_comparison = rep(lower_margin, 3)
  ),
  list_uequi.tol = list(
    EU_comparison = rep(upper_margin, 3),
    US_comparison = rep(upper_margin, 3)
  ),
  exposure = rep(exposure, 3),
  cor_mat = endpoint_corr,
  dtype = "parallel",
  nsim = 1000,
  seed = 1234,
  k = 3,
  adjust = "none"
)

joint_independent_result <- update(
  joint_result,
  cor_mat = diag(3),
  nsim = 500,
  seed = 1234
)

joint_comparison <- data.frame(
  approach = c(
    "Separate searches; joint power evaluated afterward",
    "Joint search; independent endpoints",
    "Joint search; correlated endpoints"
  ),
  n_per_arm = c(
    required_per_arm,
    joint_independent_result$n_per_arm,
    joint_result$n_per_arm
  ),
  joint_power = c(
    joint_power_at_separate$power,
    joint_independent_result$power,
    joint_result$power
  )
)

joint_comparison

The separate-search sample size is smaller because it guarantees only the individual 80% targets. Its joint power can therefore be substantially below 80%, even without an endpoint-wise adjustment, whereas the joint search targets the requested 80% trial-level power directly. This is the main advantage of joint planning: it answers the question that matters for the confirmatory trial, namely the probability that all required comparisons and endpoints succeed in the same simulated study. Positive endpoint correlation can reduce the joint sample size because endpoint successes tend to occur together, but it does not remove the need for a joint calculation.

Negative-binomial sensitivity analysis

If prior evidence suggests overdispersion, repeat the final joint Poisson sample-size calculation with the negative-binomial model. All study settings are kept the same as in joint_result; update() changes only the distribution and dispersion. The dispersion parameter controls the amount of overdispersion; larger values imply greater variability.

nb_dispersion <- 0.10

joint_negative_binomial_result <- update(
  joint_result,
  distribution = "nbinom",
  dispersion = nb_dispersion,
  nsim = 500,
  seed = 1234
)

nb_sample_size <- data.frame(
  n_per_arm = joint_negative_binomial_result$n_per_arm,
  n_total = joint_negative_binomial_result$n_total,
  achieved_power = joint_negative_binomial_result$power,
  k = joint_negative_binomial_result$k,
  adjustment = joint_negative_binomial_result$adjust,
  dispersion = nb_dispersion
)

nb_sample_size