# Empirical Likelihood Estimation for Population Pharmacokinetic Study Based on Generalized Linear Model

Authors: Fang-rong Yan; Jin-guan Lin; Yuan Huang; Jun-lin Liu; Tao Lu
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

To obtain efficient estimation of parameters is a major objective in population pharmacokinetic study. In this paper, we propose an empirical likelihood-based method to analyze the population pharmacokinetic data based on the generalized linear model. A nonparametric version of the Wilk's theorem for the limiting distributions of the empirical likelihood ratio is derived. Simulations are conducted to demonstrate the accuracy and efficiency of empirical likelihood method. An application illustrating our methods and supporting the simulation study results is presented. The results suggest that the proposed method is feasible for population pharmacokinetic data.


---

## Body

## 1. Introduction

The parameter estimation in population pharmacokinetics (PPK) is a significant problem in clinical research. It is a novel problem in pharmacokinetic (PK) study that combines classic PK models with group statistical models. The PPK parameters, including group typical values, fixed effect parameter, interindividual variation, and intraindividual variation, which are the determinant factors of drug concentration in patients, are taken into consideration. PPK is capable of quantitatively describing the effects of different factors in drug metabolism, such as pathology, physiology, and combined medication, then providing guidance on the adjustment of therapeutic regimen, thus increasing the efficacy and safety in a new drug evaluation.Many statistical models have been proposed to fit PPK parameters. The most popular analytic statistical model for PPK data is the random effects model proposed by Larid and Ware [1]. Besides, the classic compartmental model and nonlinear mixed effect model proposed by Sheiner et al. in 1977 [2] are the most commonly used statistical models for PPK. And recently, generalized linear model proposed by Salway and Wakefield in 2008 [3] has received wide attention.In fact nonlinear models are used to estimate parameters of a chosen compartmental model. This model can provide generally good results. However, a disadvantage of the nonlinear models is that it is generally more difficult to fit. From a computational point of view, one also faces the usual challenges associated with nonlinear regression such as choosing starting values, problem with convergence, nonlinear regression diagnostics, and so forth. Other nonlinear models for the analysis of PPK data see [4–6]. Salway and Wakefield [3] proposed a generalized linear model with gamma distribution to deal with PPK data. The GLM is one of the most widely used regression models for statistical analysis. The relevant theoretical investigations are rather complete and the results of parameter estimation have been the standard results in the literature and textbooks. Many methods could solve the generalized linear models, including general least square (GLS) method, quasi-likelihood method, normal theory maximum likelihood, quadratic estimating equations, extended quasi-likelihood, and general estimation equation (GEE). McCulloch [7] proposed three kinds of effective algorithm which can estimate the parameter in generalized linear model accurately. They are Monte Carlo EM (MCEM) algorithm, Monte Carlo Newton-Raphson (MCNR) algorithm, and Simulated Maximum likelihood (SML) separately. In the progress of maximum, Newton-Raphson algorithms, Fisher algorithms, EM algorithms and other basic algorithms can also be used to get the approximate solution through iteration. However, referring to multidimensional integral, it is difficult to solve a likelihood equation, and sometimes the algorithm is not convergent.In all these methods the covariance structure is an especially crucial problem. In general, we divided variance into four models following in normal distribution:(1) the constant error model, (2) the proportional error model, (3) a combined error model, (4) an alternative combined error model. Therefore, how to choose variance structure appropriately is really confusing. Wakefield et al. [8] suggested the use of integrated nested Laplace approximations for GLMs in this context. In this case, the empirical likelihood is an alternative good choice. Empirical likelihood as a nonparametric data-driven technique was first proposed by Owen [9].  It employs the likelihood function without specifically assuming a distribution for the data and incorporates side information through constraints or a prior distribution, which maximizes the efficiency of the method. What’s more, if the calculation procedure of the algorithm involves iteration, the problem of convergence may exist. Sometimes we cannot obtain convergent solution. But in principle, empirical likelihood can give estimates and confidence regions with almost any method that yields a “reasonable” set of estimating function. By “reasonable” we mean that the estimating functions have zero mean, finite variance, and either based on independent and identically distributed sampling or have higher order moments which allow use of a triangular array argument. Empirical likelihood only requires that the estimating functions have expectation zero. Hence empirical likelihood will be robust to mis-specification of the variance function, as long as the equation for the mean is correct.In the past two decades, unique and desirable properties of the empirical likelihood method have been studied by a number of researchers. These properties include, but are not limited to, range-respecting, transformation-preserving, asymmetric confidence intervals, Bartlett correctability, and better coverage probability or shorter confidence interval. Empirical likelihood has both the easy implementation of nonparametric inference method and the effectiveness of likelihood inference method. Compared with parametric reference method, it just needs a mild assumption of the existence of the matrix. When the sample size is small, the confidence interval we get by this method is superior to that by the asymptotic normal method. The shape and direction of the confidence interval confidence depend completely on the data itself. The estimation of variance is generally not required. Besides, empirical likelihood has unique superiority in parameter estimation and test of goodness of fit. An updated comprehensive overview of the empirical likelihood is available in Owen [10]. In a large range of situations, empirical likelihood has been shown to have properties analogous to a real likelihood, see Qin and Lawless [11, 12], Li [13], Jing [14], and others. A general framework of empirical likelihood based on estimation function is discussed in Qin and Lawless [12]. Wang et al. [15] incorporated the within-subject correlation structure of the repeated measurements into the empirical likelihood method. They also performed empirical likelihood inference on individual components of regression function.To the best of our knowledge, there is little literature on the empirical likelihood inference for PPK study with longitudinal data. In this paper, we develop an empirical likelihood method for inference of the parametric component in generalized linear models with population pharmacokinetic data through constructing auxiliary random vectors. The model is as follows:(1.1)E(yij(tij))=μij=Diexp(βi0+βi1tij+βi2tij),βij=βi+ηij,i=1…n,j=0,1,2,
where yij is the response variable. A few conditions needed to be assumed for the transferred additive error model to be valid to use. Assumptions are listed as follows:
(1.2)yij>0,βi0>0,βi1<0,βi2<0.
These assumptions are used to insure a ascending absorption phase and a descending elimination phase. The model is first proposed by Salway and Wakefield [3].The rest of the paper is organized as follows. In Section2, the empirical likelihood method for PPK models is proposed. The asymptotic normality of the proposed empirical estimator and asymptotic chi-square distribution of the proposed empirical likelihood ratio are derived under some regularity assumptions. In Section 3, the simulation is presented and the results also support the theory. A real data analysis is reported in Section 4. Finally, the paper concludes with some discussions in Section 5. Proofs are given in the appendix.


## 2. Methods

In recent years, numerous scholars have conducted a detailed study on PPK parameter estimation and a variety of methods appeared, for example GEE (Liang and Zeger [16]). Here, we use the method of quasi-likelihood to derive estimating functions for empirical likelihood. In the application part, we will compare the effectiveness of GEE and empirical likelihood method. Note that, since maximum likelihood and quasi-likelihood coincide for all linear exponential families, the class of estimating functions derived using the former method automatically falls within that of the latter. See Morris [17] and Nelder and Pregibon [18] Standard empirical likelihood for generalized linear models has been considered by Kolaczyk [19], based on constraints derived from the score function of the quasi-likelihood. We employ empirical likelihood and deduce as follows.Assume there aren subjects with m observations per subject in a study. From (1.1), the generalized error model can be rewritten as
(2.1)Yi=Ui+εi,
where Yi=[logyi1(ti1),logyi2(ti2),…,logyim(tim)] and Ui=[logμi1(ti1),logμi2(ti2),…,logμim(tim)] with var(εi)=∑.Therefore, change both side of (1.1) to logarithmic, and following Xue and Zhu [20], we can define an empirical log-likelihood ratio function as follows. Suppose that β is a 3*1 dimensional parameters vector to be estimated. From model (1.1), an auxiliary random vector can be defined as
(2.2)Zi(β)=XiTVi-1(β^)(Yi-Ui),
where Xi=(11×m,ti,1/ti), ti=(ti1,ti2,…,tim), and Vi(β) is an invertible covariance matrix dependent on the parameter vector β, and will be equal to cov(Yi). Vi(β) can be estimated by using the method of moments [16]. Note that E(Zi(β)) = 0 if β is the true parameters. Using this, an empirical log-likelihood-ratio function is defined as
(2.3)l~n(β)=-2sup{∑i=1nlog(npi)∣∑i=1n‍piZi(β)=0,∑i=1npi=1,pi≥0,i=1,…,n}.
By the Lagrange multiplier method, the optimal value for pi is given by
(2.4)pi=1n11+tTZi(β),
where t is Lagrange multiplierBy combining (2.3) and (2.4),   l~n(β) can be represented as
(2.5)l~n(β)=-2∑i=1nlog(1+tTZi(β)),
where t is determined by
(2.6)1n∑i=1n‍Zi(β)1+tTZi(β)=0.
Define β~ as the unique maximizer of l~n(β). We state the asymptotic properties of the empirical likelihood ratio. The following theorems indicate that β~ is asymptotically normal and the empirical log-likelihood ratio for β is a standard chi-square distribution. Some assumptions and proofs of our results are given in the appendix.Theorem 2.1.
In addition to LemmaA.1 mentioned in the appendix, one assumes that ∂2Zi(β)/∂β∂βT is continuous in a neighborhood of the true value β0. Then if ||∂2Zi(β)/∂β∂βT|| can bounded by some integrals function Z(x) in the neighborhood, then one has
(2.7)n(β~-β0)→N(0,V),
where V=[E(∂Z/∂β)T(EZZT)-1E(∂Z/∂β)]-1.Theorem 2.2.
Using the assumptions of Theorem2.1, if β is the true value of the parameter and the last component of β is a positive number, then one has
(2.8)ln(β)→χp2,
where → stands for convergence in distribution, and χp2 is a chi-squared distribution with p degrees of freedom.As a consequence of Theorem2.1, confidence regions for the parameter β can be constructed by (2.3). More precisely, for any 0≤α<1, let cα be the fraction quintile such that P(χp2≤α)=1-α. Then lα(β~)={β~:lα(β~)≤cα} constitutes a confidence region for β with asymptotically correct coverage probability 1-α.


## 3. Simulation

In this section some simulation results are reported to show the performance of the proposed empirical likelihood method. A two-stage generalized linear random effect model that incorporates a finite mixture model is considered in the simulation.At the first stage(3.1)Yi=Dexp(T*βi)+εi,
where βi=(βi0,βi1,βi2)′.At the second stage(3.2)βi=β+ηi,
where β=(β0,β1,β2)′ is the associated fixed effect.To perform the simulations, 500 datasets are generated, each consisting ofn=10, 50, 100 subjects. The time schedule {tij} is set as [1 2 4 8 10 12 16 32]/h. Dosage is set to 5 mg. The fixed effect β is set as [0.8-0.04-0.2]. The random effect is generated from a normal distribution as follows: εi~i.i.dN(0,Σ), and ηi~i.i.dN(0,Ω). For the variance matrix Σ and Ω, we choose a diagonal matrix with elements of 0.052 and 0.0052, respectively, for the reason that neither within-individual nor between-individual variation has been considered here.The conjugate gradient method used to search for the optimal empirical likelihood parameter estimates by modifying a routing is given in Press et al. [21]. For each dataset, we compute the MELE beta for fixed effects, and the 95% intervals. To illustration the efficacy of empirical likelihood, we also analysis each dataset by GEE method. The simulation results for different sample size are summarized in Table 1, based on 500 simulations.Table 1
Means of parameter estimates forβ with its standard derivation (SD) and mean length (ML) of the confidence region, when the estimation method are empirical likelihood (EL) and general estimation equation (GEE), respectively.
n
EL
GEE
Mean
CI (95% level)
Mean
CI (95% level)
10
0.8001
(0.7851, 0.8151)
0.8001
(0.7834, 0.8177)
−0.0398
(−0.0414, −0.0382)
−0.0398
(−0.0420, −0.0378)
−0.2003
(−0.2049, −0.1956)
−0.2003
(−0.2059, −0.1946)
50
0.8001
(0.7951, 0.8051)
0.8001
(0.7841, 0.8160)
−0.0401
(−0.0405, −0.0395)
−0.0401
(−0.0410, −0.0389)
−0.2000
(−0.2016, −0.1984)
−0.2000
(−0.2029, −0.1973)
100
0.8001
(0.7966, 0.8036)
0.8001
(0.7908, 0.8094)
−0.0400
(−0.0403, −0.0397)
−0.0400
(−0.0404, −0.0395)
−0.2000
(−0.2011, −0.1989)
−0.2000
(−0.2023, −0.1984)The results are listed in Table1. The “Mean” is constructed according to the 500 simulated point estimators, which was computed as the mean value of 500 simulated results of β~. The confidence intervals are listed in the columns “CI”. The lower and upper bounds of “CI” are the average of 500 simulated lower bounds, and upper bounds, respectively. From Table 1, it is fairly that, as the sample size increases, the confidence intervals become smaller, and the length of EL confidence interval shows advantage over GEE method. For the estimator β~, these show that the empirical likelihood method for the model can obtain accurate estimator and can construct the effective confidence intervals of β.


## 4. Application

We apply the empirical likelihood in the proposed generalized linear models to the theophylline data from bronchial asthma study. They were originally analyzed in Upton et al. [22] and are available from the Resource Facility for Population Kinetics. Initially this data set was used for Bayesian analysis of linear and nonlinear population models by using the Gibbs sampler. Details about the related design, methods, and medical implications of the asthma study have been described by Upton et al. [22]. Although all the individuals were scheduled to have clinical visits followed the time schedule, due to the various reasons, some individuals missed scheduled visits, which resulted in unequal number of measurements and different measurement times across individuals. It has been also used in Bayesian analysis of linear and nonlinear population models employing Gibbs sampler with normal errors. Wakefield et al. [8] suggested the generalized linear mixed model with the gamma error utilize for it. Based on this model, he had obtained a Bayesian result better than the compartment model result [8]. A Bayesian analysis of the dataset was completed by Salway and Wakefied [3] and Yan et al. [23].Here we applied our method to the data of 12 subjects who were administered an oral dose of the antiasthmatic agent theophylline, with 10 concentration measurements obtained from each individual over the subsequent 25 hours. Confidence regions for each parameters are shown in Figure1, respectively. The ranges of the confidence regions are small, which suggests a good fit for the parameters.The estimated value ofβ and its confidence region.
(a)
(b)A comparable representation with the GEE methods completed for the derived parameters can be seen in Table2. We compared the proposed empirical likelihood method with GEE method and the typically compartment model results. It concludes that these methods are comparable in terms of accuracy and efficiency. Table 2 indicates that EL-based confidence region is smaller than the GEE-based region and compartment results. Moreover, these findings are similar to the ones of Salway and Wakefield [3] and Yan et al. [23].Table 2
Comparing parameter estimates for theophylline data in terms of accuracy and efficiency. EL: empirical likelihood; GEE: general estimation equation; CP: compartment model.
EL
95% CI
GEE
95% CI
CP
95% CI
t
max
2.1676
(2.03, 2.30)
2.1676
(2.02, 2.38)
2.26
(1.73, 2.96)
c
max
8.3182
(7.53, 9.24)
8.3182
(7.52, 9.28)
9.06
(8.15, 10.1)
t
1
/
2
7.3862
(7.29, 7.97)
7.3862
(7.00, 7.96)
7.75
(7.06, 8.51)


## 5. Conclusions

This paper shows that empirical likelihood is an effective method for analyzing the PPK model. The empirical log-likelihood ratios are proven to be asymptotically chi-squared, and then the corresponding confidence regions for the parameters of interest are constructed. A simulation study is used to investigate the performance of the empirical likelihood method. The results show that the empirical likelihood could obtain accurate estimator and construct the effective confidence regions. We also analyze the theophylline data from the Upton et al. [22], and the results are similar and comparable to the ones of Wakefield et al. [8] and Yan et al. [23].At last, we have to mention the two disadvantages of this method. For one, the likelihood function in the empirical likelihood method is completely determined by the data, while avoiding the estimation of the variance, so the estimation of individual effect may not be obtained. But actually, it is an important parameter in PPK study and directly impacts doctor's recommendation of individual administration. The following two steps are advised to solve this problem. First, estimate the parameters in the model through empirical likelihood function and then calculate the random effect with the estimators.For another, in this paper, we did not consider the correlation structure of within-subject observation. However, in the current literature on empirical likelihood methods, the potential correlation structure of within-subject observation was first taken into consideration by Wang et al. [15]. They also demonstrated that it was easy to incorporate the correlation structure. It offered a promising direction to study this problem. Furthermore, we have only accomplished the complete data set at present and incomplete data will be further studied in future. From the application point of view, missing data and censoring data are worth being studied.



---

*Source: 250909-2012-12-27.xml*