Optimal designs for random effect models with correlated errors with applications in population pharmacokinetics
aa r X i v : . [ s t a t . A P ] N ov The Annals of Applied Statistics (cid:13)
Institute of Mathematical Statistics, 2010
OPTIMAL DESIGNS FOR RANDOM EFFECT MODELS WITHCORRELATED ERRORS WITH APPLICATIONSIN POPULATION PHARMACOKINETICS By Holger Dette , Andrey Pepelyshev and Tim Holland-Letz Ruhr University at Bochum, St. Petersburg State University andRuhr University at Bochum
We consider the problem of constructing optimal designs for pop-ulation pharmacokinetics which use random effect models. It is com-mon practice in the design of experiments in such studies to as-sume uncorrelated errors for each subject. In the present paper anew approach is introduced to determine efficient designs for nonlin-ear least squares estimation which addresses the problem of corre-lation between observations corresponding to the same subject. Weuse asymptotic arguments to derive optimal design densities, and thedesigns for finite sample sizes are constructed from the quantiles ofthe corresponding optimal distribution function. It is demonstratedthat compared to the optimal exact designs, whose determination isa hard numerical problem, these designs are very efficient. Alterna-tively, the designs derived from asymptotic theory could be used asstarting designs for the numerical computation of exact optimal de-signs. Several examples of linear and nonlinear models are presentedin order to illustrate the methodology. In particular, it is demon-strated that naively chosen equally spaced designs may lead to lessaccurate estimation.
1. Introduction.
The work presented in this paper is motivated by someproblems encountered in the design of experiments in a clinical trial to es-tablish the pharmacokinetics of Uzara R (cid:13) , a digitoxin related herbal diarrheamedication [based on Th¨urmann, Neff and Fleisch (2004)]. These kinds oftrials pose methodological design challenges because they require the esti- Received August 2009; revised December 2009. Supported by Deutsche Forschungsgemeinschaft (SFB 823, “Statistik nichtlinearerdynamischer Prozesse.” Teilprojekt C2). Supported in part by NIH Grant IR01GM072876:01A1 and the BMBF projectSKAVOE.
Key words and phrases.
Random effect models, nonlinear least squares estimate, cor-related observations, compartmental models, asymptotic optimal design density.
This is an electronic reprint of the original article published by theInstitute of Mathematical Statistics in
The Annals of Applied Statistics ,2010, Vol. 4, No. 3, 1430–1450. This reprint differs from the original in paginationand typographic detail. 1
H. DETTE, A. PEPELYSHEV AND T. HOLLAND-LETZ mation of global population parameters in the presence of correlated mea-surement errors. The trial in question included a number of patients eachgiven an oral application of Uzara R (cid:13) as well as, after a washout period, anintravenous application of digoxin (Lanicor R (cid:13) ), where in both cases the re-sulting serum concentration of digitoxin was measured repeatedly duringthe next 36 hours.The relation between the time and the concentration in the analysis of theUzara R (cid:13) trial can be described using the theory of one-compartment modelswith oral and, respectively, intravenous applications [Atkinson et al. (1993),Shargel (1993)]. In the intravenous case, the medication reaches the maxi-mum concentration in the blood almost immediately, and, after that, it isgradually eliminated from the body over time. Thus, the digitoxin concen-tration is modeled using the exponential elimination model η ( t, b ) = b e − b t . (1.1)In the case of the oral application, there is a gradual build up of the con-centration in the blood as the medication is absorbed through the digestivetract, while there is a simultaneous elimination process of the medicationin the blood. Therefore, the concentration function is the solution of thedifferential equation of these two parallel processes. The resulting function η ( t, b ) = b ( e − b t − e − b t )(1.2)is known as the (3 parameter) Bateman function [Garrett (1994)]. In bothmodels, the function η ( t, b ) denotes the concentration function, t is the time(in hours), b = ( b , b ) and, respectively, b = ( b , b , b ) are the vectors ofparameters. The parameters are assumed to vary between patients and theaim of the experiment is to estimate their global means (and sometimesvariances) over all patients.Measurements within the same patient are usually correlated, and weassume this correlation to be proportional to the time lag between measure-ments, which is plausible as the random errors are usually caused by tempo-rary changes in the patients physical condition. Measurements for differentpatients are assumed to be independent. In the trial at question Thuermannconsidered n = 15 (oral), respectively, n = 14 (intravenous) measurementseach on K = 18 patients. After a preliminary discussion with experts themeasurements were taken at nonoptimized time points. An approximationof the covariance of a single patient can then be expressed asΣ pop = ∂η ( t , b ) ∂b T V p ∂η ( t , b ) ∂b + V ε , (1.3)where η ( t , b ) = ( η ( t , b ) , . . . , η ( t n , b )) T denotes the vector of expected re-sponses at t , . . . , t n and V ε is the covariance matrix corresponding to this PTIMAL DESIGNS FOR RANDOM EFFECT MODELS data. This expression includes two sources of variation, the usual variation V ε caused by random errors as well as the additional variation V p due to therandom effect.Situations of this kind are rather common in the evaluation of the phar-macokinetics and the pharmacodynamics of drugs [see Buelga et al. (2005),Colombo et al. (2006), among others]. The corresponding processes are usu-ally modeled by linear or nonlinear random effects models, which try to esti-mate the population parameters, that is, the mean and the inter-individualvariability of the parameters. Under the additional assumption of a normaldistribution, the population characteristics are usually estimated by maxi-mum likelihood methods. In many cases, the likelihood cannot be evaluatedexplicitly and approximations are used to calculate the estimate. Efficientalgorithms for estimation are available for this purpose [see Aarons (1999)].Loosely speaking, under a Gaussian assumption this approach correspondsto weighted nonlinear least squares estimation. It was pointed out by sev-eral authors that the application of an appropriate design in these studiescan substantially increase the accuracy of estimation of the population pa-rameters. Usually, the construction of a good design is based on the Fisherinformation matrix which cannot be derived explicitly in pharmacokineticmodels with random effects. For this reason, many authors propose an ap-proximation of the likelihood [see, e.g., Mentr´e, Mallet and Baccar (1997),Retout, Mentr´e and Bruno (2002), Retout and Mentr´e (2003), Schmelter(2007a, 2007b), among others], which is used to derive an approximationfor the Fisher information matrix. This matrix is considered in various op-timality criteria, which have been proposed for the construction of optimaldesigns for random effect regression models.In the present paper the investigation is motivated by the following is-sues. First, the estimation of the population mean and the constructionof corresponding optimal designs for population pharmacokinetics stronglydepends on the Gaussian assumption, which is usually made for compu-tational convenience. Moreover, the maximum likelihood estimates may beinconsistent if the basic distributional assumption is violated. As a con-sequence, the derived optimal designs might be inefficient. Second, mostauthors derive the approximation for the Fisher information matrix underthe additional assumption that the random errors corresponding to the mea-surements of each individual are uncorrelated [see, e.g., Retout, Duffull andMentre (2001), Retout, Mentr´e and Bruno (2002) and Retout and Mentr´e(2003), among many others]. However, this assumption is not realistic inmany applications of population pharmacokinetics. Thus, a general conceptfor constructing optimal designs in the presence of correlated observationsis still missing. Third, even if the Gaussian assumption and the assumptionof uncorrelated errors for each subject can be justified, the numerical con-struction of the estimate and the corresponding optimal design is extremelyhard. H. DETTE, A. PEPELYSHEV AND T. HOLLAND-LETZ
In the present paper we address these issues. To derive a good design,we consider nonlinear least squares estimation in random effect regressionmodels. Note that this estimation does not require a specification of theunderlying distributions. For this estimation, we introduce a methodologywhich can be used to derive efficient or optimal designs in very generalsituations. More precisely, we embed the discrete optimal design problem ina continuous optimal design problem, where a nonlinear functional of thedesign density has to be minimized or maximized. This approach takes intoaccount the correlation dependence and yields an asymptotic optimal designdensity, which has to be determined numerically in all cases of practicalinterest. For finding the optimal density, we propose an algorithm based onpolynomial approximation. For a fixed sample size and for each individual,an exact design can be obtained from the quantiles of the correspondingoptimal distribution function. It is demonstrated by examples that thesedesigns are very efficient. Moreover, the designs, derived from the asymptoticoptimal design density, are very good starting designs for any procedure oflocal optimization for finding the exact optimal designs. To our knowledge,the proposed method is the first systematic approach to determine optimaldesigns for linear and nonlinear mixed effect models with correlated errors.The remaining part of this paper is organized as follows. In Section 2 weconsider the case of a linear random effect model and explain the basic designconcepts. In Section 3 we introduce the approach for deriving the asymptoticoptimal designs for linear regression models with correlated observations byemploying results of Bickel and Herzberg (1979). In Section 4 we presentresults for nonlinear regression models with correlated random errors andderive D -optimal designs and optimal designs for estimating the area underthe curve in the compartmental model. Finally, the Uzara R (cid:13) and Lanicor R (cid:13) trials are re-analyzed and optimal designs for the model (1.2) are determined.In particular, we show that the design proposed by the experts is ratherefficient, while a naively chosen equidistant design can yield substantiallyless accurate estimates.
2. Statement of the problem.
Consider the common random effect linearregression model Y ij = b Ti g ( t ij ) + ε ij , i = 1 , . . . , K ; j = 1 , . . . , n i , (2.1)where Y ij denotes the j th observation of the i th subject at the experimen-tal condition t ij , ε , . . . , ε K,n K are centered random variables with vari-ances depending on t , Var( ε ij ) = σ h ( t ij ) for some positive function h ( t ), g ( t ) = ( g ( t ) , . . . , g p ( t )) T is a given vector of linearly independent regressionfunctions, and b i is a p -dimensional random vector representing the individ-ual parameters of the i th subject, i = 1 , . . . , K . The explanatory variables t ij PTIMAL DESIGNS FOR RANDOM EFFECT MODELS can be chosen by the experimenter from a compact interval T . We assumethat errors ε i = ( ε i , . . . , ε in i ) for different subjects are independent, but theerrors for the same subject are correlated, that is,Cov( ε ij , ε is ) = σ h ( t ij ) h ( t is )( γr ( t ij − t is ) + (1 − γ ) δ j,s ) , (2.2)where γ ∈ [0 ,
1] is a constant, r ( t ) is a given correlation function such that r (0) = 1, and δ j,s denotes Kronecker’s symbol. Let V ε be the correspondingcovariance matrix. Assume that the individual parameters b i are drawn froma population distribution with mean β and covariance matrix V p , and theyare independent of the random variables ε i . This means that the covariancebetween two observations at the time t ij and the time t is ( j = s ) isCov( Y ij , Y is ) = g T ( t ij ) V p g ( t is ) + σ h ( t ij ) h ( t is ) γr ( t ij − t is ) , while the variance of Y ij is given by g T ( t ij ) V p g ( t ij ) + σ h ( t ij ). It was shownby Schmelter (2007b) that an optimal design necessarily advises the exper-imenter to perform observations of all subjects at the same experimentalsettings, that is, t ij = t j ( i = 1 , . . . , K , j = 1 , . . . , n ). Consequently, we definean exact design ξ = { t , . . . , t n } as an n -dimensional vector which describesthe experimental conditions for each subject. Without loss of generality, weassume that the design points are ordered, t < · · · < t n .Suppose that n observations are taken according to the design ξ . Thenthe model (2.1) for the i th subject can be written as Y i = X g b i + ε i ; i = 1 , . . . , K, (2.3)where Y i = ( Y i , . . . , Y in ) T , the matrix X g is given by X g = ( g ( t ) , . . . , g ( t n )) T ,and ε i is now a centered random variable with variance σ h ( t i ). This modelis a special case of the random-effect models discussed in Harville (1976),which are called generalized MANOVA. According to Fedorov and Hackl(1997), the (ordinary) least square estimate of β minimizes K X i =1 n X j =1 h − ( t j )( Y ij − f ( t j ) β ) . Define f ( t ) = g ( t ) /h ( t ) and X = ( f ( t ) , . . . , f ( t n )) T . Then the covariancematrix of the ordinary least squares estimate ˆ β OLS is given by D ( ˆ β OLS ) = 1 K ( X T X ) − X T ( V ε + XV p X T ) X ( X T X ) − (2.4) = 1 K (( X T X ) − X T V ε X ( X T X ) − + V p ) . Alternatively, if the covariance matrix V ε of the errors and the covariancematrix of the random effects V p were known (or can be well estimated), the H. DETTE, A. PEPELYSHEV AND T. HOLLAND-LETZ (weighted) least squares statisticˆ β WLS = 1 K K X i =1 ( X T ( V ε + XV p X T ) − X ) − X T ( V ε + XV p X T ) − Y i , (2.5)can be used to estimate the parameter β . The covariance matrix of theestimate ˆ β WLS is given by D ( ˆ β WLS ) = 1 K ( X T ( V ε + XV p X T ) − X ) − . (2.6)Since the expression (2.4) is simpler than (2.6) (which requires two dif-ferent inversions), the design methodology developed in this paper is basedon the covariance matrix of the ordinary least squares estimate. Addition-ally, we will demonstrate that the optimal designs obtained by minimizingfunctionals of the covariance matrix of the ordinary least squares estimateare also very efficient for weighted least squares estimation.We call a design optimal design if it minimizes an appropriate functionalof the covariance matrix of the least squares estimate. Since we consideroptimality criteria that are linear with respect to scalar multiplication ofthe covariance matrix, we put K = 1 without loss of generality.
3. Asymptotic optimal designs.
Although the theory of optimal designhas been discussed intensively for uncorrelated observations [see, e.g., Fe-dorov (1972), P´azman (1986) and Atkinson and Donev (1992)], less resultscan be found for dependent observations. For linear and nonlinear randomeffect models, optimal designs under the assumption of uncorrelated errorshave been investigated in Schmelter (2007a, 2007b), Mentr´e, Mallet andBaccar (1997) and Retout and Mentr´e (2003), among others. For fixed effectregression models with the presence of correlated errors, it was suggested toderive optimal designs by asymptotic considerations. Sacks and Ylvisaker(1966, 1968) have considered a fixed design space, where the number ofdesign points in this set tends to infinity. As a result, the asymptotic op-timal designs depend only on the behavior of the correlation function in aneighborhood of the point 0. In the present paper we use the approach ofBickel and Herzberg (1979) and Bickel, Herzberg and Schilling (1981), whohave considered a design interval expanding proportionally to the numberof observation points. This case is equivalent to the consideration of a fixedinterval with the correlation function depending on the sample size. To beprecise, we assume that the design space is given by an interval T and thedesign points of a sequence of designs ξ n = { t n , . . . , t nn } are generated by afunction a ( · ) in the form t jn = a (( j − / ( n − , j = 1 , . . . , n ;(3.1) PTIMAL DESIGNS FOR RANDOM EFFECT MODELS and a : [0 , → T denotes the inverse of a distribution function. Note that thefunction a ( · ) is obtained from the density of the weak limit of the sequence ξ n as n → ∞ . For example, if T = [ − , a ( u ) = (2 u −
1) cor-responds to the equally-spaced designs with distribution function a − ( x ) = x +12 and density ( a − ) ′ ( x ) = I [ − , ( x ). Furthermore, we assume that thecorrelation function r ( t ) of the errors ε i in (2.2) depends on n in the form r n ( t ) = ρ ( nt ) , (3.2)such that ρ ( t ) = o ( t ) as t → ∞ . For the numerical construction of asymptoticoptimal designs, we derive an asymptotic representation for the covariancematrix of the least squares estimate in the following lemma. For this purpose,we impose the following regularity assumptions:(C1) The functions f ( t ) , . . . , f p ( t ) are linearly independent, bounded on theinterval T and satisfy the first order Lipschitz condition: | f i ( t ) − f i ( s ) | ≤ M | t − s | and | f i ( t ) | ≤ M for all t, s ∈ T, i = 1 , . . . , p. (C2) The function a ( · ) is twice differentiable and there exists a positiveconstant M < ∞ such that for all u ∈ (0 , M ≤ a ′ ( u ) ≤ M, | a ′′ ( u ) | ≤ M. (3.3)(C3) The correlation function ρ is differentiable with bounded derivativeand satisfies ρ ′ ( t ) ≤ t .Assumption (C1) refers to the continuity of the response as a function of t and it is satisfied for most of the commonly used regression models. More-over, most of the commonly used correlation structures satisfy assumption(C3) on a compact interval. Finally, assumption (C2) refers to the charac-teristics of a design, requiring, loosely speaking, that observations shouldnot be clustered. This is quite reasonable due to ethical aspects. The follow-ing result is obtained by similar arguments as given in Bickel and Herzberg(1979) and, hence, its proof is omitted. Lemma.
Assuming that conditions (C1) , (C2) and (C3) are satisfied,then the covariance matrix of the least squares estimate have the form D ( ˆ β OLS ) = σ n ( W − ( a ) + 2 γW − ( a ) R ( a ) W − ( a )) + V p + o (cid:18) n (cid:19) , (3.4) where the matrices W ( a ) and R ( a ) are defined by W ( a ) = (cid:18)Z f i ( a ( u )) f j ( a ( u )) du (cid:19) pi,j =1 ,R ( a ) = (cid:18)Z f i ( a ( u )) f j ( a ( u )) Q ( a ′ ( u )) du (cid:19) pi,j =1 H. DETTE, A. PEPELYSHEV AND T. HOLLAND-LETZ and the function Q ( · ) is given by Q ( t ) = P ∞ j =1 ρ ( jt ) . Note that only the first term in (2.4) [and (3.4)] depends on a (asymptotic)design, and we propose to use this term for the construction of optimaldesigns. If the function a ( · ) is the inverse of a continuous distribution withdensity, say, ϕ , then a ′ ( t ) = 1 /ϕ ( t ) and, for large n , the first term of thecovariance matrix can be approximated by the matrix V ( ϕ ) /n , where the p × p matrix V is given by V ( ϕ ) := σ ( W − ( ϕ ) + 2 γW − ( ϕ ) R ( ϕ ) W − ( ϕ )) . (3.5)In (3.5) the matrices W ( ϕ ) and R ( ϕ ) have the form W ( ϕ ) = (cid:18)Z T f i ( t ) f j ( t ) ϕ ( t ) dt (cid:19) pi,j =1 ,R ( ϕ ) = (cid:18)Z T f i ( t ) f j ( t ) Q (1 /ϕ ( t )) ϕ ( t ) dt (cid:19) pi,j =1 . A density will be called an asymptotic optimal density ϕ ∗ , if it minimizes anappropriate functional of the matrix V ( ϕ ). Numerous criteria have been pro-posed in [Silvey (1980), Atkinson and Donev (1992), Pukelsheim (1993)] and,exemplarily, we consider the D - and c -optimality criteria which minimize − det V ( ϕ ) and c T V ( ϕ ) c for a given vector c ∈ R p , respectively. The applica-tion of the proposed methodology to other optimality criteria is straightfor-ward. The general procedure for constructing an efficient design minimizinga given functional of the covariance matrix of the least squares estimate isas follows:(1) Specify the correlation function ρ ( · ) in (3.2) and compute Q ( · ).(2) Compute the asymptotic optimal design density ϕ ∗ that minimizes anappropriate functional of the matrix V ( ϕ ) in (3.5).(3) Derive the exact design for a fixed sample size n by calculation of thequantiles of the distribution function Φ ∗ that corresponds to ϕ ∗ , namely, t i,n = (Φ ∗ ) − (cid:18) i − n − (cid:19) ; i = 1 , . . . , n. (3.6)The optimal density ϕ ∗ in step (2) can be determined as follows. We considerthe parametric representation of a density by a polynomial in the form ϕ ( t ) = ( p + p t + · · · + p r t r ) + R T ( p + p t + · · · + p r t r ) + dt and apply the Nelder–Mead algorithm to find the optimal density that min-imizes the specified functional of the matrix V ( ϕ ) with respect to p , . . . , p r . PTIMAL DESIGNS FOR RANDOM EFFECT MODELS One can run the algorithm for different degrees of the polynomial and differ-ent initial values and choose the density corresponding to the minimal valueof the optimality criterion. All integrals can be calculated by the Simpsonquadrature formula. We found that the minimal value of the criterion isnegligibly decreasing for polynomials of degree larger than r = 6. We alsoinvestigated the case where a density is represented in terms of rational, ex-ponential or spline functions. The results were very similar and, on the basisof our numerical experiments, we can conclude that the optimal density ϕ ∗ can be very well approximated by 6-degree polynomials.The derived designs from the asymptotic theory can be used to constructefficient designs for a given sample size as specified in step (3) of the algo-rithm. Alternatively, exact optimal designs can be determined by employingthe derived designs as initial values in a (discrete) optimization procedure.More precisely, for the determination of an exact optimal design, the aboveprocedure can be extended by the fourth step:(4) Determine an exact optimal design, that minimizes a functional of thecovariance matrix in (2.4), by using the Nelder–Mead algorithm with aninitial n -point design, which has been found in step (3).We illustrate this methodology for a quadratic regression model by the fol-lowing example. Also, in Section 4 we extend the approach for designingfor nonlinear models and investigate its performance for the compartmen-tal model. In particular, we demonstrate that the designs derived from theasymptotic consideration are very efficient compared to exact optimal de-signs. Example 1 (Quadratic regression). We illustrate the proposed approachfor constructing optimal designs for the classical quadratic regression modelwith homoscedastic errors. For this model, we have p = 3, f ( t ) = 1, f ( t ) = t , f ( t ) = t and T = [ − , ρ ( t ) = e − λt and r n ( t ) = e − λnt , λ > . (3.7)The asymptotic D -optimal densities for different choices of the parameters λ and γ are shown in Figure 1. Note that the numerically calculated optimaldesign densities are symmetrical, but we were not able to prove the symmetryof the asymptotic optimal density. We observe that the D -optimal densityconverges to the density of the uniform design, if γ → λ →
0. On theother hand, if γ → λ → ∞ , it can be seen that the asymptotic D -optimaldesign density is more concentrated at the points −
1, 0 and 1, which arethe points of the exact D -optimal design for the quadratic fixed effect modelwith uncorrelated observations [see Gaffke and Krafft (1982)]. Such behavioris natural because the errors are less correlated as γ → λ → ∞ . H. DETTE, A. PEPELYSHEV AND T. HOLLAND-LETZ
Fig. 1.
Asymptotic D -optimal design densities for least squares estimation in a randomeffect quadratic model for different choices of parameters in the covariance function (2.2) with r ( t ) defined by (3.7) . Right panel γ = 0 . ; left panel λ = 0 . . Further, we investigate the efficiency of an exact design derived from theasymptotic theory for least squares estimation, where the parameters in thecorrelation function (2.2) are given by γ = 0 . , σ = 0 . , V p = diag( σ β , σ β ) = diag(0 . , . ) . (3.8)Let ξ un be an n -point equidistant design and ξ an be an n -point design obtainedby the transformation (3.6) from the asymptotic optimal density. The pointsof the design ξ an and the exact optimal designs for least squares estimationfor λ = 1 . OLS ( ξ ) = (cid:18) det[( X T X ) − X T V ε X ( X T X ) − + V p ]det[( X T OLS X OLS ) − X T OLS V ε X OLS ( X T OLS X OLS ) − + V p ] (cid:19) /p , and for weighted least squares estimation,eff WLS ( ξ ) = (cid:18) det[ X T ( V ε + XV p X T ) − X ]det[ X T WLS ( V ε + X WLS V p X T WLS ) − X WLS ] (cid:19) /p for two designs: the design ξ an , derived from asymptotic theory, and theuniform design ξ un . Here X denotes the design matrix obtained for the design ξ under investigation, while X OLS and X WLS correspond to the optimalexact design for ordinary and weighted least squares estimation respectively.We observe that the design points concentrate in three regions containingthe points of the exact D -optimal design for a quadratic regression withuncorrelated errors. It is also noteworthy that the designs derived from theasymptotic theory are very efficient, in particular, for weighted least squaresestimation. PTIMAL DESIGNS FOR RANDOM EFFECT MODELS Fig. 2.
Left part: Various designs for ordinary and weighted least squares estimationin the random effect quadratic regression model. Exact D -optimal designs derived fromasymptotic theory: ball; exact D -optimal for ordinary least squares estimation: diamond;exact D -optimal designs for weighted least squares estimation: triangle. Right part: Effi-ciency of the designs ξ an and equidistant design ξ un for ordinary and weighted least squaresestimation. The parameters are given by (3.8) , where λ = 1 . .
4. Nonlinear random effect models.
In this section we extend the method-ology to the case of nonlinear random effect models, which have found con-siderable interest in the literature on pharmacokinetics. In this case, theresults of experiments are modeled by Y ij = η ( t j , b i ) + ε ij , i = 1 , . . . , K ; j = 1 , . . . , n. (4.1)Since the model (4.1) is nonlinear with respect to the variables b i , there isno analytical expression for the likelihood function and various approxima-tions have been considered in the literature [see Mentr´e, Mallet and Baccar(1997), Retout, Mentr´e and Bruno (2002), Retout and Mentr´e (2003), amongothers]. These approximations are used for the calculation of the maximumlikelihood estimate and the corresponding Fisher information matrix. Al-ternatively, an estimate of the population mean β can be obtained as anaverage of the nonlinear least squares estimates ˆ b i for the different individ-uals, but due to the nonlinearity of the model, an explicit representation ofthe corresponding covariance matrix cannot be derived. Following Retoutand Mentr´e (2003), we employ a first-order Taylor expansion to derive anapproximation of this covariance matrix. To be precise, we obtain (undersuitable assumptions of differentiability of the regression function) the ap-proximation η ( t, b ) ≈ η ( t, β ) + g ( t, β )( b − β ) T , (4.2)where g ( t, b ) = ∂η ( t, b ) ∂b H. DETTE, A. PEPELYSHEV AND T. HOLLAND-LETZ denotes the gradient of the regression function with respect to b . This meansthat the nonlinear model (4.1) is approximated by the linear model (4.2).For the construction of the optimal design, we assume that knowledge aboutthe parameter β is available from previous or similar experiments. This cor-responds to the concept of locally optimal designs, introduced by Chernoff(1953) in the context of fixed effect nonlinear regression models. Usually,locally optimal designs serve as benchmarks for commonly used designs andare the basis for the construction of optimal designs with respect to moresophisticated optimality criteria including the Bayesian and minimax ap-proach [see Chaloner and Verdinelli (1995) or Dette (1995)].Following the discussion in Section 3, we define the function f ( t, b ) = g ( t, b ) /h ( t ) to account for heteroscedasticity. Note that the covariance matrixof the nonlinear least squares estimate in the model (4.1) is approximatedby replacing the matrix X in model (2.1) with f ( t ) = f ( t, b ) | b = β , and themethodology described in Sections 2 and 3 can be applied to determineefficient designs. In the context of dose finding studies, it has been shownby means of a simulation study that the approximation (4.2) has sufficientaccuracy for the construction of optimal designs [see Section 5 in Dette et al.(2008) for more details].We further illustrate this concept by giving several examples for the case ofhomoscedastic errors. First, we investigate D - and c -optimal designs for therandom effect model, which has recently been studied by Atkinson (2008).Next, we re-analyze the Uzara R (cid:13) and Lanicor R (cid:13) trials introduced in Section 1.4.1. D -optimal design for a random effect compartmental model. Weconsider the random effect compartmental model with first-order absorp-tion, η ( t, b ) = b b − b ( e − b t − e − b t ) . (4.3)The model (4.3) is a special case of the Bateman function, defined in theintroduction [see Garrett (1994)], and has found considerable attention inchemical sciences, toxicology and pharmacokinetics [see, e.g., Gibaldi andPerrier (1982)]. The optimal design problem in the compartmental fixed ef-fect model has been studied by numerous authors [see, e.g., Box and Lucas(1959), Atkinson et al. (1993), Dette and O’Brien (1999), Biedermann, Detteand Pepelyshev (2004), among others], but much fewer results are availableunder the assumption of random effects. Recently optimal approximate de-signs for a random-effect compartmental model (4.3) have been determinedby Atkinson (2008), but we did not find results about exact designs for thesemodels in the presence of correlated errors. In the present paper we derivesuch designs from the asymptotic optimal design density and compare thesedesigns with the exact optimal designs. PTIMAL DESIGNS FOR RANDOM EFFECT MODELS Fig. 3.
Asymptotic D -optimal design densities for nonlinear least squares estimation inthe compartmental model (4.3) for different choices of the parameters in the covariancefunction (2.2) with r ( t ) defined by (3.7) . Left part: λ = 0 . ; right part: γ = 0 . . Note that the gradient of the function η ( t, b ) with respect to b is given by g ( t, b ) = (cid:18) b ( e − b x − e − b x ) + ( b x − b b x ) e − b x ( b − b ) , (4.4) b ( e − b x − e − b x ) + ( b x − b b x ) e − b x ( b − b ) (cid:19) T . In order to illustrate the methodology, we consider the same scenario as inAtkinson (2008) and assume that the parameters of the population and errordistributions are the following: β (0) = (1 , . T , γ = 0 . , σ = 0 . , V p = diag( σ β , σ β ) = diag(0 . , . ) , (4.5)and the design space is given by the interval T = [0 , r ( t ) in (2.2) is given by (3.7). The asymptotic D -optimal designdensities for different values of the parameters are shown in Figure 3. Weobserve that, for γ → λ →
0, the D -optimal design densities approximatethe uniform design, while, for larger values of λ or small values of γ , theasymptotic D -optimal designs put more weight at two specific regions of thedesign space. This fact corresponds to intuition, because the (approximate) D -optimal design for the model (4.3) with uncorrelated observations is atwo-point design [see, e.g., Box and Lucas (1959)].Further, we investigate the performance of the uniform and an exact de-sign derived from asymptotic theory. For this purpose we define ξ un as an n -point equidistant design { /n, /n, . . . , } and ξ an as the n -point designobtained by the transformation t j = (Φ ∗ ) − ( j/ ( n + 1)) , j = 1 , . . . , n, H. DETTE, A. PEPELYSHEV AND T. HOLLAND-LETZ where Φ ∗ denotes the distribution function corresponding to the asymptotic D -optimal design density. Note this transformation is slightly different fromthe transformation (3.1) in order to exclude the point 0 from the designpoints. Obviously, it is not reasonable to take observations t = 0 in model(4.3), because it is assumed that the drug is administered at time t = 0.The corresponding points of the exact designs are depicted in the left partof Figure 4, while the right part of the figure shows the D -efficiencies of thedifferent designs. We observe that the designs derived from the asymptotictheory have a substantially larger D -efficiency compared to the uniformdesign. For example, for n = 6, the D -efficiency of the uniform design isapproximately 50% for ordinary and weighted nonlinear least squares esti-mation, while the D -efficiency of the design ξ an is close to 90%.It is worthwhile to mention that, in nonlinear random effect models, theoptimal designs depend additionally on the mean β of the distribution of thepopulation parameters b i . Therefore, it is of interest to investigate the sensi-tivity of the designs with respect to a misspecification of this parameter. Fora study of the impact of such a misspecification on the efficiency of the result-ing designs, we consider the case n = 4 and n = 6 and the corresponding de-signs ξ a = { . , . , . , . } and ξ a = { . , . , . , . , . , . } ,respectively. In Figure 5 we display the efficiencies (cid:18) det[( X T X ) − X T V ε X ( X T X ) − + V p ]det[( X T OLS ,β X OLS ,β ) − X T OLS ,β V ε X OLS ,β ( X T OLS ,β X OLS ,β ) − + V p ] (cid:19) − /p (4.6)for different values of β . Here X denotes the design matrix obtained fromthe design ξ an under the assumption that β (0) = (1 , , T , while the matrix Fig. 4.
Left part: Various designs for ordinary and weighted nonlinear least squares es-timation in the compartmental model (4.3) . Exact D -optimal designs derived from asymp-totic theory: ball; exact D -optimal for ordinary least squares estimation: diamond; exact D -optimal designs for weighted least squares estimation: triangle. Right part: Efficiency ofthe designs ξ an and ξ un for ordinary and weighted least squares estimation. The parametersare given by (4.5) , where λ = 0 . . PTIMAL DESIGNS FOR RANDOM EFFECT MODELS Fig. 5. D -efficiencies of the designs ξ a = { . , . , . , . } (left part) and the design ξ a = { . , . , . , . , . , . } (right part) (these designs are obtained from asymp-totic density) for nonlinear least squares estimation in the random effect compartmentalmodel (4.3) , if the mean of the population distribution has been misspecified. The param-eters of the population distribution are given by (4.5) with λ = 0 . . X OLS ,β corresponds to the exact D -optimal design for nonlinear least squaresestimation for a specific β . The efficiencies are plotted in Figure 5 for therectangle [ β (0) i − σ β (0) i , β (0) i + 3 σ β (0) i ] = [0 . , . × [0 . , . D -efficiency with respect to misspecification of the population mean β over a broad range.4.2. Optimal designs for estimating the AUC.
In some bioavailabilitystudies, the aim of experiments is the estimation of the area under curveAUC = Z ∞ η ( t, β ) dt. For the compartmental model (4.3), we obtain AUC = 1 /b . It can be shownthat the locally AUC-optimal design for model (4.3) minimizes the varianceof the nonlinear least squares estimate for the parameter β . This varianceis approximately proportional to(0 , X T X ) − X T V ε X ( X T X ) − + V p )(0 , T . This expression corresponds to the c -optimality criterion, which has beendiscussed extensively in the literature for fixed effect models with uncorre-lated observations [see, e.g., Ford, Torsney and Wu (1992), Fan and Chaloner(2003) and Dette et al. (2008), among others]. The asymptotic optimal de-sign densities for estimating the area under the curve are shown in Figure6. We observe again that the optimal density is close to the uniform designdensity if λ → γ →
1. On the other hand, if λ is large or γ →
0, the H. DETTE, A. PEPELYSHEV AND T. HOLLAND-LETZ
Fig. 6.
Asymptotic optimal densities for estimating the area under the curve in the com-partmental model (4.3) for different choices of the parameters in the covariance function (2.2) with r ( t ) defined by (3.7) . Left part: λ = 0 . , right part: γ = 0 . . AUC-optimal design density has a narrow support. This fact reflects thatthe optimal design for estimating the area under the curve in the fixed effectcompartmental model with uncorrelated observations is a one-point design.In Figure 7 we show the designs derived from the asymptotic optimal designdensities, and the exact optimal designs for estimating the area under thecurve in the compartmental model. We observe that the designs ξ an , derivedfrom the asymptotic optimal design density, are very close to the exactoptimal designs for least squares estimation of the area under the curve.Moreover, the design ξ an yields a substantial improvement in efficiency com-pared to the uniform design. Again we observe that the designs derived forordinary least squares estimation also have excellent efficiencies for weightedleast squares estimation of the AUC (see the right part of Figure 7).4.3. Optimal designs for estimating the AUC in the Uzara R (cid:13) and Lanicor R (cid:13) trials. In this section we consider the optimal design problem for estimat-ing the AUC in the two examples presented in the introduction. The med-ical background of these examples was a small pilot trial of 4 patients[Th¨urmann, Neff and Fleisch (2004)], where it was observed that patientstaking the over-the-counter herbal diarrhea medication Uzara R (cid:13) (in the formof drops, i.e., oral application) showed high values in medical assays designedto measure the blood serum concentration of digitoxin, a potent treatmentagainst heart insufficiency. This is caused by the chemical similarity of thesetwo substances and can result in major complications in establishing treat-ment programs for heart insufficiency. It was thus decided to compare thepharmacokinetic properties of an oral application of Uzara R (cid:13) to the proper-ties of the usual intravenous application of a regular digitoxin medication PTIMAL DESIGNS FOR RANDOM EFFECT MODELS Fig. 7.
Left part: Various designs for estimating the area under the curve in the compart-mental model (4.3) . Designs derived from asymptotic theory: ball; exact optimal designsfor least squares estimation of the area under the curve: diamond; exact optimal designsfor weighted least squares estimation: triangle. Right part: Efficiency of the designs ξ an and ξ un for ordinary and weighted least squares estimation of the AUC. The parameters aregiven by (4.5) , where λ = 1 . . (Lanicor R (cid:13) ) on a larger sample size of 18 patients, with 15 measurementseach. The main focus of the comparison was the area under the concentra-tion curve as a measure of the total effect of an application. A preliminarydesign was proposed by experts in order to allow precise estimation of thisproperty, and the study was carried out according to this design. We willnow investigate the efficiency of this and a naively chosen equally spaceddesign compared to a design generated using the methods presented in thispaper.In order to compute the design densities, it is necessary to estimate theparameters for both of the models. To do so, we used the full 18 patient dataset, however, this could have been done as well using only the 4 patients ofthe pilot study. We estimated parameters using a combination of maximumlikelihood and least squares techniques [see Pinheiro and Bates (1995)]. Forthe intravenous part corresponding to model (1.1), we have received theestimates ˆ β = (30 , . T , ˆ V p = (cid:18) . .
189 0 . (cid:19) for the population parameters, while the estimates of the parameters in thecovariance matrix (1.3) are given by ˆ γ = 0 .
8, ˆ λ = 0 .
05, ˆ σ = 0 . H. DETTE, A. PEPELYSHEV AND T. HOLLAND-LETZ
For the oral application corresponding to model (1.2), the estimates ofthe parameters are given by ˆ β = (0 . , . , T ,ˆ V p = . . . . and ˆ γ = 0 .
8, ˆ λ = 0 .
01, ˆ σ = 0 .
2. Note that the entries in the positions (1 , , ,
1) and (3 ,
2) are not identical 0, but smaller than 10 − . Basedon a preliminary discussion with experts, the physicians decided to takemeasurements for both experiments at nearly identical (nonoptimized) timepoints 0 , . , , . , , , , , , , , , , , . , . , . , , , , , , , , , , , , respectively. Note that in the second design an additional measurement wastaken at t = − .
25, that is, before the intravenous injection. Since this pointis out of the scope of the exponential evasion model, the observation at thistime has been excluded from further considerations. Thus, the Lanicor R (cid:13) trial has n = 14 observations.For the estimated parameters, we have derived the asymptotic optimaldesign densities, which are depicted in the left and right part of Figure 8 forthe Uzara R (cid:13) and Lanicor R (cid:13) trials, respectively, where the design interval isgiven by T = [0 , . , . , . , . , . , . , . , . , . , . , . , . , . , . , Fig. 8.
Asymptotic optimal density for estimating the area under the curve in the -parameter Bateman (left part) and exponential evasion model (right part). PTIMAL DESIGNS FOR RANDOM EFFECT MODELS for the Uzara R (cid:13) trial and0 . , . , . , . , . , . , . , . , . , . , . , . , . , R (cid:13) trial. This means that the asymptotic optimal design forthe oral application is close to an equally spaced design, while the optimal de-sign for the intravenous application is much more focused on measurementsat early time points. This result is plausible in an exponential eliminationmodel. Calculating the efficiencies for ordinary least squares estimation ofthese asymptotic optimal designs compared to the exact optimal design, weobtain efficiencies 0 .
99 and 0 .
95, respectively, indicating that these designsare quite good.Upon numerical studies, we can further conclude that the constructeddesigns are robust to misspecification of the initial guess of parameters. Forexample, variations of λ (or any of the other parameters) by up to 50% yielda drop of efficiency by less than 0 . R (cid:13) ) and 97% (Lanicor R (cid:13) ), even betterthan the efficiencies for ordinary least squares estimation.We now investigate the performance of the designs which were actuallyused in the clinical trial. We found that these designs have efficiency 0 . .
92 for estimating the AUC through OLS estimation in the Uzara R (cid:13) andLanicor R (cid:13) trials, respectively. Thus, the preliminary designs, recommendedby experts, are rather efficient in both trials.Let us now investigate the performance of naively chosen equidistant de-signs: the design0 , . , . , . , . , . , . , . , . , . , . , . , . , . , R (cid:13) trial, and the design0 , . , . , . , . , . , . , . , . , . , . , . , . , R (cid:13) trial. Comparing these designs to the optimal designs, weobtain efficiencies of 0 .
97 and 0 .
41 for ordinary least squares estimation ofthe AUC in the Uzara R (cid:13) and Lanicor R (cid:13) trials, respectively.It was pointed out by a referee that it is of some interest to investigate theperformance of the optimal designs proposed in this paper for the estimationof the variance of the random effects (i.e., the parameters of the matrix V p ).These parameters are usually estimated by maximum likelihood techniques[see Retout and Mentr´e (2003), among others] and the corresponding infor-mation matrix of these estimates is of a fundamentally different structure H. DETTE, A. PEPELYSHEV AND T. HOLLAND-LETZ compared to the variance matrix of the least squares estimate. For the twooptimal designs we have calculated the D -efficiencies for estimating the di-agonal elements of the matrix V p , which are 96% and 79% in the Uzara R (cid:13) and Lanicor R (cid:13) trial, respectively. The designs actually used in the trial haveefficiency 98% and 85%, while the corresponding efficiencies of the uniformdesigns are 98% and 63%, respectively. Thus, the proposed optimal designsfor AUC-estimation also have reasonable efficiency for estimation of the co-variance matrix of the population distribution.Summarizing the discussion in this example, we conclude that the equallyspaced design in the Uzara R (cid:13) trial is very close to the optimal design deter-mined by the proposed methodology, and it is for this reason very efficient.However, the equally-spaced designs do not always have high efficiency. Inthe Lanicor R (cid:13) trial, the use of naively chosen designs yields considerably lessaccurate estimates. For this reason, the application of experimental designtechniques in the context of pharmacokinetics trials is strictly recommended. Acknowledgments.
The authors would like to thank Professor Petra Thuer-mann, who provided the data sets, and Martina Stein, who typed parts ofthe paper with considerable technical expertise. We would also like to thankthe referees and the Associate Editor for constructive comments on an earlierversion of this paper. REFERENCES
Aarons, L. (1999). Software for population pharmacokinetics and pharmacodynamics.
Clinical Pharmacokinetics Atkinson, A. C. (2008). Examples of the use of an equivalence theorem in construct-ing optimum experimental designs for random-effects nonlinear regression models.
J. Statist. Plann. Inference
Atkinson, A. C. and
Donev, A. (1992).
Optimum Experimental Designs.
ClarendonPress, Oxford.
Atkinson, A. C., Chaloner, K., Herzberg, A. M. and
Juritz, J. (1993). Optimumexperimental designs for properties of a compartmental model.
Biometrics Bickel, P. J. and
Herzberg, A. M. (1979). Robustness of design against autocorrelationin time I: Asymptotic theory, optimality for location and linear regression.
Ann. Statist. Bickel, P. J., Herzberg, A. M. and
Schilling, M. F. (1981). Robustness of designagainst autocorrelation in time II: Optimality, theoretical and numerical results for thefirst-order autoregressive process.
J. Amer. Statist. Assoc. Biedermann, S., Dette, H. and
Pepelyshev, A. (2004). Maximin optimal designs for acompartmental model. In
MODA 7—Advances in Model-Oriented Design and Analysis
Box, G. E. P. and
Lucas, H. L. (1959). Design of experiments in non-linear situations.
Biometrika Buelga, D. S., del Mar Fernandez de Gatta, M., Herrera, E. V., Dominguez-Gil, A. and
Garcia, M. J. (2005). The Bateman function revisited: A critical reeval-uation of the quantitative expressions to characterize concentrations in the one com-PTIMAL DESIGNS FOR RANDOM EFFECT MODELS partment body model as a function of time with first-order invasion and first-orderelimination. Antimicrobial Agents and Chemotherapy Chaloner, K. and
Verdinelli, I. (1995). Bayesian experimental design: A review.
Statist. Sci. Chernoff, H. (1953). Locally optimal designs for estimating parameters.
Ann. Math.Statist. Colombo, S., Buclin, T., Cavassini, M., Decosterd, L., Telenti, A., Biollaz, J. and
Csajka, C. (2006). Population pharmacokinetics of atazanavir in patients withhuman immunodeficiency virus infection.
Antimicrobial Agents and Chemotherapy Dette, H. (1995). Designing of experiments with respect to “standardized” optimalitycriteria.
J. Roy. Statist. Soc. Ser. B Dette, H. and
O’Brien, T. (1999). Optimality criteria for regression models based onpredicted variance.
Biometrika Dette, H., Bretz, F., Pepelyshev, A. and
Pinheiro, J. C. (2008). Optimal designsfor dose finding studies.
J. Amer. Statist. Assoc.
Fan, S. K. and
Chaloner, K. (2003). A geometric method for singular c -optimal designs. J. Statist. Plann. Inference
Fedorov, V. V. (1972).
Theory of Optimal Experiments . Academic Press, New York.MR0403103
Fedorov, V. V. and
Hackl, P. (1997).
Model-Oriented Design of Experiments. LectureNotes in Statist. . Springer, New York. MR1454123
Ford, I., Torsney, B. and
Wu, C. F. J. (1992). The use of canonical form in theconstruction of locally optimum designs for nonlinear problems.
J. Roy. Statist. Soc.Ser. B Gaffke, N. and
Krafft, O. (1982). Exact D -optimum designs for quadratic regression. J. Roy. Statist. Soc. Ser. B Garrett, E. R. (1994). The Bateman function revisited: A critical reevaluation of thequantitative expressions to characterize concentrations in the one compartment bodymodel as a function of time with first-order invasion and first-order elimination.
Journalof Pharmacokinetics and Biopharmaceutics Gibaldi, M. and
Perrier, D. (1982).
Parmacokinetics , 2nd ed. Dekker, New York.
Harville, D. (1976). Extension of the Gauss–Markov theorem to include the estimationof random effects.
Ann. Statist. Mentr´e, F., Mallet, A. and
Baccar, D. (1997). Optimal design in random-effectsregression models.
Biometrika P´azman, A. (1986).
Foundations of Optimum Experimental Design . D. Reidel, Dordrecht.MR0838958
Pinheiro, J. C. and
Bates, D. M. (1995). Approximations to the log-likelihood functionin the nonlinear mixed effects model.
J. Comput. Graph. Statist. Pukelsheim, F. (1993).
Optimal Design of Experiments . Wiley, New York. MR1211416
Retout, S. and
Mentr´e, F. (2003). Further developments of the Fisher information ma-trix in nonlinear mixed-effects models with evaluation in population pharmacokinetics.
Journal of Biopharmaceutical Statistics Retout, S., Duffull, S. and
Mentre, F. (2001). Development and implementation ofthe Fisher information matrix for the evaluation of population pharmakokinetic designs.
Computer Methods and Programs in Biomedicine Retout, S., Mentr´e, F. and
Bruno, R. (2002). Fisher information matrix for nonlin-ear mixed-effects models: Evaluation and application for optimal design of enoxaparinpopulation pharmacokinetics.
Stat. Med. H. DETTE, A. PEPELYSHEV AND T. HOLLAND-LETZ
Sacks, J. and
Ylvisaker, N. D. (1966). Designs for regression problems with correlatederrors.
Ann. Math. Statist. Sacks, J. and
Ylvisaker, N. D. (1968). Designs for regression problems with correlatederrors; many parameters.
Ann. Math. Statist. Schmelter, T. (2007a). Considerations on group-wise identical designs for linear mixedmodels.
J. Statist. Plann. Inference
Schmelter, T. (2007b). The optimality of single-group designs for certain mixed models.
Metrika Shargel, L. (1993).
Applied Biopharmaceutics and Pharmacokinetics . Chapman & Hall,London.
Silvey, S. D. (1980).
Optimal Design . Chapman & Hall, London. MR0606742
Th¨urmann, P., Neff, A. and
Fleisch, J. (2004). Interference of Uzara glycosides inassays of digitalis glycosides.
International Journal of Clinical Pharmacology and Ther-apeutics H. DetteFakult¨at f¨ur MathematikRuhr-Universit¨at Bochum44780 BochumGermanyE-mail: [email protected]
A. PepelyshevDepartment of MathematicsSt. Petersburg State UniversitySt. PetersburgRussiaE-mail: [email protected]