Preliminaries

A basic knowledge of R is required. To run the scripts at least version 1.2.8 (2015-07-10) of PowerTOST is required and 1.5.3 (2021-01-18) suggested. Any version of R would likely do, though the current release of PowerTOST was only tested with version 4.1.3 (2022-03-10) and later.
All scripts were run on a Xeon E3-1245v3 @ 3.40GHz (1/4 cores) 16GB RAM with R 4.3.1 on Win­dows 7 build 7601, Service Pack 1, Universal C Run­time 10.0.10240.16390.

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

sampleN.NTID() gives balanced sequences (i.e., an equal number of subjects is allocated to all sequences). Furthermore, the estimated sample size is the total number of subjects.

All examples deal with studies where the response variables likely follow a lognormal 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).

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 (in BE commonly 0.05),
2. the (potentially) narrowed BE-margin,
3. the desired (or target) power,
4. the variance (commonly expressed as a coefficient of variation), and
5. the deviation of the test from the reference treatment,

where

1. is fixed by the agency,
2. depends on the observed⁹ CV_wR (which is unknown beforehand),
3. is set by the sponsor (commonly to 0.80 – 0.90),
4. is an uncertain assumption,
5. is an uncertain assumption.

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 \(\small{CV_\textrm{w}}\) is smaller than the between-subject variability \(\small{CV_\textrm{b}}\), crossover studies are so popular.
Of note, there is no relationship between \(\small{CV_\textrm{w}}\) and \(\small{CV_\textrm{b}}\). An example are drugs which are subjected to polymorphic metabolism. For these drugs \(\small{CV_\textrm{w}\ll CV_\textrm{b}}\).
Furthermore, the drugs’ within-subject variability may be unequal (i.e., \(\small{s_{\textrm{wT}}^{2}\neq s_{\textrm{wR}}^{2}}\)). For details see the section about He­te­ro­scedasticity.

It is a prerequisite that no (unequal) carryover from one period to the next exists. Only then the comparison of treatments will be unbiased.¹² Carryover is elaborated in another article.

top of section ↩︎ previous section ↩︎

A Simple Case

We assume a CV of 0.125, a T/R-ratio of 0.975, a target a power of 0.80, and want to perform the study in a 2-se­quence 4-period full replicate study (TRTR|RTRT).

Since theta0 = 0.975,^{16 17} targetpower = 0.80, and design = "2x2x4" 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\%}\) CI). Hence, you need to specify only the CV (assuming \(\small{CV_\textrm{wT}=CV_\textrm{wR}}\)).
To shorten the output, use the argument details = FALSE.

sampleN.NTID(CV = 0.125, details = FALSE)

# 
# +++++++++++ FDA method for NTIDs ++++++++++++
#            Sample size estimation
# ---------------------------------------------
# Study design:  2x2x4 (TRTR|RTRT) 
# log-transformed data (multiplicative model)
# 1e+05 studies for each step simulated.
# 
# alpha  = 0.05, target power = 0.8
# CVw(T) = 0.125, CVw(R) = 0.125
# True ratio     = 0.975 
# ABE limits     = 0.8 ... 1.25 
# Regulatory settings: FDA 
# 
# Sample size
#  n     power
# 16   0.822780

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 stating the additional argument print = FALSE and assign the result to a data frame (here named x).

x <- sampleN.NTID(CV = 0.125, details = FALSE, print = FALSE)

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

Let’s retrieve the column names of x:

names(x)
#  [1] "Design"         "alpha"          "CVwT"           "CVwR"          
#  [5] "theta0"         "theta1"         "theta2"         "Sample size"   
#  [9] "Achieved power" "Target power"   "nlast"

Now we can access the elements of x by their names. Note that double square brackets [["foo"]] have to be used, where "foo" is the name of the element.

x[["Sample size"]]
# [1] 16
x[["Achieved power"]]
# [1] 0.82278

If you insist in accessing elements by column-numbers, use single square brackets [bar], where bar is the respective number. You could also access multiple columns, e.g., [foo:bar] or even as a vector [c(foo, bar:baz)].

x[8:9]
#   Sample size Achieved power
# 1          16        0.82278
x[c(4, 8:9)]
#    CVwR Sample size Achieved power
# 1 0.125          16        0.82278

With 16 subjects (eight per sequence) we achieve the power we desire.

What will happen if we have one dropout?

power.NTID(CV = 0.125, n = x[["Sample size"]] - 1)

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

# [1] 0.78896

Slightly below the 0.80 we desire.

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

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 twelve subjects and estimated the CV as 0.125.

A two-sided confidence interval of the CV can be obtained via the \(\small{\chi^2}\)-distribution of its error variance \(\small{\sigma^2}\) with \(\small{n-2}\) degrees of freedom. \[\eqalign{\tag{7} 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\} }\]

Let’s calculate the 95% confidence interval of the CV to get an idea.

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

# lower CL upper CL 
#  0.08717  0.22120

Surprised? Although 0.125 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 get a higher CV than a lower one in the planned study.

If we plan the study based on 0.125, we would opt for 16 subjects like in the example before (not adjusted for the dropout-rate).
If the CV will be lower, we loose power (narrower limits). But what if it will be higher? Depends. Since we have to narrow the limits less, we gain power. But how much?

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

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

#       CV   power
#  0.08717 0.75885
#  0.12500 0.82278
#  0.22120 0.85300

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. What happens when we plan the study based on the confidence interval of the CV?

m               <- seq(6, 12, 2)
d.f             <- data.frame(n.pilot = m, CV = 0.125,
                              l = NA_real_, u = NA_real_,
                              n.low = NA_integer_, n.CV = NA_integer_,
                              n.hi = NA_integer_)
for (i in 1:nrow(d.f)) {
  d.f[i, 3:4] <- CVCL(CV = 0.125, df = m[i] - 2,
                     side = "2-sided", alpha = 0.05)
  d.f[i, 5] <- sampleN.NTID(CV = d.f$l[i], details = FALSE,
                            print = FALSE)[["Sample size"]]
  d.f[i, 6] <- sampleN.NTID(CV = 0.125, details = FALSE,
                            print = FALSE)[["Sample size"]]
  d.f[i, 7] <- sampleN.NTID(CV = d.f$u[i], details = FALSE,
                            print = FALSE)[["Sample size"]]
}

d.f[, 3:4]      <- signif(d.f[, 3:4], 4)
names(d.f)[3:4] <- c("lower CL", "upper CL")
print(d.f, row.names = FALSE)

#  n.pilot    CV lower CL upper CL n.low n.CV n.hi
#        6 0.125  0.07471   0.3696    20   16   26
#        8 0.125  0.08037   0.2794    20   16   18
#       10 0.125  0.08425   0.2420    18   16   16
#       12 0.125  0.08717   0.2212    18   16   16

Small pilot studies are practically useless.

Furthermore, 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. Say, it was 0.975.

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

# [1] 0.9174

We must not forget the ratio-test comparing the drugs’ within-subject standard deviations. Explore the impact of different CVs and a relatively 2.5% lower T/R-ratio on power for the given sample size.

n       <- 16
CV      <- 0.125
ratio   <- seq(1, 2, 0.5) # variance ratio (T/R)
theta0  <- 0.975
r_const <- -log(0.90) / 0.1
comp1   <- data.frame(CVwT = CVp2CV(0.125, ratio)[, 1],
                      CVwR = CVp2CV(0.125, ratio)[, 2],
                      s.ratio.CL =  NA_real_,
                      L = NA_real_, U = NA_real_,
                      power = NA_real_, SABE = NA_real_,
                      s.ratio = NA_real_)
comp2  <- data.frame(CV = CV, theta0 = c(theta0, theta0 * 0.975),
                     power = NA_real_)
for (i in seq_along(ratio)) {
  comp1$s.ratio.CL[i] <- CV2se(comp1$CVwT[i])/CV2se(comp1$CVwR[i]) /
                         sqrt(qf(1-0.1/2, df1 = n -2 , df2 = n -2,
                                 lower.tail = FALSE))
  comp1[i, 4:5] <- round(100 * exp(c(-1, +1) *
                                   r_const * CV2se(comp1$CVwR[i])), 2)
  tmp <- suppressMessages(
           power.NTID(CV = c(comp1$CVwT[i], comp1$CVwR[i]),
                      theta0 = theta0, n = n, details = TRUE))
  comp1$power[i]   <- tmp[["p(BE)"]]
  comp1$SABE[i]    <- tmp[["p(BE-sABEc)"]]
  comp1$s.ratio[i] <- tmp[["p(BE-sratio)"]]
}
comp1 <- signif(comp1, 5)
names(comp1)[7:8] <- c("p(SABE)", "p(s-ratio)")
for (i in 1:2) {
  comp2$power[i] <- power.NTID(CV = CV, theta0 = comp2$theta0[i],
                               n = n)
}
comp2 <- signif(comp2, 5)
print(comp1, row.names = FALSE)
print(comp2, row.names = FALSE)

#     CVwT    CVwR s.ratio.CL     L      U   power p(SABE) p(s-ratio)
#  0.12500 0.12500     1.5760 87.71 114.02 0.82278 0.84869    0.95128
#  0.13704 0.11172     1.9302 88.93 112.45 0.67124 0.76335    0.82792
#  0.14452 0.10193     2.2288 89.84 111.31 0.50581 0.67983    0.66365
#     CV  theta0   power
#  0.125 0.97500 0.82278
#  0.125 0.95062 0.64342

If \(\small{CV_\textrm{wT}>CV_\textrm{wR}}\) (with \(\small{CV_\textrm{w}=\textrm{const}=0.125}\)) we are in trouble. The scaled limits get narrower and the impact of the ratio-test becomes increasingly important. With \(\small{\sigma_\textrm{wT}^2/\sigma_\textrm{wR}^2=2}\) the upper confidence limit of \(\small{s_\textrm{wT}/s_\textrm{wR}}\) is with ≈2.23 already close to the limit of 2.5.
If the T/R-ratio gets even just a little worse, the impact on power is substantial.

Essentially this leads to the murky waters of prospective power and sensitivity analyses.
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.125
theta0   <- 0.975
target   <- 0.80
minpower <- 0.70
pa       <- pa.NTID(CV = CV, theta0 = theta0,
                    targetpower = target,
                    minpower = minpower)
change.CV     <- 100*(tail(pa$paCV[["CV"]], 1) -
                      pa$plan[["CVwR"]]) /
                      pa$plan[["CVwR"]]
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 NTID
#  Design alpha  CVwT  CVwR theta0 theta1 theta2 Sample size
#   2x2x4  0.05 0.125 0.125  0.975    0.8   1.25          16
#  Achieved power Target power
#         0.82278          0.8
# 
# Power analysis
# CV, theta0 and number of subjects leading to min. acceptable power of ~0.7:
#  CV = (0.0716, 0.3059), theta0= 0.9569
#  n = 13 (power= 0.7064)
# 
#  parameter relative change
#         CV        +144.76%
#     theta0          -1.86%
#          n         -18.75%

Confirms what we have seen above. As expected, the method is extremely robust to changes in the CV. The sample size is also not very sensitive; many overrate the impact of dropouts on power. As usual in bioequivalence, changes in the T/R-ratio hurt most.

top of section ↩︎ previous section ↩︎

Heteroscedasticity

Since a comparison of \(\small{s_\textrm{wT}}\) with \(\small{s_\textrm{wR}}\) is mandatory, already pilot studies have to be performed in one of the fully replicated designs. If you are concerned about dropouts or the bioanalytical method requires large sample volumes, opt for one of the 2-sequence 3-period designs (TRT|RTR or TRR|RTT).
Since pharmaceutical technology – hopefully – improves, it is not uncommon that \(\small{s_\textrm{wT}<s_\textrm{wR}}\). If this is the case, you get an incentive in the sample size of the pivotal study; scaling the limits is based on \(\small{s_\textrm{wR}}\) but the 90% CI in the ABE-component on the (pooled) \(\small{s_\textrm{w}^{2}}\). Furthermore, the probability to pass the ratio-test increases.

\[\eqalign{\tag{8} s_\textrm{wT}^{2}&=\log_{e}(CV_\textrm{wT}^{2}+1)\\ s_\textrm{wR}^{2}&=\log_{e}(CV_\textrm{wR}^{2}+1)\\ s_\textrm{w}^{2}&=\left(s_\textrm{wT}^{2}+s_\textrm{wR}^{2}\right)/2\\ CV_\textrm{w}&=\sqrt{\exp(s_\textrm{w}^{2})-1}}\]

Note that when giving CV as a two-element vector in sampleN.NTID(), the first element has to be \(\small{CV_\textrm{wT}}\) and the second \(\small{CV_\textrm{wR}}\).

Heteroscedasticity can be challenging, esp. for low variability and \(\small{s_\textrm{wT}^2>s_\textrm{wR}^2}\).

# Cave: long runtime
CV  <- unique(sort(c(seq(0.05, 0.30, 0.001), se2CV(0.1), 0.214188)))
df1 <- data.frame(CVT = CVp2CV(CV = CV, ratio = 2)[, 1],
                  CVR = CVp2CV(CV = CV, ratio = 2)[, 2],
                  n = NA_integer_)
df2 <- data.frame(CVT = CV, CVR = CV, n = NA_integer_)
df3 <- data.frame(CVT = CVp2CV(CV = CV, ratio = 0.5)[, 1],
                  CVR = CVp2CV(CV = CV, ratio = 0.5)[, 2],
                  n = NA_integer_)
for (i in seq_along(CV)) {
  df1$n[i] <- sampleN.NTID(CV = c(df1$CVT[i], df1$CVR[i]),
                           details = FALSE,
                           print = FALSE)[["Sample size"]]
  df2$n[i] <- sampleN.NTID(CV = CV[i], details = FALSE,
                           print = FALSE)[["Sample size"]]
  df3$n[i] <- sampleN.NTID(CV = c(df3$CVT[i], df3$CVR[i]),
                           details = FALSE,
                           print = FALSE)[["Sample size"]]
}
dev.new(width = 4.5, height = 4.5)
op <- par(no.readonly = TRUE)
par(lend = 2, ljoin = 1, mar = c(4, 3.9, 0.1, 0.1), cex.axis = 0.9)
clr <- c("red", "magenta", "blue")
plot(100 * CV, df1$n, type = "n", xlim = 100 * range(CV),
     ylim = c(12, max(c(df1$n, df2$n, df3$n))),
     axes = FALSE, xlab = expression(italic(CV)[w]*" (%)"),
     ylab = "sample size")
grid(ny = NA)
abline(h = seq(12, max(c(df1$n, df2$n, df3$n)), 12), col = "lightgrey", lty = 3)
axis(1)
axis(2, at = seq(12, max(c(df1$n, df2$n, df3$n)), 12), las = 1)
lines(100 * CV, df1$n, type = "s", lwd = 2, col = clr[1])
lines(100 * CV, df2$n, type = "s", lwd = 2, col = clr[2])
lines(100 * CV, df3$n, type = "s", lwd = 2, col = clr[3])
abline(v = 100 * c(se2CV(0.1), 0.214188), lty = 2)
legend("topright", cex = 0.9, col = clr, box.lty = 0, bg = "white",
       legend = c(expression(italic(s)[wT]^2==2%*%italic(s)[wR]^2),
                  expression(italic(s)[wT]^2==italic(s)[wR]^2),
                  expression(italic(s)[wT]^2==0.5%*%italic(s)[wR]^2)),
       inset = 0.02, lwd = 2, y.intersp = 1.35)
box()
par(op)

Fig. 3 Sample sizes under homo- and heteroscedasticity (2×2×4 design).
Dashed lines at the switching CV and the implied upper cap.

For our example we require 16 subjects under homoscedasticity (\(\small{CV_\textrm{wT}\equiv CV_\textrm{wR}=12.5\%}\)) but 28 if \(\small{s_\textrm{wT}^2=2\times s_\textrm{wR}^2}\). On the other hand, if \(\small{s_\textrm{wT}^2=0.5\times s_\textrm{wR}^2}\), we require only twelve.
Starting with the implied upper cap at 21.4188% sample sizes increase because the additional condition ‘must pass ABE’ sets in.

top of section ↩︎ previous section ↩︎

Multiple Endpoints

In demonstrating bioequivalence the pharmacokinetic metrics C_max, AUC_0–t, and AUC_0–∞ are mandatory.

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

We design the study always for the worst case combination, i.e., based on the PK metric requiring the largest sample size. Let’s explore that with different CVs and T/R-ratios.

metrics <- c("Cmax", "AUCt", "AUCinf")
CV      <- c(0.125, 0.100, 0.110)
theta0  <- c(0.975, 0.970, 0.980)
d.f     <- data.frame(metric = metrics, CV = CV, theta0 = theta0,
                      n = NA_integer_, power = NA_real_)
for (i in seq_along(metrics)) {
  d.f[i, 4:5] <- sampleN.NTID(CV = CV[i], theta0 = theta0[i],
                              details = FALSE, print = FALSE)[8:9]
}
d.f$power <- signif(d.f$power, 5)
txt       <- paste0("Sample size based on ",
                    d.f$metric[d.f$n == max(d.f$n)], ".\n")
print(d.f, row.names = FALSE)
cat(txt)

#  metric    CV theta0  n   power
#    Cmax 0.125  0.975 16 0.82278
#    AUCt 0.100  0.970 18 0.80051
#  AUCinf 0.110  0.980 16 0.83411
# Sample size based on AUCt.

The PK metric with the lowest variability drives the sample size.

Let us assume for simplicity the same T/R-ratio of 0.975 for all metrics. Which are the extreme T/R-ratios (largest deviations of T from R) for C_max giving still the target power?

opt <- function(x) {
  power.NTID(theta0 = x, CV = d.f$CV[1],
             n = d.f$n[2]) - target
}
metrics <- c("Cmax", "AUCt", "AUCinf")
CV      <- c(0.125, 0.100, 0.110)
theta0  <- 0.975
target  <- 0.80
d.f     <- data.frame(metric = metrics, CV = CV,
                      n = NA_integer_, power = NA_real_)
for (i in seq_along(metrics)) {
  d.f[i, 3:4] <- sampleN.NTID(CV = CV[i], theta0 = theta0,
                              details = FALSE, print = FALSE)[8:9]
}
d.f$power <- signif(d.f$power, 5)
if (theta0 < 1) {
  res <- uniroot(opt, tol = 1e-8,
                 interval = c(0.80 + 1e-4, theta0))
} else {
  res <- uniroot(opt, tol = 1e-8,
                 interval = c(theta0, 1.25 - 1e-4))
}
res     <- unlist(res)
theta0s <- c(res[["root"]], 1/res[["root"]])
txt     <- paste0("Target power for ", metrics[1],
                  " and sample size ",
                  d.f$n[2], "\nachieved for theta0 ",
                  sprintf("%.4f", theta0s[1]), " or ",
                  sprintf("%.4f", theta0s[2]), ".\n")
print(d.f, row.names = FALSE)
cat(txt)

#  metric    CV  n   power
#    Cmax 0.125 16 0.82278
#    AUCt 0.100 18 0.84179
#  AUCinf 0.110 16 0.80470
# Target power for Cmax and sample size 18
# achieved for theta0 0.9626 or 1.0388.

That means, although we assumed for C_max the same T/R-ratio as for AUC, with the sample size of 18 required for AUC_t, for C_max it can be as low as ~0.963 or as high as ~1.039, which is an interesting side-effect.

top of section ↩︎ previous section ↩︎