Optimal Population Designs for Discrimination Between Two Nested Nonlinear Mixed Effects Models

Diseños poblacionales óptimos para discriminación entre dos modelos no lineales de efectos mixtos anidados

M. E. Castañeda L.^a,*
V. I. López-Ríos^b

^a Universidad de Antioquia, Medellín, Colombia.
^* Autor de correspondencia: maria.castaneda@udea.edu.co.
^b Profesor asociado, Escuela de Estadística, Universidad Nacional de Colombia, Medellín, Colombia.

Recepción: 21-sep-2015 Aceptación: 20-ene-2016

Abstract

Key words: Optimal Designs, Mixed Effects Model, T-Optimal Designs.

Resumen

En este artículo se considera el problema de encontrar diseños poblacionales óptimos para dicriminar entre dos modelos no lineales de efectos mixtos anidados, los cuales difieren en su matriz de covarianza intraindividual. El criterio propuesto es una generalización del criterio de T-optimalidad; para él se proporciona el respectivo teorema de equivalencia, y su aplicación se ilustra por medio de un ejemplo en farmacocinética.

Palabras clave: diseños T-óptimos, diseños óptimos, modelo de efectos mixtos.

1 Introduction

In the application of optimal design theory one of the basic assumptions is to assume that the model used to describe a given phenomenon or process is the correct model. However, in the practice there may exist several candidate models. One way of to select the most adequate model among several candidates is conducting an experiment designed so that the observations obtained allow us to discriminate between the models in the best way possible. This leads to the problem of to find the optimal experimental conditions using some optimality criterion for discriminate between competing models.

In the case of fixed effects models a commonly used criterion for discriminating between two competing homoscedastic models is the T-optimality criterion, which was proposed by [1]. Under normally distributed observations, the T-optimal design provides the most powerful F-test for the lack of fit of one model when the other is assumed to be true. The criterion has been generalized for other fixed effects models see [18, 19].

The nonlinear mixed effects models are particularly useful in longitudinal studies such as population pharmacokinetics experiments, assay analysis and studies of growth, in which a limited number of samples can be obtained from each individual. These models distinguish two classes of variation: the random variation among observations within a given individual (intra-individual) and random variation among individuals (inter-individual) [3, 4]. This separation of variability allows the estimation of population characteristics from sparse samples per individual in a set of subjects without requiring individual estimation of the parameters. Depending of the application and nature of data, different covariance structures may be considered to model the two class of variation. For example, in some situations it is common practice not to assume a particular structure for the inter-individual variation, whereas for intra-individual variation can be considered different structures, among them, the usual structure which assume independent observations with constant variance, compound symmetry and autoregressive structures, constant coefficient variation structure, variance function which depends on the conditional mean response, or combinations of these structures [7, 21, 12]. Under a nonlinear mixed effects model, a population design is defined by the number of individuals to study and the individual designs to be performed in the individuals (number of samples and the sampling times) [11]. Thus, assuming that the response function and the inter-individual covariance model have been correctly specified, it may be of interest the problem of designing an experiment in a group of individuals with sparse samples per individual for discrimination between two competing intra-individual covariance models which may be nested.

Although some approaches such as those described below have been proposed for discriminate between two nonlinear mixed effects models, these may be inappropriate or less efficient in situations involving nested models as the previously considered. Waterhouse et al. in [22] proposed the product D-optimality criterion based on the product of the determinants of the Fisher information matrices for to find designs useful for both parameter estimation and model discrimination. For nonlinear models, such designs may be less efficient for discriminate that the T-optimal designs. Vajjah and Duffull in [20] proposed a robust T-optimal design method which does not depend on a priori selection of the true model. However, in the particular case of two nested intra-individual variation models this methodology can be not applied because the T-optimal designs are based only on the fixed effects models without residual error. Kuczewski et al. in [9] proposed an extension of the T-optimality criterion for heteroscedastic models, for discriminate between two multiresponse models. The criterion is derived in the case of non-nested models and can be applied directly when all individuals are observed under the same experimental conditions. Therefore, alternative methods for discrimination in the case of nested models are required. In this work, we consider the problem of to find optimal population designs for discrimination between two nested nonlinear mixedeffects models which differ in their intra-individual covariance matrix. We propose a generalization of the T-optimality criterion for this case. Our approach can be applied to population studies for groups of individuals with different sampling scheme where each sampling scheme is a multidimensional point in a finite space of admissible sampling sequences.

2 Theoretical Background

2.1 Nonlinear Mixed Effects Model

We assume that for each individual i in a population of N individuals the number of different observations available is n. Let y_i = (y_i1,..., y_in)^T be the vector of repeated measurements for the ith individual and x_i = (x_i1,..., x_in)^T the n × 1 vector of sampling times where x_ij belongs to a finite set consisting of t different measurement times. It is assumed that the measurements made on different subjects are independent. To model the relationship between y_i and x_i we consider a nonlinear mixed effects model which may be written as hierarchical two-stage model, see [3]:

Stage 1. Intra-Individual Model

In this stage the variability among observations within a given individual (intra-individual) is modeled.

Suposse that

where β_i is a (p × 1) vector of parameters for the ith individual; f (x_i, β_i) is an n × 1 vector function, f (x_i, β_i) = (f (x_i1, β_i),..., f (x_in, β_i))^T where f is a known nonlinear function of β_i; and ε_i is the n × 1 random errors vector.

It is assumed that ε_i|β_i ∼ N_n(0, Cov (ε_i|β_i)) with

where R(x_i, β_i, σ² ,λ) is an (n × n) matrix and is called the intra-individual covariance matrix, which depends on the parameters σ² > 0 and λ ∈ Ω ⊂ R^d .

For a given individual, the matrix R takes into account the nature of intra-individual variation and may be choosen in such a way that reflects the heterogeneity of variance, and the correlation among observations, or both. For example, in the case of data from pharmacokinetic experiments and growth studies a common model is

and λ = θ. This matrix corresponds to uncorrelated errors with variance proportional to a power of the conditional mean. If the repeated observations are taken over time, a model for serial correlation can be considered, for example, the autoregressive (AR) model of order one for equally spaced data in time. Thus a model with this correlation structure and constant variance is and λ = α. Also, we can consider the case where the errors have nonconstant variance with correlation structure by the specification

For others structures of intra-individual correlation see [12].

Stage 2. Inter-Individual Model

In model (1), the variation among individuals (inter-individual) is modeled through the individual specific parameters β_i. In order to account the possible dependence of this variation on individual characteristics, a model for β_i is provided in this stage.

Suposse that

where β is a (p ×1)vector of population parameters and b_i is a (p × 1) vector of random effects associated with individual i. It is assumed that the b_i are independent and normally distributed with mean 0 and variance-covariance matrix D and that the b_i and ε_i are independent. The matrix D is called the inter-individual covariance matrix. The parameter σ² and the distinct elements of the covariance matrix D can be arranged in a single vector ψ of covariance parameters.

First-order approximation model

An approximation to the marginal distribution of y_i can be derived taking a first-order Taylor series expansion of the model (1) about E(b_i)=0. This expansion yields to the linealized model given by

Thus, under (4) the approximate marginal distribution of y_i is normal with approximate mean vector and variance-covariance matrix given by

2.2 Population Designs

For the given model, suppose that n_k independent observations are taken at point , where s is the number of distinct x_k. For example, consider a study in which n_k individuals are observed under the conditions vector x_k =(x_k1,...,x_kn), with the total number of individuals N. Then the collection of x_k and n_k, represented by

Since finding an optimal exact design is a discrete optimization problem which may be difficult from both analytical and computational points of view, the corresponding approximate design should be considered one in which the weights ω_k may be any real numbers from the interval [0,1]. Thus, the collection

is called approximate population design. The weight ω_k represents the proportion of total individuals that should be observed at the point x_k.

If r denotes the number of elements in the design region then the design ζ can be specified by the vector of weights

Under this representation, if ω_k =0 this means that the corresponding design point is not used in the experiment. The set of points x_k in the design region for which the design ω has nonzero weights ω_k is called the support set of ω and is denoted by supp(ω).

After optimization, the number of individuals in each group is obtained from the optimal weights by using n_k =N ×ω_k. This can yield noninteger number and therefore a rounding procedure is applied [14].

In what follows, we adopt the approximate locally optimal design approach and we use the approximate design ω =(ω₁,...,ω_r).

2.3 The Problem of Model Discrimination

In the case of fixed effects models, these are models that do not contain the level of random effects, one most commonly used criterion for model discrimination is T-optimality proposed by Atkinson and Fedorov in [1]. For two competing homoscedastic models this criterion is based on the assumption of one model f_t (x) = f₁(x,β₁)is the true model. The T-optimal design is a design that maximizes

where such that 0 ≤ ω_k ≤ 1 and , and f₂(x,β₂) is the rival model.

If the functions f_j(x,β _j) depend linearly on the parameter β _j, j = 1,2, then the quantity Δ(ξ) is proportional to noncentrality parameter of the χ² distribution of the residual sum of squares for the second model. The T-optimal design provides the most powerful F test for the lack of fit of the second model when the first is true. For nonlinear models this result is asymptotic.

For nonlinear mixed effects models, assuming that the response function f and the inter-individual model are correctly specified, we consider the discrimination problem between two nested intraindividual variation models. Specifically, let R₁ and R₂ be two alternative models to describe the intra-individual variability. It is assumed that R₂ is nested within R₁ in sense that both models involve the same structure R(x_i,β_i,σ² ,λ) and with respect to the parameter λ the parameter space Ω₂ of R₂ is a subset of the parameter space Ω₁ of R₁ defined by the imposition of κ equality constraints. That is, Ω₂ = {λ ∈ Ω₁ | h_τ (λ) = 0,τ = 1,...,κ} where the functions h_τ (λ) are assumed to be continuously differentiable.

The objective is to find the appropriate form of the intra-individual covariance matrix R. This can be achieved by performing an experiment designed in such a way that observations obtained allow us to discriminate between R₁ and R₂ in the best way possible. For to determine such experimental design, we propose an optimality criterion which corresponds to a generalization of the T-optimality criterion. This design provides the most powerful likelihood ratio test when the largest model is assumed to be the true model. In the next section the criterion is derived.

3 Criterion for Discrimination Between Two Intra-Individual Models

The discrimination between two nested intraindividual variation models leads to the discrimination between two nested nonlinear mixed effects models M₁ and M₂ such that the second stage is as in (3) for both models and the first stage of each model represents a different assumption about intraindividual model, specifically:

M₁: Model 1

Assuming that the approximation (4) is exact, these models can be represented as:

In order to discriminate between these models, assuming that the largest model is completely known, we propose to find the approximate design ω^*that maximizes the following generalization of Toptimality criterion over the set Ξ:

where for some known value . The design ω^* be called Twoptimal, where the letter W refers to the withinindividual variance-covariance matrix. This design is locally optimal because it depends on the values of γ₁.

For this class of nonlinear mixed effects models, this criterion may be considered as an extension of proposed criterion by [9] for groups with different designs and a single response.

The justification of criterion is follows.

Let γ⁰ denote the true value but unknown of parameter γ. To discriminate between the alternative models M₁ and M₂ we consider the likelihood ratio test for the model selection. Since M₂ is nested within M1, the testing problem can be formulated as:

where Γ2 = {γ =(β^T ,ψ^T ,λ^T )^T | γ ∈ Γ1, hτ (λ)= 0,τ = 1,...,κ}.

The following assumptions will be required:

Let y_k1,...,y_knk be the independent observations vectors of individuals with sampling times vector x_k, k = 1,...,s. Assuming that the approximation (4) is exact, y_k1,...,y_knk , are N_n( f (x_k,β),Σ(x_k,γ)) random vectors (k = 1,...,s).

For the testing problem (10), the likelihood function based on the s independent samples is

where are the maximum likelihood estimators over Γ₁ and Γ₂ respectively, the likelihood ratio test of H₀ against H₁, reject H₀ for large values of

Under H₀, the test statistic has asymptotically a central χ² distribution with κ degrees of freedom. Therefore, an approximate test of size α of H₀ is to reject H₀ if , where c_κ (α) denotes the upper 100α% point of the distribution.

In analysis of mean and covariance structure models, the function is known as the maximum likelihood discrepancy function which measures the discrepancy between the sample moments and the moments based in the model which depends on the parameter γ (see [16]). The minimizing of this function leads to the maximum likelihood estimator for ith group. Extensions of discrepancy function to more than one group is straightforward [2]. Specifically, the discrepancy function for s groups is defined as

Now, the power of this test, , is a function of the alternative parameter value γ. Given a specific value of N denoted by N₀ and a specific alternative parameter value close to Γ₂, this probability can be approximated by considering the asymptotic distribution of under a sequence of local alternatives converging to a point in Γ₂ (see [17]). It is assumed that is an interior point of Γ₁. The parameter value is identified then with . Since F(γ) in (12) is a discrepancy function, under assumptions (A1)-(A3) and regularity conditions, this asymptotic distribution is the noncentral Chi-square distribution with κ degrees of freedom and noncentrality parameter δ, which can be approximated by the value:

Since the power of test is a monotonically increasing function of the noncentrality parameter, from (14) the power is an increasing function of T_W (ζ_N) and hence can be maximized by the choice of design ζ_N .

Finally, the exact design ζ_N can be replaced by the corresponding approximate design ζ which is represented by ω. Thus we obtain the TW-criterion defined in (9).

A Necessary and Sufficient Condition for TW-Optimality

The following definition is fundamental for characterizing of Tw-optimal designs.

Definition 1. A design ω is called a regular design if the following set

is singleton, otherwise it is called singular design.

Hence, if ω is a regular design and , then is the unique solution of the equation

Theorem 1. Let ω^* be a regular design. Under the assumptions (A1)-(A3):

Proof. The proof of this theorem is similar to the proof of Theorem 1 in [9].

(i) First, we prove that the criterion Tw is a concave function. To this end, suposse ω₁, ω₂ ∈ Ξ and α ∈ [0, 1]. It is clear that Ξ is a convex set. Let ω = (1 − α)ω₁ + αω₂, then

Now, the directional derivative of Tw at ω in the direction of is any design, is given by

Since f (x, β) and Σ (x, γ) are continuous and twice continuously differentiable in Γ₁,it follows that is a continuous function at α and in Γ₂. Additionally, exists and is also continuous at α and in Γ₂. Hence, applying the Theorem 3.3 of [13], we get

Note that if Γ₂(ω)= , then

and using the definition of directional derivative

where

As ω^* is a regular Tw-optimal by assumption, we have and from (15) it fol

where ω is any design.

Since Tw(ω)is a concave function of ω, then the plasma concentration as a function of time. The the nonpositivity of the directional derivative simplest compartmental model assumed for such a at ω^*is a necessary and sufficient condition relationship is the nonlinear model given by, for the optimality of ω^*. From this fact and by (16), it follows that a necessary and sufficient condition for the optimality of ω^* is that ω^* fulfills the inequality

(ii) We assume the contrary, this mean there is a set and a scalar a such that ≤ a < 0 and φ(x ^* ,ω^*)= 0 for x ^* ∈ supp (ω^* )\{x₁,...,x_s1 }. Then

where s is the number of elements in supp(ω^*). From (15) taking ,wehave

This contradiction proves the assertion.

4 An Example

In this section we present an example to illustrate the use of the criterion proposed. This is a theorical pharmacokinetics example described by [8] in a simulation study and used by [10] in the application of methods to find optimal population designs to estimate population characteristics of the pharmacokinetics of a drug in sparse-sampling experiments. We reproduce the models and parameters values from the second study.

The pharmacokinetic studies seek to understand the process of drug absortion, distribution and elimination using for example kinetic models to describe the plasma concentration as a function of time. The simplest compartmental model assumed for such a relationship is the nonlinear model given by,

The nonlinear mixed effects model can be written as

Stage 1. (Intra-Individual Variation)

where, for the subject i, y_il represents the lth concentration measurement taken at time x_il, Cl_i is the clearance, V_i is the volume of distribution and β_i = (Cl_i, V_i)^T. The dose D = 1 is fixed for all individuals.

Stage 2. (Inter-Individual Variation)

where β =(Cl, V )^T is the mean values vector.

The two alternative models for the withinindividual covariance matrix R (xi,β_i,σ² ,λ)are:

M₁. Variance proportional to a power of mean response and uncorrelated errors

where

M2. Constant variance and uncorrelated errors

It is assumed that the variance σ² is a fixed constant equal to 0.15 for both models.

Since the model M₂ is nested within M₁, the true model is M₁ with the population parameter vector given by

where Cl₁ = 0.5, V₁ = 0.2, ψ_Cl1 = 0.01, ψV₁ = 0.0016 and θ =1.

For the alternative model M₂ the population parameter vector is

The set of possible sampling times considered is ={0.05,0.15,0.3, 0.6,1} hours after administration. We assume that three observations are available for each patient, n = 3, without replicates at an identical time point, which in a pharmacokinetic study is a sparse sampling situation. Therefore, the design region is given by the set of combinations of three sampling times from , that is, containing 10 elements. The sequences are:

x₁ =(0.05,0.15,0.3), x2 =(0.05,0.15,0.6),

x₃ =(0.05,0.15,1), x4 =(0.05, 0.3,0.6),

x₅ =(0.05,0.3, 1), x6 =(0.05, 0.6,1),

x₇ =(0.15,0.3, 0.6), x8 =(0.15,0.3, 1),

x₉ =(0.15,0.6, 1), x10 =(0.3,0.6, 1).

To find the locally TW-optimal design we used the nominal values of parameters defined previously like the local parameters and the design region . The optimal design was calculated optimizing the Tw-criterion implemented through an algorithm in R[15]. The function nlminb was used for the optimization in the design region .

The local TW-optimal design obtained is a twopoint design

Thus, if the objective is to design a new experiment with sparse sampling for discriminate between the M₁ and M₂ models then for about a 82% of the patients, blood samples should be taken at 0.05, 0.15 and 0.3 hours, and for about 18% blood samples should be taken at 0.05, 0.15 and 0.6 hours.

To check that the design obtained is Tw-optimal we use the equivalence theorem. First, we enumerate all candidate sampling sequences, that is, the elements of and calculate the sensitivity function for each sampling sequence. Then plot the sensitivity function as a function of index i. The resulting plot is shown in Figure 1. From this plot it is clear that the design ζ_* consisting of the x₁ and x₂ sequences is Tw-optimal.

5 Conclusions

A generalization of T-optimality criterion has been proposed for discriminate between two nested nonlinear mixed effects models. The first stage of each model represents a different assumption about intraindividual random variation and the second stage is the same for both models. Assuming that the response function is common for both models and it is correctly specified, we observe that the criterion development in this paper may be considered an extension of the proposed criterion by [9] for groups with different designs and a single response.

In the case of nested models an alternative criterion for discriminate between models is the DS-criterion which is appropiate when interest is in estimating a subset of s parameters. Since the rival models considered in this paper are nested this criterion may be applied. The comparison between the performances of Tw-and DS-optimal designs will be studied in future papers.

Another future work involve the study of designs with multiple objectives as the compound design. For example, the compound criteria for parameter estimation and for discrimination between models using D-optimality with Tw-optimality and Doptimality with DS-optimality.

References

[1] A. C. Atkinson, and V. V. Fedorov, "The design of experiments for discriminating between two rival models", Biometrika, vol. 62, pp. 57-70, 1975.         [ Links ]

[2] K. Bollen, "Estructural Equations with Latents Variables", John Wiley & Sons, New York, Inc, 1989.         [ Links ]

[3] M. Davidian, and D. M. Giltinan, "Nonlinear Models for Repeated Measurement Data", Chapman & Hall, London, 1995.         [ Links ]

[4] E. Demidenko, "Mixed Models: Theory and Applications", John Wiley & Sons Inc, New York, 2004.         [ Links ]

[5] V. V. Fedorov, and P. Hackl, "Model-Oriented Design of Experiments", Springer, New York, 1997.         [ Links ]

[6] R. Gagnon, and S. Leonov, "Optimal population designs for PK models with serial sampling", Journal of Biopharmaceutical Statistics, vol. 15, pp. 143-163, 2005.         [ Links ]

[7] D. Hand, and M. Crowder, "Practical Longitudinal Data Analysis", Chapman & Hall, London, 1996.         [ Links ]

[8] Y. Hashimoto, and L. Sheiner, "Designs for population pharmacodynamics: value of pharmacokinetic data and population analysis", Journal of Pharmacokinetics and Biopharmaceutics, vol. 19, pp. 333-353, 1991.         [ Links ]

[9] B. Kuczewski, B. Bogacka, and D. Uci´nski, "Optimum designs for discrimination between two nonlinear multivariate dynamic mixedeffects models", Biometrical Letters, vol. 45, pp. 1-28, 2008.         [ Links ]

[10] F. Mentré, P. Burtin, Y. Merlé, J. Bree, A. Mallet, and J. Steimer, "Sparse-sampling optimal designs in pharmacokinetics and toxicokinetics", Drug Information Journal, vol. 29, pp. 997-1019, 1995.         [ Links ]

[11] F. Mentré, A. Mallet and D. Baccar, "Optimal design in random-effects regression models", Biometrika, vol. 84, pp. 429-442, 1997.         [ Links ]

[12] J. C. Pinheiro and D. M. Bates, "Mixed-Effects Models in S and S-PLUS", Springer, New York, 2000.         [ Links ]

[13] B. N. Pshenichny, "Necessary Conditions for an Extremum", Marcel Dekker, New York, 1971.         [ Links ]

[14] F. Pukelsheim, Optimal Design of Experiments, John Wiley & Sons Inc, New York, 1993.         [ Links ]

[15] R Development Core Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, 2014, ISBN 3-900051-07-0, [Online], Available: http://www.R-project.org.         [ Links ]

[16] J. H. Steiger, A. Shapiro, and W. Browne, "On the multivariate asymptotic distribution of sequential Chi-Square statistics", Psychometrika, vol. 50, pp. 253-264, 1985.         [ Links ]

[17] T. W. Stroud, "Fixed alternatives and Wald's formulation of the noncentral asymptotic behavior of the likelihood ratio statistic", Annals of Mathematical Statistics, vol. 43, pp. 447-454, 1972.         [ Links ]

[18] D. Ucinski, and B. Bogacka, "T-optimum designs for multiresponse dynamic heteroscedastic models", MODa 7-Advances in Modeloriented Design and Analysis, eds A. Di Bucchianico, H. Laüter and H. P. Wynn, pp. 191-199, Heidelberg: Physica, 2004.         [ Links ]

[19] D. Ucinski, and B. Bogacka, "T-optimum designs for discrimination between two multiresponse dynamic models", Journal of the Royal Statistical Society: Series B, vol. 67, pp.3-18, 2005.         [ Links ]

[20] P. Vajjah, and S. Duffull, "A generalisation of T-optimality for discriminating between competing models with an application to pharmacokinetic studies", Pharmaceutical Statistics, vol. 11, pp. 503-510, 2012.         [ Links ]

[21] E. F. Vonesh, and V. M. Chinchilli, "Linear and Nonlinear Models for The Analysis of Repeated Measures", Marcel Dekker, New York, 1997.         [ Links ]

[22] T. H. Waterhouse, S. Redmann, S. B. Duffull and J. A. Eccleston, "Optimal design for model discrimination and parameter estimation for itraconazole population pharmacokinetics in cystic fibrosis patients", Journal of Pharmacokinetics and Pharmacodynamics, vol. 32, pp. 521-545, 2005.         [ Links ]