**1. Introduction**

In pharmacokinetic/pharmacodynamic studies, repeated measurements are performed on a sample of individuals (units/subjects), and responses for all experimental subjects are assumed to be described by a common structural model. This model contains both fixed and random effects to distinguish between individuals in a population, leading to a mixed effects model in which fixed effects represent fixed parameters for all individuals in the population and random effects account for individual differences. Moreover, mixed effects modeling has been shown to be useful in pharmacokinetic/pharmacodynamic studies; particularly for modeling total variation in and between individuals, the models used in pharmacokinetic/ pharmacodynamic analysis are often presented by a system of deterministic differential equations (ordinary, partial, or delay). However, the pharmacological processes in reality are always exposed to incomprehensible effects that are difficult to model, and ignoring these effects can affect the parameter estimation results and their interpretations. So, the introduction of stochastic components to deterministic models is an important tool of analysis [1] and is more appropriate to model the intra-individual variations rather than ODEs. In addition, the extension of ODEs to the SDEs makes it possible to explain the differences between the observations and the predictions by two types of noise: dynamic noise that enters through the dynamics of the system and that can result from its random fluctuations or from the shortages of model and the measurement error which is added in the case of an indirectly observed process, which may be due to a test error or to the existence of a disturbance and represent the uncorrelated part of the residual variability. In the theory, there are rich and developed resources for mixed effects models whether deterministic [2–6] or stochastic and linear or nonlinear models (see many applications of stochastic NLME models in biomedical [7–9] and in pharmacokinetic [10–12] fields).

approximating methods of the likelihood function when the integral given the random effects has no analytic solution. Then, we propose an optimization algorithm to obtain maximum likelihood estimators in order to facilitate the estimation procedure for these models. Finally, the methodology is evaluated by simulation

*Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential…*

implemented on the minimal model describing the glucose and insulin kinetics.

Consider an N-dimensional continuous and stochastic process *Yt* in the state space *E* ⊂ *<sup>N</sup>* described by the general first-order nonlinear stochastic differential

> *t* , *Y<sup>i</sup>*

*<sup>t</sup>* is defined as the solution of the SDME model Eq. (1) that exists under some conditions that we suppose satisfied [31–33] and represents the observation of individual i from M different experimental units, (*i* ¼ 1, *::*, *M*), at the

n o is assumed to verify the same model structure Eq. (1) according

∈ *B*⊆ *<sup>q</sup>* are the q-dimensional individual random parameters assumed mutually

<sup>j</sup><sup>Ψ</sup> � � depending on a population parameter <sup>Ψ</sup>; in the population approach, this parameter vector allows for a data from several subjects to be considered simulta-

Gaussian distribution; however, any other distributions may be considered contin-

The joint density function of the vector *b<sup>i</sup>* is parameterized by a q-dimensional parameter *<sup>ϑ</sup>*∈*υ*<sup>⊂</sup> *<sup>q</sup>* and a *<sup>q</sup>* � *<sup>q</sup>*-dimensional matrix *<sup>ϕ</sup>*<sup>∈</sup> <sup>Φ</sup> <sup>⊂</sup> *<sup>q</sup>*�*<sup>q</sup>* representing the covariance matrix of *b<sup>i</sup>* and specifying the parameters of the marginal distributions

unknown, they must be considered as random effects since they change between individuals; but in some cases it can be known and assumed equal to a real constant.

> 0 <sup>j</sup> *yt*0, *bi*

ð Þ*t* are the standard Brownian motions, and they are assumed mutually independent with *bj* for all 1≤*i*, *<sup>j</sup>*<sup>≤</sup> *<sup>M</sup>*. The functions *<sup>μ</sup>*ð Þ� : *<sup>E</sup>* � � <sup>Θ</sup> � *<sup>B</sup>* ! and

a strictly positive density with regard to the Lebesgue measure on E:

*y* ! *PY y*, *t* � *t*

ter which represents the same and common characteristics for all subjects, and

independent, also called random effects because they are not the same for all individuals; they change between them according to a distribution of density

<sup>0</sup> <sup>¼</sup> *<sup>y</sup><sup>i</sup>*

*<sup>t</sup>*<sup>0</sup> is the initial state of *Yt* for each subject. The process

; *θ* ∈ Θ ⊂ *<sup>p</sup>* is a p-dimensional fixed effects parame-

*<sup>l</sup>* may follow a different distribution, ð Þ *l* ¼ 1, … , *q* , and a

*<sup>b</sup><sup>i</sup>* � *<sup>i</sup>:i:d N*ð Þ *<sup>ϑ</sup>*, *<sup>ϕ</sup> :* (2)

, 1ð Þ ≤*l* ≤*q* ; the components of *ϕ* and *ϑ* represent the popula-

ð Þ*<sup>t</sup>* given *<sup>b</sup><sup>i</sup>*

0, it is not necessarily known, and when its components are

, *θ* � � and *Y<sup>i</sup> t*

, *θ* � �>0, *y*∈*E:* (3)

<sup>0</sup> ð Þ¼ *yt*0, *t*

<sup>0</sup> <*t*, has

<sup>j</sup><sup>Ψ</sup> � � of the vector *bi* could be the

<sup>0</sup>, *i* ¼ 1, … , *M*, (1)

studies on the bidimensional Ornstein-Uhlenbeck model, and then it is

*t* , *θ*, *bi* � �*dW<sup>i</sup>*

**2. Theoretical tools**

*dY<sup>i</sup>*

where *Y<sup>i</sup>*

moment *t*≥*t*

*Yi t* � �

*bi*

*PB b<sup>i</sup>*

equations of the Itô type [30]:

, *<sup>t</sup>*, *<sup>θ</sup>*, *<sup>b</sup><sup>i</sup>* � �*dt* <sup>þ</sup> <sup>Σ</sup> *<sup>Y</sup><sup>i</sup>*

*DOI: http://dx.doi.org/10.5772/intechopen.90751*

<sup>0</sup> <sup>¼</sup> *<sup>Y</sup><sup>i</sup>*

standard choice for the joint density function *PB bi*

*l*

Also, we assume that the distribution of *Y<sup>i</sup>*

*<sup>t</sup>* <sup>¼</sup> *<sup>μ</sup> <sup>Y</sup><sup>i</sup> t*

> *i* 0, and *Y<sup>i</sup>*

*<sup>t</sup>*≥<sup>0</sup>, *i* ¼ 1*::M*

to the individual parameters *b<sup>i</sup>*

neously. Each component *b<sup>i</sup>*

uous or discrete:

of the components *b<sup>i</sup>*

*W<sup>i</sup>*

**93**

tion parameters Ψ. For *Y<sup>i</sup>*

Parameter estimation in mixed effects models with SDEs, known by stochastic differential mixed effects models, is not an obvious procedure except in some cases simpler [14] because it is often difficult to write the likelihood function in its closed form. In this context, we propose a review on estimation methods of SDME models in [13, 15] and an example case that treats a generalized linear mixed models in [16] and also an example to approximate the likelihood function of an NLME with the likelihood of a linear mixed effects model in [17]. Moreover, to strengthen knowledge on estimation methods of SDME, we refer to [18, 19] that propose an example of stochastic mixed effects model with random effects log-normally distributed with a constant diffusion term.

In general, it is difficult to obtain an explicit likelihood function because the transition density of the stochastic process is often unknown or that the integral in the likelihood given the random effects cannot be computed analytically, and although the size of the random effects increases, the complexity of the problem increases also rapidly. Therefore, this requires a significant need for approximate methods to compute the transition density in an approximate closed form and for effective numerical integration methods to compute or approximate the integral in the likelihood function, using, for example, the Laplacian and Gaussian quadrature approximation [3, 8, 20] or other approaches [21, 22]. In the literature, several solutions have been proposed to approximate the transition density and have shown their effectiveness despite certain limitations. For example, the transition density could be approximated by the solution of the partial differential equations of Kolmogorov [23]; by the derivation of a Hermite expansion of closed form at the transition density [24–26] (this method has been reviewed and applied for many known stochastic processes for one-dimensional [8] and multidimensional [20] time-homogeneous SDME model); or by simulating the process to Monte-Carlointegrate the transition density [27–29]. These techniques are very useful and can solve the problem, but unfortunately, they involve intense calculations which make the problem always complicated.

In this work, we focus on two fundamental issues concerning the implementation of SDEs in NLME models. The first is how the transition density of an SDME model can be approximated when it is not known, and the second is about

*Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential… DOI: http://dx.doi.org/10.5772/intechopen.90751*

approximating methods of the likelihood function when the integral given the random effects has no analytic solution. Then, we propose an optimization algorithm to obtain maximum likelihood estimators in order to facilitate the estimation procedure for these models. Finally, the methodology is evaluated by simulation studies on the bidimensional Ornstein-Uhlenbeck model, and then it is implemented on the minimal model describing the glucose and insulin kinetics.
