R: Comparative BA-calculation for Average Bioequivalence with...

method.B {replicateBE}

R Documentation

Comparative BA-calculation for Average Bioequivalence with Expanding Limits by the EMA's 'Method B'

Description

This function performs the required calculations for the mixed (or aggregate) BE decision via Average Bioequivalence with Expanding Limits (ABEL) based on a linear mixed effects model with subjects as a random effect (‘Method B’) as specified in Annex I.

Usage

method.B(alpha = 0.05, path.in, path.out = tempdir(), file, set = "",
         ext, na = ".", sep = ",", dec = ".", logtrans = TRUE,
         regulator = "EMA", ola = FALSE, print = TRUE, details = FALSE,
         verbose = FALSE, ask = FALSE, plot.bxp = FALSE, fence = 2,
         data = NULL, option = 2)

Arguments

`alpha`	Type I Error (TIE) probability (nominal level of the test). Conventionally set to 0.05, resulting in a 100(1 – 2α) confidence interval. If `regulator = "HC"` and `alpha = 0.5` only the point estimate will be assessed (for highly variable C_max within 80.0–125.0%).
`path.in`	Path to the data file for import.
`path.out`	Path to save the result file if `print = TRUE`. You must have write-permission to the folder. For simplicity your home folder `⁠"~/"⁠` can be used. If missing, R’s standard temporary folder will be used. If a box plot of outliers should be saved (`plot.bxp = TRUE`), this path will be used as well.
`file`	Name of the dataset for import (without extension). Must be a string (i.e., enclosed in single or double quotation marks).
`set`	Name of the sheet of an Excel-file (mandatory). Must be a string (i.e., enclosed in single or double quotation marks).
`ext`	File-extension enclosed in single or double quotation marks. Acceptable are `⁠"csv"⁠` for character delimited variables (CSV) or `⁠"xls"⁠`, `⁠"xlsx"⁠` for Excel-files. The file-extension is not case-sensitive.
`na`	Character string denoting missing values. Acceptable are `⁠"NA"⁠` (not available), `⁠"ND"⁠` (not determined), `⁠"."⁠` (SAS), `⁠"Missing"⁠` (Phoenix WinNonlin), and `⁠""⁠` (Excel; empty cell). Missings will be converted to `⁠NA⁠` in the imported data. Defaults to `⁠"."⁠`.
`sep`	Variable separator in the CSV-file. Acceptable are `⁠","⁠` (comma = `⁠ASCII 44⁠`), `⁠";"⁠` (semicolon = `⁠ASCII 59⁠`), and `⁠"\t"⁠` (tabulator = `⁠ASCII 9⁠`). Defaults to `⁠","⁠`.
`dec`	Decimal separator in the CSV-file. Acceptable are `⁠"."⁠` (period = `⁠ASCII 46⁠`) or `⁠","⁠` (comma = `⁠ASCII 44⁠`). Defaults to `⁠"."⁠`.
`logtrans`	If `TRUE` (default) the raw data (provided in column `⁠PK⁠`) will be internally log-transformed and used in the calculations. If `FALSE` the already log-transformed data (provided in the column `⁠logPK⁠`) will be used in the calculations.
`regulator`	Set regulatory conditions. If `"EMA"` (default) conventional ABEL will be used. If `"HC"` Health Canada’s upper cap of scaling (~57.4%) will be applied. If `"GCC"` direct widening to 75.00–133.33% will be used if `CVwR > 30%`.
`ola`	Defaults to `FALSE`. If `TRUE` an outlier analysis based on the studentized and standardized (aka internally studentized) residuals of the model estimating `CVwR` is performed.
`print`	If `TRUE` (default), the function prints its results to a file. If `FALSE`, returns a data frame of results.
`details`	Defaults to `FALSE`. If `TRUE`, the function sends its results in full precision to a data frame.
`verbose`	Defaults to `FALSE`. If `TRUE` the model-table is send to the console. If `ola = TRUE` additional information about outliers are shown.
`ask`	Defaults to `FALSE`. If `TRUE` the user will be asked whether an already existing result file (and if outliers are found, the box plot) should be overwritten.
`plot.bxp`	Only observed if `ola = TRUE` and at least one outlier is found. If `FALSE` (default) the box plot will be shown in the graphics device. If `TRUE` the box plot will be saved in PNG format to `path.out`.
`fence`	Only observed if `ola = TRUE`. The limit for outlier detection as a multiplier of the interquartile range. Defaults to 2. Less outliers will be detected with higher values (not recommended).
`data`	Specification of one of the internal reference datasets (`⁠rds01⁠` to `⁠rds30⁠`). If given, the arguments `⁠path.in⁠`, `⁠file⁠`, `⁠set⁠`, and `⁠ext⁠` are ignored. For its use see the examples. If not given, defaults to `NULL` (i.e., import data from a file).
`option`	If `2` (default), the model will be evaluated by `⁠lme()⁠` of package `⁠nlme⁠`. The degrees of freedom of the treatment comparison will be equivalent to SAS’ `⁠DDFM=CONTAIN⁠` and Phoenix WinNonlin’s `⁠Residual⁠`. If `1` or `3`, the model will be evaluated by `⁠lmer()⁠` of package `⁠lmerTest⁠`. With `1` the degrees of freedom of the treatment comparison will be equivalent to SAS’ `⁠DDFM=SATTERTHWAITE⁠` and Phoenix WinNonlin’s `⁠Satterthwaite⁠`. `3` uses the Kenward-Roger approximation equivalent to Stata’s `⁠dfm=Kenward Roger (EIM)⁠`. If `⁠regulator = "HC"⁠`, only `1` or `3` are supported.

Details

The model for the estimation of CVwR is
lm(log(PK) ~ sequence + subject%in%sequence + period, data = data[data$treatment == "R", ])
where all effects are fixed.
The model for the treatment comparison is with the default option=2
lme(log(PK) ~ sequence + period + treatment, random = ~1|subject, data = data)
and with option=1, option=3
lmer(log(PK) ~ sequence + period + treatment + (1|subject), data = data)
where sequence, period, and treatment are fixed effects and subject(sequence) is a random effect.

Tested designs

4-period 2-sequence full replicates
⁠TRTR | RTRT⁠
⁠TRRT | RTTR⁠
⁠TTRR | RRTT⁠
2-period 4-sequence replicate
⁠TR | RT | TT | RR ⁠ (Balaam’s design)
4-period 4-sequence full replicates
⁠TRTR | RTRT | TRRT | RTTR⁠
⁠TRRT | RTTR | TTRR | RRTT⁠
3-period 2-sequence full replicates
⁠TRT | RTR⁠
⁠TRR | RTT⁠
3-period (partial) replicates
⁠TRR | RTR | RRT⁠
⁠TRR | RTR ⁠ (extra-reference design)

Data structure

Columns must have the headers subject, period, sequence, treatment, PK, and/or logPK.
Any order of columns is acceptable.
Uppercase and mixed case headers will be internally converted to lowercase headers.
- subject must be integer numbers or (any combination of) alphanumerics
  ⁠[A-Z, a-z, -, _, #, 0-9]⁠
- period must be integer numbers.
- sequence must be contained in the tested designs (numbers or e.g., ABAB are not acceptable).
- The Test treatment must be coded T and the Reference R.

Value

Prints results to a file if argument print = TRUE (default).
If argument print = FALSE, returns a data.frame with the elements:

`Design`	e.g., TRTR\|RTRT
`Method`	B-`option` (`1`, `2`, or `3`)
`n`	total number of subjects
`nTT`	number of subjects with two treatments of `T` (full replicates only)
`nRR`	number of subjects with two treatments of `R`
`Sub/seq`	number of subjects per sequence
`Miss/seq`	if the design is unbalanced, number of missings per sequence
`Miss/per`	if the design is incomplete, number of missings per period
`alpha`	nominal level of the test
`DF`	degrees of freedom of the treatment comparison
`CVwT(%)`	intra-subject coefficient of variation of the test treatment (full replicates only)
`CVwR(%)`	intra-subject coefficient of variation of the reference treatment
`swT`	intra-subject standard deviation of the test treatment (full replicates only)
`swR`	intra-subject standard deviation of the reference treatment
`sw.ratio`	ratio of intra-subject deviations of `T` and `R` (full replicates only)
`sw.ratio.CL`	upper confidence limit of `sw.ratio` (full replicates only)

If reference-scaling is applicable (i.e., CVwR(%) >30):

L(%) lower expanded limit of the acceptance range (AR)

U(%) upper expanded limit of the acceptance range (AR)
If reference-scaling is not applicable (i.e., ≤30):

BE.lo(%) lower limit of the conventional AR (⁠ 80⁠)

BE.hi(%) upper limit of the conventional AR (⁠125⁠)

`CL.lo(%)`	lower confidence limit of the treatment comparison
`CL.hi(%)`	upper confidence limit of the treatment comparison
`PE(%)`	point estimate of the treatment comparison (aka GMR)
`CI`	assessment whether the 100(1 – 2α) CI lies entirely within the acceptance range (`⁠pass\|fail⁠`)
`GMR`	assessment whether the PE lies entirely within the GMR-restriction 80.00--125.00% (`⁠pass\|fail⁠`)
`BE`	mixed (aggregate) assessment whether the study demonstrates bioequivalence (`⁠pass\|fail⁠`)
`log.half-width`	half-width of the confidence interval in log-scale

If ola = TRUE and at least one studentized outlier was detected:

`outlier`	outlying subject(s)
`CVwR.rec(%)`	intra-subject coefficient of variation of `R`; recalculated after exclusion of outlier(s)
`swR.rec`	intra-subject standard deviation of the reference treatment after exclusion of outlier(s)
`sw.ratio.rec`	ratio of intra-subjectstandard deviations of `T` and `R` after exclusion of outlier(s); full replicates only
`sw.ratio.rec.CL`	upper confidence limit of `sw.ratio.rec` (full replicates only)

If reference-scaling is applicable (i.e., CVwR.rec(%) >30):

L.rec(%) recalculated lower expanded limit of the AR

U.rec(%) recalculated upper expanded limit of the AR
If reference-scaling is not applicable (i.e., CVwR.rec(%) ≤30):

BE.rec.lo(%) lower limit of the conventional AR (⁠ 80⁠)

BE.rec.hi(%) upper limit of the conventional AR (⁠125⁠)

`CI.rec`	assessment whether the 100(1–2α) CI lies entirely within the new acceptance range (`⁠pass\|fail⁠`)
`GMR.rec`	assessment whether the PE lies entirely within the GMR-restriction 80.00--125.00% (`⁠pass\|fail⁠`)
`BE.rec`	mixed (aggregate) assessment whether the study demonstrates bioequivalence (`⁠pass\|fail⁠`)

Warning

Files may contain a commentary header. If reading from a CSV-file, each line of the commentary header must start with ⁠"# "⁠ (hashmark space = ⁠ASCII 35 ASCII 32⁠). If reading from an Excel-file all lines preceding the column headers are treated as a comment.

Clarification

The ‘ASCII line chart’ in the result file gives the confidence limits with filled black squares and the point estimate as a white rhombus. If a confidence limit exceeds the maximum possible expansion limit, it is shown as a triangle. Expanded limits are given as double vertical lines. Unscaled limits, the GMR restriction, and 100% are given with single vertical lines. The ‘resolution’ is approximatelly 0.5% and therefore, not all symbols might be shown. The CI and PE take presedence over the limits and the expanded limits over unscaled ones.

Disclaimer

Program offered for Use without any Guarantees and Absolutely No Warranty. No Liability is accepted for any Loss and Risk to Public Health Resulting from Use of this R-Code.

Note

The EMA’s model specified as ‘Method B’ in Annex I assumes 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 CV_wT and CV_wR can be estimated). Hence, amongst biostatisticians it is called the “crippled model” because the replicative nature of the study is ignored.
The method for calculating the degrees of freedom is not specified in the SAS code provided by the EMA in Annex I. Hence, the default in PROC MIXED, namely DDFM=CONTAIN is applied.
For incomplete data (i.e., missing periods) Satterthwaite’s approximation of the degrees of freedom (option = 1) or Kenward-Roger (option = 3) might be better choices – if stated as such in the statistical analysis plan.
The half-width of the confidence interval in log-scale allows a comparison of methods (B v.s. A) or options (2 v.s. 1). A higher value might point towards a more conservative decision. Quoting the Q&A-document:
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.
In the provided reference datasets – with one exception – the conclusion of BE (based on the mixed CI and GMR criteria) agrees between ‘Method A’ and ‘Method B’. However, for the highly incomplete dataset 14 ‘Method A’ was liberal (passing by ANOVA but failing by the mixed effects model).

Reference-scaling is acceptable for C_max (immediate release products) and C_max,ss, C_τ,ss, and _partialAUC (modified release products). However, quoting the BE guideline:
The applicant should justify that the calculated intra-subject variability is a reliable estimate and that it is not the result of outliers.
Quoting the Q&A on the Revised EMA Bioequivalence Guideline:
... 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.

The EMA’s method of reference-scaling for highly variable drugs / drug products is currently recommended in other jurisdictions as well (e.g., the WHO; ASEAN States, Australia, Belarus, Brazil, Chile, Egypt, the Eurasian Economic Union, the East African Community, New Zealand, the Russian Federation).

Health Canada’s variant of ABEL (upper cap of scaling ~57.4% limiting the expansion at 67.7–150.0%) is only approximate because a mixed-effects model would be required.

In a pilot phase the WHO accepted reference-scaling for AUC (4-period full replicate studies are mandatory in order to assess the variability associated with each product). It was an open issue how this assessment should be done. In Population Bioequivalence (PBE) and Individual Bioequivalence (IBE) the s_wT/s_wR ratio was assessed and similar variability was concluded for a ratio within 0.667–1.500. However, the power of comparing variabilities in a study designed to demonstrate ABE is low. This was one of the reasons why PBE and IBE were not implemented in regulatory practice. An alternative approach is given in the FDA’s draft ANDA guidance. Variabilities are considered comparable if the upper confidence limit of σ_wT/σ_wR is less than or equal to 2.5.
In 2021 the requirement of comparing variabilities was lifted by the WHO.

Author(s)

Helmut Schütz, Michael Tomashevskiy, Detlew Labes

References

European Medicines Agency, Committee for Medicinal Products for Human Use. Guideline on the Investigation of Bioequivalence. CPMP/EWP/QWP/1401/98 Rev. 1/ Corr **. London. 20 January 2010. online

European Generic and Biosimilar Medicines Association. 3^rd EGA Symposium on Bioequivalence. Questions and Answers on the Revised EMA Bioequivalence Guideline. London. 1 June 2010. online

European Medicines Agency, Committee for Medicinal Products for Human Use. Questions & Answers: positions on specific questions addressed to the Pharmacokinetics Working Party (PKWP). EMA/618604/2008 Rev. 13. London. 19 November 2015. online

European Medicines Agency. Clinical pharmacology and pharmacokinetics: questions and answers. 3.1 Which statistical method for the analysis of a bioequivalence study does the Agency recommend? Annex I. EMA/582648/2016. London. 21 September 2016. online

Executive Board of the Health Ministers’ Council for GCC States. The GCC Guidelines for Bioequivalence. Version 3.0. May 2021. online

Health Canada. Guidance Document. Conduct and Analysis of Comparative Bioavailability Studies. Ottawa. 2018/06/08. online

European Medicines Agency, Committee for Medicinal Products for Human Use. Guideline on the pharmacokinetic and clinical evaluation of modified release dosage forms. EMA/CPMP/EWP/280/96 Corr1. London. 20 November 2014. online

World Health Organization, Prequalification Team: medicines. Guidance Document: Application of reference-scaled criteria for AUC in bioequivalence studies conducted for submission to PQTm. Geneva. 22 November 2018. online

World Health Organization. Application of reference-scaled criteria for AUC in bioequivalence studies conducted for submission to PQT/MED. Geneva. 02 July 2021. online

U.S. Food and Drug Administration, Center for Drug Evaluation and Research. Draft Guidance for Industry. Bioequivalence Studies with Pharmacokinetic Endpoints for Drugs Submitted Under an ANDA. August 2021. download

Examples


# Importing from a CSV-file, using most of the defaults: variable
# separator colon, decimal separator period, no outlier-analyis,
# print to file.
# Note: You must adapt the path-variables. The example reads from
# the data provided by the library. Write-permissions must be granted
# for 'path.out' in order to save the result file. Here the default
# (R's temporary folder) is used. If you don't know where it is,
# type tempdir() in the console.
path.in <- paste0(find.package("replicateBE"), "/extdata/")
method.B(path.in = path.in, file = "DS", set = "01", ext = "csv")
# Should result in:
#   CVwT               :  35.16%
#   swT                :   0.34138
#   CVwR               :  46.96% (reference-scaling applicable)
#   swR                :   0.44645
#   Expanded limits    :  71.23% ... 140.40% [100exp(±0.760·swR)]
#   swT / swR          :   0.7647 (similar variabilities of T and R)
#   sw-ratio (upper CL):   0.9324 (comparable variabilities of T and R)
#   Confidence interval: 107.17% ... 124.97%  pass
#   Point estimate     : 115.73%              pass
#   Mixed (CI & PE)    :                      pass
#
# Internal reference dataset 01 used and results to R's temporary
# folder. Additional outlier-analyis and box plot saved as PNG.
method.B(ola = TRUE, plot.bxp = TRUE, data = rds01)
# Should give the same as above. Additionally:
#   Recalculation due to presence of 2 outliers (subj. 45|52)
#   CVwR (outl. excl.) :  32.16% (reference-scaling applicable)
#   swR  (recalc.)     :   0.31374
#   Expanded limits    :  78.79% ... 126.93% [100exp(±0.760·swR)]
#   swT / swR (recalc.):   1.0881 (similar variabilities of T and R)
#   sw-ratio (upper CL):   1.3282 (comparable variabilities of T and R)
#   Confidence interval: pass
#   Point estimate     : pass
#   Mixed (CI & PE)    : pass
#
# Same dataset. Show information about outliers and the model-table.
method.B(ola = TRUE, print = FALSE, verbose = TRUE, data = rds01)
# data.frame of results (full precision) shown in the console.
x <- method.B(ola = TRUE, print = FALSE, details = TRUE, data = rds01)
print(x, row.names = FALSE)
# Compare Method B with Method A for all reference datasets.

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 A and B agree.\n")
  } else {
    if (A == "fail" & B == "pass") {
      cat("Method A is conservative.\n")
    } else {
      cat("Method B is conservative.\n")
    }
  }
}
# should give
#   rds01: A pass, B pass - Methods A and B agree.
#   ...
#   rds14: A pass, B fail - Method B is conservative.
#   ...

# Health Canada: Only the PE of Cmax has to lie within 80.0-125.0%
# (i.e., no CI is required). With alpha = 0.5 the CI is practically
# supressed (zero width) and ignored in the assessment.
x    <- method.B(alpha = 0.5, regulator = "HC", option = 1,
                 data = rds03, print = FALSE, details = TRUE)[19:20]
x[1] <- round(x[1], 1) # only one decimal place for HC
print(x, row.names = FALSE)
# Should result in:
# PE(%)  GMR
# 124.5 pass

[Package replicateBE version 1.1.3 Index]

`method.A`	evaluation by a fixed effects model (ANOVA)
`ABE`	evaluation for conventional (unscaled) Average Bioequivalence

`L(%)`	lower expanded limit of the acceptance range (AR)
`U(%)`	upper expanded limit of the acceptance range (AR)

`BE.lo(%)`	lower limit of the conventional AR (`⁠ 80⁠`)
`BE.hi(%)`	upper limit of the conventional AR (`⁠125⁠`)

`L.rec(%)`	recalculated lower expanded limit of the AR
`U.rec(%)`	recalculated upper expanded limit of the AR

`BE.rec.lo(%)`	lower limit of the conventional AR (`⁠ 80⁠`)
`BE.rec.hi(%)`	upper limit of the conventional AR (`⁠125⁠`)