Notes on the methods
Estimation of intra-subject variability
The EMA proposed a
linear model of log-transformed PK
responses of the reference treatment
sequence, subject(sequence),
period
where all effects are fixed. Estimated via function lm() of
package stats:
modCV <- lm(log(PK) ~ sequence + subject%in%sequence + period,
data = data[data$treatment = "R", ])In full replicate designs the same model is run with
data = data[data$treatment = "T", ] for informational
purposes.
Special conditions for the sample size in three period full replicate designs:
The question raised asks if it is possible to use a design where subjects are randomised to receive treatments in the order of TRT or RTR.
The CHMP bioequivalence guideline requires that at least 12 patients are needed to provide data for a bioequivalence study to be considered valid, and to estimate all the key parameters. Therefore, if a 3-period replicate design, where treatments are given in the order TRT or RTR, is to be used to justify widening of a confidence interval for Cmax then it is considered that at least 12 patients would need to provide data from the RTR arm. This implies a study with at least 24 patients in total would be required if equal number of subjects are allocated to the 2 treatment sequences.
— Q&A document14
If less than twelve subjects are eligible in sequence
RTR of a TRT | RTR design (and in analogy in
sequence TRR of a TRR | RTT design), the user
is notified about the ‘uncertain’ estimate of CVwR.
However, in a sufficiently powered study such a case is extremely unlikely.
Let us explore the 95% confidence interval of the CV:
# Estimate sample sizes of full replicate designs (theta0 0.90,
# target power 0.80) and CI of the CV with package PowerTOST
CV <- 0.30
design <- c("2x2x4", "2x2x3") # 4- and 3-period full replicate designs
res <- data.frame(design = rep(design, 2), n = c(rep(NA, 2), rep(12, 2)),
df = NA, CV = CV, lower = NA, upper = NA)
# 95% confidence interval of the CV
for (i in 1:nrow(res)) {
if (is.na(res$n[i])) {
res$n[i] <- PowerTOST::sampleN.scABEL(CV = CV,
design = res$design[i],
details = FALSE,
print = FALSE)[["Sample size"]]
}
if (i > 2 ) {
n.1 <- res$n[i-2] / 2 # no dropouts in one sequence
n.2 <- 12 # only 12 eligible subjects in the other
res$n[i] <- n.1 + n.2
}
if (res$design[i] == "2x2x4") {
res$df[i] <- 3 * res$n[i] - 4
} else {
res$df[i] <- 2 * res$n[i] - 3
}
res[i, 5:6] <- round(PowerTOST::CVCL(CV = CV,
df = res$df[i],
side = "2-sided",
alpha = 0.05), 3)
}
print(res, row.names = FALSE)
# design n df CV lower upper
# 2x2x4 34 98 0.3 0.262 0.351
# 2x2x3 50 97 0.3 0.262 0.352
# 2x2x4 29 83 0.3 0.259 0.357
# 2x2x3 37 71 0.3 0.256 0.362
# Rows 1-2: Sample sizes for target power
# Rows 3-4: Only 12 eligible subjects in one sequenceIt is incomprehensible why a four period full replicate design is considered by the EMA to give a more ‘reliable’ estimate of the CV than a three period full replicate design.
Model structure
The EMA’s models assume equal [sic] intra-subject variances of Test and Reference (like in 2×2×2 trials) – even if proven false in one of the full replicate designs (were both CVwT and CVwR can be estimated). Hence, among biostatisticians they are called ‘crippled models’ because the replicative nature of the study is ignored.
It should be noted that by default ANOVA() produces a
‘Type I’ table, i.e., all tests are performed against the
residual error. Starting with version 1.1.0 ‘Type III’ is employed in
order to obtain the correct test for carryover, where subject
has to be tested against subject(sequence).15
library(replicateBE) # attach the library to run example scripts
method.A(data = rds01, print = FALSE, verbose=TRUE)
#
# Data set DS01: Method A by lm()
# ───────────────────────────────────
# Type III Analysis of Variance Table
#
# Response: log(PK)
# Df Sum Sq Mean Sq F value Pr(>F)
# sequence 1 0.0077 0.007652 0.00268 0.9588496
# period 3 0.6984 0.232784 1.45494 0.2278285
# treatment 1 1.7681 1.768098 11.05095 0.0010405
# sequence:subject 75 214.1296 2.855061 17.84467 < 2.22e-16
# Residuals 217 34.7190 0.159995
#
# treatment T – R:
# Estimate Std. Error t value Pr(>|t|)
# 0.14547400 0.04650870 3.12788000 0.00200215
# 217 Degrees of FreedomIn ‘Type I’ the F value of sequence would be \(\small{\frac{0.007652}{0.159995}=0.04783}\) with \(\small{p(F_{q=0.04783,\,df_1=1,\,df_2=217})\sim 0.827}\), whereas in ‘Type III’ it is \(\small{\frac{0.007652}{2.855061}=0.00268}\) with \(\small{p(F_{q=0.00268,\,df_1=1,\,df_2=75})\sim 0.959}\).
The nested structure subject(sequence) of the methods leads to an over-specified model.16 The
simple model
sequence, subject, period,
treatment
gives identical estimates of the residual variance and the treatment
effect and hence, its confidence interval.
The same holds true for the EMA’s model to estimate
CVwR. The simple model
subject, period
gives an identical estimate of the residual variance.
Reference-scaling is acceptable for Cmax (immediate release products: BE Guideline) and Cmax, Cmax,ss, Cτ,ss, partialAUC (modified release products17). The intention to expand the limits has to be stated in the protocol and – contrary to the FDA’s RSABE – a clinical justification provided.
Those HVDP for which a wider difference in Cmax is considered clinically irrelevant based on a sound clinical justification can be assessed with a widened acceptance range. The request for widened interval must be prospectively specified in the protocol.
— CHMP Bioequivalence Guideline
BE limits, PE restriction, rounding issues
The limits can be expanded based on CVwR.
- \(\small{CV_\textrm{wR}\leq
30\%}\)
Lower cap, i.e., no scaling, conventional limits:
\(\small{\left\{{L,U}\right\} = 80.00 - 125.00\%}\) - \(\small{30\%<CV_\textrm{wR}\leq
50\%}\) (EMA), \(\small{30\%<CV_\textrm{wR}\leq
\sim57.4\%}\) (Health Canada)
Expanded limits based on \(s_\textrm{wR}\):
\(\small{\left\{L,\,U\right\} = 100\cdot \exp(\mp 0.760\cdot{s_\textrm{wR}})}\) - \(\small{CV_\textrm{wR}>50\%}\)
(EMA), \(\small{CV_\textrm{wR}>\sim57.4\%}\)
Upper cap, i.e., applying \(\small{s^*_\textrm{wR}=\sqrt{\log_{e}(CV_\textrm{cap}^{2}+1)}}\) in the expansion formula or
\(\small{\left\{L,U\right\} = 69.84 - 143.19\%}\) (EMA)
\(\small{\left\{L,U\right\} = 66.7 - 150.0\%}\) (Health Canada).
The special case recommended by the Gulf Cooperation Council (Bahrain, Kuwait, Oman, Qatar, Saudi Arabia, the United Arab Emirates) is:
- \(\small{CV_\textrm{wR}\leq
30\%}\)
Lower cap, i.e., no scaling, conventional limits:
\(\small{\left\{{L,U}\right\} = 80.00 - 125.00\%}\)
- \(\small{CV_\textrm{wR}>30\%}\)
Widened limits:
\(\small{\left\{{L,U}\right\} = 75.00 - 133.33\%}\)
In reference-scaling a so-called mixed (a.k.a. aggregate) criterion is applied. In order to pass ABEL,
- the 90% confidence interval has to lie entirely within the acceptance range \(\small{\left\{L,U\right\}}\) and
- the point estimate (PE) has to lie within 80.00 – 125.00%.
To avoid discontinuities due to double rounding, expanded limits are calculated in full numeric precision and only the confidence interval is rounded according to the guidelines.
# Calculate limits with package PowerTOST
CV <- c(30, 50, 57.382)
df <- data.frame(CV = CV,
reg1 = "EMA", L1 = NA, U1 = NA, cap1 = "",
reg2 = "HC", L2 = NA, U2 = NA, cap2 = "",
reg3 = "GCC", L3 = NA, U3 = NA, cap3 = "",
stringsAsFactors = FALSE)
for (i in seq_along(CV)) {
df[i, 3:4] <- sprintf("%.3f",
PowerTOST::scABEL(CV[i]/100,
regulator = df$reg1[i])*100)
df[i, 7:8] <- sprintf("%.1f",
PowerTOST::scABEL(CV[i]/100,
regulator = df$reg2[i])*100)
df[i, 11:12] <- sprintf("%.3f",
PowerTOST::scABEL(CV[i]/100,
regulator = df$reg3[i])*100)
}
df$cap3[df$CV <= 30] <- df$cap2[df$CV <= 30] <- df$cap1[df$CV <= 30] <- "lower"
df$cap1[df$CV >= 50 & df$reg1 == "EMA"] <- "upper"
df$cap2[df$CV >= 57.382 & df$reg2 == "HC"] <- "upper"
names(df) <- c("CV(%)", rep(c("reg", "L(%)", "U(%)", "cap"), 3))
print(df, row.names = FALSE)
# CV(%) reg L(%) U(%) cap reg L(%) U(%) cap reg L(%) U(%) cap
# 30.000 EMA 80.000 125.000 lower HC 80.0 125.0 lower GCC 80.000 125.000 lower
# 50.000 EMA 69.837 143.191 upper HC 69.8 143.2 GCC 75.000 133.333
# 57.382 EMA 69.837 143.191 upper HC 66.7 150.0 upper GCC 75.000 133.333
Discrete limits resulting from rounding (EMA)
Degrees of freedom, comparison of methods
The SAS code provided by the
EMA in the
Q&A document does not
specify how the degrees of freedom should be calculated in Method B.
Hence, the default in PROC MIXED, namely
DDFM=CONTAIN is applied, i.e.,
method.B(..., option = 2). For incomplete data (missing
periods) Satterthwaite’s approximation of the degrees of freedom,
i.e., method.B(..., option = 1) or Kenward-Roger
method.B(..., option = 3) might be a better choice – if
stated as such in the
SAP.
For background about approximations in different software packages see
the electronic Supplementary
Material of Schütz et al.
The EMA seemingly prefers Method A:
A simple linear mixed model, which assumes identical within-subject variability (Method B), may be acceptable as long as results obtained with the two methods do not lead to different regulatory decisions. However, in borderline cases […] additional analysis using Method A might be required.
— Q&A document (January 2011 and later revisions)
The half-width of the confidence interval in log-scale allows a
comparison of methods (B v.s. A) where a higher value
might point towards a more conservative decision.18 In the
provided example datasets – with one exception – the conclusion of
BE (based on the mixed criterion)
agrees between Method A and Method B.
However, for the highly incomplete dataset 14 Method A is
liberal (passing by
ANOVA but failing by the
random effects model):
# Compare Method B acc. to the GL with Method A for all reference data sets.
ds <- substr(grep("rds", unname(unlist(data(package = "replicateBE"))),
value = TRUE), start = 1, stop = 5)
for (i in seq_along(ds)) {
A <- method.A(print = FALSE, details = TRUE, data = eval(parse(text = ds[i])))$BE
B <- method.B(print = FALSE, details = TRUE, data = eval(parse(text = ds[i])))$BE
r <- paste0("A ", A, ", B ", B, " - ")
cat(paste0(ds[i], ":"), r)
if (A == B) {
cat("Methods agree.\n")
} else {
if (A == "fail" & B == "pass") {
cat("Method A is conservative.\n")
} else {
cat("Method B is conservative.\n")
}
}
}
# rds01: A pass, B pass - Methods agree.
# rds02: A pass, B pass - Methods agree.
# rds03: A pass, B pass - Methods agree.
# rds04: A fail, B fail - Methods agree.
# rds05: A pass, B pass - Methods agree.
# rds06: A pass, B pass - Methods agree.
# rds07: A pass, B pass - Methods agree.
# rds08: A pass, B pass - Methods agree.
# rds09: A pass, B pass - Methods agree.
# rds10: A pass, B pass - Methods agree.
# rds11: A pass, B pass - Methods agree.
# rds12: A fail, B fail - Methods agree.
# rds13: A fail, B fail - Methods agree.
# rds14: A pass, B fail - Method B is conservative.
# rds15: A fail, B fail - Methods agree.
# rds16: A fail, B fail - Methods agree.
# rds17: A fail, B fail - Methods agree.
# rds18: A fail, B fail - Methods agree.
# rds19: A fail, B fail - Methods agree.
# rds20: A fail, B fail - Methods agree.
# rds21: A fail, B fail - Methods agree.
# rds22: A pass, B pass - Methods agree.
# rds23: A pass, B pass - Methods agree.
# rds24: A pass, B pass - Methods agree.
# rds25: A pass, B pass - Methods agree.
# rds26: A fail, B fail - Methods agree.
# rds27: A pass, B pass - Methods agree.
# rds28: A pass, B pass - Methods agree.
# rds29: A pass, B pass - Methods agree.
# rds30: A fail, B fail - Methods agree.Exploring dataset 14:
A <- method.A(print = FALSE, details = TRUE, data = rds14)
B1 <- method.B(print = FALSE, details = TRUE, data = rds14, option = 1)
B2 <- method.B(print = FALSE, details = TRUE, data = rds14) # apply default option
B3 <- method.B(print = FALSE, details = TRUE, data = rds14, option = 3)
# Rounding of CI according to the GL
A[17:21] <- round(A[17:21], 2) # all effects fixed
B1[17:21] <- round(B1[17:21], 2) # Satterthwaite's df
B2[17:21] <- round(B2[17:21], 2) # df acc. to Q&A
B3[17:21] <- round(B3[17:21], 2) # Kenward-Roger df
cs <- c(2, 10, 17:25)
df <- rbind(A[cs], B1[cs], B2[cs], B3[cs])
names(df)[c(1, 3:6, 11)] <- c("Meth.", "L(%)", "U(%)",
"CL.lo(%)", "CL.hi(%)", "hw")
df[, c(2, 11)] <- signif(df[, c(2, 11)], 5)
print(df[order(df$hw, df$BE, decreasing = c(FALSE, TRUE),
method = "radix"), ], row.names = FALSE)
# Meth. DF L(%) U(%) CL.lo(%) CL.hi(%) PE(%) CI GMR BE hw
# A 197.44 69.84 143.19 69.21 121.27 91.62 fail pass fail 0.28043
# A 192.00 69.84 143.19 69.21 121.28 91.62 fail pass fail 0.28046
# A 195.99 69.84 143.19 69.21 121.28 91.62 fail pass fail 0.28052
# A 192.00 69.84 143.19 69.99 123.17 92.85 pass pass pass 0.28261All variants of Method B are more conservative than Method A. Before
rounding the confidence interval, option = 2 with 192
degrees of freedom would be more conservative (lower
CL 69.21029) than
option = 1 with 197.44 degrees of freedom (lower
CL 69.21286). Given the
incompleteness of this data set (four missings in period 2, twelve in
period 3, and 19 in period 4), Satterthwaite’s or Kenward-Roger degrees
of freedom are probably the better choice.
For detailed comparisons between methods based on simulations see the electronic Supplementary Material of Schütz et al.
Outlier analysis
It is an open issue how outliers should be handled.
The applicant should justify that the calculated intra-subject variability is a reliable estimate and that it is not the result of outliers.
— CHMP Bioequivalence Guideline
The author ‘suggested’ box plots as a mere joke [sic] at the EGA/EMA symposium, being aware of their nonparametric nature and the EMA’s reluctance towards robust methods. Alas, this joke was included in the Q&A document.
[…] a study could be acceptable if the bioequivalence requirements are met both including the outlier subject (using the scaled average bioequivalence approach and the within-subject CV with this subject) and after exclusion of the outlier (using the within-subject CV without this subject).
An outlier test is not an expectation of the medicines agencies but outliers could be shown by a box plot. This would allow the medicines agencies to compare the data between them.
— EGA/EMA Q&A document
With the additional argument ola = TRUE in
method.A() and method.B() an outlier analysis
is performed, where the default fence = 2.19
Results differ slightly depending on software’s algorithms to calculate the median and quartiles. Example with the nine types implemented in R:
# Compare different types with some random data
set.seed(123456)
x <- rnorm(48)
p <- c(25, 50, 75)/100
q <- matrix(data = "", nrow = 9, ncol = 4,
dimnames = list(paste("type =", 1:9),
c("1st quart.", "median", "3rd quart.",
"software / default")))
for (i in 1:9) {
q[i, 1:3] <- sprintf("%.5f", quantile(x, prob = p, type = i))
}
q[c(2, 4, 6:8), 4] <- c("SAS, Stata", "SciPy", "Phoenix, Minitab, SPSS",
"R, S, MATLAB, Octave, Excel", "Maple")
print(as.data.frame(q))
# 1st quart. median 3rd quart. software / default
# type = 1 -0.75397 0.05761 0.96782
# type = 2 -0.74451 0.06358 0.96984 SAS, Stata
# type = 3 -0.75397 0.05761 0.96782
# type = 4 -0.75397 0.05761 0.96782 SciPy
# type = 5 -0.74451 0.06358 0.96984
# type = 6 -0.74924 0.06358 0.97084 Phoenix, Minitab, SPSS
# type = 7 -0.73978 0.06358 0.96883 R, S, MATLAB, Octave, Excel
# type = 8 -0.74609 0.06358 0.97017 Maple
# type = 9 -0.74569 0.06358 0.97009Note the differences even in the medians.
- Box plots of studentized20 and standarized21 model residuals are constructed.22
- Potential outliers are flagged based on the argument
fenceprovided by the user. - With the additional argument
verbose = TRUEdetailed information is shown in the console.
Example for the reference dataset 01:
Outlier analysis
(externally) studentized residuals
Limits (2×IQR whiskers): -1.717435, 1.877877
Outliers:
subject sequence stud.res
45 RTRT -6.656940
52 RTRT 3.453122
standarized (internally studentized) residuals
Limits (2×IQR whiskers): -1.69433, 1.845333
Outliers:
subject sequence stand.res
45 RTRT -5.246293
52 RTRT 3.214663If based on studentized residuals outliers are detected, additionally to the expanded limits based on the complete reference data, tighter limits are calculated based on CVwR after exclusion of outliers and BE assessed with the new limits. Note that standardized residuals are given for informational purposes only and not used for exclusion of outliers.
Output for the reference dataset 01 (re-ordered for clarity):
CVwR : 46.96% (reference-scaling applicable)
swR : 0.44645
Expanded limits : 71.23% ... 140.40% [100exp(±0.760·swR)]
Assessment based on original CVwR 46.96%
────────────────────────────────────────
Confidence interval: 107.11% ... 124.89% pass
Point estimate : 115.66% pass
Mixed (CI & PE) : pass
╟────────┼─────────────────────┼───────■────────◊─────────■───────────────╢
Outlier fence : 2×IQR of studentized residuals.
Recalculation due to presence of 2 outliers (subj. 45|52)
─────────────────────────────────────────────────────────
CVwR (outl. excl.) : 32.16% (reference-scaling applicable)
swR (recalculated) : 0.31374
Expanded limits : 78.79% ... 126.93% [100exp(±0.760·swR)]
Assessment based on recalculated CVwR 32.16%
────────────────────────────────────────────
Confidence interval: pass
Point estimate : pass
Mixed (CI & PE) : pass
╟┼─────────────────────┼───────■────────◊─────────■─╢Note that the PE and its
CI are not affected since the
entire data are used and therefore, these values not reported in the
second analysis (only the conclusion of the assessment).
The ‘ASCII line plot’ is given for informational purposes since its
resolution is only ~0.5%. The filled squares ■ are the
lower and upper 90% confidence limits, the rhombus ◊ the
point estimate, the vertical lines │ at 100% and the
PE restriction (80.00 – 125.00%),
and the double vertical lines ║ the expanded limits. The
PE and
CI take presedence over other
symbols. In this case the upper limit of the
PE restriction is not visible.
Since both analyses arrive at the same conclusion, the study should be
acceptable according to the
Q&A document.
An assessment of outliers is not required for Brazil’s ANVISA.
Health Canada23 accepts exclusion of outliers (irrespective of the treatment) if
- the procedure is stated in the protocol and followed before the results of the analysis are summarised into confidence intervals (i.e., regardless of whether results meet the standard),
- no more than 5% of the subjects may be considered to be outliers, unless there are 20 or fewer subjects, in which case only one subject may be removed,
- it is recommended that a simple outlier test such as a studentised residual being greater than 3 be used, and
- the subject in question should be identified as an outlier for all parameters, for either the test or reference product, upon which the bioequivalence decision is to be based.
Note that assessing of outliers of the test treatment is currently not
supported. In order to comply with conditions #1 and #3 use (here for
the internal dataset 01 with a fence narrower than the default):
B1 <- method.B(ola = TRUE, verbose = TRUE, print = FALSE,
details = FALSE, fence = 1.5, data = rds01,
option = 1)You will get
Outlier analysis
(externally) studentized residuals
Limits (1.5×IQR whiskers): -1.631514, 1.553557
Outliers:
subject sequence stud.res
41 RTRT 1.877877
45 RTRT -6.656940
46 TRTR -1.717435
52 RTRT 3.453122
standarized (internally studentized) residuals
Limits (1.5×IQR whiskers): -1.612749, 1.53832
Outliers:
subject sequence stand.res
41 RTRT 1.845333
45 RTRT -5.246293
46 TRTR -1.694330
52 RTRT 3.214663There are 73 subjects in the dataset with two treatments of the reference. Hence, while observing condition #2, we are allowed to delete 73×5%~3 outliers. Delete subject 45 and 52 with residuals outside ±3 from the dataset. Run the analysis on the reduced dataset as usual.
It must be mentioned that Health Canada’s approach might increase the Type I Error.24