Multivariate Assessment for Bioequivalence Based on the Correlation of Random Effect

Introduction

Bioequivalence tests are conducted to assess equivalence in bioavailability between a test and reference drug product.¹ When they are pharmaceutically comparable based on the active ingredient, formulation, dosage strength, and route of administration, and their pharmacokinetic profiles are not significantly different, the products are therapeutically equivalent.^1–4 Bioequivalence is commonly evaluated using the area under the concentration-time curve (AUC) and the peak concentration (C_max), which are assumed to be log-normally distributed.⁵ If the 90% confidence interval (CI) for the geometric mean ratio (GMR) of a test product to a reference product falls within the bioequivalence limit of 0.8 to 1.25 for both AUC and C_max, the test product is considered to be bioequivalent.⁵ Bioequivalence tests have also been conducted in comparative pharmacokinetic studies and drug–drug interaction studies, where similarity in pharmacokinetic profiles is of interest.

In general, noncompartmental pharmacokinetic analysis has been considered as the standard approach for bioequivalence testing.⁶ However, because of certain limitations of noncompartmental analysis including sufficient samples at time points post dose, observed sampling time points for AUC_0–t/AUC_inf > 0.8, systemic circulation related to the exposure of site of action and linear pharmacokinetics, the model-based approach may be more appropriate in specific situations.⁶ The model-based approach is individual-based, separately fitted data and is able to evaluate the pure metric of rate reflecting the individual absorption rate constants (Ka).⁶ As a result, model-based approaches can be beneficial for drugs with complex pharmacokinetics including nonlinear pharmacokinetics, non-elimination from the sampling compartment and endogenous homeostatic feedback mechanism, which may not be applicable to the assumption of noncompartmental analysis.⁶ Furthermore, biosimilars and topical dermatological drugs are preferable candidates to apply the model-based approach due to the comparison of rate and extent of exposure at the site of action.⁶.

Using a linear mixed model (LMM) with period, sequence, and treatment as fixed effects and subject as a random effect, and fitting this model with log(AUC) or log(C_max) as a response, the GMR, the exponent of the estimate for the treatment effect, can be obtained for bioequivalence testing.¹ Although, 80–125% for the ratio of the product averages was suggested as the bioequivalence criteria in general, in the case of highly variable drugs with an expected within-subject variability of 30% or greater in the pharmacokinetic measures C_max or AUC (ie), concerns about the bioequivalence limit have been raised.⁷ This is because the high variability with a wide 90% CI for the GMR makes it difficult to be included in the standard bioequivalence limit.⁷

To address the issues with highly variable drugs, regulatory agencies generally agree that direct expansion of the bioequivalence limits, enlargement of the sample size, or special study designs including pharmacokinetic evaluation at steady state or replicated design are required.^3,7,8 For broadening of the bioequivalence limits, the acceptance range for C_max can be widened to a maximum of 0.6984–1.4319 according to within-subject CV of C_max.³ However, these proposed solutions for drugs with high intrasubject variability gave rise to safety or efficacy concerns and ethical issues regarding unnecessarily large numbers of subjects.⁷

It is well known that AUC and C_max are strongly correlated because they are two interrelated parameters in the same subject that represent the extent and rate of the exposure to a drug, respectively.^9,10 However, usually the 90% CIs of GMR between a test and reference product have been estimated separately for AUC and C_max assuming that these endpoints are independent.⁵ This could lead to wider 90% CIs, thereby increasing the probability that a test product fails to meet bioequivalence criteria, with respect to the reference product, when in fact, it is bioequivalent.^4,7

A multivariate approach to address bioequivalence has been applied in the last few decades. Chinchilli and Elswick suggested a multivariate assessment based on the acceptable region of the reference profile and another region of the test profile.¹¹ Quan et al developed a simultaneous equivalence assessment and an intersection union test for testing on multiple endpoints, which were developed by Berger and Hsu, and Wang et al.^12–14 In addition, another Bayesian modeling of multivariate average bioequivalence was introduced by Ghosh and Gonen.⁹ Most methods are likelihood approaches that consider all endpoints, including AUC and C_max, come from a multivariate normal distribution, resulting in correlated parameter estimates and test statistics.⁴ Lee et al developed multivariate models by assuming correlations among random effects for different responses.¹⁵ They showed several examples that have multiple correlated responses from the same subject and compared the correlated model to independent models.¹⁶

The objective of the present work is to apply multivariate analysis to the bioequivalence test in practice and to compare the performances of this model and separate models. To this end, three datasets that have been described previously by Choi et al and Yi et al were analyzed.^17–19 In particular, two datasets by Yi et al are interesting case studies because of the 90% CI of the GMR being out of the equivalence limit.^17,18 More precisely, one is a comparative pharmacokinetic study that has a 90% CI of the GMR for C_max slightly beyond the limit, and the other is a drug–drug interaction study in which the product consists of highly variable drugs with large within-subject variability C_max. Using the data shared by the authors, our goal was to build a multivariate model for two co-primary endpoints, AUC and C_max, which incorporated the correlation between responses by multivariate random effects. We sought to compare both 90% CIs for the GMR of C_max using the multivariate model and separate models.

The remainder of this paper is organized as follows. First, we describe the three datasets being considered and remind the reader of the bioequivalence test problem and conventional LMM. This is followed by an explanation of the multivariate HGLM used to test bioequivalence in detail. In the following section, the results of different analyses for the three datasets are reported and discussed before concluding.

Materials and Methods

Pharmacokinetic Study Data

Our analysis was based on three pharmacokinetic datasets that have been described previously in detail by Choi et al and Yi et al.^17–19 The names of the patients were anonymized, and the patients provided informed consent.

The fixed-dose combination (naproxen and esomeprazole) data were obtained from a comparative pharmacokinetic study of a newly developed (HCP1004) and marketed (VIMOVO^®) fixed-dose combination of naproxen (500mg) and esomeprazole (20mg); it was an open-label, randomized, single dose, 2-treatment 2-sequence crossover study in 66 healthy male subjects.¹⁹ In each period, a reference (VIMOVO^®) or a test (HCP1004) drug was administered orally, and blood samples for PK analysis were obtained for 72 h post dose. The C_max and AUC_0-t of naproxen and esomeprazole were estimated, to compare the pharmacokinetic parameters between the two drugs using two separate LMMs (1) with the period, sequence, and treatment as fixed effects and subject as random effects.

The pharmacokinetic results of tramadol were obtained from a comparative pharmacokinetic study between an extended-release (ER) and immediate-release (IR) fixed-dose combination (tramadol, 37.5 mg/acetaminophen 325 mg); it was an open-label, multiple dose, randomized, 2-sequence, 2-period crossover study in 12 healthy male subjects.¹⁷ All subjects were orally administered both formulations: IR formulations for 4 days followed by ER formulations for 4 days, or vice versa. The C_max and AUC_0–12 at steady state (on day 4) were calculated to compare the pharmacokinetic profiles of ER and IR formulations using the two separate LMMs (1).¹⁷

The fimasartan data were obtained from a drug–drug interaction study between fimasartan and amlodipine; it was an open-label, multiple-dose, 1-sequence 2-period crossover study in 19 healthy male subjects.¹⁸ All subjects received fimasartan 120mg alone for 7 days in period I. With a 5-day washout period, they were administered amlodipine 10mg and fimasartan 120mg for 7 days in period II.¹⁸ The effect of the co-administration of fimasartan and amlodipine at steady-state was investigated by the C_max and AUC_0–24 on day 7 using two separate LMMs (1) without period and sequence effects.

Linear Mixed Models and Multivariate Hierarchical Generalized Linear Models

Let be the repeatedly measured outcome from subject i in period j with a crossover design (i = 1, …, m; j = 1, …, q). We then considered a linear mixed model (LMM) with period, sequence, and treatment as fixed effects and subject as random effects.

(1)

where for the test and for the reference, and are covariates indicating period j and the sequence of subject i, and is the subject-specific random effect.⁴ The between (inter)-subject variance and within (intra)-subject variance are presented by variances and, respectively. The standard error of the treatment effect determines the width of the 90% CI for GMR of a test to the reference, and it depends largely on the within-subject variance .⁴ The within-subject coefficient of variation is defined as

(2)

When we applied model (1) for log(AUC) or log(C_max) as a response, then , the estimate for the GMR of the test to the reference, can be calculated for bioequivalence assessment.⁴

If there are two or more coprimary endpoints that are correlated, such as C_max and AUC, in clinical pharmacology data, multivariate HGLMs suggested by Lee et al can be applied.¹⁵ Let and be jth two outcomes from subject i (i = 1, …, m; j = 1, …, q), and υ_1i and υ_2i be the subject-specific random effects that are correlated.⁴ It is assumed that and are conditionally independent given and ⁴ Then, the following bivariate HGLM can be proposed.⁴

(3)

where , the vector of two subject random effects, follows a multivariate normal distribution with variance-covariance matrix composed of , , and . The correlation coefficient of subject random effects satisfies .⁴ Correlations between and within a same subject i are explained through the correlations between subject-specific random effects and .⁴ If , then two outcomes and from the same subject are independent and this model is equivalent to two conventional separate LMMs (1).⁴ If , then it is equivalent to two conventional separate LMMs (1).⁴ These two LMMs are nested in the bivariate HGLM, and the likelihood ratio (LR) test can be used to select a better model.⁴ All statistical analyses were conducted using the R software (version 3.6.1) with the mdhglm package and the R codes are provided as Supplement 1.

Results

Fixed-Dose Combination (Naproxen and Esomeprazole) Data

In the original study of Choi et al, the C_max and AUC_0-t were separately analyzed by the LMM (1) with the period, sequence, and treatment as fixed effects and subject as random effects assuming independence between two responses for each ingredient (naproxen or esomeprazole).¹⁹ As a result of the conventional analysis (using separate two LMMs), the GMRs and 90% CIs (HCP1004 to VIMOVO^®) of the C_max and AUC_0–72 of naproxen were 0.993 (0.936, 1.054) and 0.996 (0.984, 1.009), respectively. For esomeprazole, the GMRs (with 90% CI) for C_max and AUC_0–12 were 0.985 (0.825, 1.177) and 1.038 (0.915, 1.178), respectively (Table 1).

Table 1 Summaries of Multivariate (Bivariate HGLM) and Conventional Analysis (Separate LMMs) for the Fixed-Dose Combination (Naproxen and Esomeprazole) Data

We also performed multivariate analysis (bivariate HGLM) to assess the bioequivalence between the two products. For naproxen, estimates and standard errors for the treatment effects in both C_max and AUC_0–72 by the bivariate HGLM (3) are the same as those of the separate LMMs; therefore, GMRs with 90% CIs remained the same for both C_max and AUC_0–72. In the case of esomeprazole, the standard error of the treatment effect for C_max by the bivariate HGLM (3) was 0.095, under 0.108, that by the separate LMMs, and this leads to a narrower 90% CI of the GMR for C_max ((0.843, 1.152) vs (0.825, 1.177)) while remaining the same for AUC_0–12.

Based on the restricted likelihood, LR was calculated to test the statistical significance of the correlation coefficient ρ. Because the independent model (separate two LMMs) is nested in the correlated model (bivariate HGLMs), the LR test for is used to select a better model between two separate LMMs and the bivariate HGLM. It rejects the null hypothesis because of the deviance difference ; thus, the LR test selects the bivariate HGLM as a better model.

Figure S1a and b in the Supplement 2 show that there are positive correlations between subject random effects of each LMM for log(C_max) and log(AUC_0-t) in naproxen and esomeprazole (r = 0.286 and r = 0.883, respectively). For naproxen in the fixed-dose combination data, the correlation coefficient r = 0.286 is relatively small and the reported within-subject coefficient of variation () in a previous study was less than 30%; therefore, the estimated results using bivariate HGLM were unchanged.²⁰ However, in the case of esomeprazole, the correlation coefficient (r = 0.883) and reported in a previous study were all relatively high () and the results from bivariate HGLMs showed obvious differences.²⁰

Tramadol Data

In the original study by Yi et al, the C_max and AUC_0–12 at steady state were separately analyzed by the LMM (1) with the period, sequence, and treatment as fixed effects and subject as random effects assuming independence between two responses.¹⁷ The GMR of the extended-release (ER) to immediate-release (IR) formulation for AUC_0–12 was 0.947 with a 90% CI (0.910, 0.986). However, the C_max was 0.866 with a 90% CI (0.797, 0.941), which was slightly lower than the bioequivalence criteria. Moreover, Figure S1c in the Supplement 2 shows strong positive correlations between subject random effects for log(C_max) and those for log(AUC_0–12) in the tramadol results (r = 0.966) with relatively low within-subject coefficient of variation ( = 0.106 as maximum) as reported in a previous study. We can expect different results to come from the bivariate HGLMs.²¹

The results of the two conventional separate LMMs and bivariate HGLMs are summarized in Table 2. The LR test for rejects the null hypothesis because the deviance difference between separate LMMs and the bivariate HGLM is ; thus, the LR test selects the bivariate HGLM. The standard error for the treatment effect in C_max by the bivariate HGLM (3) was 0.040, under 0.046, that by the separate LMMs, and this led to a narrower 90% CI of the GMR for C_max. As a result of the bivariate HGLM, both the 90% CI of GMR for C_max and AUC_0–12 were included in the bioequivalence limit (0.8, 1.25), and the ER formulation was equivalent to the IR formulation.

Table 2 Summaries of Multivariate (Bivariate HGLM) and Conventional Analysis (Separate LMMs) for the Tramadol Data

Fimasartan Data

Yi et al used conventional independent LMMs (1) without sequence and period effects for the analysis of C_max and AUC_0–24 at steady state.¹⁸ The GMR and 90% CI of co-administration of fimasartan with amlodipine to fimasartan alone for C_max were 1.096 and 0.761–1.579, respectively, and those for AUC_0–24 were 1.163 and 1.009–1.341, respectively. These findings indicate that fimasartan is a highly variable drug with large within-subject variability for C_max () and moderate within-subject variance for AUC (), similar to that reported in previous studies.^22,23 As shown in the tramadol data, the fimasartan data have strong positive correlations (r = 0.832) between subject random effects for log(C_max) and those for log(AUC_0–24) in Figure S1d. Therefore, there are differences between the results obtained using two separate LMMs and a bivariate HGLM.

We performed multivariate analysis using bivariate HGLM (3) without sequence and period effects on fimasartan data. A summary of the results of the two methods is shown in Table 3. The LR test for H₀: ρ = 0 rejects the null hypothesis because the deviance difference between separate LMMs and the bivariate HGLM is ; thus, the LR test selects the bivariate HGLM. The bivariate HGLM yields smaller standard errors than conventional LMMs for treatment effects, which results in narrower 90% confidence intervals.

Table 3 Summaries of Multivariate (Bivariate HGLM) and Conventional Analysis (Separate LMMs) for the Fimasartan Data

Discussion

In the assessment of bioequivalence, a valid statistical evaluation of pharmacokinetic endpoints representing bioavailability, including AUC and C_max, is required to guarantee the safety and efficacy of the therapeutic effect.⁹ While AUC is the most accepted measure to assess bioequivalence, no single parameter can be concluded to be universally superior to the other endpoints.⁹ In general, the extent of absorption is evaluated by AUC, and C_max is often used for the rate of absorption.²⁴ Because they are highly correlated, C_max is not enough to reflect the pure rate of absorption, and alternative markers such as the ratio of C_max/AUC or time to achieve maximum concentration (T_max) were proposed instead of C_max.^2,9,24 Although these markers were relatively independent of AUC, their usefulness is limited due to the low power caused by their discrete nature or the assumption of normality for random effects.^25,26 Essentially, bioequivalence is a multivariate endpoint from the same biological sample measurements, and correlation between primary endpoints has to be considered for clearly appropriate analysis.^9,27 However, currently accepted practice usually involves univariate bioequivalence analysis.⁹ Although several average bioequivalence (ABE) studies have been generated, there are a limited number of multivariate analyses.^9,12,25,26

In the present study, we adapted an advanced bioequivalence test model based on multivariate analysis for pharmacokinetic studies and compared the results from a multivariate HGLM and two conventional independent LMMs in real datasets. Because pharmacokinetic parameters such as AUC and C_max are highly correlated, these correlations need to be incorporated into the multivariate model, which deals with multiple endpoints simultaneously. It is clear that the multivariate HGLM, including the correlation structure between the two responses, leads to smaller standard errors for estimates of treatment effects. This is because when point estimates of fixed effects do not change and total variations remain the same, the additional parameter in the structured variance-covariance matrix increases the proportion of the explained variance by the model with respect to the total variance. In other words, the correlation parameter fed into the model improved the power of explanation for intrasubject variations; thus, it was observed that the GMRs did not change, but 90% confidence intervals were prone to shrink (Figure 1).

Figure 1 Geometric mean ratio (GMR) and 90% confidence interval (CI) using separated two linear mixed models (LMMs) and bivariate hierarchical generalized linear model (HGLM).

According to regulatory agencies, including the US Food and Drug Administration (FDA), the ABE approach for evaluating differences in exposure between test and reference formulations using both AUC and C_max has been suggested, which requires equivalence between the population means of two formulations.^7,9,28,29 However, traditional bioequivalence studies based on this ABE assessment do not exclude the possibility of differences in exposure at an individual level by intrasubject variability or subject-by-formulation effects.²⁹ In a previous study, this intrasubject variability may have played a crucial role in the variation of drug exposure rather than formulation-dependent variations.²⁹ Despite the benefit for intrasubject variability, due to the limitation of tolerability in healthy subjects and sensitivity of differences in exposure, single-dose bioequivalence studies are still preferred for regulatory acceptance in general.³⁰ The width of the 90% CI is proportional to the drug variability, especially intrasubject variability in a crossover design, and inversely proportional to the number of subjects.⁷ As a result, proper and reduced estimation of the 90% CI of primary endpoints using a multivariate model could be attributed to data-specific bioequivalence analysis. Because this study was not conducted on pure bioequivalence datasets between same ingredients and formulations, further investigation on the intrasubject variability from partial replicate design would be required.

In conclusion, the proposed multivariate bioequivalence approach using a multivariate HGLM helps to reduce the 90% CI by adding correlation parameters and thus has an advantage in evaluating the bioequivalence of highly variable drugs or biosimilars with large intrasubject variability.

Conclusions

In this study, we applied multivariate analysis (a multivariate HGLM) to the bioequivalence test in practice and compared the performances of the multivariate model and the conventional separate models. For real data analysis, previously studied fixed-dose combinations (naproxen and esomeprazole), tramadol and fimasartan datasets were introduced, and we developed a multivariate model for two endpoints, AUC and C_max, in which the correlations between responses were incorporated by multivariate normally distributed random effects. As a result, the 90% CIs for the GMR of C_max decreased, and both AUC and C_max showed bioequivalence in the tramadol data. Because of its highly variable nature, the fimasartan data did not meet the standard bioequivalence criteria (0.8, 1.25). However, it is clear that multivariate HGLM performs well with the bioequivalence test for multiple endpoints. This method is easily applicable to more general cases and consequently will help resolve a number of issues in bioequivalence evaluation without expanding the bioequivalence limit or requiring a change in study design to a replicated design. However, further investigation of other datasets is required to validate these findings.