"De-Shrinking" EBEs: The Solution for Bayesian Therapeutic Drug Monitoring.

BACKGROUND: Therapeutic drug monitoring (TDM) aims at individualising a dosage regimen and is increasingly being performed by estimating individual pharmacokinetic parameters via empirical Bayes estimates (EBEs). However, EBEs suffer from shrinkage that makes them biased. This bias is a weakness for TDM and probably a barrier to the acceptance of drug dosage adjustments by prescribers.
OBJECTIVE: The aim of this article is to propose a methodology that allows a correction of EBE shrinkage and an improvement in their precision.
METHODS: As EBEs are defined, they can be seen as a special case of ridge estimators depending on a parameter usually denoted λ. After a bias correction depending on λ, we chose λ so that the individual pharmacokinetic estimations have minimal imprecision. Our estimate is by construction always better than EBE with respect to bias (i.e. shrinkage) and precision.
RESULTS: We illustrate the performance of this approach with two different drugs: iohexol and isavuconazole. Depending on the patient's actual pharmacokinetic parameter values, the improvement given by our approach ranged from 0 to 100%.
CONCLUSION: This innovative methodology is promising since, to the best of our knowledge, no other individual shrinkage correction has been proposed.

Key Points

Introduction

Therapeutic drug monitoring (TDM) is used for some drugs to control drug exposure by choosing an appropriate dosage regimen (i.e. drug dosage and frequency of administration) [1].

Classical TDM checks whether or not a concentration is within a concentration range if measured at a predefined time and considers the exposure (i.e. the area under the curve) to be abnormal when the concentration is outside this range. It is very simple to use in practice. However, this classical TDM assumes a perfect correlation between the measured concentration and the associated exposure, which is not the case for all patients. Indeed, two patients could have two very similar concentrations for two samples collected at the same time and still have different pharmacokinetic profiles and hence different exposures. Interpreting a measured concentration without considering the physio-pathological characteristics of the patients on the one hand and by referring to the same range of concentrations for all patients on the other hand can lead to inappropriate dose adjustments. The use of a population pharmacokinetic (PPK) model to analyse the concentrations allows a better evaluation of the patient's exposure than classic TDM, through an estimation of their own pharmacokinetic parameters [2, 3, 4]. These estimates are generally the empirical Bayes estimates (EBEs), computed by assuming that all the PPK model parameters are known without any imprecision [5, 6].

Unfortunately, EBEs can suffer from shrinkage [7]. Shrinkage is an average deviation of the individual pharmacokinetic (IPK) parameter estimate toward the population mean. Because shrinkage is a kind of bias (i.e. a systematic deviation), it should be possible to remove it. For example, Lavielle and Riba [8] proposed a method to limit the shrinkage when validating a PPK model. In the case of TDM, the problem is not an average but an individual deviation, so-called individual shrinkage. This individual shrinkage is specific to each patient (i.e. to their own pharmacokinetic parameter values), to the time when the blood samples are collected, and to the population's parameter values; this is why it is not usually removed [9].

Instead, an average shrinkage is provided for each pharmacokinetic parameter after a PPK model analysis. Usually this shrinkage is expressed as a percentage. When this average shrinkage is large, the general recommendation is to discard or at least not trust the IPK parameters calculated for a patient from only a few concentrations [9, 10]. Thus, a large shrinkage can be a limitation of TDM based on the estimation of the IPK parameters, but it is not the only one. Indeed, it has been shown [11-14] that it is not reasonable to expect a precise estimation of the pharmacokinetic parameters for a patient with only a few concentrations obtained from the patient, even using an optimal sampling strategy [15, 16, 17]. Even without any shrinkage, several combinations of pharmacokinetic parameters can lead to the drug concentrations observed in a patient. One way to quantify this imprecision is thus the length of the interval, or more generally, the volume of the prediction region of all the possible values of pharmacokinetic parameters that are consistent with the observed concentrations. This imprecision is intrinsic to the way the IPK parameters are estimated; in other words, this imprecision is intrinsic to the EBEs. Using the EBEs to estimate the IPK parameters cannot therefore lead to any significant improvements.

The purpose of this work is to propose an estimate of IPK parameters with no shrinkage and with less imprecision than the EBEs.

Materials and Methods

Models and Notations

When a PPK analysis is performed, a model describing the evolution of concentrations over time is used. This model is often written as in Eq. (1). where is the jth concentration measured in individual i at the time , is a vector containing the individual pharmacokinetic parameters for the ith individual, and g are two known functions, is the jth residual error for individual i at the time assumed to be independent of , h is a known function usually assumed to be , is a known matrix of covariates of the ith individual, is a fixed effect vector, is the vector of the random effects involved in inter-individual variability for the ith individual, is the variance–covariance matrix of the s, and is the vector of the parameters involved in the residual error model.

Data analysis with a PPK model as described by Eq. (1) provides estimates of , and , respectively denoted by , and . Of course, these estimates are imprecise. However, we will not take into account this uncertainty. In other words, herein, we consider these estimates to be fixed.

Once the model has been created and validated, it can be used daily for hospitalised patients. So, let us consider the concentrations observed at the time for a hospitalised patient with covariates for whom TDM is requested. TDM should be based on the patient’s IPK parameters , which are unknown when the blood sample is performed. However, a prediction of can be computed using both population parameters of the model and the patient’s concentrations .

A classical predictor of is the EBE obtained as , where (Eq. 2)

In a more detailed manner (Eq. 3),

The first term of the previous equation is a least squares term. The second term corresponds to a penalty. It penalises values deviating far from the population average.

In Eq. (3), when increases, decreases, whereas tends towards a fixed value. In this case, the patient’s concentrations govern the prediction and the shrinkage should be low. By contrast, when is low, dominates and the population parameters dominate the prediction, i.e. the is closer to the population mean than it should be. In this case, one can expect to be “shrunk” toward . Equation (3) can be interpreted as a special case of a criterion giving a ridge estimator, defined as (Eq. 4):

Actually, .

The only difference between Eqs. (3) and (4) is the presence of in . If the predictions of were good, should be close to , which is unknown. To circumvent this difficulty (i.e. being unknown), we evaluate the distance between and using simulations by proceeding as follows:

A set of K vectors are simulated where for Then, for each , are computed.

For each , are simulated where for j = 1,…,n.

For a given , one computes from Eq. (4) by using instead of .

In other words, for a vector containing covariate values for a hospitalised patient, it is possible to simulate K different sets of IPK parameters (line 1). For each set of IPK parameters, it is possible to simulate concentrations at the same sampling times as those of the hospitalised patient (line 2). From these concentrations , it is possible to estimate IPK parameters thanks to a Newton–Raphson algorithm (Eq. 4) for a given (line 3). Therefore, for K simulated patients, we have actual IPK parameters and predicted IPK parameters (line 3).

Once this is obtained, one can compute = . We can see whether is a good predictor of by drawing as a function of or equivalently as a function of If it is a good predictor, the points should be concentrated around the first bisector line. Unfortunately, there may exist a systematic departure from this line, signalling a bias. This is why we look for a function , which minimises the distance between the vectors and (Eq. 5):

Depending on the space where the functions are to be found, there should exist at least one function so that the minimum is reached. There exist many possible choices for . As any continuous functions can be locally approximated by polynomials, we chose as the set of the local polynomial functions. This polynomial function can be estimated, for example, by locally estimated scatterplot smoothing (LOESS) regression. Next, we look for the value , so that (Eq. 6):

Because , Eq. (6) can be solved by looking for lambda in an adaptive discretisation of ]0; ]. For the practical implementation, can be chosen close to 10.

At the end of these two optimisation steps, the best prediction with respect to the mean square error of their actual η for all patients having the same characteristics (i.e. same covariates values) and the same sampling times, is given by . Obviously, this predictor is not exact: several values from the same population correspond to a given value of obtained from the hospitalised patient (blue double arrow in Fig.1). The length of this range can be used as a prediction interval of and so the uncertainty of the estimation. A more traditional way of evaluating this uncertainty is to compute empirically a 95% prediction interval defined as the 2.5 and 97.5% percentiles of the s corresponding to the obtained from the patient of interest. The way to build these intervals (regression quantile methods) is not discussed in this article.

Fig. 1

The y-axis represents the that have been used to simulate the concentrations; the x-axis contains the obtained from the procedure described in section 2.1. The red line is the first bisector line. For a patient P, having a clearance individual random effect, denoted by , far below or above zero or, equivalently, an individual pharmacokinetic parameter far from the typical value, there is a systematic error. The systematic error is represented on the plot by the distance of the green curve from the red curve. The green curve represents . The blue band is the 95% prediction interval. Its width, for a given (i.e., for a given patient) gives an evaluation of the imprecision of our method

These ideas can be illustrated by the graph as shown in Fig. 1. To simplify, we concentrate on the component of linked to the clearance of the drug. We denote the kth simulated η for the clearance and we restrict as a function of clearance only.

The graph in Fig. 1 can be disturbing because intuitively, we would tend to represent the predicted values as a function of the simulated/true values. However, the aim is to minimise the difference between the true value of and its prediction (i.e. between and ). This presents itself as a classic nonlinear regression problem, in which we minimise the vertical distance between a point and its predictor.

A crude but global evaluation of the improvement that can be obtained using our prediction instead of the usual EBE can be computed by Eq. (7):

This non-negative quantity evaluates how much closer our prediction is to the actual than is the EBE (=1) given by . When this improvement is close to 0, the EBE and our predictor are close. As a consequence, our predictor does not improve the estimation of the individual pharmacokinetic parameter. On the contrary, any positive improvement shows that our predictor is better than the EBE. This indicator is a global evaluation that does not give any information on any specific patient located in a specific area of the green curve in Fig. 1.

We also provide estimates of the shrinkage for the EBE and for the estimator that we propose and determine the prediction of the shrinkage based on the method proposed by Combes et al. [14].

Finally, from this prediction, we deduce a prediction for IPK parameters, as (Eq. 8):

Once lambda and gamma ) are defined, they are used to improve predictions made for a given patient (Eq. 8). In a more detailed manner, when we have concentration(s) available for a patient P, simulations allow to compute and . Next, Eq. (4) is solved using the concentrations observed in the patient and with computed at the previous step instead of

From this computation, we obtain and therefore , which is the de-shrunk estimate of its clearance individual random effect. By taking into account the imprecision observed during the simulations, it is also possible to define the prediction interval from which the individual pharmacokinetic parameter is estimated.

For example, we estimate a for patient P equal to 0.2. Figure 1 shows that the true median value for the individual pharmacokinetic parameter is 0.18 and the prediction interval is between 0.05 and 0.29.

Simulations Study

Datasets were simulated using R version 3.5.2.

This methodology was applied to two different models from the literature. The performance evaluation of the proposed methodology was performed on 5000 subjects simulated by the models described in the corresponding articles.

Usually, it is necessary to perform K simulations for each subject; a subject being defined by its sampling times and its covariate values A. Since we simulated patients with exactly the same sample times and covariate values, K = 5000 simulations could be used for all simulated patients.

For the rest of this work, we focus on clearance that controls exposure and use only one concentration to estimate individual clearance as per hospital practice.

For each drug tested, the sampling times were chosen in line with the model itself (i.e. sampling times used for the model construction) and according to a clinical practicability allowing a later clinical use of this methodology.

First Model

The first model used relates to the pharmacokinetics of iohexol in dogs [18]. This molecule allows establishment of the glomerular filtration rate (GFR). As it is an exogenous molecule eliminated mainly by glomerular filtration, its clearance represents GFR. The model was built from iohexol pharmacokinetic data in 49 dogs. Each kinetic profile includes five points (i.e. 5, 15, 60, 90, and 180 min).

The PPK model is a two-compartment model parameterised with clearances and volumes. This implies the function of the first equation. The function was chosen as .

Part of the clearance variability is explained by two covariates: the health status and the plasma creatinine concentration. Health status is a dichotomous covariate defined on the basis of the International Renal Interest Society (IRIS). The plasma creatinine concentration is a continuous covariate.

This results in the h function of Eq. (1) being such that

The indicator function = 1 for a diseased individual and 0 otherwise (Table 1).

Table 1

Population pharmacokinetic parameters of iohexol published by Baklouti et al. [18]

	Value
Parameters
θ^_Cl	0.00432 L/min/kg
θ^_V_c	0.163 L/kg
θ^_V_p	0.0584 L/kg
θ^_Q	0.00336 L/min/kg
Covariates
θ^_1	− 0.345
θ^_2	− 0.517 dL/mg
Standard deviation of η=η_Cl,η_V_p,η_V_c
ω^_Cl	0.207
ω^_V_c	0.249
ω^_V_p	0.213
Residual error
Proportional	6.12%

The variance–covariance matrix of is assumed to be a diagonal matrix.

The dose chosen to perform the simulation is equal to 1600 mg. This corresponds to the average dose that was used to build the PPK model in the study by Baklouti et al. [18]. Similarly, we performed the simulations with sick dogs with a plasma creatinine level of 2.25 mg/dL that roughly corresponds to the average concentration reported in diseased animals [18].

We tested this methodology when a single sample is available at 180 min, corresponding to the best sampling time as identified in the study by Baklouti et al. [18].

Second Model

The second model relates the pharmacokinetics of isavuconazole [19]. This molecule is used in the treatment of invasive fungal infections. The model was built from 471 concentrations obtained from 79 patients.

The PPK model has two compartments. The error model for residual variability is a proportional error. Two covariates are included in the model, body mass index (BMI) and sex. BMI explains a part of peripheral distribution volume variability, and sex explains a part of clearance variability.

This results in the h function of Eq. (1) being such that

The indicator function = 1 for a man and 0 for a woman (Table 2).

Table 2

Population pharmacokinetics model of isavuconazole published by Wu et al. [19]

Parameters	Value
θ^_Cl	4.28 L/h
θ^_V_c	57.6 L
θ^_V_p	468 L
θ^_Q	37.4 L/h
Covariates
θ^_1(sex on clearance)	2.73 L/h
θ^_2(BMI on Vp)	39.7 L
Standard deviation of η=η_Cl,η_V_p
ω^_Cl	0.461
ω^_V_p	0.444
Residual error
Proportional	20.5%

BMI body mass index, Vp peripheral volume of distribution

The dose in the simulation was 200 mg, which corresponds to the standard dose administered in hospital. The covariate value for the BMI was 25.4 kg/m2, corresponding to the average value found in the study [19]. The mean clearance selected was the clearance for men because the model was constructed with a majority of men.

The sampling time selected to perform simulations was 195 h because this is when steady state is reached and TDM can be proposed [20].

Results

The results for simulated and predicted random effects on clearances are represented in Fig. 2. Figure 2a and b show the results for iohexol, and the Fig. 2c and d are for isavuconazole. On all plots, the vertical axis corresponds to . The horizontal axis represents the EBEs (i.e.) on the left plots (a and c) and on the right plots (b and d). Therefore, the left plots correspond to the results before shrinkage correction, and the right plots show the results obtained after our corrections. For the two tested models, the value is different from 1 (iohexol ; isavuconazole ), so are different from the and therefore from the EBEs. The difference between the median curve (obtained using LOESS regression) in green and the first bisector line in red represents the shrinkage.

Fig. 2

Graphical representation of as a function of (i.e. empirical Bayes estimate for clearance) (left) and of as a function of (right). The red curves represent the first bisector, the green curve represents the regression line, and the blue curves represent the prediction interval. After correction, the point clouds are refocused on the first bisector

Figure 2a and c shows that estimations of the IPK parameters provided by EBEs are not good: they are biased (i.e. the two clouds of points are not centred on the first bisector) and imprecise (i.e. the two clouds are thick).

Figure 2a (iohexol) shows that EBEs overestimated the IPK parameters, particularly for the lowest IPK parameter values. In this case, elimination clearance is overestimated and an unjustified dose increase is applied. Conversely, for isavuconazole, elimination clearance is underestimated for the lowest values, leading to an unjustified dose reduction. Thus, for drugs with a narrow therapeutic range, this kind of mistake can be detrimental to the patient; this means a therapeutic failure or, conversely, increased toxicity.

Surprisingly, it can be seen that, for these two experimental designs (covariates, sampling times, PPK model), shrinkage does not affect the high IPK parameter values. However, these results could be radically different for another experimental design and evolve in an unpredictable way. As a consequence, for each experimental design (covariates, sampling times, PPK model), it is necessary to repeat the shrinkage correction.

The second important information observed in these plots is that the IPK parameter estimates are imprecise. Indeed, one can see that a large number of actual IPK parameters are consistent with a single estimated IPK parameter, i.e. the cloud of points is thick. If the role of is to reduce the thickness of the cloud of points, it seems that this decrease is modest in these two examples. Indeed, there is an incompressible limit to imprecision, mainly governed by the number of available concentrations for the patient. Moreover, the imprecision varies according to the values of η: the further η is from the population average, the greater the imprecision.

Figure 3 shows the improvement obtained with our method (iohexol in Fig. 3a and isavuconazole in Fig. 3b), with respect to EBEs, for different values of departure from median population clearance (i.e. The improvement is maximal for patients with extreme actual clearance (lowest and highest values). Conversely, our method does not improve the clearance estimation for patients with actual clearances close to the median population clearance; this is logical. Finally, one would expect this improvement to be the same for extreme clearance values (small and large). This is true for isavuconazole but surprisingly is not the case for iohexol. To date, we have no explanation for this phenomenon, and it most likely depends on the experimental design used. The overall improvement, as defined by Eq. (7), is 21.2% for iohexol and 15.5% for isavuconazole.

Fig. 3

Graphical representation of the evolution of “Improv” of iohexol and isavuconazole. Improvement is much more important for values far from the average than for values close to the mean. The improvement of the estimates following the correction of the empirical Bayes estimate is therefore not always the same. So, this methodology does not provide a great correction for patients who have a pharmacokinetic profile close to the average. However, for pharmacokinetic profiles far from average pharmacokinetic profiles (i.e. patients for whom therapeutic drug monitoring is usually indicated), this methodology greatly improves results

Table 3 presents the results of shrinkage (for the EBE or the proposed estimator) and the prediction of shrinkage provided by the Combes et al. [14] method. The proposed estimator has a lower shrinkage than the EBE, whatever the application (i.e. the PPK model, the sampling times, etc.). Surprisingly, the method proposed by Combes et al. [14] does not provide a good prediction of EBE shrinkage, especially for isavuconazole.

Table 3

Shrinkage estimates calculated from empirical Bayes estimate and the proposed estimator.

	Empirical Bayes estimate (%)	Proposed estimator (%)	Shrinkagea (%)
Iohexol	40	20	24
Isavuconazole	57	49	9.3

aShrinkage of the proposed estimator computed using the Combes et al. [14] method

Discussion

We proposed a method to improve the prediction of IPK parameters using the measured concentration for one patient. The EBE is the traditional predictor; however, this individual prediction suffers from an “attraction” toward the typical population value and moves the prediction away from the actual value of the IPK parameters. This “distance” is called shrinkage [7]. It is usually computed over both different individuals (i.e. with different IPK parameters) and different sampling times. This number gives information on the average departure from the actual value (i.e. actual clearance, volume of distribution, etc.), but it cannot be used as it stands to correct the prediction of the IPK parameters for a particular patient [21, 22].

In contrast, the method we propose allows a specific shrinkage correction for each patient by considering the sampling times, the PPK model, and the patient’s covariate values. Moreover, for each patient, the level of shrinkage correction depends on the proximity of the patient’s pharmacokinetic parameters to the population's pharmacokinetic parameters. Indeed, we have shown that the correction was at its maximum for extreme values. This information is particularly important as, in hospital practice, when TDM is requested, the pharmacokinetic parameters of the patients are quite extreme values (i.e. far from the mean population values). This is the case for example for patients in the intensive care unit and those with cystic fibrosis or obesity, among others [23-25]. Thus, this approach allows us to expect a significant improvement in the estimates of the IPK parameters. In addition to TDM, a better estimate of EBEs could help to improve PPK models during their construction and validation steps.

As stated in Eq. (6), was chosen to minimise the shrinkage on clearance. Thus, nothing guarantees a shrinkage decrease of other individual pharmacokinetic parameters. This is a not problem because, for drugs with a linear kinetic, clearance drives exposure and then TDM. We estimated this pharmacokinetic parameter with only one concentration both to mimic hospital practice and to guarantee a high shrinkage to be removed. This choice does not exclude the application of the same kind of methodology for all the other pharmacokinetic parameters (i.e. to choose a and function specific to the pharmacokinetic parameter of interest).

In this work, the method of Combes et al. [18] poorly predicted EBE shrinkage. We consider that this is mainly due to the assumed asymptotic framework. It is clear that, with only one observation per individual, the variance of can be quite different from the corresponding Bayesian Fisher information matrix obtained by expanding at its mode. Actually, the quality of approximation of by a Gaussian distribution is difficult to anticipate.

One of the limits to our methodology is the need for a validated PPK model. Indeed, all these results have been obtained conditional on . In other words, an important assumption behind the method is that these estimations are precise enough to be used as if they were the actual value of the parameters. This hopefully happens for estimates obtained for a validated model built with a large number of individuals. If not, a model misspecification could lead to biased predictions. In this case, the shrinkage correction could be negligible compared with the bias. As such, the impact of a model misspecification on the final individual parameter estimate should be investigated.

However, although validation of PPK models for TDM is beyond the scope of this paper, the implicit assumption is that the chosen model properly describes the kinetic profile of the hospitalised patient. From our experience, this is far from the case. An easy way to identify a poor PPK model is to look at the residual terms. The most often used residual error models are the proportional error model (i.e. ) and in some cases a combined error model (i.e. ). A large (e.g. > 40%) indicates a poor PPK model and thus precludes any TDM based on IPK parameter estimations.

To conclude, this innovative methodology is promising since, to the best of our knowledge, no individual shrinkage correction has been proposed and it consistently provided better results than EBEs.

This study provides a methodology to suppress the individual shrinkage of the usual individual estimates of the pharmacokinetic parameters (empirical Bayes estimate).

This method decreases the estimate's imprecision.

The use of this methodology is illustrated for the estimations of the individual pharmacokinetic parameters of isavuconazole and iohexol.