Preliminaries

A basic knowledge of R is required. To run the scripts at least version 1.4.8 (2019-08-29) of PowerTOST is required and at least version 1.5.3 (2021-01-18) suggested. Any version of R would likely do, though the cur­rent release of PowerTOST was only tested with R 4.3.3 (2024-02-29) and later. All scripts were run on a Xeon E3-1245v3 @ 3.40GHz (1/4 cores) 16GB RAM with R 4.4.1 on Windows 7 build 7601, Service Pack 1.

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.scABEL() 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 – and not the number of subjects per sequence like in some other software packages.

All 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).

top of section ↩︎ 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 (potentially expanded) BE-margin,
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,

where

1. is fixed by the agency,
2. depends on the observed⁴⁰ \(\small{CV_\text{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. estimated, post hoc, a posteriori) power is not only futile but plain nonsense.⁴² For details see another article.

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

Power is a hypersurface depending on \(\small{\theta_0}\), \(\small{n}\), and \(\small{CV}\). Don’t try to imagine four-dimensional space – it may fry your brain! We can only show spatial projections at a certain \(\small{CV\text{-}}\)value.

Fig. 4 Spatial projection of power for \(\small{\theta_{1,2}=0.7215-1.3859,\,\log_{e}\theta_{1,2}=\mp 0.32636}\), \(\small{n=22-40}\), and \(\small{CV_\text{wT}=CV_\text{wR}=0.45}\). Four period full replicate design for the EMA’s ABEL.

Inspecting the power-surface for a given \(\small{n}\), we see the maximum at \(\small{\theta_0=1}\). Due to the symmetry in \(\small{\log_{e}}\)-scale we require the same sample size for any given \(\small{\theta_0}\) and \(\small{1/\theta_0}\). Furthermore, at the expanded limits of the acceptance range (\(\small{\theta_1=0.7215,\theta_2=1.3859}\)) power \(\small{\leq 0.05}\), demonstrating that the Type I Error is controlled.⁴⁴
For any given \(\small{\theta_0}\) power increases with \(\small{n}\) (asympotically approaching \(\small{1}\)).
Intersections with the horizontal power-planes (here \(\small{0.7,\,0.8,\,0.9}\)) show the U-shaped dependency of \(\small{n}\) on \(\small{\theta_0}\).

Note that the power-surface is ‘flatter’ than in the corresponding example of ABE.

top of section ↩︎ previous section ↩︎

A Simple Case

We assume a CV of 0.45, a T/R-ratio of 0.90, a target a power of 0.80, and want to perform the study in a two sequence four period full replicate study (TRTR|RTRT or TRRT|RTTR or TTRR|RRTT) for the EMA’s ABEL.

Since theta0 = 0.90,⁴⁷ targetpower = 0.80, and regulator = "EMA" 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). Hence, you need to specify only the CV (assuming \(\small{CV_\text{wT}=CV_\text{wR}}\)) and design = "2x2x4".
To shorten the output, use the argument details = FALSE.

sampleN.scABEL(CV = 0.45, design = "2x2x4", details = FALSE)

# 
# +++++++++++ scaled (widened) ABEL +++++++++++
#             Sample size estimation
#    (simulation based on ANOVA evaluation)
# ---------------------------------------------
# Study design: 2x2x4 (4 period full replicate)
# log-transformed data (multiplicative model)
# 1e+05 studies for each step simulated.
# 
# alpha  = 0.05, target power = 0.8
# CVw(T) = 0.45; CVw(R) = 0.45
# True ratio = 0.9
# ABE limits / PE constraint = 0.8 ... 1.25 
# Regulatory settings: EMA 
# 
# Sample size
#  n     power
# 28   0.8112

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 x).

x <- sampleN.scABEL(CV = 0.45, design = "2x2x4",
                    details = FALSE, print = FALSE)

Let’s retrieve the column names of x:

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

Now we can access the elements of x by their names. Note that double square brackets and double or single quotes ([["…"]], [['…']]) have to be used.

x[["Sample size"]]
# [1] 28
x[['Achieved power']]
# [1] 0.81116

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[8:9]
#   Sample size Achieved power
# 1          28        0.81116

With 28 subjects (14 per sequence) we achieve the power we desire.

What happens if we have one dropout?

power.scABEL(CV = 0.45, design = "2x2x4",
             n = x[["Sample size"]] - 1)

# Unbalanced design. n(i)=14/13 assumed.

# [1] 0.79848

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 16 subjects and estimated the CV as 0.45.

The \(\small{\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. \[\eqalign{\tag{6} s_\text{w}^2&=\log_{e}(CV_\text{w}^2+1)\\ L=\frac{(n-1)\,s_\text{w}^2}{\chi_{\alpha/2,\,n-2}^{2}}&\leq\sigma_\text{w}^2\leq\frac{(n-1)\,s_\text{w}^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  <- 16 # pilot study
ci <- CVCL(CV = 0.45, df = m - 2,
           side = "2-sided", alpha = 0.05)
ci <- c(ci, mean(ci))
names(ci)[3] <- "midpoint"
signif(ci, 4)

# lower CL upper CL midpoint 
#   0.3223   0.7629   0.5426

Surprised? Although 0.45 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 that we will face a higher CV than a lower one in the planned study.

If we plan the study based on 0.45, we would opt for 28 subjects like in the examples before (not adjusted for the dropout-rate).
If the CV will be lower, we loose power (less expansion). But what if it will be higher? Depends. Since we may expand more, we gain power. However, if we cross the upper cap of scaling (50% for the EMA), we will loose power. But how much?

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

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

#      CV   power
#  0.3223 0.73551
#  0.4500 0.81116
#  0.5426 0.79868
#  0.7629 0.60158

Might hurt.

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(12, 30, 6)
df              <- data.frame(n.pilot = m, CV = 0.45,
                              l = NA, u = NA,
                              n.low = NA, n.CV = NA, n.hi = NA)
for (i in 1:nrow(df)) {
  df[i, 3:4] <- CVCL(CV = 0.45, df = m[i] - 2,
                     side = "2-sided",
                     alpha = 0.05)
  df[i, 5] <- sampleN.scABEL(CV = df$l[i], design = "2x2x4",
                             details = FALSE,
                             print = FALSE)[["Sample size"]]
  df[i, 6] <- sampleN.scABEL(CV = 0.45, design = "2x2x4",
                             details = FALSE,
                             print = FALSE)[["Sample size"]]
  df[i, 7] <- sampleN.scABEL(CV = df$u[i], design = "2x2x4",
                             details = FALSE,
                             print = FALSE)[["Sample size"]]
}

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

#  n.pilot   CV lower CL upper CL n.low n.CV n.hi
#       12 0.45   0.3069   0.8744    36   28   56
#       18 0.45   0.3282   0.7300    34   28   42
#       24 0.45   0.3415   0.6685    34   28   38
#       30 0.45   0.3509   0.6334    34   28   34

One leading generic company has an internal rule to perform pilot studies of HVDPs in a four period full replicate design and at least 24 subjects. Makes sense.

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.90.

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

# [1] 0.7515

Explore the impact of a relatively 5% lower CV (less expansion) and a relatively 5% lower T/R-ratio on power for the given sample size.

n      <- 28
CV     <- 0.45
theta0 <- 0.90
comp1  <- data.frame(CV = c(CV, CV*0.95),
                     power = NA)
comp2  <- data.frame(theta0 = c(theta0, theta0*0.95),
                     power = NA)
for (i in 1:2) {
  comp1$power[i] <- power.scABEL(CV = comp1$CV[i],
                                 theta0 = theta0,
                                 design = "2x2x4",
                                 n = n)
}
comp1$power <- signif(comp1$power, 5)
for (i in 1:2) {
  comp2$power[i] <- power.scABEL(CV = CV,
                                 theta0 = comp2$theta0[i],
                                 design = "2x2x4",
                                 n = n)
}
comp2$power <- signif(comp2$power, 5)
print(comp1, row.names = FALSE)
print(comp2, row.names = FALSE)

#      CV   power
#  0.4500 0.81116
#  0.4275 0.80095
#  theta0   power
#   0.900 0.81116
#   0.855 0.61952

---

Interlude 2

    

Fig. 4 Power curves for n = 28 (2×2×4 design, EMA ABEL).

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.90 and 1.1111). Contrary to ABE, power is maintained, unless we cross the upper scaling limit, where additionally the PE-constraint becomes increasingly important.

<nitpick>

• A common flaw in protocols is the phrase
        »The sample size was calculated [sic] based on a T/R-ratio of 0.90 – 1.10«
  If you assume a deviation of 10% of the test from the reference and are not sure about its direction (lower or higher 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.10 = 0.9091.
  Given, sometimes you will need a higher sample size with the lower T/R-ratio. There’s no free lunch.

CV      <- 0.45
delta   <- 0.10 # direction unknown
design  <- "2x2x4"
theta0s <- c(1 - delta, 1 / (1 + delta),
             1 + delta, 1 / (1 - delta))
n       <- sampleN.scABEL(CV = CV, theta0 = 1 - delta,
                          design = design,
                          details = FALSE,
                          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.scABEL(CV = CV,
                                 theta0 = comp1$theta0[i],
                                 design = design, n = n)
}
n       <- sampleN.scABEL(CV = CV, theta0 = 1 + delta,
                          design = design,
                          details = FALSE,
                          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.scABEL(CV = CV,
                                 theta0 = comp2$theta0[i],
                                 design = design, 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.45 0.9000  TRUE 28 0.8112
#  0.45 0.9091 FALSE 28 0.8388
#  0.45 1.1000 FALSE 28 0.8397
#  0.45 1.1110 FALSE 28 0.8101
#    CV theta0  base  n  power
#  0.45 0.9000 FALSE 26 0.7846
#  0.45 0.9091 FALSE 26 0.8149
#  0.45 1.1000  TRUE 26 0.8140
#  0.45 1.1110 FALSE 26 0.7837

</nitpick>

---

Essentially this leads to the murky waters of prospective sensitivity analyses, 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.45
theta0   <- 0.90
target   <- 0.80
minpower <- 0.70
pa       <- pa.scABE(CV = CV, theta0 = theta0,
                     targetpower = target,
                     minpower = minpower,
                     design = "2x2x4")
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] <- "rel. change"
print(pa, plotit = FALSE)
print(comp, row.names = FALSE)

# Sample size plan scABE (EMA/ABEL)
#  Design alpha CVwT CVwR theta0 theta1 theta2 Sample size
#   2x2x4  0.05 0.45 0.45    0.9    0.8   1.25          28
#  Achieved power Target power
#         0.81116          0.8
# 
# Power analysis
# CV, theta0 and number of subjects leading to min. acceptable power of ~0.7:
#  CV= 0.6629, theta0= 0.8719
#  n = 22 (power= 0.7185)
# 
#  parameter rel. change
#         CV     +47.32%
#     theta0      -3.12%
#          n     -21.43%

Confirms what we have seen above. As expect the method is robust to changes of the CV. The sample size is also not very sensitive; many overrate the impact of dropouts on power.

top of section ↩︎ previous section ↩︎

Heteroscedasticity

As we have seen already above, for an ANOVA we have to assume homoscedasticity.

I recommend to perform pilot studies in one of the fully replicated designs. When you are concerned about drop­outs in a four period design or the bioanalytical method requires large sample volumes, opt for one the two sequence three period designs (TRT|RTR or TRR|RTT).
Contrary to the partial replicate design (TRR|RTR|RRT) we get estimates of both \(\small{CV_\text{wT}}\) and \(\small{CV_\text{wR}}\). Since pharmaceutical technology improves, it is not uncommon that \(\small{CV_\text{wT}<CV_\text{wR}}\). If this is the case, we get an incentive in the sample size \(\small{n}\) of the planned pivotal study (expanding the limits is based on \(\small{CV_\text{wR}}\) but the 90% CI on the – pooled – \(\small{s_\text{w}^{2}}\)).

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

Say, we performed two pilot studies. In the partial replicate we estimated the \(\small{CV_\text{wR}}\) with 0.45 (estimation of \(\small{CV_\text{wT}}\) is not possible). In the full replicate we estimated \(\small{CV_\text{wR}}\) with 0.484 and \(\small{CV_\text{wT}}\) with 0.414. Note that the \(\small{CV_\text{w}}\) according to \(\small{(7)}\) is 0.45 as well. How will that impact the sample size \(\small{n}\) of the pivotal 4-period full replicate design?

comp <- data.frame(pilot = c("TRR|RTR|RRT", "TRT|RTR"),
                   CVwT = c(0.45, 0.414),
                   CVwR = c(0.45, 0.484),
                   CVw = NA,
                   n = NA, power = NA)
for (i in 1:nrow(comp)) {
  comp[i, 4]   <- signif(
                    mse2CV((CV2mse(comp$CVwT[i]) +
                            CV2mse(comp$CVwR[i])) / 2), 3)
  comp[i, 5:6] <- sampleN.scABEL(CV = c(comp$CVwT[i], comp$CVwR[i]),
                                 design = "2x2x4", details = FALSE,
                                 print = FALSE)[8:9]
}
print(comp, row.names = FALSE)

#        pilot  CVwT  CVwR  CVw  n   power
#  TRR|RTR|RRT 0.450 0.450 0.45 28 0.81116
#      TRT|RTR 0.414 0.484 0.45 24 0.80193

Since bioanalytics drives study costs to a great extent, we may save ~14%.

Note that when you give CV as two-element vector, the first element has to be \(\small{CV_\text{wT}}\) and the second \(\small{CV_\text{wR}}\).

Although according to the guidelines it is not required to estimate \(\small{CV_\text{wT}}\), its value is ‘nice to know’. Sometimes studies fail only due to the large \(\small{CV_\text{wR}}\) thus inflating the confidence interval. In such a case you have at least ammunation to start an argument.

Even if you plan the pivotal study in a partial replicate design (why on earth?) knowing both \(\small{CV_\text{wT}}\) and \(\small{CV_\text{wR}}\) is useful.

comp <- data.frame(pilot = c("TRR|RTR|RRT", "TRT|RTR"),
                   CVwT = c(0.45, 0.414),
                   CVwR = c(0.45, 0.484),
                   CVw = NA,
                   n = NA, power = NA)
for (i in 1:nrow(comp)) {
  comp[i, 4]   <- signif(
                    mse2CV((CV2mse(comp$CVwT[i]) +
                            CV2mse(comp$CVwR[i])) / 2), 3)
  comp[i, 5:6] <- sampleN.scABEL(CV = c(comp$CVwT[i], comp$CVwR[i]),
                                 design = "2x3x3", details = FALSE,
                                 print = FALSE)[8:9]
}
print(comp, row.names = FALSE)

#        pilot  CVwT  CVwR  CVw  n   power
#  TRR|RTR|RRT 0.450 0.450 0.45 39 0.80588
#      TRT|RTR 0.414 0.484 0.45 36 0.80973

Again, a smaller sample size is possible and we may still save ~7.7%.

Note that sampleN.scABEL() is inaccurate for the partial replicate design if and only if \(\small{CV_\text{wT}>CV_\text{wR}}\). Let’s reverse the values and compare the results with subject simulations by the function sampleN.scABEL.sdsims().

comp <- data.frame(method = c("key", "subj"), n = NA,
                   power = NA, rel.speed = NA)
st   <- proc.time()[[3]]
comp[1, 2:3] <- sampleN.scABEL(CV = c(0.484, 0.414),
                               design = "2x3x3",
                              details = FALSE,
                              print = FALSE)[8:9]
et   <- proc.time()[[3]]
comp$rel.speed[1] <- et - st
st   <- proc.time()[[3]]
comp[2, 2:3] <- sampleN.scABEL.sdsims(CV = c(0.484, 0.414),
                                      design = "2x3x3",
                                      details = FALSE,
                                      print = FALSE)[8:9]
et   <- proc.time()[[3]]
comp$rel.speed[2] <- et - st
comp$rel.speed    <- signif(comp$rel.speed /
                            comp$rel.speed[1], 3)
comp[1] <- c("\u2018key statistics\u2019",
             "subject simulations")
print(comp, row.names = FALSE)

#               method  n   power rel.speed
#     ‘key statistics’ 45 0.80212       1.0
#  subject simulations 48 0.80938      84.6

In such a case, always use the function sampleN.scABEL.sdsims().

top of section ↩︎ previous section ↩︎

Multiple Endpoints

For demonstrating 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 \(\alpha\), we are protected by the intersection-union principle.^{53 54} We design the study always for the worst case combination, i.e., based on the PK metric requiring the largest sample size. In most jurisdictions wider BE limits are acceptable only for C_max. Let’s explore that with different CVs and T/R-ratios.

metrics <- c("Cmax", "AUC")
methods <- c("ABEL", "ABE")
CV      <- c(0.45, 0.35)
design  <- "2x2x4"
theta0  <- c(0.90, 0.925)
df      <- data.frame(metric = metrics, method = methods,
                      CV = CV, theta0 = theta0, n = NA,
                      power = NA)
df[1, 5:6] <- sampleN.scABEL(CV = CV[1], theta0 = theta0[1],
                             design = design, details = FALSE,
                             print = FALSE)[8:9]
df[2, 5:6] <- sampleN.TOST(CV = CV[2], theta0 = theta0[2],
                           design = design, print = FALSE)[7:8]
df$power <- signif(df$power, 5)
txt      <- paste0("Sample size is driven by ",
                   df$metric[df$n == max(df$n)], ".\n")
print(df, row.names = FALSE); cat(txt)

#  metric method   CV theta0  n   power
#    Cmax   ABEL 0.45  0.900 28 0.81116
#     AUC    ABE 0.35  0.925 36 0.81604
# Sample size is driven by AUC.

Commonly the PK metric evaluated by ABE drives the sample size. Consequently, with 36 subjects the study will be ‘overpowered’ (~0.886) for C_max.

Let us assume the same T/R-ratios for both 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.scABEL(theta0 = x, CV = df$CV[1],
               design = design,
               n = df$n[2]) - target
}
metrics <- c("Cmax", "AUC")
methods <- c("ABEL", "ABE")
CV      <- c(0.45, 0.30)
theta0  <- c(0.90, 0.90)
design  <- "2x2x4" # must be a replicate
target  <- 0.80
if (!design %in% c("2x2x4", "2x3x3", "2x2x3"))
  stop("design = \"", design, "\" is not supported.")
df      <- data.frame(metric = metrics, method = methods,
                      CV = CV, theta0 = theta0, n = NA,
                      power = NA)
df[1, 5:6] <- sampleN.scABEL(CV = CV[1], theta0 = theta0[1],
                             design = design, details = FALSE,
                             print = FALSE)[8:9]
df[2, 5:6] <- sampleN.TOST(CV = CV[1], theta0 = theta0[2],
                           design = design, print = FALSE)[7:8]
df$power <- signif(df$power, 5)
if (theta0[1] < 1) {
  res <- uniroot(opt, tol = 1e-8,
                 interval = c(0.80 + 1e-4, theta0[1]))
} else {
  res <- uniroot(opt, tol = 1e-8,
                 interval = c(theta0[1], 1.25 - 1e-4))
}
res     <- unlist(res)
theta0s <- c(res[["root"]], 1/res[["root"]])
txt     <- paste0("Target power for ", metrics[1],
                  " and sample size ",
                  df$n[2], "\nachieved for theta0 ",
                  sprintf("%.4f", theta0s[1]), " or ",
                  sprintf("%.4f", theta0s[2]), ".\n")
print(df, row.names = FALSE); cat(txt)

#  metric method   CV theta0  n   power
#    Cmax   ABEL 0.45    0.9 28 0.81116
#     AUC    ABE 0.30    0.9 84 0.80569
# Target power for Cmax and sample size 84
# achieved for theta0 0.8340 or 1.1990.

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

For HVD(P)s Health Canada accepts ABEL only for AUC, whereas for C_max the T/R-ratio has to lie within 80.0 – 125.0% (i.e., without assessing its CI).

metrics <- c("Cmax", "AUC")
methods <- c("PE", "ABE(L)")
alpha   <- c(0.50, 0.05)
CV      <- c(0.65, 0.35)
theta0  <- c(0.90, 0.90)
design  <- "2x2x4" # must be a replicate
target  <- 0.90
if (!design %in% c("2x2x4", "2x3x3", "2x2x3"))
  stop("design = \"", design, "\" is not supported.")
df      <- data.frame(design = design, metric = metrics, method = methods,
                      CV = CV, theta0 = theta0, n = NA, power = NA)
# alpha = 0.5 leads to a zero-width confidence
# interval, i.e., only the T/R-ratio is relevant
df[1, 6] <- sampleN.TOST(alpha = alpha[1],
                         CV = CV[1],
                         theta0 = theta0[1],
                         design = design,
                         targetpower = target,
                         print = FALSE)[["Sample size"]]
if (df[1, 6] < 12) df[1, 6] <- 12 # acc. to the guidance
df[1, 7] <- power.TOST(alpha = alpha[1],
                       CV = CV[1],
                       theta0 = theta0[1],
                       design = design,
                       n = df[1, 6])
df[2, 6] <- sampleN.scABEL(alpha = alpha[2],
                           CV = CV[2],
                           theta0 = theta0[2],
                           design = design,
                           regulator = "HC",
                           targetpower = target,
                           details = FALSE,
                           print = FALSE)[["Sample size"]]
if (df[2, 6] < 12) df[2, 6] <- 12 # acc. to the guidance
df[2, 7] <- power.scABEL(alpha = alpha[2],
                         CV = CV[2],
                         theta0 = theta0[2],
                         regulator = "HC",
                         design = design,
                         n = df[2, 6])
driving  <- which(df$n == max(df$n))
easy     <- which(df$n != max(df$n))
if (driving == 1) {
  easy.pwr <- power.TOST(alpha = alpha[easy],
                         CV = CV[easy],
                         theta0 = theta0[easy],
                         design = design,
                         n = df$n[driving])
} else {
  easy.pwr <- power.scABEL(alpha = alpha[easy],
                           CV = CV[easy],
                           theta0 = theta0[easy],
                           regulator = "HC",
                           design = design,
                           n = df$n[driving])
}
df$power <- signif(df$power, 5)
names(df)[5] <- "T/R"
txt      <- paste0(sprintf("Sample size is driven by %s.\n", df$metric[driving]),
                   sprintf("With a sample size of %.0f power for %s = %.5f.\n",
                           df$n[driving], metrics[easy], easy.pwr))
print(df, row.names = FALSE); cat(txt)

#  design metric method   CV T/R  n   power
#   2x2x4   Cmax     PE 0.65 0.9 42 0.90058
#   2x2x4    AUC ABE(L) 0.35 0.9 50 0.90897
# Sample size is driven by AUC.
# With a sample size of 50 power for Cmax = 0.91980.

Even for a CV of C_max which is substantially higher than the one of AUC, the latter still drives the sample size.

top of section ↩︎ previous section ↩︎