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 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 (Nelsen 2006).
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 <- 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. The Poisson and negative-binomial planning framework follows published count-outcome equivalence methods (Chang et al. 2017; Zhu 2017).
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
#> 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.
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
#> 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.800For 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. 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
#> n_per_arm n_total achieved_power k adjustment
#> 1 1443 4329 0.804 3 noneFor 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 (Nelsen 2006). 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 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
#> approach n_per_arm joint_power
#> 1 Separate searches; joint power evaluated afterward 903 0.437
#> 2 Joint search; independent endpoints 1477 0.802
#> 3 Joint search; correlated endpoints 1443 0.804The 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.
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
#> n_per_arm n_total achieved_power k adjustment dispersion
#> 1 1482 4446 0.8 3 none 0.1