D-optimal experimental designs for a mixture of Adair models

Chemometrics and Intelligent Laboratory Systems 154 (2016) 150–161

Contents lists available at ScienceDirect

Chemometrics and Intelligent Laboratory Systems journal homepage: www.elsevier.com/locate/chemolab

Review

D-optimal experimental designs for a mixture of Adair models J. López-Fidalgo1 , M.M. Rodríguez-Hernández* University of Castilla-La Mancha, Department of Mathematics, Institute of Mathematics applied to Science and Engineering, Avda. Camilo José Cela 3, 13071 Ciudad Real, Spain

A R T I C L E

I N F O

Article history: Received 4 March 2015 Received in revised form 9 March 2016 Accepted 12 March 2016 Available online 28 March 2016 Keywords: Adair model Covariance matrix D-optimality EM algorithm Information matrix Mixture of distributions

A B S T R A C T In a Pharmacokinetic reaction a ligand may be bound to different types of macro-molecules, each with different number of binding sites. They are frequently involved in certain diseases diagnostics. Adair equation is used very often to model the reactions of biological macro-molecules with a ligand. This equation relates the saturation ratio with the free ligand concentration when the equilibrium is reached and it depends on some association constants of the chemical reaction, which have to be estimated. The main problem considered in this paper is the computation of optimal experimental designs for a mixture of Adair models when different types of macromolecules are mixed in the experiment. The main contribution of this work is obtaining the Fisher Information Matrix for a model with a mixture of probability distributions. Since this is not anymore in the Exponential family the expectation cannot be obtained analytically. Then the computation of optimal designs through the information matrix cannot be done with traditional methods. In data analysis when one has the data this expectation can be computed empirically from the data. But in experimental design data are not available when the experiment is being scheduled. Assuming nominal values of the parameters, as is usually done for nonlinear models, simulations were performed for each point in a suitable discretized design space. The number of simulations and the sample size used in each simulation were empirically tuned for both. A sensitivity analysis was performed for different possible true values of the parameters. Since this meant an important computational burden fractional designs were used to cover a reasonable neighborhood of the nominal values of the parameters. In order to display the results in a friendly way for the practitioner, “safe” neighborhoods of the optimal designs are provided. © 2016 Elsevier B.V. All rights reserved.

Contents Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Designs and estimates of a mixture of Adair models 2.1. Initial and free lingand relationship for a mixture of Adair models . . . . . . . . . . . . 2.2. Optimal designs for a mixture of Adair models 3. Models with a mixture of distributions. The EM algorithm . . . . . . . . . 3.1. Mixture of Adair models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. Illustrative example 4.1. Comparison of the Information Matrix and the covariance matrix . . . 4.2. Sensitivity of the design with respect to the true values of the parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5. Discussion . . . . . . . . . . . . . . . . . . . . . . . . Conflict of interests Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1. 2.

. . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

* Corresponding author. Tel.: +34 926 295300. E-mail addresses: jesus.lopezfi[email protected] (J. López-Fidalgo), [email protected] (M. Rodríguez-Hernández). 1

Tel.: +34 926 295212; fax: +34 926 295361.

http://dx.doi.org/10.1016/j.chemolab.2016.03.021 0169-7439/© 2016 Elsevier B.V. All rights reserved.

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

151 152 152 153 153 155 155 157 159 160 160 160 161

J. López-Fidalgo, M. Rodríguez-Hernández / Chemometrics and Intelligent Laboratory Systems 154 (2016) 150–161

1. Introduction Adair’s equation is used to model reactions of biological macromolecules with a ligand. This equation, which depends on the association constants of the chemical reaction, relates the saturation rate and the concentration of free ligand when the chemical equilibrium is reached. Monomers and dimers are chemical compounds of much interest in the fields of biochemistry and medicine, where they are involved in the diagnosis of certain diseases. Frequently one ligand binds to different types of macromolecules. A well known example of a mixture of two macro-molecules with different binding sites is the binding of oxygen to myoglobin (monomer) and hemoglobin (tetramer). This paper provides optimal experimental designs for a mixture of two Adair models depending on the free ligand. Since the free ligand is a random variable that is not under the control of the experimenter it is necessary to find a relationship between the free ligand and another variable that will be controlled experimentally. This variable is the initial ligand set at the beginning of the experiment. Then, the initial ligand is chosen in such a way the final free ligand will be optimal from the point of view of the experimental design theory. The Adair model is used at the moment the chemical reaction of the macro-molecules and ligands reaches the equilibrium. These reactions are usually reversible. Adair [1] proposed a model of how oxygen molecules bind to the binding sites of hemoglobin. For more details see [12]. In this paper we consider an experiment consisting of a container with an initial ligand concentration, [L0 ], in it. Then a semipermeable dialysis bag with concentrations of two types of macro-molecules (0) (1) (proteins) with N binding sites, [MN ] and [MN ] is introduced into the container. They could be two types of macro-molecules, two types of proteins or a mixture of a macro-molecule and a protein. Hereafter we will refer to macro-molecules for simplicity. Fig. 1 (left) shows the initial situation of the two macro-molecules just before unions occur. The initial ligand molecules begin to enter the dialysis bag, but the two types of macro-molecules cannot leave the bag. The equilibrium is reached when there is the same unbound ligand concentration inside and outside the bag, called free ligand, [L]. Fig. 1 (right) displays the situation, in a naïf way, for both types of macro-molecules when the equilibrium is reached. From top to bottom we see both types of bound macro-molecules, free ligand and the two types of unbound macro-molecules. In general, if the macro-molecules have N binding sites the ligand can bind them through 0, 1, 2, . . . , N sites. The concentration of the two types, s = 0, 1, of macro-molecules bound to the (s)i ligand will be denoted by [MN ]i = 1, 2, . . . , N while the unbound (s)0 macro-molecules will be denoted by [MN ]. The models for s = 0, 1 can be formulated as:

y(s) =

(s) (s) (s) K [L] + 2K1 K2 [L]2 +  1 (s) (s) (s) N 1 + K1 [L] + K1 K2 [L]2

... +

(s) NK1 (s) + K1

+ ...

...

(s) KN [L]N  (s) KN [L]N

...

151

An optimal design will provide both, good estimates and savings in the use of experimental resources [5]. The theory of optimal experimental design provides tools for computing good designs from different points of view. A model relating uncorrelated observations (responses), y1 , . . . , yn with some explanatory variables with the corresponding values [L]1 , . . . , [L]n that are under the control of the experimenter is stated by giving a probability density function (pdf),

f (yi , [L]i ; h), i = 1, . . . , n,

where h is the vector of parameters to be estimated. The collection of values [L]1 , . . . , [L]n from a design space w are usually called exact experimental design of size n. Some of these values may be repeated, e.g. there are k different values, say [L]1 , . . . , [L]k without loss of generality. This suggests the definition of a finite discrete probability distribution assigning to each point a weight proportional to the number of replicates. Using this idea Kiefer [8] introduced the concept of approximate design, n, as any probability measure defined on a design space w, where the design points may be chosen. The Fisher Information Matrix (FIM) associated to an approximate design n is, 

 M(n; h) =

I([L]; h)n(d[L]), I([L]; h) = E w

−∂ 2 L(h) ∂ h2

 ,

(2)

where L(h) is the log-likelihood function. For non-linear models this matrix depends on the parameters. Caratheodory’s theorem allows restricting to finite designs,  n=

[L]1 ... [L]k p1 ... pk

 ,

where n([L]i ) = pi define the mass probability function. A realizable design has to be “approximated” as ni ≈ n × pi replicates for the i-th experiment. D-optimality is one of the most popular optimality criteria. A D-optimal design maximizes the determinant of the Information Matrix, det[M(n; h)]. The general equivalence theorem establishes that a design, n*, is D-optimal if   tr M−1 (n∗ , h)I([L]; h) ≤ m,

[L] ∈ w,

(3)

+ e. (1)

For more details on the procedure see Section 2 of [9]. In order to estimate the parameters of the model different values of the initial ligand concentration are tried in different experiments. When the equilibrium is reached, the ligand that is out of the bag is measured, [L]. Inside the bag, the saturation of each of the two types of macro-molecules, yi , is also measured from a sample. The aim of this paper is to compute optimal experimental designs for the saturation ratio for a mixture of two types of macro-molecules in the following cases: • one binding site (two monomers), • one and two binding sites (one monomer and one dimer), • two binding sites (two dimers).

where m is the number of parameters. The left term is usually referred as the sensitivity function. Moreover, the equality is reached at the design points. For more details on optimal design see Section 2 of [9] and the references provided there. The main issue in this paper is that a mixture of distributions is not in the Exponential family. Thus, it is not possible to find an analytic expression for the FIM since the expectation in Eq. (2) cannot be computed analytically. It needs to be approximated in some way. In this paper it has been approximated through simulations at different points of the design space, once some nominal values of the parameters are established. As is well known (see e.g [3], [11] or [10]), computing Maximimum Likelihood Estimates (MLE) directly for models with mixture of distributions becomes rather tedious for many families of distributions. This problem may be solved by using the EM algorithm for missing data.

152

J. López-Fidalgo, M. Rodríguez-Hernández / Chemometrics and Intelligent Laboratory Systems 154 (2016) 150–161

Fig. 1. The initial phase of the experiment (left) and the experiment at equilibrium (right) for the mixture of two types of macro-molecules with just one binding site.

A general introduction of optimal experimental design has been given in this section. The estimation of the parameters with the EM algorithm will be considered in Section 3. Optimal designs will be computed in Section 2.2 through simulation. Estimation of parameters will be made in Section 3.1. The specific case of the mixture between a ligand and two macro-molecules with one and two binding sites will be illustrated with a example in Section 4. Efficiencies of designs around the optimal design were computed and the neighborhood with efficiencies over 95% was shown with some graphics. In Section 4.1 the parameters of the mixture will be estimated using the EM algorithm for a number of simulations. It will be shown that the empirical covariance matrix obtained from the estimates is very similar to the inverse of the information matrix for the optimal design divided by the appropriate sample size. Since the information matrix depends on the nominal values chosen for the parameters with the advice of the experimenter, locally optimal designs will be considered (see e.g [4]). A robustness analysis will be performed with respect to the choice of the nominal values of the parameters in Section 4.2. Due to the huge number of possible combinations of possible true values of the parameters a fractional factorial design will be used to reduce the number of simulations. Finally, the results and the conclusions are discussed in Section 5. Some ideas are also provided for further research.

2. Designs and estimates of a mixture of Adair models 2.1. Initial and free lingand relationship for a mixture of Adair models If an analytic relationship between the free ligand [L] and the initial ligand [L0 ] introduced in the experiment is known it can be used to transform the optimal design on [L] into a design on [L0 ], which can be controlled experimentally. A ligand that can bind to two types (0) of macro-molecules with one and two binding sites, say [M1 ] or (1) [M2 ] is considered here. The reversible reaction that takes place for the monomer is

inside and outside the bag. Thus, the so called association constant for the monomer is (0)1

K (0) =

[M1

]

(0)0 [M1 ][L]

,

and for the dimer, (1)1

(1)

K1 =

[M2

(1)2

]

(0)0

(0)1

+ L  M1

,

and the reaction carried out by the dimer is (1)0

M2

(1)1

+ L  M2

(1)2

+ L  M2

(1)0 [M2 ][L]

]

(1)1 [M2 ][L]

.

(5)

(0)

The initial macro-molecule concentration of the monomer, [M1 ], is the sum of the concentrations of the free macro-molecule and the bound macro-molecule, (0)

(0)0

[M1 ] = [M1

(0)1

] + [M1

].

(6) (1)

The initial concentration of protein with two binding sites, [M2 ], (1)0 is the sum of concentrations of the free protein, [M2 ], the protein (1)1 with one site bound, [M2 ], and the protein with two sites bound, (1)2 [M2 ], (1)

(1)0

[M2 ] = [M2

(1)1

] + [M2

(1)2

] + [M2

].

(7)

Let [L0 ] be the initial ligand concentration introduced into the container at the beginning of the experiment. Part of this ligand is bound to the macro-molecules and part remains free. Once the chemical equilibrium is reached the same free ligand concentration is inside and outside the bag, [L], that is, (0)1

(1)1

] + [M2

(1)2

] + 2[M2

].

(8)

Solving the equation system of Eqs. (4), (5), (6) and (7), plugging it into (8) the solution is obtained as a function of the free ligand concentration, the initial macro-molecule concentrations and the association constants,

.

When the equilibrium is reached, part of the initial ligand is (0)1 (1)1 bound to some parts of the macro-molecule, [M1 ], [M2 ] and (1)2 [M2 ]. The rest of the concentration of the macro-molecule is not (0)0 (1)0 bound, [M1 ], and [M2 ]. Then, [L] is the free ligand concentration

[M2

(1)

, K2 =

[L0 ] = 2[L] + [M1 M1

(4)

 [L0 ] = [L]

(0)

K (0) [M1 ] 1+

K (0) [L]

(1)

+

2 + K1

  ⎞  (1) (1) (1) (1) 2K2 [L]2 + [M2 ] + 2[L] 1 + K2 [M2 ] ⎠. (9)   (1) (1) 1 + K1 [L] 1 + K2 [L]

J. López-Fidalgo, M. Rodríguez-Hernández / Chemometrics and Intelligent Laboratory Systems 154 (2016) 150–161

The cases of mixture of two monomers and mixture of two dimers are developed in a similar way. In particular, the initial ligand concentration for the mixture of two Adair monomers models is (0)

[L0 ] = 2[L] +

(1)

K (0) [M1 ][L] 1+

+

K (0) [L]

K (1) [M1 ][L] 1 + K (1) [L]

,

(10)

⎛

(0) K1



(0) (0) 2K2 [L][M2 ]

2.2. Optimal designs for a mixture of Adair models The general model for any number of binding sites was derived in [9], Section 2.1. Assuming normality in the observations, the probability density function of the mixture of the two macromolecules is

fy (y, [L]; h) = (1 − p)

+p



Ly (h) = logfy (y, [L]i ; h)



(0) [M2 ]

+   (0) (0) 1 + K1 [L] 1 + K2 [L]   ⎞  (1) (1) (1) (1) (1) 2 + K1 2K2 [L]2 +[M2 ]+2[L] 1 + K2 [M2 ] ⎠. (11)   + (1) (1) 1+K1 [L] 1 + K2 [L]

[L0 ] = [L] ⎝

differentiating twice the log-likelihood function with respect to h and computing the expectation of with the distribution given by fy . The log-likelihood function is,

= log{(1 − p)

and for the mixture of two Adair dimer models,

1 2ps02

1 2ps12

 exp − 

exp −

1  2s02

y − g0 ([L]; h

(0)

2



)

 2 1  (1) , ([L]; h ) y − g 1 2s12

(12)

where g0 ([L]; h(0) ) and g1 ([L]; h(1) ) are the Adair models for the mean for the first and the second population respectively. The vector of parameters to be estimated is h = (p, h(0) , h(1) , s02 , s12 )T . It is a nonlinear model, and [L] is again the concentration of the free ligand once equilibrium in the experiment is reached. As mentioned above the models considered are the monomer or the dimer. The monomer has only one binding site while the dimer has two binding sites for the ligand. This binding yields a substantial change in the protein, which is transmitted to the rest of the binding sites. The FIM for the mixture model, not taking into account the uncertainty of the free ligand concentration, is computed by

153

+p

1 2ps12

1 2ps02 

exp −

 exp −

2 1  y − g0 ([L]; h(0) ) 2 2s0

 2 1  (1) }. ([L]; h ) y − g 1 2s12



(13)

For this model the expectation cannot be computed analytically. It will be computed through simulations at each value of [L]. It will be detailed in the Section 4. The mixture of two Adair models depends on the free ligand, which is not completely under the control of the experimenter. Just the initial ligand can be designed. Thus, it would be enough to plug the free ligand in an equation giving the initial ligand. This equation has already been calculated for some marginal cases in [9]. In particular, for macromolecules with two binding sites the formula is quite large and once it is plugged in a model with mixtures, it becomes too large to be used in practice. This is the reason why the optimal designs for the mixture will be calculated as a function of free ligand in the original Adair model. Then the locally D-optimal design will be n∗ = arg max |M(n; h)|. n

The design points, [L]∗ , cannot be controlled experimentally, since the free ligand is a random variable. Plugging [L]∗ into Eqs. (9), (10) or (11) an initial concentration of the ligand can be obtained for the tree cases of study. The design obtained in this way will not be D-optimal, but it is expected to have high D-efficiency, as it happened for the simplest cases where there was no mixture [9]. It will be called quasi D-optimal. In the case of mixture of two Adair monomer models it will be seen that the quasi D-optimal design has a single point. In the mixture of one monomer and one dimer the quasi D-optimal design will have two support points. Finally, for a mixture of two dimers the quasi D-optimal design will have three support points.

3. Models with a mixture of distributions. The EM algorithm Let y be the random variable that takes different values of observations coming from two mixed populations and let z be a Bernoulli random variable with parameter p. The value z = 0 means the observation comes from the first population and z = 1 means it comes from the second one. The conditional probability distribution of y given z = 0 (or 1) will have a pdf f0 (y; h(0) ) (or f1 (y; h(1) )). Let h = (p, h(0) , h(1) , s02 , s12 )T be the vector of all the parameters. For simplicity of notation hereafter the parameters not involved in the distribution will be omitted, e.g. fs (y) ≡ fs (y; h(s) , ss2 ), s = 0, 1. The marginal pdf of y is then fy (y) ≡ fy (y; h) = (1 − p)f0 (y) + pf1 (y). For each of the points [L]1 , . . . , [L]k some observations are taken independently. Let yi = (yi1 , . . . , yini )T be the vector of observations taken at [L]i , i = 1, . . . , k. The log-likelihood function is then

Ly (h) =

ni k  

logfy (yij ).

(14)

i=1 j=1

Maximizing this function is usually unaffordable as mentioned above. In order to compute approximate values of the MLEs of the parameters (say hˆ MLE ) the EM algorithm is frequently used in practice. This algorithm considers first the expectation (E step) of the logarithm of the joint distribution and then maximizes (M step) it. Specifically, let w = (y, z), where y is observed and the value of z is unknown (missing). For

154

J. López-Fidalgo, M. Rodríguez-Hernández / Chemometrics and Intelligent Laboratory Systems 154 (2016) 150–161

each value of y, say yij , the probability that this observation comes from the population z = 1, which is the same as the conditional expected value of z given yij , is defined by using the Bayes theorem,

uij = Ez|y (z|y = yij ) = P(z = 1|y = yij ) = fz|y (1|yij ) =

pf1 (yij ) , i = 1, . . . , k, j = 1, . . . , ni . (1 − p)f0 (yij ) + pf1 (yij )

This value depends on the initial parameters. The probability that the observation yij comes from population z = 0 is then P(z = 0|y = yij ) = 1 − uij . The conditional pdf of y given z can be written as fy|z (yij |zij ) = (1 − zij )f0 (yij ) + zij f1 (yij ), i = 1, . . . , k, j = 1, . . . , ni . Therefore, the joint pdf of w can be expressed in two possible ways: fw (yij , zij ) = fy|z (yij |zij )fz (zij ) = (1 − zij )(1 − p)f0 (yij ) + zij pf1 (yij )   = fz|y (zij |yij )fy (yij ) = (1 − zij )(1 − uij ) + zij uij fy (yij ). Thus, the log-likelihood function for the joint distribution is then

Lw (h) =

ni k  

((1 − zij )log(1 − p) + zij logp + (1 − zij )logf0 (yij ) + zij logf1 (yij ))

i=1 j=1

= Ly (h) +

ni k   

 (1 − zij )log(1 − uij ) + zij log(uij ) .

i=1 j=1

Taking the conditional expectation given the observed values of y (E step) then:   L∗y (h) = E Lw (h)|yij , i = 1, . . . , k, =

ni k  

j = 1, . . . , ni ,

((1 − uij )log(1 − p) + uij logp + (1 − uij )logf0 (yij ) + uij logf1 (yij ))

i=1 j=1

= Ly (h) +

ni k   

 (1 − uij )log(1 − uij ) + uij loguij .

(15)

i=1 j=1

Maximizing this function in order to obtain an estimate (say hˆ EM ) is again too difficult (M step). The EM algorithm uses this idea to approximate the MLEs. It consists on approximating the values of uij first and then computing the maximum of Eq. (15). Specifically, the EM algorithm may be formulated as follows: (0)

1. Given an initial estimation of the parameters, say h0 , all the probabilities uij , i = 1, . . . , k, j = 1, . . . , ni , are computed.

(t−1)

2. At step t − 1 of the algorithm there is an estimation of the parameters, which is used to calculate the probabilities u11 the new estimates are computed as follows:

(t−1)

, . . . , ukn

k

. Then,

hˆ (t) = arg max E (Lw (h)) h

= arg max h

nk  k  

h(t)

(t−1)

) log(1 − p) +uij

(t−1)

log p + (1 − uij

(t−1)

) log f0 (yij ) + uij

 log f1 (yij ) .

i=1 j=1

In particular, pˆ (t) =

hˆ (t) = arg max

(t−1)

(1 − uij

1 n

k

nk  k  

i=1

nk j=1

(t−1)

(1 − uij

(t−1)

uij

and

(t−1)

) log f0 (yij ) + uij

 log f1 (yij ) .

i=1 j=1

(t−1)

The last optimization procedure usually yields the ordinary MLEs with the corresponding weighted means through 1 − uij

for

(t−1) uij

population 0 and for population 1, i = 1, . . . , k, j = 1, . . . , ni . 3. The procedure stops whenever all the estimates obtained at one step do not differ from those obtained at the previous step in more than some previously defined quantities. It is well known that the convergence of the algorithm for some distributions is rather sensitive to the initial estimates of the parameters. Therefore, it is quite important to have good information about the possible values of the parameters.

J. López-Fidalgo, M. Rodríguez-Hernández / Chemometrics and Intelligent Laboratory Systems 154 (2016) 150–161

155

3.1. Mixture of Adair models The EM algorithm for this model starts assuming some initial values of the parameters and they will be re-estimated at each step. At step t the algorithm may be formulated as follows: (t)

1. Given a previous estimate, say h(t−1) , the proba bilities uij , i = 1, . . . , k, j = 1, . . . , ni are computed. It is assumed that ni experiments are performed for a specific concentration [L]i . 2. The parameters are estimated as follows. For the proportion p we have the usual estimate: k ni 1  (t) pˆ(t) = uij . ni i=1 j=1

The estimation of the parameters of the Adair monomer and dimer models is shown below. For other mixtures the formulae are similar. The estimates of the association constants are then updated as follows, k Kˆ (0)(t) =

i=1

k i=1

k (1)(t) Kˆ 1

i=1

=

k i=1

k i=1

(1)(t) = Kˆ 2

ni

(t)

(1−uij )yij [L]i j=1 (1+Kˆ (0)(t−1) [L] )3 i

ni

(t)

(1−uij )(1−yij )[L]2i

,

j=1 (1+Kˆ (0)(t−1) [L] )3 i

ni

uij (1+2Kˆ 2 (t)

(1)(t−1)

[L]i )yij [L]i

j=1 (1+Kˆ (1)(t−1) (1+Kˆ (0)(t−1) [L] )[L] )2 i i 1 2

ni

uij (1+2Kˆ 2 (t)

(1)(t−1)

[L]i )2 [L]2i

,

j=1 2(1+Kˆ (0)(t−1) [L] (1+Kˆ (1)(t−1) [L] ))3 i i 1 2

ni j=1

(t) (1)(t−1) (1)(t−1) (1)(t−1) uij (2+Kˆ 1 [L]i )(2yij +Kˆ 1 [L]i (−1+2yij +2Kˆ 2 yij [L]i ))[L]2i (1)(t−1) (1)(t−1) 2(1+Kˆ [L] (1+Kˆ [L] ))3 i

1

k i=1

ni

(t)

2 (1)(t−1)

(1)(t−1)

i

uij (2+k1 [L]i )k1 [L]4i j=1 (1+Kˆ (1)(t−1) [L] (1+Kˆ (1)(t−1) [L] ))3 i i 1 2

,

and the estimates of the variances are k 2(t) sˆ 0 =

i=1

 (t) (1 − u ) yij − j=1 ij

ni

k i=1

k 2(t) sˆ 1 =

i=1

ni j=1

(t)

uij

ni j=1

(1 −

 yij −

Kˆ (0)(t−1) [L]i 1+Kˆ (0)(t−1) [L] (t) uij )

2 i

,

(1)(t−1) (1)(t−1) ˆ (1)(t−1) Kˆ 1 K2 [L]i +2Kˆ 1 [L]2i

(1)(t−1) (1)(t−1) ˆ (1)(t−1) K2 2(1+Kˆ 1 [L]i +Kˆ 1 [L]2i )

k

i=1

ni

j=1

(t)

uij

2

.

3. The procedure stops whenever all the estimates obtained at one step do not differ from those obtained in the previous step in more than some previously defined quantities d = (d1 , d2 , d3 . . . ). 4. Illustrative example The mean value of the proportion of the number of occupied sites will be modeled by,

g([L]; h) = (1 − p)g0 ([L]; h(0) ) + pg1 ([L]; h(1) ),

with the vector of parameters h(0) and h(1) corresponding to the first and second Adair models. Normality in the observations from each population is assumed. Thus, the logarithm of the likelihood function is given by Eq. (13). For this model there is not an analytic expression for the expectation involved in the information matrix Eq. (2). The expectation is being approximated through the procedure described below. It depends on the parameters and then some nominal values need to be assigned.

Using data this matrix can be computed empirically, but in the time of designing the experiment there are not data. Thus, some approximation has to be used instead. The main problem for that is that fy (y, [L]i ; h)2 appears in the denominator making the usual approximation with the Taylor expansion rather inefficient. Some possibilities, as approximating the global matrix or just the function fy (y, [L]i ; h), were carried out until the second order, but the results were not satisfactory at all. Most of the available possibilities were local approximations both in the possible values of y and in the parameters. This was not very promising from a robustness point of view. As mentioned above for nonlinear models the FIM depends on the parameters. Usually, nominal values are considered for them in order to obtain locally optimal designs. With the nominal values of the parameters a simulation can be preformed and the results can be used to estimate the expected value. This is the main contribution of this work. Then a FIM can be approximated for each point. Taking

156

J. López-Fidalgo, M. Rodríguez-Hernández / Chemometrics and Intelligent Laboratory Systems 154 (2016) 150–161

into account that the information matrix is additive it can be constructed for any possible experimental design. Here it is an explanation of how this was performed. In particular, a value of the saturate rate, y, for a particular value [L], is generated in the following way: 1. Nominal values of the parameters, h = (p, h(0) , h(1) , s02 , s12 )T , are chosen, say h0 . In this paper p will be assumed known. In any case, at this point there are specific values for all the parameters in order to perform the simulations. 2. Then, a value, z (0 or 1), is generated from the Bernoulli distribution b(p). 3. Independently a value from the standard Gaussian distribution will be generated, say e. 4. The value of y is then built using the values z and e as follows: • If z = 0 then the observation generated is y = g0 ([L]; h(0) ) + s 0 e. • If z = 1 then the observation generated is y = g1 ([L]; h(1) ) + s 1 e. 5. This process is repeated a number of times, say S, and then the expectation is approximated in the usual way,   S −1  ∂ 2 log fy (y; [L], h) I([L]; h) ≈ . S ∂ h2

(16)

Fig. 2. Determinant of the information matrix of one-point designs against this optimal point (X-axis). The big dot corresponds to the D-optimal design.

The nominal values of the parameters were taken as K(0) = 0.1, K(1) = 0.3, s02 = (0.02)2 and s12 = (0.015)2 . The D-optimal design obtained is a one–point design at [L]∗ = 5.8, n∗ =



5.8 1



.

r=1

6. The design space has to be discretized in a feasible number of points {[L]1 , . . . , [L]k }. Thus, after this process there will be an approximated FIM for each point, 

 I([L]1 ; h0 ), . . . , I([L]k ; h0 ) .

The FIM for any experimental design in the discretized space is then

M(n; h0 ) =

h  ni I([L]i ; h0 ). n

(17)

i=1

A number of fixed values of zs and es are generalized. The same will be used for building the values of ys for each [L]i . At each point 10,000 simulations are made for each of the three cases of mixtures studied. The proportion of the concentration of the two types of macromolecules is assumed known, p = 0.75, for the three cases. Then, the remaining parameters, h(0) , h(1) , s02 and s12 , have to be estimated. A neighborhood around the optimal designs will be checked identifying the region with efficiencies greater than 95% (Figs. 2, 3 and 4). The efficiency of a design n will be computed as follows:  ef fD (n) =

|M(n, h0 )| |M(|n∗h , h0 )|

Fig. 5 (left) shows that the equivalence theorem is verified for this design and therefore it is optimum. The maximum is 4, the number of unknown parameters to be estimated. Fig. 2 shows the determinant of the information matrix of one-point designs. For 4.6 < [L] < 7.3 the efficiency is larger than 95%. This gives an idea of the robustness of the optimal design versus misspecifications of the nominal values of the parameters. Case 2. Mixture of a monomer and a dimer The means are now g0 ([L]; h(0) ) = g1 ([L]; h(1) ) =

K (0) [L] 1 + K (0) [L]

and

(1)

(1)

(1)

K1 [L] + 2K1 K2 [L]2 2(1 +

(1) K1 [L]

+

(1) (1) K1 K2 [L]2 )

.

The values of the monomer K(0) = 0.5, s02 = (0.015)2 and the (1) (1) values of the dimer K1 = 0.1, K1 = 2, s12 = (0.02)2 were assumed

 m1 ,

0

where m is the number of parameters and n∗h is the D-optimal design 0 obtained for nominal values of the parameters h0 . Case 1. Mixture of two monomers Let the means for the two monomers be modeled as g0 ([L]; h(0) ) =

K (0) [L] K (1) [L] , g1 ([L]; h(1) ) = . 1 + K (0) [L] 1 + K (1) [L]

Fig. 3. Graph representing the first (X-axis) and the second (Y-axis) points of designs with a D-efficiency greater than 95% for a mixture of two monomers. The big dot corresponds to the D-optimal design.

J. López-Fidalgo, M. Rodríguez-Hernández / Chemometrics and Intelligent Laboratory Systems 154 (2016) 150–161

157

Fig. 4. Graph representing the second (X-axis) and the third (Y-axis) points of designs with a D-efficiency greater than 95% for a mixture of two dimers in three cases for the first point: [L]1 = 0.9 (left), [L]2 = 1 (center), [L]3 = 1.1 (right). The big dot corresponds to the D-optimal design.

as nominal values of the parameters. The D-optimal design obtained is a two-point design at [L]∗1 = 1.3 and [L]∗2 = 4.9,

n˜ ∗ =



4.9 1.3 0.5775 0.4225



.

The equivalence theorem was checked numerically for this design and it was verified (Fig. 5, center), so it is actually D-optimal. The maximum is 5, the number of unknown parameters to be estimated. Fig. 3 shows the first (X-axis) and the second (Y-axis) points of designs with D-efficiencies above 95%. This gives an idea of the robustness of the optimal design versus misspecifications of the nominal values of the parameters.

Case 3. Mixture of two dimers The means are now

g0 ([L]; h(0) ) = g1 ([L]; h(1) ) =

4.1. Comparison of the Information Matrix and the covariance matrix

(0) (0) (0) K1 [L] + 2K1 K2 [L]2 (0) (0) (0) 2(1 + K1 [L] + K1 K2 [L]2 ) (1) (1) (1) K1 [L] + 2K1 K2 [L]2 (1) (1) (1) 2(1 + K1 [L] + K1 K2 [L]2 )

and

.

The nominal values of the parameters used for computations (0) (0) (1) (1) were K1 = 0.1, K2 = 2, K1 = 0.5, K2 = 0.5, s02 = (0.015)2 and s12 = (0.02)2 . The D-optimal design obtained has three points in its support: [L]∗1 = 1, [L]∗2 = 2.2 and [L]∗3 = 4.8,

n˜ ∗ =



4.8 1 2.2 0.5352 0.1115 0.3533

The equivalence theorem was checked numerically for this design and it was verified (Fig. 5, right), so it is actually D-optimal. The maximum is the number of unknown parameters to be estimated, 6 in this case. Fig. 4 shows three graphs for concentrations of free ligand at the first point of the design [L]1 = 0.9 (left), [L]2 = 1 (center), [L]3 = 1.1 (right), where the X-axis and Y-axis show the second and third points of the values of the concentrations of free ligand of the design such that the D-efficiency is greater than 0.95. These calculations were performed for values of the second point from 1.1 to 4.7 and the third point from 4.3 to 6. The graphs of Figs. 2, 3 and 4 provide the experimenter the loss of efficiency when the design moves away from the optimal design, which can be used to provide appropriate nominal values. All these calculations were made with the same nominal values, changing the support points and optimizing the weights.



.

Optimal designs were computed through the information matrix, assuming there is a good approximation of the covariance matrix with the inverse of the FIM divided by the sample size. In this section an empirical verification of this approximation is being analyzed. The EM algorithm will be used to compute the MLEs in a number of simulations in the three cases of study. The stopping rule for the EM algorithm was set as 10 −4 for all the parameters. As expected, for small samples the matrices and the determinants differ significantly. For computing the covariance matrix n observations are generated in ten thousand simulations. At each point [L]i of the optimal  design ni ≈ n × pi observations were obtained, with ki=1 ni = n. The ni observations, yij , come from the first population with a proportion of 25% and from the second with a proportion of 75%. The generated observations are used for estimation of the parameters, which is carried out using the EM algorithm given in Section 3.1. Then the matrices are compared in the three cases.

Fig. 5. Sensitivity functions (Y-axis) for the free ligand concentrations of the design space (X-axis) showing that the equivalence theorem is satisfied for the mixture of two monomers (left), the mixture of one monomer and one dimer (center) and the mixture of two dimers (right).

158

J. López-Fidalgo, M. Rodríguez-Hernández / Chemometrics and Intelligent Laboratory Systems 154 (2016) 150–161

Table 1 Mean, bias and variance of the MLEs. The nominal values of the parameters were taken as K(0) = 0.1,K(1) = 0.3,s02 = (0.02)2 and s12 = (0.015)2 .

ˆ E[h] ˆ |h − h| ˆ Var[h]

Kˆ (0)

Kˆ (1)

sˆ 02

sˆ 12

0.10 3.03 10 −5 3.04 10 −6

0.30 7.82 10 −6 4.96 10 −6

3.84 10 −4 1.52 10 −5 1.31 10 −8

2.22 10 −4 2.17 10 −6 1.37 10 −9

ˆ E[h] ˆ |h − h| ˆ Var[h]

Case 1. Mixture of two monomers The D-optimal design corresponds to the concentration of the free ligand [L]* = 5.8 and the inverse of the FIM calculated from the simulated observations divided by n = 100 is, ⎛

3.08 ⎜ 1.23 ⎜ ⎝ 1.33 −4.78

10−6 10−29 10−9 10−31

1.23 10−29 1.33 4.96 10−6 6.24 6.24 10−30 1.44 9.43 10−11 −2.41

⎞ 10−9 −4.78 10−31 10−30 9.43 10−11 ⎟ ⎟, 10−8 −2.41 10−31 ⎠ 10−31 1.26 10−9

with determinant 2.7910 −28 . The empirical covariance matrix from the simulations is ⎛

3.04 ⎜ −1.94 ⎜ ⎝ 3.21 −1.82

10−6 −1.94 10−8 10−8 4.96 10−6 10−10 5.62 10−10 10−10 −9.57 10−10

3.21 5.62 1.31 1.46

⎞ 10−10 −1.82 10−10 −10 −10 ⎟ 10 −9.57 10 ⎟, 10−8 1.46 10−11 ⎠ 10−11 1.37 10−9

with determinant 2.7210 −28 . It is observed that the variances of the parameters which are in the main diagonal of the inverse of the information matrix and of the above two matrices have similar values while the covariance parameters almost vanished. The value of the determinant is also very similar. Table 1 shows the MLEs obtained by using the EM algorithm. The means of the parameter estimates are quite similar to the nominal values of the parameters, i.e. the biases, as well as the variances, are rather small. Case 2. Mixture of a monomer and a dimer In this case the optimal design has two points, [L]∗1 = 1.3 and [L]∗2 = 4.9, with weights 0.577 and 0.422 respectively. The inverse of the FIM divided by n = 100 was not so similar to the covariance matrix. The sample size was increased to 300 observations and now the inverse of the FIM divided by 300 resulted very similar to the covariance matrix. The FIM was ⎛

1.58 10−5 ⎜ −5.96 10−7 ⎜ ⎜ 1.35 10−5 ⎜ ⎝ 7.93 10−9 −4.47 10−9

−5.96 10−7 3.95 10−5 −5.49 10−4 −1.15 10−8 9.41 10−9

1.35 10−5 7.93 10−9 −4.47 10−9 8.05 10−4 −1.15 10−8 9.41 10−9 1.66 10−2 2.53 10−7 −2.05 10−7 2.53 10−7 1.68 10−9 −9.76 10−11 −2.05 10−7 −9.76 10−11 1.42 10−9

Table 2 Mean, bias and variance of the MLEs. The nominal values of the parameters were taken (1) (1) as K(0) = 0.5,K1 = 0.1,K2 = 2,s02 = (0.015)2 and s12 = (0.02)2 .

⎞ ⎟ ⎟ ⎟, ⎟ ⎠

Kˆ (0)

(1) Kˆ 1

(1) Kˆ 2

sˆ 02

sˆ 12

0.50 3.44 10 −5 1.59 10 −5

0.10 1.36 10 −3 3.35 10 −5

1.97 2.04 10 −2 1.42 10 −2

2.22 10 −4 2.92 10 −6 1.55 10 −9

3.96 10 −4 3.04 10 −6 1.53 10 −9

The EM algorithm stopped when all the estimates obtained at one step did not differ from those obtained at the previous step in more than di = 10 −3 then the determinants of the covariance matrices were quite different from that obtained in the inverse of the matrix information. Choosing di = 10 −4 the covariance matrix and its determinant were very similar to the inverse of the information matrix and its determinant. As mentioned above a suitable sample size in this case was n = 300, i.e. 173 observations at [L]1 = 1.3 and 127 at [L]2 = 4.9. Table 2 shows some properties of the MLEs obtained by means of the EM algorithm. The biases are small and the variances minimal. Case 3. Mixture of two dimers The optimal design has three points, [L]∗1 = 1, [L]∗2 = 2.2 and ∗ [L]3 = 4.8, whose weights are 0.535, 0.111 and 0.353 respectively. The inverse of the information matrix divided by n was quite different from the covariance matrix for n = 100 and n = 300. The inverse of the FIM divided by n = 1000 was ⎛

2.96 10−5

⎜ −6.22 10−4 ⎜ ⎜ ⎜ 2.64 10−6 ⎜ ⎜ −3.13 10−6 ⎜ ⎜ ⎝ 3.45 10−8 2.42 10−9

−6.22 10−4

2.64 10−6

−3.13 10−6

3.45 10−8

2.42 10−9

1.33 10−2

5.61 10−5

7.01 10−5

−7.24 10−7

−7.87 10−8

−5.61 10−5

4.37 10−5

−4.54 10−5

4.59 10−9

−2.38 10−8

7.01 10−5

−4.54 10−5

5.71 10−5

−2.59 10−10

2.81 10−8

−7.24 10−7

4.59 10−9

−2.59 10−10

6.97 10−10

−2.97 10−12

2.81 10−8

−2.97 10−12

5.29 10−10

−7.87 10−8 −2.38 10−8

⎞ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎠

(19) with determinant 1.01610 −36 . The empirical covariance matrix was ⎛

2.31 10−5

⎜ −4.84 10−4 ⎜ ⎜ ⎜ −5.64 10−7 ⎜ ⎜ −3.59 10−7 ⎜ ⎜ ⎝ 2.88 10−8 1.47 10−9

−4.84 10−4 −5.64 10−7 1.04 10−2

9.18 10−6

−6

−5

9.18 10

4.69 10

−1.39 10−6 −4.97 10−5 −6.09 10

−7

−5.80 10

−8

5.80 10

−9

−3.71 10

−8

3.59 10−7

2.88 10−8

−1.39 10−6 −6.09 10−7 −4.97 10

−5

6.27 10−5 −2.44 10 4.52 10

−9

−8

5.80 10

−9

−2.44 10−9 6.40 10

−10

−2.31 10

−11

1.47 10−9 −5.80 10−8 −3.71 10

−8

4.52 10−8 −2.31 10−11

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

6.04 10−10

with determinant 1.02910 −36 . Table 3 displays some properties of the MLEs.

(18)

with determinant 4.1110 −31 . The empirical covariance matrix was

⎛

1.59 10−5 ⎜ ⎜ −5.48 10−7 ⎜ ⎜ 1.26 10−5 ⎜ ⎜ ⎝ 6.72 10−9 −3.52 10−9

−5.48 10−7 1.26 10−5

6.72 10−9

3.35 10−5 −6.82 10−4 −9.72 10−9 −6.82 10−4 1.42 10−2

1.94 10−7

−9.72 10−9 1.94 10−7

1.55 10−9

7.80 10−9 −1.57 10−7 −9.80 10−11

with determinant 4.5610 −31 .

−3.52 10−9

⎞

⎟ 7.80 10−9 ⎟ ⎟ −1.57 10−7 ⎟ ⎟ ⎟ −9.80 10−11 ⎠ 1.53 10−9

Remark 1. It is observed empirically that suitable sample sizes are those ensuring that there are at least about 30 observations at each design point for each of the normal distributions in the mixture. This result was observed after a number of simulations with a great computational burden. The sample size was tuned doing simulations for the designs used in practice increasing the replications for each point until the empirical covariance matrix was similar to the inverse of the information matrix (divided by the sample size). After the computation of the optimal design it was checked that the approximation worked reasonably. For example, for the Case 2 a sample size of n = 300 meant 173 observations at the first design point. Since there is a mixture of two populations in a proportion of 25% and 75% for

J. López-Fidalgo, M. Rodríguez-Hernández / Chemometrics and Intelligent Laboratory Systems 154 (2016) 150–161

159

Table 3 Mean, bias and variance of the estimation of the parameters obtained by the EM algorithm. The nominal values of the parameters were (0) (0) (1) (1) taken as K1 = 0.1,K2 = 2,K1 = 0.5,K2 = 0.5,s02 = (0.015)2 and s12 = (0.02)2 .

ˆ E[h] ˆ |h − h] ˆ Var[h]

(0) Kˆ 1

(0) Kˆ 2

(1) Kˆ 1

(1) Kˆ 2

sˆ 02

sˆ 12

0.10 1.64 10 −3 2.31 10 −5

1.97 2.97 10 −2 1.04 10 −2

0.49 4.26 10 −4 4.69 10 −5

0.50 5.95 10 −4 6.27 10 −5

2.24 10 −4 1.37 10 −7 6.40 10 −10

3.98 10 −4 1.06 10 −6 6.04 10 −10

both populations the 173 observations were split in about 43 and 130 observations respectively. Similarly, for the second point with 127 observations divided in 32 and 95 for the two populations. The number of simulations was for a good approximation of the expectation was also tuned in a similar way. It is also important to note that the stopping bound chosen in the EM algorithm was di = 10 −4 . A bit higher values still provide good estimates of the parameters but not good approximations of the the covariance matrix. For the three cases considered the variances of MLEs, which are in the main diagonal of the covariance matrix had similar values and the covariances were near zero. The values of the determinants of the inverse of the FIM (divided by n) and the covariance matrix were also very similar.

4.2. Sensitivity of the design with respect to the true values of the parameters

Case 1. Mixture of two monomers Let h0 ± e = (K (0) ± eK (0) , K (1) ± eK (1) , s02 ± es 2 , s12 ± es 2 ) be a neigh0

Case 2. Mixture of one monomer and one dimer (1) The neighborhood of h0 is given by h0 ± e = (K (0) ± eK (0) , K1 ± (1) 2 2 eK (1) , K2 ± eK (1) , s0 ± es 2 , s1 ± es 2 ), where eK (0) = 0.05, eK (1) = 0.05, 1

 ef fh∗ (nh0 ) =

|M(nh0 , h∗ )| |M(nh∗ , h∗ )|

 m1

.

This efficiency will be studied for values of h* in a neighborhood of h0 . In the usual theory there is an explicit expression of the FIM and the computations of the optimal designs and the efficiencies, although time consuming, are reasonable for a number of possible values of h*. In this framework the FIM for each particular design point, and for each particular values of h*, has to be computed doing simulations. This means that computing all these efficiencies takes a considerable time. For instance, choosing just three possible values of each parameter (the nominal, one smaller and one larger values) would mean 3m experiments. In the first case there are m = 4 parameters, then 81 efficiencies will be needed, which is a lot, but reasonable. For the second case with m = 5 parameters, 243 would be needed and in the third case with m = 6 parameters, 729 will be needed. In order to reduce those numbers in a meaningful way fractional factorial designs of size 8 are being considered for the three cases. The number of simulations performed in all cases were 10,000.

0

2

1

1

eK (1) = 0.05, es 2 = 10−4 and es 2 = 10−4 . The 8 cases following the 0

2

As mentioned above the designs and the information matrix depend on the choice of the nominal values of the parameters. A sensitivity analysis has been conducted in order to check how the efficiency changes if the nominal values, h0 , are wrongly chosen. This will be made basically in the same way Optimal Experimental Design theory does. This means precisely checking how risky would be choosing wrongly the nominal values. As a matter of fact, this analysis helps the practitioner to choose appropriate nominal values. In particular, the true value, say h*, is unknown, but it is expected to be near the guessed value, h0 . Let nh0 be the optimal design computed with the nominal values of the parameters and let nh∗ be the optimal design computed with a true possible value of the parameters. Then the efficiency of the design nh0 will be obtained assuming in the calculation of the FIM that h* is the true value and thus nh∗ is the true optimal design,

1

borhood with eK (0) = 0.05, eK (1) = 0.05, es 2 = 10−4 and es 2 = 10−4 . 0 1 This can be seen as a factorial design of 16 points. A fractional design (0) (1) 2 2 of 8 points was chosen with the rule K K s0 s1 = I. Table 4 shows the efficiencies for the eight cases with the typical notation of the signs meaning addition or subtraction from the nominal values of the parameters. The efficiencies are greater than 0.96, which means the design is rather robust. For example, if the experimenter gives a wrong value of the parameter to calculate the optimal design, the true value may be assumed to be in a neighborhood of the nominal value used to compute the optimal design. In the neighborhood considered here there is only a loss of efficiency of at most 4%.

(1)

1

(1)

(1)

(1)

criteria K (0) K1 K2 s02 s12 = I and K (0) K1 K2 = I were considered. The efficiencies are approximately over 0.80 as shown in Table 5, which means that the efficiency is safe enough. In this case the loss of efficiency is under 20%, which is important, but not too risky. It is observed that varying the values of some parameters designs with two and three points are obtained. In particular, the third and the seventh deliver three-point designs. The efficiency in these cases was higher than for the other experiments whose designs were two points. It is noticeable that if the signs of the association constants are the same the optimal designs obtained have similar points and weights. For instance, the first and the fifth are quite similar as well as the second and the sixth. Case 3. Mixture of two Adair dimers (0) For the last case the neighborhood of h0 is taken by h0 ±e = (K1 ± (0) (1) (1) eK (0) , K2 ±eK (0) , K1 ±eK (1) , K2 ±eK (1) , s02 ±es 2 , s12 ±es 2 ), where eK (0) = 1

2

1

0

2

1

1

0.05, eK (0) = 0.05, eK (1) = 0.05, eK (1) = 0.05, es 2 = 10−4 and es 2 = 2

1

2

0

(0)

1

(0)

(1)

(1)

10−4 . There are 8 cases matching the criteria K1 K2 K1 K2 s02 s12 = (0) (0) (1) (1) (0) (0) (1) (1) I, K1 K2 K1 K2 s02 = I and K1 K2 K1 K2 s12 = I. The efficiencies, shown in Table 6, have high values for the experiments three, five,

Table 4 Relative efficiencies of the optimal design for h0 as function of the possible true values, h* for the first case. h*

nh∗

ef fh∗ (h0 )

(−, −, −, −) (+, +, −, −) (+, −, +, −) (−, +, +, −) (+, −, −, +) (−, +, −, +) (−, −, +, +) (+, +, +, +)

8.9 4.4 5.2 7.6 5.2 7.6 8.9 4.4

0.9608 0.9808 0.9960 0.9861 0.9959 0.9861 0.9608 0.9808

160

J. López-Fidalgo, M. Rodríguez-Hernández / Chemometrics and Intelligent Laboratory Systems 154 (2016) 150–161

Table 5 Optimal design point (weights) and relative efficiencies of the optimal design for h0 as function of the possible true values h* for de second case. Numbers in parentheses are the weights of the optimal design points. h* (+, −, −, −, −) (−, +, −, −, −) (−, −, +, −, −) (+, +, +, −, −) (+, −, −, +, +) (−, +, −, +, +) (−, −, +, +, +) (+, +, +, +, +)

1.5(0.550) 0.8(0.407) 1.6(0.430) 1(0.497) 1.4(0.556) 0.8(0.347) 1.6(0.374) 0.9(0.488)

nh∗

ef fh∗ (h0 )

3.9(0.450) 2.9(0.593) 2.9(0.321) 3.3(0.503) 3.4(0.444) 3.2(0.653) 2.4(0.344) 3.7(0.512)

0.8663 0.8084 0.8973 0.8504 0.8737 0.8676 0.9209 0.8598

7.5(0.249)

6.6(0.282)

seven and eight over 0.96. In all cases the efficiencies are greater than 0.89. Thus, giving a wrong value h0 and assuming this neighborhood for the true value there is not more loss of efficiency than 11%. This is just to let the experimenter know the risk he or she is assuming. Then, increasing 11% the sample size may be wanted. Again for different combinations of the nominal values of the parameters designs with two and three points are obtained. In the third and fifth experiments the designs have two support points.

5. Discussion For any non-linear model the Information Matrix depends on the nominal values of the parameters. If the distribution considered is in the Exponential Family, e.g. a Normal distribution, there is an explicit expression of the Information Matrix. Then there are mainly three different ways to tackle the problem of the unknown parameters. On the on hand a sequential design, if possible, improves the nominal values of the parameters once a new observation is available. A second possibility is a Bayesian approach considering a prior distribution of the parameters. A third way is just to choose nominal values of the parameters once for all and design the experiment using them. These are the so called locally optimal designs. In this case a sensitivity analysis needs to be performed in order to check the impact of a wrong choice of the nominal values on the efficiency of the optimal design. This third approach has been used in this paper with an important remark. Since a mixture of two Normal distributions does not belong to the exponential family an analytic expression of the Information Matrix is not available anymore. Then we perform simulations to approximate it at different points of a discretized design space. In order to choose a right sample size and a sufficient number of simulations a study was performed to determine both empirically. In practical terms the optimal design depends here on the nominal values of the parameters in the same way it depends in an ordinary Optimal Design problem with the Exponential family. Of course, the other two approaches can also be considered under this new philosophy and can be implemented in further research.

Table 6 Optimal design point (weights) and relative efficiencies of the optimal design for h0 as function of the possible true values h* for the second case. Numbers in parentheses are the weights of the optimal design points. h* (−, −, −, −, +, +) (+, +, −, −, +, +) (+, −, +, −, +, +) (−, +, +, −, +, +) (+, −, −, +, +, +) (−, +, −, +, +, +) (−, −, +, +, +, +) (+, +, +, +, +, +)

3.3(0.325) 2.3(0.201) 4.1(0.479) 3.5(0.362) 3.8(0.441) 3.7(0.385) 4.3(0.444) 1.3(0.440)

p

nh∗

0 0.01 0.02 0.03 0.05 0.1 0.3 0.5 0.7 0.9 1

2(1) 3(1) 3.1(1) 1.7(0.9203) 1.2(0.5248) 1.4(0.6056) 1.2(0.5352) 1.2(0.5402) 1.3(0.5896) 1.3(0.5762) 1.2(0.5)

5(0.0797) 5.1(0.4752) 5(0.3944) 5.1(0.4648) 5(0.4598) 4.8(0.4104) 4.9(0.4238) 4(0.5)

Mixing two macro-molecules with one and two binding sites and a ligand there is a mixture of two Adair’s models. Although there are two random variables, the response y and the explanatory variable [L], in practice the independent variable may be considered without random error since it is easy to measure and has a very small error. D-optimal designs were computed for improving the estimates of the parameters. Since the concentration of free ligand [L] is not controllable experimentally a transformation to initial ligand concentration, [L0 ], which is under the control of the experimenter, was made. Three real cases were studied through simulations and estimates of the parameters were obtained. Then the means, variances and biases were obtained empirically. The covariance matrix was compared to the inverse of the information matrix to verify empirically the asymptotic convergence. Sample sizes and number of simulations were adjusted to obtain good approximations. For the nominal values of the unknown parameters, v0 , v1 and for the association constants of the Adair models, different p values were checked in order to see its influence in the optimal designs. The extreme cases correspond to no mixture, the monomer for p = 0 and the dimer for p = 1. If 0 < p < 1 there is an actual mixture of two possible Adair models. By varying the values of p in the mixture, it is observed that the design points were varying smoothly from p = 1 to p = 0.3, as shown in Table 7. Then there is a jump. The weight of the first point is increasing rapidly when p decreases from p = 0.02 to p = 0 as well as the first point is approaching the optimal value for the monomer. The optimal designs seem to depend more on changes of the association constants and the variances as it can be guessed from Tables 4, 5 and 6. These computations were performed with the same Gaussian errors, as detailed in Section 2.2, in order to make them more comparable. A robustness study was performed in order to see the impact of choosing wrongly the nominal values of the parameters. The efficiencies were over 0.80 in all cases meaning a desirable situation. The number of points of the optimal design depends on the nominal values of the parameters. The stopping bound for the EM algorithm was set to obtain good approximations. A possible future research may consist in finding optimal designs for the case of heteroscedasticity in the mixture as well as considering other models of interest with a macro-molecule and various ligands that can be bound to it. Conflict of interests

ef fh∗ (h0 )

nh ∗ 1.2(0.375) 0.7(0.373) 1(0.521) 1.1(0.384) 0.9(0.558) 1.2(0.411) 1.2(0.464) 0.9(0.184)

Table 7 Optimal design point (weights) for different values of p for the mixture of one monomer and one dimer.

5.6(0.300) 4.1(0.574) 5.8(0.252) 6.4(0.204) 7.4(0.092) 4.2(0.624)

0.9289 0.8949 0.9818 0.9158 0.9651 0.9281 0.9838 0.9783

There is no conflict of interest. Acknowledgments This work has been supported by Ministerio de Educación y Ciencia and Fondos FEDER MTM2010–20774–C03–01, Junta de Comunidades de Castilla la ManchaPEII10–0291–1850 and Ministerio de

J. López-Fidalgo, M. Rodríguez-Hernández / Chemometrics and Intelligent Laboratory Systems 154 (2016) 150–161

Economía y Competitividad and fondos FEDER MTM2013-47879-C21-P. The authors want to thank Dr. Burguillo for his helpful comments about the chemical problem and the referees for the review and appropriate suggestions.

References [1] G.S. Adair, The hemoglobin system: VI. The oxygen dissociation curve of hemoglobin, J.Biol. Chem. 63 (1925) 529–545. [3] J. López-Fidalgo, M.M. Rodríguez-Hernández, Experimental designs for the Adair model, Chemometr. Intell. Lab. Syst. 138 (2014) 133–141.

161

[4] V.V. Fedorov, Theory of Optimal Experiments, Academic press, New York, 1972. [5] J. Kiefer, General equivalence theory for optimum designs (approximate theory), Ann. Stat. 2 (1974) 849–879. [8] J.M. Marin, K. Mengersen, Ch.P. Robert, Bayesian modelling and inference on mixtures of distributions, Handbook of Statistics, Elsevier-Sciences 2005, 25 [9] H. Chernoff, Locally optimal designs for estimating parameters, Ann. Math. Stat. 24 (1953) 586–602. [10] A.C. Atkinson, A.N. Donev, R.D. Tobias, Optimum experimental designs, with SAS, Oxford Universiy Press 2007. [11] M. Imhof, J. López-Fidalgo, W.K. Wong, Efficiencies of optimal approximate designs for small samples, Statista Neerlandica 55 (2001) 301–315. [12] J. Kiefer, Optimum experimental designs, J. R. Stat. Soc. Ser. B 21 (1959) 272–319.