sampleSize: Sample Size Calculation for Bioequivalence and Multi-Endpoint...

In SimTOST: Sample Size Estimation for Bio-Equivalence Trials Through Simulation

Sample Size Calculation for Bioequivalence and Multi-Endpoint Studies

Description

Computes the required sample size to achieve a target power in studies with multiple endpoints and treatment arms. The function employs modified root-finding algorithms to estimate sample size while accounting for correlation structures, variance assumptions, and equivalence bounds across endpoints. It is particularly useful for bioequivalence trials and multi-arm studies with complex endpoint structures.

Usage

sampleSize(
  distribution = c("norm", "lnorm", "pois", "nbinom"),
  mu_list = NULL,
  varcov_list = NA,
  sigma_list = NA,
  cor_mat = NA,
  sigmaB = NA,
  rate_list = NULL,
  exposure = 1,
  dispersion = 0.1,
  Eper = c(0, 0),
  Eco = c(0, 0),
  rho = 0,
  TAR = rep(1, length(mu_list)),
  arm_names = NA,
  ynames_list = NA,
  type_y = NA,
  list_comparator = NA,
  list_y_comparator = NA,
  power = 0.8,
  alpha = 0.05,
  lequi.tol = NA,
  uequi.tol = NA,
  list_lequi.tol = NA,
  list_uequi.tol = NA,
  dtype = "parallel",
  ctype = "ROM",
  vareq = TRUE,
  k = NA,
  adjust = "no",
  dropout = NA,
  nsim = 5000,
  seed = 1234,
  ncores = 1,
  optimization_method = "fast",
  lower = 2,
  upper = 500,
  step.power = 6,
  step.up = TRUE,
  pos.side = FALSE,
  maxiter = 1000,
  verbose = FALSE,
  keep_sim_data = FALSE,
  .warn_redundant_bon = TRUE
)

Arguments

distribution	Outcome distribution. Choose the R distribution names "norm", "lnorm", "pois", or "nbinom". Matching is case-insensitive. The longer labels such as "normal", "lognormal", and "poisson" are accepted for compatibility; results store the R-style labels above.
mu_list	Named list of arithmetic means per treatment arm. Each element is a vector representing expected outcomes for all endpoints in that arm.
varcov_list	List of variance-covariance matrices, where each element corresponds to a comparator. Each matrix has dimensions: number of endpoints × number of endpoints.
sigma_list	List of standard deviation vectors, where each element corresponds to a comparator and contains one standard deviation per endpoint.
cor_mat	Matrix specifying the correlation structure between endpoints. For continuous outcomes it is used with sigma_list to calculate varcov_list; for count outcomes it is used by the joint count engine.
sigmaB	Numeric. Between-subject standard deviation parameter for the continuous 2×2 design; for count outcomes it is the log-rate standard deviation in the 2×2 kernel.
rate_list	Named list of equal-length endpoint-rate vectors, one per count-outcome arm.
exposure	Exposure per subject for count outcomes. Supply a scalar or endpoint vector shared by arms, or a named list of arm-specific values.
dispersion	Positive negative-binomial dispersion parameter. The per-subject negative-binomial size is 1 / dispersion; parallel-arm totals use size n / dispersion. Supply a scalar or endpoint vector shared by arms, or a named list of arm-specific values.
Eper	Optional numeric vector of length 2 specifying period effects. For count outcomes these are log-rate effects applied to periods 1 and 2.
Eco	Optional numeric vector of length 2 specifying carry-over effects in the order reference carry-over and treatment carry-over. For count outcomes these are log-rate effects in period 2.
rho	Numeric. Correlation parameter applied uniformly across all endpoint pairs. Used with sigma_list to compute varcov_list when cor_mat or varcov_list are not provided.
TAR	Numeric vector specifying treatment allocation rates per arm. The order must match arm_names. Defaults to equal allocation across arms if not provided.
arm_names	Optional character vector of treatment names. If not supplied, names are derived from mu_list.
ynames_list	Optional list of vectors specifying endpoint names per arm. If names are missing, arbitrary names are assigned based on order.
type_y	Integer vector indicating endpoint types: 1 for co-primary endpoints, 2 for secondary endpoints.
list_comparator	List of comparators. Each element must be a vector of length 2 in the form c(test, reference). The first arm is treated as the treatment arm and the second arm as the reference arm throughout the simulation and estimand calculations.
list_y_comparator	List of endpoint sets per comparator. Each element is a vector containing endpoint names to compare. If not provided, all endpoints common to both comparator arms are used. For count outcomes, the selected endpoints define the count multiplicity and effective k; joint count analyses require the same endpoint set for every comparison.
power	Numeric. Target power (default = 0.8).
alpha	Numeric. Significance level (default = 0.05).
lequi.tol	Numeric. Lower equivalence bounds (e.g., -0.5) applied uniformly across all endpoints and comparators.
uequi.tol	Numeric. Upper equivalence bounds (e.g., 0.5) applied uniformly across all endpoints and comparators.
list_lequi.tol	List of numeric vectors specifying lower equivalence bounds per comparator.
list_uequi.tol	List of numeric vectors specifying upper equivalence bounds per comparator.
dtype	Character. Trial design: "parallel" (default) for parallel-group or "2x2" for crossover (only for 2-arm studies).
ctype	Character. Continuous-outcome test type: "DOM" (Difference of Means) or "ROM" (Ratio of Means). For Poisson and negative-binomial outcomes, "RR" (event-rate ratio) is used; an unavailable value triggers a warning and the applicable default is used.
vareq	Logical. Assumes equal variances across arms if TRUE (default = FALSE).
k	Integer vector. Minimum number of successful endpoints required for global bioequivalence per comparator. Defaults to all endpoints per comparator.
adjust	Character. Alpha adjustment method: "k" (K-fold), "bon" (Bonferroni across all selected endpoints), "sid" (Sidak), "t" (Mielke's strong \(k\)-out-of-\(m\) adjustment using alpha / (m - k + 1); legacy "pc" aliases are accepted), "no" (none), or "seq" (sequential).
dropout	Numeric vector specifying dropout proportion per arm.
nsim	Integer. Number of simulated studies (default = 5000).
seed	Integer. Seed for reproducibility.
ncores	Integer. Number of processing cores for parallel computation. Defaults to 1. Set to NA for automatic detection (ncores - 1). For count outcomes, Monte Carlo trials are split into independent seeded chunks and evaluated by the count C++ kernel on each worker.
optimization_method	Character. Sample size optimization method: "fast" (default, root-finding algorithm) or "step-by-step".
lower	Integer. Minimum sample size for search range (default = 2).
upper	Integer. Maximum sample size for the search range (default = 500). For count outcomes, this is the maximum number of subjects per arm; the plotted and returned total sample size is this value multiplied by the number of trial arms.
step.power	Numeric. Initial step size for sample size search, defined as 2^step.power. Used when optimization_method = "fast".
step.up	Logical. If TRUE (default), search increments upward from lower; if FALSE, decrements downward from upper. Used when optimization_method = "fast".
pos.side	Logical. If TRUE, finds the smallest integer i closest to the root such that f(i) > 0. Used when optimization_method = "fast".
maxiter	Integer. Maximum iterations allowed for sample size estimation (default = 1000). Used when optimization_method = "fast".
verbose	Logical. If TRUE, prints progress and messages during execution (default = FALSE).
keep_sim_data	Logical. If TRUE, retain model-scale observations for every simulated trial in sim_data for distribution diagnostics. Defaults to FALSE because retained data can be large.
.warn_redundant_bon	Logical. If TRUE, warn when a requested multiplicity adjustment is redundant or uncalibrated for the selected endpoint decision.

Details

The common planning arguments are power, alpha, list_comparator, list_lequi.tol, list_uequi.tol, k, adjust, dtype, dropout, nsim, seed, lower, and upper. Use the following distribution-specific arguments in addition to those common arguments:

Normal and Log Normal

Supply mu_list, sigma_list or varcov_list; use cor_mat or rho for endpoint dependence. The ctype argument selects DOM or ROM testing.

Poisson and Negative Binomial

Supply rate_list, list_comparator, and comparator-specific list_lequi.tol and list_uequi.tol. Use exposure and, for negative-binomial outcomes, dispersion; both may be scalar, endpoint-specific, or named arm-specific lists. Continuous-outcome arguments are ignored.

For count outcomes, optimization_method = "fast" brackets the first sample size whose simulated power reaches the target and refines the bracket by integer bisection. The "step-by-step" option remains available when a complete candidate-by-candidate power table is preferred. The fast method assumes the usual approximately monotone power curve and uses the same seed at each candidate to reduce simulation noise. The effective endpoint count is comparator-specific: when list_y_comparator is omitted, only endpoints present in both arms are tested; when it is supplied, only the listed endpoints are tested. k is validated against that comparator-specific count and oversized values are capped with a warning. Formal endpoint-wise adjustment is unnecessary when all selected endpoints are required (k = m), although requested Bonferroni or Sidak adjustment remains available with a warning. For k < m, adjust = "no" is explicitly reported as an uncalibrated choice. For a k-of-m decision, adjust = "t" applies Mielke's strong \(k\)-out-of-\(m\) calibration alpha / (m - k + 1). The legacy adjust = "pc" label is accepted as an alias. For continuous and count outcomes, type_y is used with adjust = "seq"; named endpoint vectors are aligned to the selected comparator endpoints. Count analyses use the same primary-gate and secondary-family decision rule as the continuous kernels. The unified function returns primary class simss for all outcome distributions. Count results retain countss as a secondary compatibility class. Use summary() and plot() to inspect the result.

Value

A list containing:

response

Array summarizing simulation results, including estimated sample sizes, achieved power, and confidence intervals.

table.iter

Data frame showing estimated sample sizes and calculated power at each iteration. For count outcomes, one row is retained for every evaluated candidate.

table.test

Data frame containing test results for all simulated trials. For count outcomes, this contains complete-trial, comparator, and endpoint decision indicators for each simulated trial and candidate; the count kernel returns aggregate decision counts rather than raw endpoint-level test statistics.

param.u

Original input parameters.

param

Final adjusted parameters used in sample size calculation.

param.d

Trial design parameters used in the simulation.

sim_data

Optional long-format simulated observations, returned when keep_sim_data = TRUE.

References

Schuirmann, D. J. (1987). A comparison of the Two One-Sided Tests procedure and the Power approach for assessing the equivalence of average bioavailability. Journal of Pharmacokinetics and Biopharmaceutics, 15(6), 657-680. doi:10.1007/BF01068419

Mielke, J., Jones, B., Jilma, B., & König, F. (2018). Sample size for multiple hypothesis testing in biosimilar development. Statistics in Biopharmaceutical Research, 10(1), 39-49. doi:10.1080/19466315.2017.1371071

Berger, R. L., & Hsu, J. C. (1996). Bioequivalence trials, intersection-union tests, and equivalence confidence sets. Statistical Science, 283-302.

Sozu, T., Sugimoto, T., Hamasaki, T., & Evans, S. R. (2015). "Sample Size Determination in Clinical Trials with Multiple Endpoints." SpringerBriefs in Statistics. doi:10.1007/978-3-319-22005-5

Examples

mu_list <- list(SB2 = c(AUCinf = 38703, AUClast = 36862, Cmax = 127.0),
                EUREF = c(AUCinf = 39360, AUClast = 37022, Cmax = 126.2),
                USREF = c(AUCinf = 39270, AUClast = 37368, Cmax = 129.2))

sigma_list <- list(SB2 = c(AUCinf = 11114, AUClast = 9133, Cmax = 16.9),
                   EUREF = c(AUCinf = 12332, AUClast = 9398, Cmax = 17.9),
                   USREF = c(AUCinf = 10064, AUClast = 8332, Cmax = 18.8))

# Equivalent boundaries
lequi.tol <- c(AUCinf = 0.8, AUClast = 0.8, Cmax = 0.8)
uequi.tol <- c(AUCinf = 1.25, AUClast = 1.25, Cmax = 1.25)

# Arms to be compared
list_comparator <- list(EMA = c("SB2", "EUREF"),
                        FDA = c("SB2", "USREF"))

# Endpoints to be compared
list_y_comparator <- list(EMA = c("AUCinf", "Cmax"),
                          FDA = c("AUClast", "Cmax"))

# Equivalence boundaries for each comparison
lequi_lower <- c(AUCinf = 0.80, AUClast = 0.80, Cmax = 0.80)
lequi_upper <- c(AUCinf = 1.25, AUClast = 1.25, Cmax = 1.25)

# Run the simulation
sampleSize(power = 0.9, alpha = 0.05, mu_list = mu_list,
           sigma_list = sigma_list, list_comparator = list_comparator,
           list_y_comparator = list_y_comparator,
           list_lequi.tol = list("EMA" = lequi_lower, "FDA" = lequi_lower),
           list_uequi.tol = list("EMA" = lequi_upper, "FDA" = lequi_upper),
           adjust = "no", dtype = "parallel", ctype = "ROM", vareq = FALSE,
           distribution = "lnorm", ncores = 1, nsim = 50, seed = 1234)

# The same entry point for a two-arm Poisson count-rate calculation:
sampleSize(power = 0.80, distribution = "Poisson",
           rate_list = list(TEST = 0.21, REF = 0.20),
           list_comparator = list(TEST_vs_REF = c("TEST", "REF")),
           list_lequi.tol = list(TEST_vs_REF = 0.80),
           list_uequi.tol = list(TEST_vs_REF = 1.25),
           exposure = 10, lower = 20, upper = 500,
           nsim = 100, seed = 1234)

SimTOST documentation built on Oct. 9, 2026, 5:07 p.m.