Preliminaries

A basic knowledge of R is required. To run the scripts at least version 1.4.3 (2016-11-01) of PowerTOST is suggested. Any version of R would likely do, though the current release of PowerTOST was only tested with version 4.0.5 (2021-03-31) and later.
All scripts were run on a Xeon E3-1245v3 @ 3.40GHz (1/4 cores) 16GB RAM with 64 bit R 4.1.3 on Windows 7.

Note that in all functions of PowerTOST the arguments (say, the assumed T/R-ratio theta0, the BE-limits (theta1, theta2), the assumed coefficient of variation CV, etc.) have to be given as ratios and not in percent.

sampleN.TOST() gives balanced sequences (i.e., the same number of subjects is allocated to each of the sequences). Furthermore, the estimated sample size is the total number of subjects (not subjects per sequence – like in some other software packages).

Higher-Order Designs with three and four treatments are supported.

"3x3" denotes the Latin Square (ABC|BCA|CAB).

"3x6x3" denotes the 6-sequence Williams’ design (ABC|ACB|BAC|BCA|CAB|CBA).

"4x4" denotes the Latin Square (ABCD|BCDA|CDAB|DABC) or
any of the possible Williams’ designs with four periods (ADBC|BACD|CBDA|DCAB, ADCB|BCDA|CABD|DBAC, ACDB|BDCA|CBAD|DABC, ACBD|BADC|CDAB|DBCA, ABDC|BCAD|CDBA|DACB, ABCD|BDAC|CADB|DCBA).

df are the degrees of freedom, where n is the sample size.

“Suppose we have a bioequivalence study with three treatments – A, B, and C – and the objective of the study is to make pairwise comparisons among the treatments. Suppose further that treatment C is different in kind from A and B, so that the assumption of homogeneous variance among the three treatments is questionable. One way to do the analyses, under normality assumptions, is Two at a Time – e.g., to test hypotheses about A and B, use only the data from A and B. Another way is All at Once – include the data from all three treatments in a single analysis, making pairwise comparisons within this analysis. If the assumption of homogeneous variance is correct, the All at Once approach will provide more d.f. for estimating the common variance, resulting in increased power. If the variance of C differs from that of A and B, the All at Once approach may have reduced power or an inflated type I error rate, depending on the direction of the difference in variances.

Which one you should use is discussed further down in the sub-sections

• Two at a Time
• All at Once

of Multiple Treatments.

Most examples deal with studies where the response variables likely follow a log-normal distribution, i.e., we assume a multiplicative model (ratios instead of differences). We work with \(\small{\log_{e}}\)-transformed data in order to allow analysis by the t-test (requiring differences). This is the default in most functions of PowerTOST and hence, the argument logscale = TRUE does not need to be specified.

previous section ↩︎

Terminology

It may sound picky but ‘sample size calculation’ (as used in most guidelines and alas, in some publications and textbooks) is sloppy terminology. In order to get prospective power (and hence, a sample size), we need five values:

1. The level of the test \(\small{\alpha}\) (in BE commonly 0.05),
2. the BE-margins (commonly 0.80 – 1.25),
3. the desired (or target) power \(\small{\pi}\),
4. the variance (commonly expressed as a coefficient of variation), and
5. the deviation of the test from the reference treatment.

1 – 2 are fixed by the agency,
3 is set by the sponsor (commonly to 0.80 – 0.90), and
4 – 5 are just (uncertain!) assumptions.

In other words, obtaining a sample size is not an exact calculation like \(\small{2\times2=4}\) but always just an estimation.

“Power Calculation – A guess masquerading as mathematics.

Of note, it is extremely unlikely that all assumptions will be exactly realized in a particular study. Hence, calculating retrospective (a.k.a. post hoc, a posteriori) power is not only futile but plain nonsense.⁶

Since generally the within-subject variability is lower than the between-subject variability, crossover studies are so popular. The efficiency of a crossover study compared to a parallel study is given by \(\small{\frac{\sigma_\textrm{intra}^2\;+\,\sigma_\textrm{inter}^2}{0.5\,\times\,\sigma_\textrm{intra}^2}}\). If, say, \(\small{\sigma_\textrm{intra}^2=0.5\times\sigma_\textrm{inter}^2}\) in a paralled study we need six times as many subjects than in a crossover to obtain the same power. On the other hand, in a crossover we have two measurements per subject, which makes the parallel study approximately three times more costly.

Of note, there is no relationship between \(\small{CV_\textrm{intra}}\) and \(\small{CV_\textrm{inter}}\).
An example are drugs which are subjected to polymorphic metabolism. For these drugs \(\small{CV_\textrm{intra}\ll CV_\textrm{inter}}\). On the other hand, some HVD(P)s show \(\small{CV_\textrm{intra}>CV_\textrm{inter}}\).

However, it is a prerequisite that no carryover from one period to the next exists. Only then the comparison of treatments will be unbiased. For details see another article.⁷

Subjects have to be in the same physiological state⁸ throughout the study – guaranteed by sufficiently long wash­out phases. Crossover studies can not only be performed in healthy volunteers but also in patients with a stable disease (e.g., asthma). Studies in patients with an instable disease (e.g., in oncology) must be performed in a pa­ral­lel design (covered in another article).
If crossovers are not feasible (e.g., for drugs with a very long half life), studies could be performed in a parallel design as well.

previous section ↩︎

Planned Evaluation

First of all we have have to consider the type of study, its purpose, and how we intend to evaluate it.

1. A pilot study with candidates vs a reference.
2. A pivotal study of a test treatment vs references from different regions,
   a pivotal study of test vs reference in fasting and fed state,
   a pivotal study of ODTs administered with and without water.
3. A pivotal study with candidates vs a reference.

Let’s discuss the cases from above and whether we can evaluate the study by the common \(\small{\alpha=0.05}\) (i.e., with a \(\small{90\%}\) confidence interval) or have to adjust for multiplicity (\(\small{\alpha<0.05 \to\; >90\%}\) CI).

1. Say, we have two candidate treatments (\(\small{\textrm{C}_1}\), \(\small{\textrm{C}_2}\)) and one reference (\(\small{\textrm{R}}\)).
   We would select the candidate with \(\small{\min \left \{\left|\log_{e}\theta_{\textrm{C}_1}\right|,\left|\log_{e}\theta_{\textrm{C}_2}\right|\right\}}\) for the pivotal study because it is the most promising one. If \(\small{\min \left \{\left|\log_{e}\theta_{\textrm{C}_1}\right|\sim\left|\log_{e}\theta_{\textrm{C}_2}\right|\right\}}\), we would select the candidate with lower within-subject variance.
   No multiplicity correction has to be applied.
2. In the ‘Two at a Time’ approach we obtain not only point estimates of the pairwise comparisons (like in the ‘All at Once’ approach) but also their within-subject variances (e.g., in the comparisons \(\small{\textrm{T}}\) vs \(\small{\textrm{R}_1}\) and \(\small{\textrm{T}}\) vs \(\small{\textrm{R}_2}\)).¹⁴
   In both approaches no multiplicity correction has to be applied.¹⁵
3. Sometimes the pilot was indecisive (very similar T/R ratios and CVs) and companies are wary to select one based on ‘gut feelings’ and include both in the pivotal study.
   Submit the ‘better’ one to the authority and stop developing the other.
   Opinion split amongst statisticians:

1. Since only one product will be marketed, this approach does not increase the patient’s risk (the \(\small{90\%}\) CI is sufficient).
2. The company gets to two chances to show BE and selects one product based on the outcome. This inflates the Type I Error and a Bonferroni adjustment is recommended.¹⁶

A Simple Case

We assume a CV of 0.25, a T/R-ratio of 0.95, and target a power of 0.80.

Since theta0 = 0.95 and targetpower = 0.80 are defaults of the function, we don’t have to give them explicitely. As usual in bioequivalence, alpha = 0.05 is employed (we will assess the study by a \(\small{100\,(1-2\,\alpha)=90\%}\) confidence interval) and the BE-limits are theta1 = 0.80 and theta2 = 1.25. Since they are also defaults of the function, we don’t have to give them as well. Hence, we need to specify only the CV and the design.

sampleN.TOST(CV = 0.25, design = "3x3") # All at Once

# 
# +++++++++++ Equivalence test - TOST +++++++++++
#             Sample size estimation
# -----------------------------------------------
# Study design: 3x3 crossover 
# log-transformed data (multiplicative model)
# 
# alpha = 0.05, target power = 0.8
# BE margins = 0.8 ... 1.25 
# True ratio = 0.95,  CV = 0.25
# 
# Sample size (total)
#  n     power
# 27   0.803494

sampleN.TOST(CV = 0.25, design = "2x2") # Two at a Time

# 
# +++++++++++ Equivalence test - TOST +++++++++++
#             Sample size estimation
# -----------------------------------------------
# Study design: 2x2 crossover 
# log-transformed data (multiplicative model)
# 
# alpha = 0.05, target power = 0.8
# BE margins = 0.8 ... 1.25 
# True ratio = 0.95,  CV = 0.25
# 
# Sample size (total)
#  n     power
# 28   0.807439

Don’t get confused by Study design: 2x2 crossover in the output. It only means that we will evaluate the study by the ‘Two at a Time’ approach as IBDs. In the ‘All at Once’ approach we have more degrees of freedom (50) than in the ‘Two at a Time’ approach (26). An example to get the correct sample size is given in the Interlude further down.

Sometimes we are not interested in the entire output and want to use only a part of the results in subsequent calculations. We can suppress the output by the argument print = FALSE and assign the result to a data frame (here named x).

x <- sampleN.TOST(CV = 0.25, design = "3x3",
                  print = FALSE)

Although you could access the elements by the number of the column(s), I don’t recommend that, since in various functions these numbers are different and hence, difficult to remember.

Let’s retrieve the column names of x:

names(x)
# [1] "Design"         "alpha"          "CV"            
# [4] "theta0"         "theta1"         "theta2"        
# [7] "Sample size"    "Achieved power" "Target power"

Now we can access the elements of x by their names. Note that double square brackets [[…]] have to be used.

x[["Sample size"]]
# [1] 27
x[["Achieved power"]]
# [1] 0.8034938

Although you could access the elements by the number of the column(s), I don’t recommend that, since in various other functions of PowerTOST these numbers are different and hence, difficult to remember. If you insist in accessing elements by column-numbers, use single square brackets […].

x[7:8]
#   Sample size Achieved power
# 1          27      0.8034938

With 27 subjects (nine per sequence) we achieve the power we desire.

What happens if we have one dropout?

power.TOST(CV = 0.25, design = "3x3",
           n = x[["Sample size"]] - 1)

# Unbalanced design. n(i)=9/9/8 assumed.

# [1] 0.7868869

Slightly below the 0.80 we desire.

get.sequences <- function(design) {
  known <- known.designs()[, c(2, 5)]
  ns    <- as.integer(known[known$design == design,
                            "steps"])
  return(ns)
}
make.equal <- function(n, ns) {
  return(as.integer(ns * (n %/% ns + as.logical(n %% ns))))
}
design <- "3x3"
CV     <- 0.25
ns     <- get.sequences(design)
n.TaT  <- sampleN.TOST(CV = CV, design = "2x2",
                       print = FALSE)[["Sample size"]]
n.AaO  <- sampleN.TOST(CV = CV, design = design,
                       print = FALSE)[["Sample size"]]
n      <- make.equal(n = n.TaT, ns = ns)
df     <- data.frame(design = design,
                     evaluation = c("\u2018All at Once\u2019",
                                    "\u2018Two at a Time\u2019"),
                     n = c(n.AaO, n), power = NA)
df[1, 4] <- power.TOST(CV = CV, design = design,
                       n = n.AaO)
df[2, 4] <- power.TOST(CV = CV, design = "2x2",
                       n = n)
print(df, row.names = FALSE)

#  design      evaluation  n     power
#     3x3   ‘All at Once’ 27 0.8034938
#     3x3 ‘Two at a Time’ 30 0.8342518

Since dropouts are common, it makes sense to include / dose more subjects in order to end up with a number of eligible subjects which is not lower than our initial estimate.

Let us explore that in the next section.

top of section ↩︎ previous section ↩︎

Lost in Assumptions

As stated already in the Introduction, the CV and the T/R-ratio are only assumptions. Whatever their origin might be (literature, previous studies) they carry some degree of uncertainty. Hence, believing²⁴ that they are the true ones may be risky.
Some statisticians call that the ‘Carved-in-Stone’ approach.

Say, we performed a pilot study in 12 subjects and estimated the CV as 0.25.

The \(\alpha\) confidence interval of the CV is obtained via the \(\small{\chi^2}\)-distribution of its error variance \(\small{\sigma^2}\) with \(\small{n-2}\) degrees of freedom. \[\begin{matrix}\tag{6} s^2=\log_{e}(CV^2+1)\\ L=\frac{(n-1)\,s^2}{\chi_{\alpha/2,\,n-2}^{2}}\leq\sigma^2\leq\frac{(n-1)\,s^2}{\chi_{1-\alpha/2,\,n-2}^{2}}=U\\ \left\{lower\;CL,\;upper\;CL\right\}=\left\{\sqrt{\exp(L)-1},\sqrt{\exp(U)-1}\right\} \end{matrix}\] Let’s calculate the 95% confidence interval of the CV to get an idea.

m  <- 12 # pilot study
ci <- CVCL(CV = 0.25, df = get.df(design, m),
           side = "2-sided", alpha = 0.05)
signif(ci, 4)

# lower CL upper CL 
#   0.1901   0.3671

Surprised? Although 0.25 is the best estimate for planning the next study, there is no guarantee that we will get exactly the same outcome. Since the \(\small{\chi^2}\)-distribution is skewed to the right, it is more likely to face a higher CV than a lower one in the pivotal study.

If we plan the study based on 0.25, we would opt for 27 subjects like in the examples before (not adjusted for the dropout-rate).
If the CV will be lower, we gain power. Fine.
But what if it will be higher? Of course, we will loose power. But how much?

Let’s explore what might happen at the confidence limits of the CV.

m    <- 12
ci   <- CVCL(CV = 0.25, get.df(design, m),
             side = "2-sided", alpha = 0.05)
n    <- 27
comp <- data.frame(CV = c(ci[["lower CL"]], 0.25,
                          ci[["upper CL"]]),
                   power = NA)
for (i in 1:nrow(comp)) {
  comp$power[i] <- power.TOST(CV = comp$CV[i],
                              design = design,
                              n = n)
}
comp[, 1] <- signif(comp[, 1], 4)
comp[, 2] <- signif(comp[, 2], 6)
print(comp, row.names = FALSE)

#      CV    power
#  0.1901 0.951526
#  0.2500 0.803494
#  0.3671 0.417982

Ouch!

What can we do? The larger the previous study was, the larger the degrees of freedom and hence, the narrower the confidence interval of the CV. In simple terms: The estimate is more certain. On the other hand, it also means that very small pilot studies are practically useless.

m               <- seq(6, 16, 2)
df              <- data.frame(n = m, CV = 0.25,
                              l = NA, u = NA)
for (i in 1:nrow(df)) {
  df[i, 3:4] <- CVCL(CV = 0.25,
                     df = get.df(design, m[i]),
                     side = "2-sided",
                     alpha = 0.05)
}
df[, 3:4]      <- signif(df[, 3:4], 4)
names(df)[3:4] <- c("lower CL", "upper CL")
print(df, row.names = FALSE)

#   n   CV lower CL upper CL
#   6 0.25   0.1675   0.4992
#   8 0.25   0.1779   0.4238
#  10 0.25   0.1849   0.3883
#  12 0.25   0.1901   0.3671
#  14 0.25   0.1940   0.3528
#  16 0.25   0.1973   0.3425

A Bayesian method is implemented in PowerTOST, which takes the uncertainty of estimates into account. We can explore the uncertainty of the CV or the T/R-ratio and of both.

m      <- 12   # sample size of pilot
CV     <- 0.25
theta0 <- 0.95
df     <- expsampleN.TOST(CV = CV, theta0 = theta0,
                          targetpower = 0.80,
                          design = design,
                          prior.parm = list(
                            m = m, design = design),
                          prior.type = "CV",
                          details = FALSE)

# 
# ++++++++++++ Equivalence test - TOST ++++++++++++
#        Sample size est. with uncertain CV
# -------------------------------------------------
# Study design:  3x3 crossover 
# log-transformed data (multiplicative model)
# 
# alpha = 0.05, target power = 0.8
# BE margins = 0.8 ... 1.25 
# Ratio = 0.95
# CV = 0.25 with 20 df
# 
# Sample size (ntotal)
#  n   exp. power
# 30   0.804060

Not that bad. 11% more subjects required.

What about an uncertain T/R-ratio?

m      <- 12
CV     <- 0.25
theta0 <- 0.95
df     <- expsampleN.TOST(CV = CV, theta0 = theta0,
                          targetpower = 0.80,
                          design = design,
                          prior.parm = list(
                            m = m, design = design),
                          prior.type = "theta0",
                          details = FALSE)

# 
# ++++++++++++ Equivalence test - TOST ++++++++++++
#      Sample size est. with uncertain theta0
# -------------------------------------------------
# Study design:  3x3 crossover 
# log-transformed data (multiplicative model)
# 
# alpha = 0.05, target power = 0.8
# BE margins = 0.8 ... 1.25 
# Ratio = 0.95
# CV = 0.25
# 
# Sample size (ntotal)
#  n   exp. power
# 183   0.881163

Terrible! What a nightmare.

That should not take us by surprise. We don’t know where the true T/R-ratio lies but we can calculate the lower 95% confidence limit of the pilot study’s point estimate to get an idea about a worst case.

m      <- 12
CV     <- 0.25
pe     <- 0.95
ci     <- round(CI.BE(CV = CV, pe = 0.95,
                      n = m, design = design), 4)
if (pe <= 1) {
  cl <- ci[["lower"]]
} else {
  cl <- ci[["upper"]]
}
print(cl)

# [1] 0.7988

Exlore the impact of a relatively 5% higher CV and a relatively 5% lower T/R-ratio on power for a given sample size.

n      <- 27
CV     <- 0.25
theta0 <- 0.95
comp1  <- data.frame(CV = c(CV, CV*1.05),
                     power = NA)
comp2  <- data.frame(theta0 = c(theta0, theta0*0.95),
                     power = NA)
for (i in 1:2) {
  comp1$power[i] <- power.TOST(CV = comp1$CV[i],
                               theta0 = theta0,
                               design = design,
                               n = n)
}
comp1$power <- signif(comp1$power, 6)
for (i in 1:2) {
  comp2$power[i] <- power.TOST(CV = CV,
                               theta0 = comp2$theta0[i],
                               design = design,
                               n = n)
}
comp2$power <- signif(comp2$power, 6)
print(comp1, row.names = FALSE)
print(comp2, row.names = FALSE)

#      CV    power
#  0.2500 0.803494
#  0.2625 0.765085
#  theta0    power
#  0.9500 0.803494
#  0.9025 0.550795

Interlude 2

Fig. 1 Power curves for n = 27 (3×3 design evaluated ‘All at Once’).

Note the log-scale of the x-axis. It demonstrates that power curves are symmetrical around 1 (\(\small{\log_{e}(1)=0}\), where \(\small{\log_{e}(\theta_2)=\left|\log_{e}(\theta_1)\right|}\)) and we will achieve the same power for \(\small{\theta_0}\) and \(\small{1/\theta_0}\) (e.g., for 0.95 and 1.0526). Furthermore, all curves intersect 0.05 at \(\small{\theta_1}\) and \(\small{\theta_2}\), which means that for a true \(\small{\theta_0}\) at one of the limits \(\small{\alpha}\) is strictly controlled.

<nitpick>

• A common flaw in protocols is the phrase
        »The sample size was calculated [sic] based on a T/R-ratio of 0.95 – 1.05…«
  If you assume a deviation of 5% of the test from the reference and are not sure about its direction (lower or higer than 1), always use the lower T/R-ratio. If you would use the upper T/R-ratio, power would be only preserved down to 1/1.05 = 0.9524.
  Given, sometimes you will need a higher sample size with the lower T/R-ratio. There’s no free lunch.

CV      <- 0.25
delta   <- 0.05 # direction unknown
theta0s <- c(1 - delta, 1 / (1 + delta),
             1 + delta, 1 / (1 - delta))
n       <- sampleN.TOST(CV = CV, theta0 = 1 - delta,
                        design = design,
                        print = FALSE)[["Sample size"]]
comp1   <- data.frame(CV = CV, theta0 = theta0s,
                      base = c(TRUE, rep(FALSE, 3)),
                      n = n, power = NA)
for (i in 1:nrow(comp1)) {
  comp1$power[i] <- power.TOST(CV = CV,
                               theta0 = comp1$theta0[i],
                               n = n,
                               design = design)
}
n       <- sampleN.TOST(CV = CV, theta0 = 1 + delta,
                        design = design,
                        print = FALSE)[["Sample size"]]
comp2   <- data.frame(CV = CV, theta0 = theta0s,
                      base = c(FALSE, FALSE, TRUE, FALSE),
                      n = n, power = NA)
for (i in 1:nrow(comp2)) {
  comp2$power[i] <- power.TOST(CV = CV,
                               design = design,
                               theta0 = comp2$theta0[i],
                               n = n)
}
comp1[, c(2, 5)] <- signif(comp1[, c(2, 5)] , 4)
comp2[, c(2, 5)] <- signif(comp2[, c(2, 5)] , 4)
print(comp1, row.names = FALSE)
print(comp2, row.names = FALSE)

#    CV theta0  base  n  power
#  0.25 0.9500  TRUE 27 0.8035
#  0.25 0.9524 FALSE 27 0.8124
#  0.25 1.0500 FALSE 27 0.8124
#  0.25 1.0530 FALSE 27 0.8035
#    CV theta0  base  n  power
#  0.25 0.9500 FALSE 27 0.8035
#  0.25 0.9524 FALSE 27 0.8124
#  0.25 1.0500  TRUE 27 0.8124
#  0.25 1.0530 FALSE 27 0.8035

</nitpick>

Since power is much more sensitive to the T/R-ratio than to the CV, the results obtained with the Bayesian method should be clear now.

Essentially this leads to the murky waters of prospective power analysis, which is covered in another article.
An appetizer to show the maximum deviations (CV, T/R-ratio and decreased sample size due to dropouts) which give still a minimum acceptable power of ≥ 0.70:

CV       <- 0.25
theta0   <- 0.95
target   <- 0.80
minpower <- 0.70
pa       <- pa.ABE(CV = CV, theta0 = theta0,
                   targetpower = target,
                   minpower = minpower,
                   design = design)
change.CV     <- 100*(tail(pa$paCV[["CV"]], 1) -
                      pa$plan$CV) / pa$plan$CV
change.theta0 <- 100*(head(pa$paGMR$theta0, 1) -
                      pa$plan$theta0) /
                      pa$plan$theta0
change.n      <- 100*(tail(pa$paN[["N"]], 1) -
                      pa$plan[["Sample size"]]) /
                      pa$plan[["Sample size"]]
comp     <- data.frame(parameter = c("CV", "theta0", "n"),
                       change = c(change.CV,
                                  change.theta0,
                                  change.n))
comp$change <- sprintf("%+.2f%%", comp$change)
names(comp)[2] <- "relative change"
print(pa, plotit = FALSE)
print(comp, row.names = FALSE)

# Sample size plan ABE
#  Design alpha   CV theta0 theta1 theta2 Sample size
#     3x3  0.05 0.25   0.95    0.8   1.25          27
#  Achieved power
#       0.8034938
# 
# Power analysis
# CV, theta0 and number of subjects leading to min. acceptable power of ~0.7:
#  CV= 0.2827, theta0= 0.9276
#  n = 22 (power= 0.7106)
# 
#  parameter relative change
#         CV         +13.09%
#     theta0          -2.36%
#          n         -18.52%

Confirms what we have seen above. Interesting that the sample size is the least sensitive one. Many overrate the impact of dropouts on power.

If you didn’t stop reading in desperation (understandable), explore both uncertain CV and T/R-ratio:

m      <- 12
CV     <- 0.25
theta0 <- 0.95
df     <- expsampleN.TOST(CV = CV, theta0 = theta0,
                          targetpower = 0.80,
                          design = design,
                          prior.parm = list(
                            m = m, design = design),
                          prior.type = "both",
                          details = FALSE)

# 
# ++++++++++++ Equivalence test - TOST ++++++++++++
#   Sample size est. with uncertain CV and theta0
# -------------------------------------------------
# Study design:  3x3 crossover 
# log-transformed data (multiplicative model)
# 
# alpha = 0.05, target power = 0.8
# BE margins = 0.8 ... 1.25 
# Ratio = 0.95 with 20 df
# CV = 0.25 with 20 df
# 
# Sample size (ntotal)
#  n   exp. power
# 198   0.866206

I don’t suggest that you should use it in practice. AFAIK, not even the main author of this function (Benjamin Lang) does. However, it serves educational purposes to show that it is not that easy and why even properly planned studies might fail.

An alternative to assuming a T/R-ratio is based on ‘Statistical Assurance’.²⁵ This concept uses the distribution of T/R-ratios and assumes an uncertainty parameter \(\small{\sigma_\textrm{u}}\). A natural assumption is \(\small{\sigma_\textrm{u}=1-\theta_0}\), i.e., for the commonly applied T/R-ratio of 0.95 one can use the argument sem = 0.05 in the function expsampleN.TOST(), where the argument theta0 has to be fixed at 1.

CV        <- 0.25
theta0    <- 0.95
target    <- 0.80
sigma.u   <- 1 - theta0
comp      <- data.frame(theta0 = theta0,
                        n.1 = NA, power = NA,
                        sigma.u = sigma.u,
                        n.2 = NA, assurance = NA)
comp[2:3] <- sampleN.TOST(CV = CV,
                          targetpower = target,
                          theta0 = theta0,
                          design = design,
                          details = FALSE,
                          print = FALSE)[7:8]
comp[5:6] <- expsampleN.TOST(CV = CV,
                             theta0 = 1, # fixed!
                             targetpower = target,
                             design = design,
                             prior.type = "theta0",
                             prior.parm =
                               list(sem = sigma.u),
                             details = FALSE,
                             print = FALSE)[9:10]
names(comp)[c(2, 5)]  <- rep("n", 2)
print(signif(comp, 6), row.names = FALSE)

#  theta0  n    power sigma.u  n assurance
#    0.95 27 0.803494    0.05 27  0.813331

top of section ↩︎ previous section ↩︎

Multiple Endpoints

It is not unusual that equivalence of more than one endpoint has to be demonstrated. In bioequivalence the pharmacokinetic metrics C_max and AUC_0–t are mandatory (in some jurisdictions like the FDA additionally AUC_0–∞).

We don’t have to worry about multiplicity issues (inflated Type I Error) since if all tests must pass at level \(\small{\alpha}\), we are protected by the intersection-union principle.^{26 27}

We design the study always for the worst case combination, i.e., based on the PK metric requiring the largest sample size. In some jurisdictions wider BE limits for C_max are acceptable. Let’s explore that with different CVs and T/R-ratios.

metrics <- c("Cmax", "AUCt", "AUCinf")
CV      <- c(0.25, 0.19, 0.20)
theta0  <- c(0.95, 1.04, 1.06)
theta1  <- c(0.75, 0.80, 0.80)
theta2  <- 1 / theta1
target  <- 0.80
df      <- data.frame(metric = metrics,
                      theta1 = theta1, theta2 = theta2,
                      CV = CV, theta0 = theta0, n = NA)
for (i in 1:nrow(df)) {
  df$n[i] <- sampleN.TOST(CV = CV[i],
                          theta0 = theta0[i],
                          theta1 = theta1[i],
                          theta2 = theta2[i],
                          targetpower = target,
                          design = design,
                          print = FALSE)[["Sample size"]]
}
df$theta1 <- sprintf("%.4f", df$theta1)
df$theta2 <- sprintf("%.4f", df$theta2)
txt       <- paste0("Sample size based on ",
                     df$metric[df$n == max(df$n)], ".\n")
print(df, row.names = FALSE)
cat(txt)

#  metric theta1 theta2   CV theta0  n
#    Cmax 0.7500 1.3333 0.25   0.95 18
#    AUCt 0.8000 1.2500 0.19   1.04 15
#  AUCinf 0.8000 1.2500 0.20   1.06 21
# Sample size based on AUCinf.

Even if we assume the same T/R-ratio for two PK metrics, we will get a wider margin for the one with lower variability.

Let’s continue with the conditions of our previous examples, this time assuming that the CV and T/R-ratio were applicable for C_max. As common in PK, the CV of AUC is lower, say, only 0.20. That means, the study is ‘overpowered’ for the assumed T/R-ratio of AUC.

Which are the extreme T/R-ratios (largest deviations of T from R) giving still the target power?

opt <- function(x) {
  power.TOST(theta0 = x, CV = df$CV[2],
             theta1 = theta1, theta2 = theta2,
             n = df$n[1], design = design) - target
}
metrics <- c("Cmax", "AUC")
CV      <- c(0.25, 0.20) # Cmax, AUC
theta0  <- 0.95          # both metrics
theta1  <- 0.80
theta2  <- 1.25
target  <- 0.80
df      <- data.frame(metric = metrics, theta0 = theta0,
                      CV = CV, n = NA, power = NA)
for (i in 1:nrow(df)) {
  df[i, 4:5] <- sampleN.TOST(CV = CV[i],
                             theta0 = theta0,
                             theta1 = theta1, theta2 = theta2,
                             targetpower = target,
                             design = design, print = FALSE)[7:8]
}
df$power <- signif(df$power, 6)
if (theta0 < 1) {
  res <- uniroot(opt, tol = 1e-8,
                 interval = c(theta1 + 1e-4, theta0))
} else {
  res <- uniroot(opt, tol = 1e-8,
                 interval = c(theta0, theta2 - 1e-4))
}
res     <- unlist(res)
theta0s <- c(res[["root"]], 1/res[["root"]])
txt     <- paste0("Target power for ", metrics[2],
                 " and sample size ",
                 df$n[1], "\nachieved for theta0 ",
                 sprintf("%.4f", theta0s[1]), " or ",
                 sprintf("%.4f", theta0s[2]), ".\n")
print(df, row.names = FALSE)
cat(txt)

#  metric theta0   CV  n    power
#    Cmax   0.95 0.25 27 0.803494
#     AUC   0.95 0.20 18 0.808949
# Target power for AUC and sample size 27
# achieved for theta0 0.9164 or 1.0912.

That means, although we assumed for AUC the same T/R-ratio as for C_max – with the sample size of 27 required for C_max – for AUC it can be as low as ~0.92 or as high as ~1.09, which is a soothing side-effect.

Furthermore, sometimes we have less data of AUC than of C_max (samples at the end of the profile missing or unreliable \(\small{\widehat{\lambda}_\textrm{z}}\) in some subjects and therefore, less data of AUC_0–∞ than of AUC_0–t). Again, it will not hurt because for the originally assumed T/R-ratio we need only 18 subjects.

Since – as a one-point metric – C_max is inherently more variable than AUC, Health Canada does not require assessment of its confidence interval – only the point estimate has to lie within 80.0 – 125.0%. We can explore that by setting alpha = 0.5 for it.²⁸

metrics <- c("Cmax", "AUC")
CV      <- c(0.60, 0.20)
theta0  <- 0.95
target  <- 0.80
alpha   <- c(0.50, 0.05)
df      <- data.frame(metric = metrics, CV = CV, theta0 = theta0,
                      alpha = alpha, n = NA)
for (i in 1:nrow(df)) {
  df$n[i] <- sampleN.TOST(alpha = alpha[i],
                          CV = CV[i],
                          theta0 = theta0,
                          targetpower = target,
                          design = design,
                          print = FALSE)[["Sample size"]]
}
txt       <- paste0("Sample size based on ",
                     df$metric[df$n == max(df$n)], ".\n")
print(df, row.names = FALSE)
cat(txt)

#  metric  CV theta0 alpha  n
#    Cmax 0.6   0.95  0.50 24
#     AUC 0.2   0.95  0.05 18
# Sample size based on Cmax.

Only with a relatively high CV (compared to the one of AUC) you may face a situation where you have to base the sample size on C_max.

top of section ↩︎ previous section ↩︎

Multiple Treatments

First of all we have have to consider the purpose of the study and how we will evaluate it (see above).

The well-known Bon­ferroni-adjustment can be used \[\begin{equation}\tag{6} \alpha_\textrm{adj}=\alpha\,/\,k \end{equation}\] where \(\small{k}\) is the number of simultaneous tests. The Type I Error \(\small{TIE}\) will always be controlled because \[\begin{equation}\tag{7} TIE=1-{(1-\alpha_\textrm{adj})}^k \end{equation}\] For \(\small{\alpha=0.05}\) and \(\small{k=2}\) we get \(\small{0.049375<\alpha}\).

Of course, all comparisons in the study should be done by \(\small{100\,(1-2\,\alpha_\textrm{adj})}\) confidence intervals.

Extending the basic example and using the argument alpha:

k         <- 2        # comparisons
alpha.adj <- 0.05 / k # Bonferroni
sampleN.TOST(alpha = alpha.adj, CV = 0.25,
             theta0 = 0.95, targetpower = 0.80,
             design = design)

# 
# +++++++++++ Equivalence test - TOST +++++++++++
#             Sample size estimation
# -----------------------------------------------
# Study design: 3x3 crossover 
# log-transformed data (multiplicative model)
# 
# alpha = 0.025, target power = 0.8
# BE margins = 0.8 ... 1.25 
# True ratio = 0.95,  CV = 0.25
# 
# Sample size (total)
#  n     power
# 36   0.827798

The sample size increased by ~33% from the 27 subjects required in a single comparison.

With an increasing number of comparisons, the Bon­ferroni-adjustment quickly becomes overly conservative. Another issue is that the tests are not strictly independent (at least when we compare different tests to the same reference treatment some – unknown – correlation exists). More powerful methods (Holm, Hoch­berg, …) are out of scope of this article.

Three treatments intended for evaluation ‘All at Once’ or ‘Two at a Time’. Equivalence is seeked for both (say, a tablet and a capsule vs a reference).

CV <- 0.25
df <- data.frame(design = c("3x3", "3x6x3",
                            "2x2x2"),
                 n = NA, power = NA)
for (i in 1:nrow(df)) {
  df[i, 2:3] <- sampleN.TOST(alpha = 0.05 / 2,
                             CV = CV,
                             design = df$design[i],
                             print = FALSE)[7:8]
}
df$power <- signif(df$power, 4)
print(df, row.names = FALSE)

#  design  n  power
#     3x3 36 0.8278
#   3x6x3 36 0.8278
#   2x2x2 36 0.8161

Four treatments (e.g., Test fasting and fed state, Reference fasting and fed state). Two comparisons are confirmatory (Test vs Reference fasting, Test vs Reference fed) and two are exploratory (Test fed vs fasting, Reference fed vs fasting). We have to adjust only for the former since the latter are only ‘nice to know’.

CV <- 0.25
df <- data.frame(design = c("4x4", "2x2x2"),
                 n = NA, power = NA)
for (i in 1:nrow(df)) {
  df[i, 2:3] <- sampleN.TOST(alpha = 0.05 / 2,
                             CV = CV,
                             design = df$design[i],
                             print = FALSE)[7:8]
}
df$power <- signif(df$power, 4)
print(df, row.names = FALSE)

#  design  n  power
#     4x4 36 0.8315
#   2x2x2 36 0.8161

top of section ↩︎ previous section ↩︎

NTIDs

So far we employed the common (and hence, default) BE-limits theta1 = 0.80 and theta2 = 1.25. In some jurisdictions tighter limits of 90.00 – 111.11% have to be used.²⁹

Generally NTIDs show a low within-subject variability (though the between-subject CV can be much higher – these drugs require quite often dose-titration).

You have to provide only the lower BE-limit theta1 (the upper one will be automatically calculated). In the examples I used a lower CV of 0.125.

sampleN.TOST(CV = 0.125, theta1 = 0.90,
             design = design)

# 
# +++++++++++ Equivalence test - TOST +++++++++++
#             Sample size estimation
# -----------------------------------------------
# Study design: 3x3 crossover 
# log-transformed data (multiplicative model)
# 
# alpha = 0.05, target power = 0.8
# BE margins = 0.9 ... 1.111111 
# True ratio = 0.95,  CV = 0.125
# 
# Sample size (total)
#  n     power
# 69   0.813991

Not so nice. If you would exhaust the capacity of the clinical site, you may consider a replicate design.

sampleN.TOST(CV = 0.125, theta1 = 0.90, theta2 = 1.12,
             design = design)

# 
# +++++++++++ Equivalence test - TOST +++++++++++
#             Sample size estimation
# -----------------------------------------------
# Study design: 3x3 crossover 
# log-transformed data (multiplicative model)
# 
# alpha = 0.05, target power = 0.8
# BE margins = 0.9 ... 1.12 
# True ratio = 0.95,  CV = 0.125
# 
# Sample size (total)
#  n     power
# 69   0.813991

get.sequences <- function(design) {
  known <- known.designs()[, c(2, 5)]
  ns    <- as.integer(known[known$design == design,
                            "steps"])
  return(ns)
}
make.equal <- function(n, ns) {
  return(as.integer(ns * (n %/% ns + as.logical(n %% ns))))
}
design  <- "3x3"
target  <- 0.80
theta0  <- 0.975
CV.min  <- 0.075
CV.max  <- 0.20
CV.step <- 1000
CV      <- seq(CV.min, CV.max, length.out = CV.step)
ns      <- get.sequences(design)
res     <- data.frame(CV = CV, n.EMA = NA_integer_,
                      n.HC = NA_integer_, less = FALSE)
for (i in 1:nrow(res)) {
  n.TaT <- sampleN.TOST(CV = CV[i], theta0 = theta0,
                        targetpower = target, design = "2x2",
                        theta1 = 0.90, theta2 = 1/0.90,
                        print = FALSE)[["Sample size"]]
  res$n.EMA[i] <- make.equal(n = n.TaT, ns = ns)
  n.TaT <- sampleN.TOST(CV = CV[i], theta0 = theta0,
                        targetpower = target, design = "2x2",
                        theta1 = 0.90, theta2 = 1.12,
                        print = FALSE)[["Sample size"]]
  res$n.HC[i]  <- make.equal(n = n.TaT, ns = ns)
  if (res$n.HC[i] < res$n.EMA[i]) res$less[i] = TRUE
}
cat("design                      :", design,
    "(intended for \u2018Two at a Time\u2019 approach)",
    "\ntarget power                :", target,
    "\ntheta0                      :", theta0,
    "\nCV                          :", CV.min, "\u2013", CV.max,
    "\nSample size for HC < common :",
    sprintf("%.2f%%", 100 * sum(res$less) / CV.step),
    "of cases.\n")

# design                      : 3x3 (intended for ‘Two at a Time’ approach) 
# target power                : 0.8 
# theta0                      : 0.975 
# CV                          : 0.075 – 0.2 
# Sample size for HC < common : 12.90% of cases.

The FDA requires for NTIDs tighter batch-release specifications (±5% instead of ±10%). Let’s hope that your product complies and the T/R-ratio will be closer to 1:

sampleN.TOST(CV = 0.125,
             theta0 = 0.975, # ‘better’ T/R-ratio
             theta1 = 0.90,
             design = design)

# 
# +++++++++++ Equivalence test - TOST +++++++++++
#             Sample size estimation
# -----------------------------------------------
# Study design: 3x3 crossover 
# log-transformed data (multiplicative model)
# 
# alpha = 0.05, target power = 0.8
# BE margins = 0.9 ... 1.111111 
# True ratio = 0.975,  CV = 0.125
# 
# Sample size (total)
#  n     power
# 33   0.820725

Substantially lower sample and doable.

top of section ↩︎ previous section ↩︎

Difference of Means

Sometimes we are interested in assessing differences of responses instead of their ratios. In such a case we have to set logscale = FALSE. The limits theta1 and theta2 can be expressed in the following ways:

• As a difference of means relative to the same (underlying) reference mean.
• In units of the difference of means.

Note that in the former case the units of CV, and theta0 need also to be given relative to the reference mean (specified as ratio).

Let’s estimate the sample size for an equivalence trial of two blood pressure lowering drugs assessing the difference in means of untransformed data (raw, linear scale). In this setup everything has to be given with the same units (i.e., here the assumed difference –5 mm Hg and the lower / upper limits –15 mm Hg / +15 mm Hg systolic blood pressure). Furthermore, we assume a CV of 25 mm Hg.

sampleN.TOST(CV = 20, theta0 = -5,
             theta1 = -15, theta2 = +15,
             design = design,
             logscale = FALSE)

# 
# +++++++++++ Equivalence test - TOST +++++++++++
#             Sample size estimation
# -----------------------------------------------
# Study design: 3x3 crossover 
# untransformed data (additive model)
# 
# alpha = 0.05, target power = 0.8
# BE margins = -15 ... 15 
# True diff. = -5,  CV = 20
# 
# Sample size (total)
#  n     power
# 51   0.805426

Sometimes in the literature we find not the CV but the standard deviation of the difference. Say, it is given with 36 mm Hg. We have to convert it to a CV. In all implemented Higher-Order Crossover designs it’s \(\small{SD/\sqrt{2}}\).

SD.delta <- 36
sampleN.TOST(CV = SD.delta / sqrt(2), theta0 = -5,
             theta1 = -15, theta2 = +15,
             design = design,
             logscale = FALSE)

# 
# +++++++++++ Equivalence test - TOST +++++++++++
#             Sample size estimation
# -----------------------------------------------
# Study design: 3x3 crossover 
# untransformed data (additive model)
# 
# alpha = 0.05, target power = 0.8
# BE margins = -15 ... 15 
# True diff. = -5,  CV = 25.45584
# 
# Sample size (total)
#  n     power
# 81   0.800354

Note that other software packages (e.g., PASS, nQuery, StudySize, …) require the standard deviation of the difference as input.

previous section ↩︎

Ratio of Means

It would be a fundamental flaw if responses are assumed to a follow a normal distribution and not – like above – assess their difference and at the end instead giving an estimate of the difference, calculate the ratio of the test- and reference-means.

In such a case Fieller’s (‘fiducial’) confidence interval³⁰ has to be employed.
Note that in the functions sampleN.RatioF()and power.RatioF() alpha = 0.025 is the default, since it is intended for studies with clinical endpoints, where the 95% confidence interval is usually used for equivalence testing.³¹ Only the ‘Two at a Time’ approach is implemented, i.e., you have to specify design = "2x2". Sometimes we have to round the sample size up (see the interlude). Note that we need additionally the between-subject CV (here arbitrarily twice the CV).

planned  <- "3x3"
CV       <- 0.25
CVb      <- 2 * CV
get.sequences <- function(design) {
  known <- known.designs()[, c(2, 5)]
  ns    <- as.integer(known[known$design == design,
                            "steps"])
  return(ns)
}
make.equal <- function(n, ns) {
  return(as.integer(ns * (n %/% ns + as.logical(n %% ns))))
}
ns <- get.sequences(planned)
n  <- sampleN.RatioF(CV = CV, CVb = CVb, design = "2x2",
                     print = FALSE)[["Sample size"]]
n  <- make.equal(n = n, ns = ns)
df <- data.frame(Design = planned, Eval = "2x2",
                 alpha = 0.025, CV = CV, CVb = CVb, n = n)
df$power <- suppressMessages(
              power.RatioF(CV = CV, CVb = CVb,
                           design = "2x2", n = n))
names(df)[6:7] <- c("Sample size", "Achieved power")
print(df, row.names = FALSE)

#  Design Eval alpha   CV CVb Sample size Achieved power
#     3x3  2x2 0.025 0.25 0.5          42      0.8021525

previous section ↩︎