Methodological Assumptions in SimTOST

Purpose and scope

SimTOST estimates sample size by repeatedly simulating complete studies and applying the planned equivalence procedure to each simulated study. The same framework is available through sampleSize() for planning and simPower() for evaluating power at a fixed sample size. The distribution argument selects the outcome model using one of the documented choices: "norm", "lnorm", "pois", or "nbinom".

This vignette describes the assumptions that should be considered when using these models. It is methodological guidance, not a substitute for a prespecified statistical analysis plan.

Assumptions shared by the simulation framework

The following assumptions apply to all outcome distributions unless explicitly changed by the design inputs:

• treatment arms, endpoint definitions, exposure periods, and equivalence margins are specified before the simulation;
• the supplied means, rates, variances, correlations, and dispersion values represent the intended study population;
• observations are generated according to the selected marginal outcome model;
• dropout is represented through the supplied arm- or sequence-specific dropout proportions and is non-informative with respect to the outcome;
• the selected design (dtype = "parallel" or dtype = "2x2") matches the planned trial;
• the number of simulations is sufficiently large for the desired Monte Carlo precision; and
• the equivalence margins, endpoint rule (k), and multiplicity adjustment are the actual decision rules that will be used in the trial.

For a fixed number of simulations, the reported power is a Monte Carlo estimate. Its confidence interval describes simulation uncertainty; it does not describe uncertainty in the assumed rates, means, variances, or correlations.

Continuous outcomes

Normal outcomes

For distribution = "norm", the simulated endpoint vector is multivariate normal within each treatment arm. The supplied mu_list contains arm- and endpoint-specific means. Variability is supplied through varcov_list, or is constructed from sigma_list and cor_mat (or the common correlation rho).

The Normal model assumes:

• continuous outcomes on the analysis scale;
• approximately symmetric marginal distributions;
• a common covariance structure within each arm for the simulated study;
• a valid, positive-definite covariance matrix; and
• no unmodelled skewness or heavy tails that would materially alter the TOST statistic.

The DOM test uses additive equivalence margins. The ROM test is appropriate when the scientific question concerns a ratio of means and the outcome scale supports that interpretation.

Log-normal outcomes

For distribution = "lnorm", the simulated outcomes are positive and right-skewed. The ratio-of-means procedure is applied on the scale specified by the package implementation and the supplied means and standard deviations must be interpreted consistently with that implementation.

The Log Normal model assumes:

• strictly positive outcomes;
• approximately normal log-transformed outcomes;
• endpoint dependence represented through the supplied covariance or correlation inputs; and
• a ratio-based equivalence question rather than an additive one.

The Log Normal model should not be used for outcomes with structural zeros or negative values without an explicit transformation and a corresponding reconsideration of the estimand.

distribution = "lnorm" with ctype = "DOM" is rejected because the current implementation converts arithmetic means and covariances to the log scale and applies a DOM test. This has a ratio interpretation on the original scale, not an additive difference interpretation. Use distribution = "lnorm", ctype = "ROM" for an arithmetic mean ratio, or use the normal distribution for an additive mean difference.

Count outcomes

For count outcomes, rate_list supplies the event rate per unit exposure. For a parallel arm $a$ and endpoint $j$, the expected aggregate count is

$$\mathrm{E}(Y_{aj}) = n_a e_{aj}\lambda_{aj},$$

where $n_a$ is the number of participants, $e_{aj}$ is exposure, and $\lambda_{aj}$ is the event rate. Exposure can be scalar, endpoint-specific, or arm-specific. Count equivalence is assessed through a rate ratio and log-rate-ratio TOST.

Poisson outcomes

For distribution = "pois", each marginal count follows a Poisson model:

$$Y_{aj} \sim \operatorname{Poisson}(n_a e_{aj}\lambda_{aj}).$$

The Poisson assumption implies that the variance equals the mean. This is appropriate only when additional heterogeneity, clustering, and exposure variation are negligible or have already been incorporated into the model. Overdispersion caused by unobserved subject heterogeneity or recurrent-event dependence can make a Poisson analysis anticonservative.

Negative-binomial outcomes

For distribution = "nbinom", let $\phi_{aj} > 0$ denote the per-subject dispersion parameter. A subject-level count with mean

$$\mu^{(\mathrm{subj})}_{aj} = e_{aj}\lambda_{aj}$$

has variance

$$\operatorname{Var}(Y^{(\mathrm{subj})}_{aj}) = \mu^{(\mathrm{subj})}_{aj} + \phi_{aj}\left(\mu^{(\mathrm{subj})}_{aj}\right)^2,$$

which is the negative-binomial parameterization with $\operatorname{size}=1/\phi_{aj}$. In a parallel design, the aggregate count over $n_a$ independent subjects has mean

$$\mu_{aj} = n_a e_{aj}\lambda_{aj}$$

and the implementation uses $\operatorname{size}=n_a/\phi_{aj}$. Equivalently, its variance is

$$\operatorname{Var}(Y_{aj}) = n_a\left\{\mu^{(\mathrm{subj})}_{aj} + \phi_{aj}\left(\mu^{(\mathrm{subj})}_{aj}\right)^2\right\}.$$

Thus, larger values of dispersion produce more variation beyond the Poisson variance, while values closer to zero give behavior closer to the Poisson model. This parallel-design scaling of the size parameter is important: using $1/\phi_{aj}$ for an aggregate count would incorrectly make the overdispersion increase with sample size. The dispersion parameter is positive and should preferably be based on historical data, pilot data, or a clinically justified sensitivity range.

The Negative Binomial model assumes that the selected mean-dispersion relationship adequately represents overdispersion in the planned study.

Dependence between endpoints

Continuous outcomes

For continuous outcomes, endpoint dependence is represented directly through the covariance matrices in varcov_list, or through sigma_list together with cor_mat or rho. The covariance matrices must be compatible with the endpoint means and standard deviations and must be positive definite.

Count outcomes and the Gaussian copula

For joint count simulations, cor_mat is interpreted as the correlation matrix of latent Gaussian variables. For each arm, SimTOST:

• draws a multivariate standard-normal latent vector with correlation cor_mat;
• transforms each latent value to a uniform value using the standard normal CDF; and
• transforms each uniform value through the Poisson or negative-binomial inverse CDF for the specified arm and endpoint.

This is a Gaussian-copula construction. It preserves the selected marginal count distributions while inducing dependence between endpoints. The entries of cor_mat are not Pearson correlations of the observed counts. Observed count correlations also depend on rates, exposure, and dispersion, and discreteness means that the raw-count correlation need not equal the supplied latent correlation.

What cor_mat means in practice

For endpoints $j=1,\ldots,m$, the user-supplied matrix $\boldsymbol{R}=\texttt{cor\_mat}$ is the correlation matrix of a latent standard-normal vector using a Gaussian-copula construction (Nelsen 2006):

$$\boldsymbol{Z}_a \sim N_m(\boldsymbol{0},\boldsymbol{R})$$

for arm $a$. Each component is then transformed using the probability integral transform and the inverse marginal distribution function:

$$U_{aj}=\Phi(Z_{aj}), \qquad Y_{aj}=F^{-1}_{aj}(U_{aj}),$$

where $F_{aj}$ is the specified Poisson or negative-binomial marginal distribution. Thus, cor_mat[1, 2] = 0.8 means that endpoints 1 and 2 have latent Gaussian correlation 0.8 before they are transformed into counts. It does not mean that their observed event counts will have Pearson correlation 0.8. This Gaussian-copula construction is a standard way to generate dependent non-normal outcomes.

The matrix has the following interpretation:

• diagonal elements must be 1;
• off-diagonal elements describe dependence between endpoints within the same arm or participant-level simulation unit;
• the matrix does not describe dependence between test and reference arms;
• it does not describe dependence between treatment comparisons; and
• it does not replace the subject-level pairing used by the 2x2 crossover analysis.

An identity matrix, diag(m), gives independent latent endpoint simulations. Positive correlations generally increase the probability that endpoint tests pass or fail together, which can materially change the probability of meeting the k-endpoint rule. The matrix should therefore be based on historical or pilot information, or varied in sensitivity analyses. It must be symmetric, positive definite, and have unit diagonal.

The latent vectors are generated independently between arms. In a three-arm study, however, the same simulated test-arm outcomes are used in the test-versus-reference comparisons. Thus, joint power reflects both endpoint dependence and the shared test arm. An identity matrix, diag(m), represents latent endpoint independence. In a 2x2 crossover, the same copula is applied to endpoint-specific subject effects and to the period-specific count innovations within a participant. This preserves the intended endpoint dependence while retaining the within-participant pairing.

Design-specific assumptions

Parallel design

For a parallel design, observations are generated independently between participants and treatment arms, apart from endpoint dependence specified by the model. Allocation and dropout determine the number of analyzable participants in each arm.

2x2 crossover design

Normal and log-normal outcomes

For continuous outcomes, dtype = "2x2" represents a balanced two-sequence, two-period crossover. The two treatment sequences are reference–test (RT) and test–reference (TR), and n is interpreted as the number of subjects per sequence before dropout. On the analysis scale, the data-generating model can be written for subject i, period p, and endpoint j as follows. The two-period, two-sequence design and log-scale bioequivalence analysis are consistent with regulatory bioequivalence guidance. The subject, period, treatment, sequence, and carry-over components follow the standard AB/BA formulation described by Chow and Liu (Chow and Liu 2008).

$$Y_{ipj} = \mu_{a(i,p),j} + E_p + C_{a(i,p),s(i)} + b_i + \varepsilon_{ipj},$$

Here, $Y_{ipj}$ is the observed outcome for subject $i$, period $p$, and endpoint $j$. The other terms are:

• $\mu_{a(i,p),j}$ is the mean for endpoint $j$ under the treatment received in period $p$, where $a(i,p)$ identifies the treatment assigned to subject $i$ in that period;
• $E_p$ is the fixed effect of period $p$, common to the subjects observed in that period;
• $C_{a(i,p),s(i)}$ is the carry-over contribution associated with the treatment received in the preceding period and the sequence $s(i)$. It is zero in period 1 because there is no preceding treatment;
• $b_i$ is a subject-specific random effect shared by both periods; and
• $\varepsilon_{ipj}$ is the residual error. The collection of residuals across periods and endpoints has within-subject covariance matrix SigmaW.

In the package, the model components are supplied through the following arguments. Eper = c(E_1, E_2) supplies the period effects $E_p$, and Eco = c(reference_carryover, treatment_carryover) supplies the carry-over values. In the RT sequence, treatment in period 2 receives the reference carry-over effect; in the TR sequence, reference in period 2 receives the treatment carry-over effect. sigmaB supplies the standard deviation $\sigma_B$ of the subject effect $b_i$ on the selected analysis scale. Because the same $b_i$ is used in both periods, it induces dependence between the two observations from one subject. SigmaW supplies the residual within-subject variances and endpoint correlations. Thus, sigmaB controls between-subject heterogeneity, whereas SigmaW controls residual within-subject dependence; they represent different sources of variability.

For continuous outcomes, treatment-specific means and standard deviations are supplied through mu_list and sigma_list; SigmaW can be supplied directly or constructed from sigma_list and cor_mat. For count outcomes, the corresponding inputs are rate_list, exposure, and dispersion, as described below. The sequence structure and number of subjects per sequence are determined by dtype = "2x2" and n, respectively. These arguments define the data-generating model; the estimand is then obtained from the treatment contrast specified by ctype and the corresponding analysis kernel.

For distribution = "norm" with ctype = "DOM", the simulated outcomes are generated directly on the supplied continuous analysis scale. The equivalence test compares the treatment–reference difference with additive equivalence limits. The crossover kernel estimates treatment and reference means by averaging the corresponding sequence-period means and uses the within-subject variation in the TOST standard error.

For distribution = "lnorm", the supplied arithmetic means and standard deviations are converted to the corresponding log-scale means and covariance matrix before simulation. The equivalence limits are also transformed by the log function, and the DOM kernel is then applied on the log scale. Consequently, the resulting decision has a ratio interpretation on the original scale. A Log-Normal outcome must be strictly positive; structural zeros or negative values require a different model or a prespecified transformation.

The continuous 2x2 implementation assumes that:

• both sequences contain subjects observed in both periods after the supplied sequence-level dropout adjustment;
• treatment, period, and carry-over effects are additive on the analysis scale, with no treatment-by-period interaction unless represented through the supplied inputs;
• the within-subject covariance matrix is valid and positive definite;
• the subject effect is adequately represented by a normal random effect and is independent of the residual errors; and
• dropout is non-informative and reduces the available subjects by sequence; period-specific missing observations are not modelled by this continuous kernel.

Under balanced sequences and no differential carry-over, the treatment effect is identified after averaging the two sequence-specific estimates and the period effect cancels. Non-zero carry-over values change the simulated means and therefore change the estimand being evaluated; they should be prespecified and subjected to sensitivity analysis.

Endpoint dependence for these continuous crossover simulations is supplied directly through SigmaW or through sigma_list and cor_mat when the covariance matrix is constructed. This differs from the count implementation, where cor_mat is used as a latent Gaussian-copula correlation. Non-zero carry-over values change the simulated outcomes and therefore the estimand being evaluated; they should be prespecified and examined in sensitivity analyses.

Poisson and negative-binomial outcomes

For a 2x2 crossover design, the count extension uses a paired log-link rate analysis. Each complete participant contributes one count under the test treatment and one count under the reference treatment. Conditional on the subject effect, the count-generating model is

$$Y_{isp} \sim \begin{cases} \operatorname{Poisson}(e_{isp}\lambda_{isp}), &\text{Poisson},\\ \operatorname{NegBin}(e_{isp}\lambda_{isp},\phi_{isp}), &\text{Negative Binomial}, \end{cases}$$

with

$$\log(\lambda_{isp}) = \log(\lambda_{i}) + E_{p} + C_{isp} + b_{is}, \qquad b_{is}\sim N(0,\sigma_B^2).$$

Here, $Y_{isp}$ is the count observed for subject $i$, sequence $s$, and period $p$, with conditional mean $e_{isp}\lambda_{isp}$. Thus, exposure is the follow-up-time offset and $\lambda_{isp}$ is the event rate per unit of exposure. The Poisson model has variance equal to its conditional mean, whereas the negative-binomial model permits extra-Poisson variation. For negative-binomial outcomes, the package uses the size-parameter convention $\operatorname{size}=1/\phi$; larger dispersion therefore produces greater variability. The common crossover parameters and their package arguments were defined above; count-specific baseline rates are supplied through rate_list. These rate-ratio and Poisson/negative-binomial equivalence assumptions follow published methods for count-outcome equivalence trials (Chang et al. 2017; Zhu 2017).

Using the common crossover construction described above, the analysis forms a log-rate contrast within each complete participant, averages the contrasts within each sequence, and then averages the two sequence estimates. This removes the period effect under a balanced 2x2 design and retains the specified carry-over correction, exposure offset, sampling variation, and negative- binomial dispersion in the simulated standard error.

The within-subject contrast removes the subject random intercept from the treatment effect. Consequently, sigmaB affects the simulated paired counts but should not be interpreted as an additional treatment-effect variance in the paired contrast. This is a conditional, paired rate-ratio analysis; it is not a marginal population-average model and it does not estimate an independent carry-over coefficient. If substantial carry-over is expected, the design and estimand should be reconsidered because correction relies on the supplied Eco values.

Sequence-specific dropout is applied before the paired analysis. Participants missing either period are excluded from the within-subject contrast. The method therefore assumes non-informative dropout and requires enough complete participants in both sequences. A crossover design should not be used merely as a computational substitute for a parallel design.

For multiple endpoints in either design, the endpoint-specific analyses are combined using k, while the Gaussian copula determines the joint simulated success event. Thus, cor_mat is used by both the joint parallel engine and the multi-endpoint 2x2 engine. The supplied matrix remains a latent Gaussian correlation, not necessarily the Pearson correlation of the observed counts.

Practical sensitivity analyses

Because power depends on assumptions that are rarely known exactly, planning should examine sensitivity to:

• event rates or continuous-outcome means;
• exposure and dropout;
• standard deviations or covariance matrices;
• negative-binomial dispersion;
• endpoint correlations or the latent count copula correlation;
• the number of required endpoints k; and
• the multiplicity adjustment.

The selected scenario should be justified in the statistical analysis plan, and the final reported sample size should account for Monte Carlo uncertainty and any operational inflation required by the study.

Chang, Yu-Wei, Yi Tsong, and Zhigen Zhao. 2017. “Sample Size Determination for a Three-Arm Equivalence Trial of Poisson and Negative Binomial Responses.” Journal of Biopharmaceutical Statistics 27 (2): 239–56. https://doi.org/10.1080/10543406.2016.1269787.

Chow, Shein-Chung, and Jen-pei Liu. 2008. Design and Analysis of Bioavailability and Bioequivalence Studies. 3rd ed. Chapman; Hall/CRC.

Nelsen, Roger B. 2006. An Introduction to Copulas. 2nd ed. Springer.

Zhu, Haiyuan. 2017. “Sample Size Calculation for Comparing Two Poisson or Negative Binomial Rates in Non-Inferiority or Equivalence Trials.” Statistics in Biopharmaceutical Research 9 (1): 107–15. https://doi.org/10.1080/19466315.2016.1225594.