Update vignette explanations and references

Author: JoMu-DW

vignettes/continuous_crossover.Rmd

@@ -18,7 +18,11 @@ doc.cache <- T #for cran; change to F
 ```
 
 
-In the example below, we illustrate the use of `SimTOST` for 2x2 cross-over trials. As a first step, we load the package.
+In the example below, we illustrate the use of `SimTOST` for 2x2 crossover
+trials. This randomized two-period, two-sequence design and its standard
+subject/period/treatment model are described by Chow and Liu
+[@chow_liu_bioequivalence_2008]. The design is also used in bioequivalence
+studies [@CHMP2010]. As a first step, we load the package.
 
 ```{r, echo = T, message=F}
 library(SimTOST)

vignettes/count.Rmd

@@ -85,7 +85,9 @@ approximately `2 + 0.50 * 2^2 = 4`. Larger dispersion means more variation
 between participants and usually less information per participant, which can
 increase the required sample size. The dispersion should be chosen from pilot
 data or prior knowledge and examined in a sensitivity analysis; it is not a
-replacement for the event rate or the exposure.
+replacement for the event rate or the exposure. Poisson and negative-binomial
+rate equivalence planning is described by Chang et al. and Zhu
+[@chang_sample_2017; @zhu_sample_2017].
 
 The remaining planning inputs have the same meaning as for continuous
 outcomes: `power`, `alpha`, `nsim`, `seed`, and `dtype` specify the target
@@ -132,7 +134,8 @@ discrete event totals and reports a finite-sample Monte Carlo estimate.
 For this example, we utilize a two-endpoint, two-comparator count outcome
 study, which will serve as the foundation for the remainder of this vignette.
 The supplied endpoint correlation is a latent Gaussian-copula
-correlation, not the expected Pearson correlation of the observed counts.
+correlation, not the expected Pearson correlation of the observed counts
+[@nelsen2006].
 
 ```{r poisson-inputs}
 count_corr <- matrix(c(1, 0.5, 0.5, 1), nrow = 2,
@@ -182,7 +185,8 @@ The empirical rate distributions should be centered near the values in
 plot_distribution(poisson_diagnostic, estimand = "correlation",
                   arms = c("TEST", "REF"))
 ```
-The dashed line is the latent correlation supplied through `count_corr`. The
+The dashed line is the latent correlation supplied through `count_corr`
+[@nelsen2006]. The
 observed Pearson correlation can differ because the latent normal variables
 are transformed into discrete counts. The plot checks the direction and
 approximate magnitude of the generated dependence; `endpoint_corr` is the
@@ -315,4 +319,3 @@ settings but evaluates paired treatment-period data. The values above are
 illustrative; in an actual study, `sigmaB`, period effects, carry-over effects,
 and sequence-specific dropout should be justified from prior data or a
 sensitivity analysis.
-

vignettes/count_parallel_3A3E.Rmd

@@ -34,7 +34,8 @@ 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.
+the `cor_mat` argument. This dependence is generated with a Gaussian-copula
+construction [@nelsen2006].
 
 ## Study assumptions
 
@@ -66,6 +67,8 @@ 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
 
@@ -246,8 +249,8 @@ 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. This is therefore a joint count simulation rather than a maximum of
-separate comparison-specific searches.
+model [@nelsen2006]. This is therefore a joint count simulation rather than a
+maximum of separate comparison-specific searches.
 
 ## Why joint planning is important
 

vignettes/intropkg.Rmd

@@ -38,6 +38,8 @@ and
 $$H_{02}: \mu_T - \mu_R \geq U
 \qquad\text{versus}\qquad
 H_{12}: \mu_T - \mu_R < U.$$
+The TOST formulation follows the standard equivalence-testing framework
+[@sozu_sample_2015; @schuirmann_comparison_1987].
 
 Equivalence is established if and only if both $H_{01}$ and $H_{02}$ are
 rejected at the chosen significance level $\alpha$. Equivalently, the
@@ -63,6 +65,10 @@ For pharmacokinetic (PK) outcomes, such as the area under the curve (AUC) and ma
 \text{Logarithmic Variance} = \log\left(1 + {\text{CV}^2}\right)
 \]
 
+The use of the logarithmic scale for bioequivalence analysis is described in
+standard sample-size and regulatory guidance [@sozu_sample_2015; @CHMP2010].
+The equation itself follows from the definition of a log-normal distribution.
+
 Equivalence limits must also be specified on the log scale to align with the transformed data.
 
 
@@ -79,6 +85,8 @@ and
 $$H_{02}: \frac{\mu_T}{\mu_R} \geq U
 \qquad\text{versus}\qquad
 H_{12}: \frac{\mu_T}{\mu_R} < U.$$
+The ratio formulation is the corresponding TOST on the ratio-of-means scale
+[@sozu_sample_2015; @schuirmann_comparison_1987].
 
 Thus, equivalence requires rejection of both $H_{01}$ and $H_{02}$, or
 equivalently:
@@ -105,16 +113,16 @@ $$
 The logarithmic mean is then calculated as:
 
 $$\text{Logarithmic Mean} = \log(\text{Arithmetic Mean}) - \frac{1}{2}(\text{Logarithmic Variance})$$
-
+These transformations use the standard log-normal parameterization.
 ## Count outcomes: event-rate ratio (RR)
 
 For count outcomes, such as the number of events observed during a specified
 exposure period, SimTOST uses the event-rate ratio rather than the continuous
 outcome ratio of means. The event-rate ratio is
 
 $$RR = \frac{\lambda_T}{\lambda_R},$$
-
-where $\lambda_T$ and $\lambda_R$ are the event rates in the test and
+The event-rate ratio is the estimand used for count-outcome equivalence
+comparisons [@chang_sample_2017; @zhu_sample_2017], where $\lambda_T$ and $\lambda_R$ are the event rates in the test and
 reference groups. For lower and upper rate-ratio limits $L$ and $U$, the
 hypotheses are:
 
@@ -158,7 +166,9 @@ counts. The resulting Pearson correlation between the observed discrete
 counts is generally different and depends on the marginal count
 distributions. For normally distributed continuous outcomes, the latent
 correlation and Pearson correlation coincide because the transformation is
-linear [@chang_sample_2017].
+linear. This follows directly from the definition of Pearson correlation: a
+common linear transformation preserves the correlation, apart from a possible
+change in sign if the transformation reverses orientation.
 
 ## Regulatory requirements
 When evaluating bioequivalence, certain statistical and methodological requirements must be adhered to, as outlined in the European Medicines Agency's bioequivalence guidelines [@CHMP2010]. These requirements ensure that the test and reference products meet predefined criteria for equivalence in terms of PK parameters. The key considerations are summarized below:
@@ -185,7 +195,7 @@ A decision must be made as to whether it is desirable to
 
 ## Testing multiple co-primary endpoints
 
-When a trial defines multiple co-primary endpoints, equivalence must be demonstrated for all of them to claim overall treatment equivalence. In this setting, each endpoint is tested separately at the usual significance level ($\alpha$), and equivalence is established only if all individual tests are statistically significant. Because conclusions require rejecting all null hypotheses, a formal multiplicity adjustment is not needed to control the Type I error rate. [@cpmp_points_2002] However, as the number of co-primary endpoints ($K$) increases, the likelihood of failing to meet equivalence on at least one endpoint also rises, resulting in a higher Type II error rate (i.e., a greater risk of incorrectly concluding non-equivalence) [@mielke_sample_2018]
+When a trial defines multiple co-primary endpoints, equivalence must be demonstrated for all of them to claim overall treatment equivalence. In this setting, each endpoint is tested separately at the usual significance level ($\alpha$), and equivalence is established only if all individual tests are statistically significant. Because conclusions require rejecting all null hypotheses, the decision is an intersection--union test and a formal endpoint-wise multiplicity adjustment is not needed to control the Type I error rate [@berger1996; @hamasaki_coprimary_2018]. However, as the number of co-primary endpoints ($K$) increases, the likelihood of failing to meet equivalence on at least one endpoint also rises, resulting in a higher Type II error rate (i.e., a greater risk of incorrectly concluding non-equivalence) [@mielke_sample_2018]
 
 This has several important implications:
 
@@ -197,43 +207,51 @@ This has several important implications:
 
 ## Testing multiple primary endpoints {#multiple-primary}
 
-When a trial aims to establish equivalence for at least $k$ primary endpoints, adjustments are necessary to account for the increased risk of Type I error due to multiple hypothesis testing [@sozu_continuous_2015]. Without such adjustments, the likelihood of incorrectly concluding equivalence for at least one endpoint increases as the number of endpoints grows.
+When a trial aims to establish equivalence for at least $k$ of $m$ primary
+endpoints, the decision is a partial-conjunction test. When $k < m$, an
+endpoint-wise adjustment is generally needed to control the Type I error across
+the possible successful endpoint subsets [@benjamini2008; @mielke_sample_2018].
+The appropriate adjustment depends on the prespecified $k$-out-of-$m$ rule.
 
 For example, if a study includes $m = 3$ independent primary endpoints and uses a significance level of $\alpha = 5%$ for each test, the overall probability of falsely concluding equivalence for at least one endpoint is:
 
 $$ 1 – (1-\alpha)^m  = 1 - (1-0.05)^3 = 0.1426. $$
+This is the probability of at least one false rejection under independent
+tests [@dunn1961; @sidak1967].
 This means that the overall probability of making any false positive error, also known as the **family-wise error rate (FWER)**, increases to approximately 14%.
 
 To address this issue, adjustments to the significance level are necessary for multiple endpoint comparisons for which various methods have been proposed. In SimTOST, the following approaches are included:
 
 ### Bonferroni correction
-The most common and easiest procedure for multiplicity adjustment to control the FWER is the Bonferroni method. Each hypothesis is tested at level
+The most common and easiest procedure for multiplicity adjustment to control the FWER is the Bonferroni method [@bonferroni1936]. Each hypothesis is tested at level
 
 $$\alpha_{bon}= \alpha/m$$
 
+
 where $m$ is the total number of tests. Although simple, this method is highly conservative, particularly when tests are correlated, as it assumes all tests are independent. This conservativeness remains pronounced even for $k=1$, where only one of the $m$ hypotheses needs to be rejected. [@mielke_sample_2018]
 
 In the [sampleSize()](../reference/sampleSize.html) function, the Bonferroni correction can be applied by setting `adjust = "bon"`.
 
 ### Sidak correction
 
-The Sidak correction is an alternative method for controlling the FWER. Like the Bonferroni correction, it assumes that tests are independent. However, the Sidak correction accounts for the joint probability of all tests being non-significant, making it mathematically less conservative than the Bonferroni method. The adjusted significance level is calculated as:
+The Sidak correction [@sidak1967] is an alternative method for controlling the FWER. Like the Bonferroni correction, it assumes that tests are independent. However, the Sidak correction accounts for the joint probability of all tests being non-significant, making it mathematically less conservative than the Bonferroni method. The adjusted significance level is calculated as:
 
 $$\alpha_{sid}= 1-(1-\alpha)^ {1/m}$$
 
+
 The Sidak correction can be implemented by specifying `adjust = "sid"` in the [sampleSize()](../reference/sampleSize.html) function.
 
-### K adjustment
+### K adjustment: Mielke's weak k-out-of-m correction
 
-This correction explicitly accounts for the scenario where equivalence is required for only $k$ out of $m$ endpoints. Unlike the Bonferroni and Sidak corrections, which assume that all  $m$ tests contribute equally to the overall Type I error rate, the *k*-adjustment directly incorporates the number of endpoints ($k$) required for equivalence into the adjustment. The adjusted significance level is calculated as:
+The K adjustment [@mielke_sample_2018] explicitly accounts for the scenario where equivalence is required for only $k$ out of $m$ endpoints. Unlike the Bonferroni and Sidak corrections, which assume that all  $m$ tests contribute equally to the overall Type I error rate, the *k*-adjustment directly incorporates the number of endpoints ($k$) required for equivalence into the adjustment. The adjusted significance level is calculated as:
 
 $$\alpha_k= \frac{k*\alpha}{m}$$
 
 where $k$ is the number of endpoints required for equivalence, and $m$ is the total number of endpoints evaluated.
 
 ### T adjustment: Mielke's strong k-out-of-m correction
 
-The `t` adjustment described below is the strong k-out-of-m procedure for
+The `t` adjustment  [@mielke_sample_2018] described below is the strong k-out-of-m procedure for
 partial-conjunction decisions and is distinct from the weak `k` adjustment.
 
 When the decision rule requires equivalence for at least `k` of `m` endpoints,
@@ -245,8 +263,7 @@ $$\alpha_{t}= \frac{\alpha}{m-k+1}$$
 
 where `m` is the number of selected endpoints and `k` is the minimum number
 that must demonstrate equivalence. This adjustment controls
-the Type I error rate for the partial-conjunction decision, including cases in
-which only some endpoints satisfy the null hypothesis [@mielke_sample_2018].
+the Type I error rate for the partial-conjunction decision, including cases in which only some endpoints satisfy the null hypothesis [@mielke_sample_2018].
 
 ## Hierarchical testing of multiple endpoints
 Hierarchical testing is an approach to multiple endpoint testing where endpoints are tested in a predefined order, typically based on their clinical or regulatory importance. A fallback testing strategy is applied, allowing sequential hypothesis testing. If a hypothesis earlier in the sequence fails to be rejected, testing stops, and subsequent hypotheses are not evaluated. [@chowdhry_finding_2024]

vignettes/methodological_assumptions.Rmd

@@ -1,6 +1,8 @@
 ---
 title: "Methodological Assumptions in SimTOST"
 output: rmarkdown::html_vignette
+bibliography: 'references.bib'
+link-citations: yes
 vignette: >
   %\VignetteIndexEntry{Methodological Assumptions in SimTOST}
   %\VignetteEngine{knitr::rmarkdown}
@@ -203,13 +205,14 @@ latent correlation.
 
 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
+standard-normal vector using a Gaussian-copula construction [@nelsen2006]:
 
 \[
 \boldsymbol{Z}_a \sim N_m(\boldsymbol{0},\boldsymbol{R})
 \]
 
-for arm \(a\). Each component is transformed as
+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
@@ -220,7 +223,8 @@ 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.
+0.8. This Gaussian-copula construction is a standard way to generate
+dependent non-normal outcomes.
 
 The matrix has the following interpretation:
 
@@ -266,7 +270,11 @@ 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
+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_liu_bioequivalence_2008].
 
 \[
   Y_{ipj} = \mu_{a(i,p),j} + E_p + C_{a(i,p),s(i)} + b_i +
@@ -332,7 +340,9 @@ where `cor_mat` is used as a latent Gaussian-copula correlation.
 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. The count-generating
-model is, conditionally on the subject effect,
+model is, conditionally on the subject effect. The rate-ratio estimand and
+Poisson/negative-binomial framework is supported by published count-outcome
+equivalence methods [@chang_sample_2017; @zhu_sample_2017].
 
 \[
 Y_{isp} \sim \begin{cases}
@@ -354,7 +364,9 @@ period effects, `Eco = c(C_R,C_T)` contains reference- and treatment-
 carry-over effects, and `sigmaB` is the standard deviation of the subject
 random intercept on the log-rate scale. The negative-binomial dispersion is
 used in the subject-period count generation through the package's
-size-parameter convention \(\operatorname{size}=1/\phi\).
+size-parameter convention \(\operatorname{size}=1/\phi\). This is the
+standard mean--dispersion parameterization used for planning negative-binomial
+rate comparisons.
 
 The analysis forms a log-rate contrast within each complete participant,
 averages the contrasts within each sequence, and then averages the two