Bioequivalence Calculations

How the calculation works

be_assess() (and be_compare()) are thin verbs over a single calculation path, the coordinator be_fit_models(), which runs these stages in order:

1. be_dataset() – validate the data, resolve the value column, drop excluded rows, detect the subject/sequence/period columns, set the reference as the first factor level, and carry any covariates.
2. be_design() – classify the design (parallel, 2x2x2 crossover, full or partial replicate) and report which frameworks it can support.
3. be_fit_model_single() – fit the average-BE model (this is the only place models are fit), dispatching on the model type to be_fit_model_lmer(), be_fit_model_nlme(), be_fit_model_gls(), or be_fit_model_anova().
4. be_extract_param() – extract the geometric means, the GMR and its confidence interval, and the within-formulation standard deviations.
5. be_table() – apply the regulatory limits or criterion and decide pass/fail (using the constants from be_regulator()), or, for the "descriptive" framework, report the estimates without a decision.

The model type is selected automatically from the design and the regulator; pass model_type to be_assess() to override it. A parallel study (one measurement per subject) is analyzed by a fixed-effects ANOVA, while a crossover or replicate study (repeated measures per subject) uses the mixed model. The FDA reference-scaled family always uses intra-subject contrasts ("isc") regardless of design: each subject’s mean log test minus mean log reference, from the subjects with every administration their sequence calls for, averaged with equal weight per sequence and with n - m degrees of freedom (n subjects, m sequences), as the FDA progesterone and warfarin guidances specify. Their within-subject variability comes from each subject’s difference between its two administrations of a formulation, centered within sequence (be_within_var(model_type = "isc")); the expanding-limits frameworks use the EMA ANOVA estimator. The two agree for a TRTR/RTRT full replicate. The unscaled 80.00-125.00% criterion of the FDA frameworks (RSABE below the switching variability, NTID, and HVNTID) uses the guidance’s mixed-effects model, with correlated treatment random effects per subject, a residual variance per treatment, and Satterthwaite degrees of freedom; with complete data it gives the same ratio and interval as the intra-subject contrasts.

The table also reports, on the measurement scale, the reference and test geometric means with their 90% confidence intervals (gm_reference/gm_test and their _lower/_upper bounds) and the endpoint units (taken from the PKNCAresults object; NA when not provided, and a column because they differ between C_max and AUC).

Optional packages

The bioequivalence calculations use lme4, lmerTest, and emmeans to fit the mixed model and estimate least-squares means. They are Suggested packages, so install them if they are not already present:

install.packages(c("lme4", "lmerTest", "emmeans"))

Simulating a 2x2 crossover study

A classic average-bioequivalence study uses a two-sequence, two-period crossover design: each subject receives both the test (T) and reference (R) formulations, in an order determined by their randomized sequence (RT or TR). We simulate concentration-time data from a one-compartment model with first-order absorption, with between-subject variability on clearance and a small difference between formulations.

simulate_crossover <- function(nsub = 24, seed = 20251201) {
  set.seed(seed)
  times <- c(0, 0.25, 0.5, 1, 1.5, 2, 3, 4, 6, 8, 12, 16, 24)
  sequence <- rep(c("RT", "TR"), length.out = nsub)
  subj_cl <- stats::rnorm(nsub, mean = log(5), sd = 0.25)
  dose <- 100
  v <- 50
  ka <- 1.1
  rows <- list()
  for (i in seq_len(nsub)) {
    forms <- if (sequence[i] == "RT") c("R", "T") else c("T", "R")
    for (p in 1:2) {
      form <- forms[p]
      # The test formulation has slightly higher exposure (lower clearance)
      cl <- exp(subj_cl[i] + ifelse(form == "T", -0.05, 0) + stats::rnorm(1, sd = 0.06))
      ke <- cl / v
      conc <- (dose * ka / (v * (ka - ke))) * (exp(-ke * times) - exp(-ka * times))
      conc <- conc * exp(stats::rnorm(length(times), sd = 0.05))
      conc[times == 0] <- 0
      rows[[length(rows) + 1]] <-
        data.frame(
          subject = i, sequence = sequence[i], period = p, formulation = form,
          time = times, conc = conc
        )
    }
  }
  do.call(rbind, rows)
}

conc_data <- simulate_crossover()
head(conc_data)
#>   subject sequence period formulation time      conc
#> 1       1       RT      1           R 0.00 0.0000000
#> 2       1       RT      1           R 0.25 0.5054207
#> 3       1       RT      1           R 0.50 0.8116205
#> 4       1       RT      1           R 1.00 1.1486459
#> 5       1       RT      1           R 1.50 1.5238739
#> 6       1       RT      1           R 2.00 1.5782127

Noncompartmental analysis

The crossover structure is described in the grouping formula. The grouping variables are the sequence, period, and formulation, with the subject as the final grouping variable.

dose_data <- conc_data[conc_data$time == 0, c("subject", "sequence", "period", "formulation")]
dose_data$dose <- 100
dose_data$time <- 0

o_conc <- PKNCAconc(conc_data, conc ~ time | sequence + period + formulation + subject)
o_dose <- PKNCAdose(dose_data, dose ~ time | sequence + period + formulation + subject)

intervals <-
  data.frame(start = 0, end = Inf, cmax = TRUE, auclast = TRUE, aucinf.obs = TRUE)
o_data <- PKNCAdata(o_conc, o_dose, intervals = intervals)

o_results <- pk.nca(o_data)

The NCA parameters are available with summary() and as.data.frame() as usual:

knitr::kable(head(as.data.frame(o_results)), digits = 2)

Bioequivalence assessment

be_assess() takes the PKNCAresults object, the column that identifies the formulation (reference_col), the value of that column that is the reference (reference_value), and the regulatory framework. This 2x2 study is not replicated, so only average bioequivalence (regulator = "ABE") is feasible.

be <- be_assess(
  o_results,
  endpoints = c("cmax", "auclast", "aucinf.obs"),
  reference_col = "formulation",
  reference_value = "R",
  regulator = "ABE"
)
#> Warning: Units were not found in the data; omitting the `units` column. Supply
#> units via a PKNCAresults object or `PPSTRESU`/`PPORRESU` columns.

knitr::kable(as.data.frame(be), digits = 2)

Interpreting the results

For each endpoint and test formulation the table reports:

• units: the measurement units of the endpoint (e.g. ng/mL for C_max, hr*ng/mL for AUC), taken from the PKNCAresults object,
• gm_reference, gm_test (with _lower/_upper bounds): the reference and test geometric means and their 90% confidence intervals, on the measurement scale (in units),
• gmr_percent: the geometric mean ratio (test / reference), as a percentage,
• ci_lower, ci_upper: the bounds of the confidence interval for the GMR (90% by default, controlled by alpha),
• cvwr_percent, swr: the within-reference coefficient of variation and standard deviation (used by the reference-scaling frameworks; NA for unscaled ABE),
• limit_lower, limit_upper: the acceptance limits, and criterion for the FDA reference-scaled bound, and
• pass: the regulatory decision.

Two formulations are declared bioequivalent under ABE when the 90% confidence interval for the geometric mean ratio falls entirely within 80.00%-125.00% for the primary endpoints (commonly AUC and C_max).

For ABE the average-BE model fit for each endpoint is

\[\log(\text{value}) = \text{sequence} + \text{period} + \text{formulation} + (1 \mid \text{subject})\]

fit with lmerTest::lmer() so that Satterthwaite degrees of freedom are available, and the confidence interval is obtained by exponentiating the difference in least-squares means between the test and reference formulations.

Reference scaling and regulatory frameworks

Highly variable drugs (within-subject coefficient of variation above 30%) are difficult to show bioequivalent against the fixed 80.00-125.00% limits, so regulators allow the limits to be scaled by the within-reference variability, estimated from a replicate design in which the reference is given more than once. PKNCA implements the major frameworks through be_assess():

• ABE – average bioequivalence, the unscaled 80.00-125.00% criterion.
• EMA, Health Canada (HC), and GCC – average bioequivalence with expanding limits (ABEL); the limits widen with the within-reference standard deviation, with framework-specific caps.
• FDA RSABE – reference-scaled average bioequivalence using the linearized Howe/Hyslop criterion.
• FDA NTID and HVNTID – the (highly variable) narrow therapeutic index drug frameworks, which add a constraint on the ratio of within-subject standard deviations.

All regulatory constants are built into PKNCA (see be_regulator()); no external package is required.

simulate_replicate <- function(nsub = 24, seed = 20251015, cv_wr = 0.45, cv_wt = 0.40) {
  set.seed(seed)
  seqs <- rep(c("TRTR", "RTRT"), length.out = nsub)
  swR <- sqrt(log(cv_wr^2 + 1))
  swT <- sqrt(log(cv_wt^2 + 1))
  subj <- stats::rnorm(nsub, sd = 0.4)
  rows <- list()
  for (i in seq_len(nsub)) {
    trt <- strsplit(seqs[i], "")[[1]]
    for (p in seq_along(trt)) {
      sw <- if (trt[p] == "R") swR else swT
      eff <- if (trt[p] == "T") log(1.05) else 0
      rows[[length(rows) + 1]] <- data.frame(
        subject = i, sequence = seqs[i], period = p, treatment = trt[p],
        PPTESTCD = "auclast",
        PPORRES = exp(log(100) + eff + subj[i] + stats::rnorm(1, sd = sw))
      )
    }
  }
  do.call(rbind, rows)
}

rep_data <- simulate_replicate()

ema <- be_assess(rep_data, reference_col = "treatment", reference_value = "R",
                 endpoints = "auclast", regulator = "EMA")
#> Warning: Units were not found in the data; omitting the `units` column. Supply
#> units via a PKNCAresults object or `PPSTRESU`/`PPORRESU` columns.
knitr::kable(as.data.frame(ema), digits = 3)

fda <- be_assess(rep_data, reference_col = "treatment", reference_value = "R",
                 endpoints = "auclast", regulator = "FDA")
#> Warning: Units were not found in the data; omitting the `units` column. Supply
#> units via a PKNCAresults object or `PPSTRESU`/`PPORRESU` columns.
knitr::kable(as.data.frame(fda)[, c("gmr_percent", "ci_lower", "ci_upper",
                                    "cvwr_percent", "swr", "criterion", "pass")],
             digits = 4)

cmp <- be_compare(rep_data, reference_col = "treatment", reference_value = "R",
                  endpoints = "auclast",
                  regulators = c("EMA", "HC", "GCC", "FDA", "NTID", "HVNTID"))
#> Warning: Units were not found in the data; omitting the `units` column. Supply
#> units via a PKNCAresults object or `PPSTRESU`/`PPORRESU` columns.
summary(cmp)
#>   endpoint  EMA   HC  GCC  FDA NTID HVNTID
#>  auclast:T TRUE TRUE TRUE TRUE TRUE   TRUE
#> 
#> Caption: Pass/fail by endpoint and regulator.

Unequal variances, covariates, and descriptive comparisons

Three options extend the same calculation path beyond a regulatory bioequivalence decision.

• heteroscedastic = TRUE estimates a separate residual variance for each treatment. Crossover and replicate studies use nlme::lme() with varIdent(~ 1 | treatment), and parallel studies use nlme::gls() with the same variance structure and Satterthwaite degrees of freedom. The lmer model cannot do this, because lme4::lmer() fits a single residual variance, so heteroscedastic = TRUE with model_type = "lmer" is an error. Without reference scaling the nlme model accepts non-replicated crossovers and several test formulations; with reference scaling it needs a fully replicated design.
• covariates names columns added to every model as additive fixed effects. The least-squares means are then at the mean of numeric covariates. The intra-subject contrasts and the FDA mixed model used by the FDA frameworks and the within-subject variances used for reference scaling do not use covariates.
• regulator = "descriptive" estimates the geometric means, their ratio, and its confidence interval exactly as "ABE" does, but applies no acceptance limits and reports no pass/fail. It is meant for comparisons that are estimated rather than judged, such as food effect and drug-drug interaction studies. The caption states that no regulatory decision was applied.

As an example, a fed/fasted food-effect crossover, where the fed state is more variable than the fasted state, is summarized descriptively with treatment-specific variances and body weight as a covariate.

simulate_food_effect <- function(nsub = 20, seed = 20260927) {
  set.seed(seed)
  sequence <- rep(c("FastedFed", "FedFasted"), length.out = nsub)
  weight <- round(stats::rnorm(nsub, mean = 75, sd = 12))
  subj <- stats::rnorm(nsub, sd = 0.3)
  rows <- list()
  for (i in seq_len(nsub)) {
    states <- if (sequence[i] == "FastedFed") c("Fasted", "Fed") else c("Fed", "Fasted")
    sw <- ifelse(states == "Fed", 0.35, 0.15)
    rows[[i]] <- data.frame(
      subject = i, sequence = sequence[i], period = 1:2, state = states,
      weight = weight[i], PPTESTCD = "cmax", PPORRESU = "ng/mL",
      PPORRES = exp(log(50) + ifelse(states == "Fed", log(1.4), 0) +
                      0.01 * (weight[i] - 75) + subj[i] + stats::rnorm(2, sd = sw))
    )
  }
  do.call(rbind, rows)
}

food <- be_assess(
  simulate_food_effect(),
  reference_col = "state", reference_value = "Fasted", endpoints = "cmax",
  regulator = "descriptive", heteroscedastic = TRUE, covariates = "weight"
)
knitr::kable(as.data.frame(food), digits = 2)

The table has no limit_lower, limit_upper, criterion, or pass columns. Its caption records the model and states that no decision was applied:

attr(food, "caption")
#> [1] "Descriptive treatment comparison (90% CI). Geometric means and the geometric mean ratio are least-squares means from a mixed-effects model with treatment-specific residual variances (nlme::lme); the 90% CI comes from exponentiated least-squares-mean differences (emmeans). No regulatory decision was applied: there are no acceptance limits and no pass/fail. The model includes weight as additive covariate(s); least-squares means are at the mean of numeric covariates and averaged over the levels of factor covariates."

Reporting checklist

The assessment table supports common bioequivalence reporting requirements:

1. Framework and model: the regulator and model_type columns record the framework and the average-BE model used, and the table’s caption attribute gives a one-sentence methods statement (model, confidence-interval method, within-subject variability, and the decision rule).
2. %GMR, confidence interval, within-subject CV, limits/criterion, and the pass/fail decision: one row per endpoint and test formulation in as.data.frame(be_assess(...)).
3. Cross-regulator comparison: summary(be_compare(...)) gives an endpoint-by-regulator pass/fail grid.
4. Individual data: available from as.data.frame(o_results) for listings and figures.