count_parallel_3A3E.RmdThis 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 strong
k-out-of-m (adjust = "t"), or Šidák adjustment to the
endpoint-wise one-sided significance level. For a
k-of-m decision, the t-adjustment uses
alpha / (m - k + 1) [@mielke_sample_2018]. The joint count kernel
can also simulate correlated endpoints through the cor_mat
argument.
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
# Five units of follow-up provide a realistic event-count scale for this
# illustrative example and allow the target power to be reached efficiently.
exposure <- 5The 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.
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.
fixed_power <- do.call(rbind, lapply(names(comparisons), function(comp) {
arms <- comparisons[[comp]]
do.call(rbind, lapply(names(rate_list[[arms[1]]]), function(endpoint) {
result <- simPower(
n = 100,
distribution = "pois",
rate_list = setNames(list(
setNames(rate_list[[arms[1]]][[endpoint]], endpoint),
setNames(rate_list[[arms[2]]][[endpoint]], endpoint)
), arms),
list_comparator = setNames(list(arms), comp),
list_lequi.tol = setNames(list(lower_margin), comp),
list_uequi.tol = setNames(list(upper_margin), comp),
exposure = exposure,
dtype = "parallel",
nsim = 1000,
seed = 1234
)
data.frame(comparison = comp, endpoint = endpoint,
power = result$power,
power_LCI = result$power_LCI,
power_UCI = result$power_UCI)
}))
}))
fixed_power
#> comparison endpoint power power_LCI power_UCI
#> 1 EU_comparison Exacerbations 0.000 0.0000000000 0.004770729
#> 2 EU_comparison Hospitalizations 0.000 0.0000000000 0.004770729
#> 3 EU_comparison RescueEvents 0.000 0.0000000000 0.004770729
#> 4 US_comparison Exacerbations 0.002 0.0003464932 0.008032515
#> 5 US_comparison Hospitalizations 0.000 0.0000000000 0.004770729
#> 6 US_comparison RescueEvents 0.000 0.0000000000 0.004770729We now estimate the smallest sample size per arm that reaches 80% power for each individual comparison and endpoint.
sample_size_results <- do.call(rbind, lapply(names(comparisons), function(comp) {
arms <- comparisons[[comp]]
do.call(rbind, lapply(names(rate_list[[arms[1]]]), function(endpoint) {
result <- sampleSize(
power = 0.80,
distribution = "pois",
rate_list = setNames(list(
setNames(rate_list[[arms[1]]][[endpoint]], endpoint),
setNames(rate_list[[arms[2]]][[endpoint]], endpoint)
), arms),
list_comparator = setNames(list(arms), comp),
list_lequi.tol = setNames(list(lower_margin), comp),
list_uequi.tol = setNames(list(upper_margin), comp),
exposure = exposure,
dtype = "parallel",
nsim = 1000,
seed = 1234,
lower = 10,
upper = 2000
)
data.frame(comparison = comp, endpoint = endpoint,
n_per_arm = result$n_per_arm,
n_total_for_pair = result$n_total,
achieved_power = result$power)
}))
}))
sample_size_results
#> comparison endpoint n_per_arm n_total_for_pair achieved_power
#> 1 EU_comparison Exacerbations 428 856 0.806
#> 2 EU_comparison Hospitalizations 522 1044 0.810
#> 3 EU_comparison RescueEvents 846 1692 0.800
#> 4 US_comparison Exacerbations 793 1586 0.815
#> 5 US_comparison Hospitalizations 903 1806 0.803
#> 6 US_comparison RescueEvents 846 1692 0.800
stopifnot(all(sample_size_results$achieved_power >= 0.80))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)
#> required_per_arm required_total
#> 903 2709The 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.
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, while the Bonferroni option adjusts
across all six comparison-endpoint tests. The returned sample size is
therefore 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 = "bonferroni"
)
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
#> n_per_arm n_total achieved_power k adjustment
#> 1 1881 5643 0.8 3 bonferroni
stopifnot(all(joint_sample_size$achieved_power >= 0.80))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. This is therefore a joint count simulation rather than a maximum of separate comparison-specific searches.
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 .
The following comparison evaluates the separate result under the
actual joint criterion. The joint calculations use Bonferroni adjustment
across the six comparison-endpoint tests and require all three endpoints
for both comparison families (k = 3). 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 = "bonferroni"
)
joint_independent_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 = diag(3),
dtype = "parallel",
nsim = 500,
seed = 1234,
lower = 10,
upper = 3000,
k = 3,
adjust = "bonferroni"
)
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
#> approach n_per_arm joint_power
#> 1 Separate searches; joint power evaluated afterward 903 0.159
#> 2 Joint search; independent endpoints 1916 0.800
#> 3 Joint search; correlated endpoints 1881 0.800The separate-search sample size is smaller because it guarantees only the individual 80% targets. In this example, its joint power is approximately 5% after the six-test Bonferroni adjustment, whereas the joint searches achieve the requested 80% trial-level power. 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.
If prior evidence suggests overdispersion, repeat the calculations
with the negative-binomial model. The dispersion parameter
controls the amount of overdispersion; larger values imply greater
variability.
nb_sensitivity <- simPower(
n = required_per_arm,
distribution = "nbinom",
rate_list = list(
TEST = rate_list$TEST[["Exacerbations"]],
EU_REF = rate_list$EU_REF[["Exacerbations"]]
),
list_comparator = list(TEST_vs_EU = c("TEST", "EU_REF")),
list_lequi.tol = list(TEST_vs_EU = lower_margin),
list_uequi.tol = list(TEST_vs_EU = upper_margin),
exposure = exposure,
dtype = "parallel",
dispersion = 0.10,
nsim = 300,
seed = 1234
)
nb_sensitivity
#> Fixed-sample-size power
#> Distribution: nbinom
#> Sample size: 903
#> Power: 0.2633 [0.2152, 0.3177]This vignette demonstrates both a conservative and a joint endpoint
workflow for a three-arm, three-endpoint count study using the current
API. In the separate-search workflow, endpoint and comparison
calculations are performed independently, whereas the joint workflow
evaluates the required success event in the same simulated trial. The
distributional and dependence assumptions are described in the companion
methodological_assumptions.Rmd
vignette.
For a confirmatory three-arm/three-endpoint trial, the statistical analysis plan should specify the familywise error strategy and joint success criterion.