Bioequivalence and Bioavailability Forum

Hi Hötzi,

it is very interesting. I looked at Chow&Liu's 2009 version of the holy scriptures and they also seem not to mention RTR/RRT/TRR.

I can easily imagine that some innovator had a presub meeting with FDA (about which not much is in the public domain) and then the idea caught on in a rather un-public way and later got the nod from EMA. It is speculation of course but such things would often happen at the initiative in the private sector behind the curtains.
I don't recall the details of dossiers I assessed, I am fairly sure that I (incompetently, of course) assessed BE trials based on RTR/RRT/TRR designs from 2005 and onwards, but I am not 100% sure. Had I not sniffed all that glue, memories of the past would likely have been less foggy.

It might also be that one of the consortia like PQRI could have looked into the matter and released a white paper. I have no idea, but what a mystery
My money is on Walter Hauck, if you can get hold of him, he will tell you where it all came from.


Update: Kamal Midha referred to it as ABE3 in a presentation from 2006.
https://wayback.archive-it.org/7993/20170405065959/https://www.fda.gov/ohrms/dockets/ac/06/slides/2006-4241s2_3_files/frame.htm

Is there an easy way to download that thing, by the way????

Dear Helmut,

❝ If you know any publication (preferrably prior to 2010), please let me know.

here ya go

Dear Helmut,
in this post I mentioned the oldest study with partial replicate design (sequences RRT,RTR,TRR) I know about. The study was conducted at AstraZeneca R&D Lund, Lund, Sweden in 1999.

PDF (opening page linked to Statistical Review).

Best regards,
zizou

Dear all,

while I still hope for a reliable source (i.e., published in a peer-reviewed journal) of the partial replicate, an interlude about sample sizes and more.

Of course, I understand that people prefer three periods over four (less blood volume, lower chance of drop­outs). László Tóthfalusi mentioned (IIRC, »Dissolution Testing, Bio­avail­abi­lity & Bio­equivalence Conference«, Budapest 2006) that power depends on the number of administrations and therefore, if \(\small{N}\) is the sample size of a 2×2×2 crossover, for a 4-period replicate one would need \(\small{n=N/2}\) subjects and for 3-period replicates \(\small{n=3N/4}\). However, this is only approximate because the degrees of freedom are different (see below). Currently, this approximation is implemented – sorry, Dave and Yung-jin – in FARTSSIE and the R-package bear.

   CV design  type    df   n  df   power admins appr match dev part.full
 0.20  2x2x2 conv.   N-2  20  18 0.83468     40                         
 0.20  2x2x4  full 3*n-4  10  26 0.84331     40   10   yes             
 0.20  2x2x3  full 2*n-3  14  25 0.81793     42   16    no   +         
 0.20  2x3x3 part. 2*n-3  15  27 0.84401     45   15   yes             +
 0.25  2x2x2 conv.   N-2  28  26 0.80744     56                         
 0.25  2x2x4  full 3*n-4  14  38 0.81399     56   14   yes             
 0.25  2x2x3  full 2*n-3  22  41 0.83198     66   22   yes             
 0.25  2x3x3 part. 2*n-3  21  39 0.81434     63   21   yes             –
 0.30  2x2x2 conv.   N-2  40  38 0.81585     80                         
 0.30  2x2x4  full 3*n-4  20  56 0.82024     80   20   yes             
 0.30  2x2x3  full 2*n-3  30  57 0.82040     90   30   yes             
 0.30  2x3x3 part. 2*n-3  30  57 0.82040     90   30   yes             =
 0.35  2x2x2 conv.   N-2  52  50 0.80747    104                         
 0.35  2x2x4  full 3*n-4  26  74 0.81090    104   26   yes             
 0.35  2x2x3  full 2*n-3  38  73 0.80082    114   40    no   +         
 0.35  2x3x3 part. 2*n-3  39  75 0.81099    117   39   yes             +
 0.40  2x2x2 conv.   N-2  66  64 0.80525    132                         
 0.40  2x2x4  full 3*n-4  34  98 0.81934    136   34   yes             
 0.40  2x2x3  full 2*n-3  50  97 0.81189    150   50   yes             
 0.40  2x3x3 part. 2*n-3  51  99 0.81940    153   51   yes             +
 0.45  2x2x2 conv.   N-2  82  80 0.80691    164                         
 0.45  2x2x4  full 3*n-4  42 122 0.81823    168   42   yes             
 0.45  2x2x3  full 2*n-3  62 121 0.81220    186   62   yes             
 0.45  2x3x3 part. 2*n-3  63 123 0.81826    189   63   yes             +
 0.50  2x2x2 conv.   N-2  98  96 0.80322    196                         
 0.50  2x2x4  full 3*n-4  50 146 0.81281    200   50   yes             
 0.50  2x2x3  full 2*n-3  74 145 0.80768    222   74   yes             
 0.50  2x3x3 part. 2*n-3  75 147 0.81283    225   75   yes             +
 0.55  2x2x2 conv.   N-2 116 114 0.80386    232                         
 0.55  2x2x4  full 3*n-4  58 170 0.80539    232   58   yes             
 0.55  2x2x3  full 2*n-3  86 169 0.80089    258   88    no   +         
 0.55  2x3x3 part. 2*n-3  87 171 0.80541    261   87   yes             +
 0.60  2x2x2 conv.   N-2 134 132 0.80173    268                         
 0.60  2x2x4  full 3*n-4  68 200 0.80879    272   68   yes             
 0.60  2x2x3  full 2*n-3 100 197 0.80111    300  102    no   +         
 0.60  2x3x3 part. 2*n-3 102 201 0.80880    306  102   yes             +

The approximation works well for the 4-period replicate but not so for the 3-period replicates. Hence, I suggest PowerTOST. In most cases the partial replicate requires more subjects than the 3-period full replicate.

In a nutshell:

 Partial replicate (TRR|RTR|RRT) 

1. Generally slightly larger sample sizes than full replicates. 〰️
   Since same degrees of freedom, the sample size is a multiple of three instead of two.
   
2. Estimation of s_wR / CV_wR. ✅
   
3. Estimation of s_wT / CV_wT not possible. ❌
   
4. The EMA’s “all fixed effects” model for ABE. ✅
   
5. Sometimes convergence issues with the FDA’s covariance structure of the mixed effects model for ABE.* Changing from FA0(2) to FA0(1) or CSH generally – not always! – helps. Should be specified in the SAP. Rarely convergence is not achieved with any specification. ?
   
6. The EMA’s ABEL (CV_wR >30%). ✅
   
7. The FDA’s RSABE (s_wR ≥0.294). ✅

 3-period full replicates (TRT|RTR or TRR|RTT) 

1. Generally slightly smaller sample sizes than partial replicate. 〰️
   Since same degrees of freedom, the sample size is a multiple of two instead of three.
   
2. Estimation of s_wR / CV_wR. ✅
   
3. Estimation of s_wT / CV_wT. ✅
   
4. The EMA’s “all fixed effects” model for ABE. ✅
   
5. The FDA’s covariance structure of the mixed effects model for ABE. ✅
   
6. The EMA’s ABEL (CV_wR >30%). ✅
   
7. The FDA’s RSABE (s_wR ≥0.294). ✅

#1: Not a big deal; if sample sizes are different, the difference is pretty small.
#3: In a pilot study CV_wT is not only nice to know but useful. If CV_wT < CV_wR you get an incentive in planning the pivotal study (scaling depends on CV_wR but the BE assessment on the pooled CV_w). If CV_wT > CV_wR you can account for that and increase the sample size accordingly. If the pilot study was performed in the partial replicate design you have to assume that CV_wT = CV_wR. If CV_wT < CV_wR you waste money. If CV_wT > CV_wR your study will be underpowered. For a pivotal study you can approach the WHO for reference-scaling of AUC (cause a 4-period full replicate is recommended in the guidance).
#5: Can be a show stopper in the partial replicate design. There is a slight risk (ABE for the FDA using FA0(2) acc. to the guidance and no convergence) that the PK metric in question cannot be assessed for ABE at all.

 Given all that, I don’t see any justification
for using the partial replicate design. 

If you know one, please let me know.

---

• The model tries to estimate \(\small{s_\textrm{wT}^2}\) though the test was only administered once. Hence, the model is over-spe­ci­fied for this design. Even if convergence is achieved, the estimate is nonsense and differs (contrary to \(\small{s_\textrm{wR}^2}\)) between software packages. Only the total variance \(\small{s_\textrm{T}^2=s_\textrm{wT}^2+s_\textrm{bT}^2}\) could be estimated but not done in the model.

---

R-script for the table:

library(PowerTOST)
balance <- function(x, y) {
  return(y * (x %/% y + as.logical(x %% y)))
}
designs <- known.designs()[c(3, 8, 7, 10), 2:3]
type    <- c("conv.", "full", "full", "part.")
eval(parse(text = designs$df[1]))
CV  <- seq(0.2, 0.6, 0.05)
res <- data.frame(CV = rep(CV, each = 4), design = designs[, 1], type = type,
                  df.1 = designs[, 2], n = NA, df.2 = NA, power = NA,
                  admins = NA, appr = "", match = "", dev = "", part.full = "")
for (j in 1:nrow(res)) {
  tmp <- sampleN.TOST(CV = res$CV[j], design = res$design[j],
                      print = FALSE, details = FALSE)
  res$n[j]     <- tmp[["Sample size"]]
  res$power[j] <- signif(tmp[["Achieved power"]], 5)
  n            <- res$n[j]
  nseq         <- as.numeric(substr(res$design[j], 3, 3))
  nper         <- as.numeric(substr(res$design[j], 5, 5))
  if (res$design[j] == "2x2x2") N <- n
  if (!res$design[j] == "2x2x2") {
    if (nper == 4) res$appr[j] <- balance(N/2, nseq)
    if (nper == 3) res$appr[j] <- balance(3*N/4, nseq)
    if (res$appr[j] == res$n[j]) {
      res$match[j] <- "yes"
    } else {
      res$match[j] <- "no"
      ifelse (res$appr[j] < res$n[j], res$dev[j] <- "\u2013", res$dev[j] <- "+")
    }
    if (res$design[j] == "2x2x3") n.3full <- n
    if (res$design[j] == "2x3x3") {
      if (n < n.3full)  res$part.full[j] <- "\u2013"
      if (n == n.3full) res$part.full[j] <- "="
      if (n > n.3full)  res$part.full[j] <- "+"
    }
  }
  res$admins[j] <- nper*n
  df            <- designs$df[which(designs$design == res$design[j])]
  res$df.2[j]   <- eval(parse(text = df))
}
names(res)[c(4, 6)] <- rep("df", 2)
res[which(res[, 4] == "n-2"), 4] <- "N-2"
print(res, row.names = FALSE)

R-script to simulate 1,000 studies; CV 20–60%, T/R-ratio 0.85–0.95, power 80–90%:

library(PowerTOST)
set.seed(1234567)
sims   <- 1000
CV     <- runif(sims, 0.20, 0.60)
theta0 <- runif(sims, 0.85, 0.95)
target <- runif(sims, 0.80, 0.90)
res    <- data.frame(f = rep(NA, sims), p = NA)
for (j in 1:sims) {
  res$f[j] <- sampleN.TOST(CV = CV[j], theta0 = theta0[j], design = "2x2x3",
                           targetpower = target[j], print = FALSE)[["Sample size"]]
  res$p[j] <- sampleN.TOST(CV = CV[j], theta0 = theta0[j], design = "2x3x3",
                           targetpower = target[j], print = FALSE)[["Sample size"]]
}
bp <- boxplot(res, plot = FALSE)
plot(c(0.5, 2.5), range(bp$stats[1, ], bp$stats[5, ]), type = "n", axes = FALSE,
     xlab = "3-period replicate designs", ylab = "sample size",
     main = paste(sims, "studies simulated for ABE"), cex.main = 1.2, font.main = 1)
abline(h = axTicks(2, log = FALSE), col = "lightgrey", lty = 3)
axis(1, at = 1:2, labels = c("full", "partial"), tick = FALSE)
axis(2, las = 1)
bxp(bp, las = 1, outline = FALSE, boxwex = 0.4, add = TRUE, ann = FALSE,
    boxfill = "bisque", medcol = "blue", axes = FALSE, frame.plot = TRUE)
text(rep(1.25, 5), bp$stats[, 1], labels = bp$stats[, 1], pos = 4)
text(rep(2.25, 5), bp$stats[, 2], labels = bp$stats[, 2], pos = 4)

Dear all,

I simulated 500 data sets in the partial replicate design with \(\small{s_\textrm{wT}^2=s_\textrm{wR}^2=0.086\: (CV_\textrm{w}\approx 29.97\%),}\) \(\small{s_\textrm{bT}^2=s_\textrm{bR}^2=0.172\: (CV_\textrm{b}\approx 43.32\%),}\) \(\small{\rho=1},\) \(\small{\theta_0=1},\) i.e., no subject-by-formulation interaction. With \(\small{n=24}\) subjects 82.3797% power to demonstrate ABE. Evaluation in Phoenix/WinNonlin 8.1 with the FDA’s covariance structure FA0(2). Singularity tolerance and convergence criterion 1E-12 (instead of 1E-10), maximum iterations 250 (instead of 50).
If you want to try it in SAS or any other software: The data sets in CSV-format.

In 403 (80.6%) of the data sets PHX issued at least one warning.
In 56 (11.2%) of the data sets PHX threw this:

Negative final variance component. Consider omitting this VC structure.

Well roared lion. The model reached for the stars (namely \(\small{s_\textrm{wT}^2}\)). We know that we can get only the total variability (i.e., \(\small{s_\textrm{T}^2=s_\textrm{wT}^2+s_\textrm{bT}^2}\)) like in a parallel design where the within-subject variance (as well as the between variance) is not accessible as well. Amazingly in the other cases PHX got an estimate though it’s nonsense, of course.

In 340 (68%) of the data sets I was told:

Model may be over-specified. A simpler model could be tried.

Oh yes! »May be« is an euphemism. Actually it is. Not only in some of the data sets but in all.

In 333 (66.6%) of the data sets PHX threw this:

Newton's algorithm converged with modified Hessian. Output is suspect.

How I love to be told that results are suspect. Will assessors love that as well? Well, in PHX it’s hidden in the ‘Core Output’ and the ‘Warnings and Errors’. In SAS in a log-file… Will it be shown in your fancy output? Not necessarily.

79.6% of the data sets passed BE. That’s only slightly lower than expected and likely due to the small number of simulations (10,000 are running).
How close are the estimates to the targets?

                s²wR     s²bR     s²T       S×F      PE
─────────────────────────────────────────────────────────
target         0.08600  0.17200  0.25800  0        100.00
mean estimate  0.08733  0.16923  0.22655  0.01090   98.48
%RE            +1.54%   –1.61%   –12.19%   –       –1.52%

Fine with me, except \(\small{s_\textrm{T}^2}\) which is not directly estimated (see this lengthy thread for a promising alternative) but as the sum of two doubtful estimates (\(\small{s_\textrm{wT}^2,s_\textrm{bT}^2}\)).

We see also that the optimizer is fine in estimating CV_wR but desperate with CV_wT (only 437 values). The target was 29.9677% for both.

       min    QI     med    QIII   max
───────────────────────────────────────
CVwR  14.25  26.75  30.11  33.20  46.13
CVwT   1.33  18.29  23.23  28.05  45.33

I evaluated the data sets with other covariance structures as well. Seems that FA0(1) is the winner.

Convergence          FA0(2)      FA0(1)        CS
─────────────────────────────────────────────────────
Achieved          160 (32.0%)  500 (100%)  500 (100%)
Modified Hessian  340 (68.0%)     –           –

Warnings                       FA0(2)      FA0(1)      CS
─────────────────────────────────────────────────────────────
Modified Hessian            333 (66.6%)     –          –
Negative variance component  56 (11.2%)  25 (5.0%)  56 (11.2%)
Both                         14 ( 2.8%)     –          –

As long as we achieve convergence, it doesn’t matter (we have seen nasty data sets in the past, where FA0(2) didn’t converge). Perhaps as long as the data set is balanced and/or does not contain ‘outliers’, all is good. At the end of the day we are interested in the 90% CI. I compared the results obtained with FA0(1) and CS to the guidances’ FA0(2). Up to the 4^th decimal (rounded to percent, i.e., 6–7 significant digits) the CI was identical in all cases. Only when I looked at the 5^th decimal for both covariance structures, 1/500 differed (the CI was wider). Since all guidelines require rounding to the 2^nd decimal, that’s not relevant.

I’m not a friend of the EMA’s ‘all effects fixed’ model because it assumes identical variances of T and R (which has be shown to be wrong in many full replicate studies). But, of course, no issues with convergence in this simple linear model.

My original simulation code contained a stupid error (THX to Detlew for detecting it!) which lead to an extreme S×F-interaction. Example of one data set where the optimizer was in deep trouble. The default maximum iterations in PHX/WNL are 50. I got:

max.iter    s²wR     %RE  -2REML LL   AIC    BIC       df         90% CI
────────────────────────────────────────────────────────────────────────────
    50    0.084393  –1.87  39.368   61.368  85.455  22.10798  82.212–111.084
   250    0.085094  –1.05  39.345   61.345  85.431  22.18344  82.223–111.070
 1,250    0.085271  –0.85  39.339   61.339  85.425  22.20523  82.225–111.066
 6,250    0.085309  –0.80  39.338   61.338  85.424  22.20991  82.226–111.066
31,250    0.085317  –0.79  39.338   61.338  85.424  22.21032  82.226–111.066

A lament in all cases:

Failed to converge in allocated number of iterations. Output is suspect.

Note that the degrees of freedom increase with the number of iterations and hence, the CI narrows. Now I understand why Health Canada requires that the optimizer’s constraints are stated in the SAP.

Welcome to the hell of mixed effects modeling. ?

---

Edit 1: Results of a large data set (10,000 simulations, 20.5 MB in CSV-format). 80.95% passed BE in all setups. Relevant estimates were identical and pretty close to the targets:

                s²wR      PE
──────────────────────────────
target         0.08600  100.00
mean estimate  0.08584   99.97
%RE            –0.19%   –0.03%

Convergence        FA0(2)  FA0(1)   CS
───────────────────────────────────────
Achieved          30.14%   99.97%  100%
Modified Hessian  69.83%     –      –
> max. iter.       0.03%    0.03%   –

Warnings                    FA0(2)  FA0(1)   CS
─────────────────────────────────────────────────
Modified Hessian            68.75%    –      –
Negative variance component  9.01%  3.82%  11.15%
Both                         2.22%  0.06%    –

Given all that, I would opt for FA0(1).

---

Edit 2: I manipulated the small data sets. Removed subject 24 to make the study imbalanced, removed the last period (T) of subject 23 to make it incomplete. Multiplied T of subject 1 with 5–10 to mimic an ‘outlier’. Only 15.4% of studies passed. Yep, even a single outlier might be the killer. 381 warnings by FA0(2), 2 by FA0(1), and 7 by CS.