**2. Theoretical tools**

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

*Numerical Modeling and Computer Simulation*

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

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

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

[10–12] fields).

with a constant diffusion term.

the problem always complicated.

**92**

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 equations of the Itô type [30]:

$$d\mathbf{Y}\_t^i = \mu\left(\mathbf{Y}\_t^i, t, \theta, b^i\right)dt + \Sigma\left(\mathbf{Y}\_t^i, \theta, b^i\right)dW\_t^i, \quad \mathbf{Y}\_0^i = \mathbf{y}\_0^i, \quad i = 1, \ldots, M,\tag{1}$$

where *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 moment *t*≥*t i* 0, and *Y<sup>i</sup>* <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 *Yi t* � � *<sup>t</sup>*≥<sup>0</sup>, *i* ¼ 1*::M* n o is assumed to verify the same model structure Eq. (1) according to the individual parameters *b<sup>i</sup>* ; *θ* ∈ Θ ⊂ *<sup>p</sup>* is a p-dimensional fixed effects parameter which represents the same and common characteristics for all subjects, and *bi* ∈ *B*⊆ *<sup>q</sup>* are the q-dimensional individual random parameters assumed mutually independent, also called random effects because they are not the same for all individuals; they change between them according to a distribution of density *PB b<sup>i</sup>* <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 simultaneously. Each component *b<sup>i</sup> <sup>l</sup>* may follow a different distribution, ð Þ *l* ¼ 1, … , *q* , and a standard choice for the joint density function *PB bi* <sup>j</sup><sup>Ψ</sup> � � of the vector *bi* could be the Gaussian distribution; however, any other distributions may be considered continuous or discrete:

$$b^i \sim i.i.d.\ N(\theta, \phi). \tag{2}$$

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 of the components *b<sup>i</sup> l* , 1ð Þ ≤*l* ≤*q* ; the components of *ϕ* and *ϑ* represent the population parameters Ψ. For *Y<sup>i</sup>* 0, it is not necessarily known, and when its components are 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. Also, we assume that the distribution of *Y<sup>i</sup>* ð Þ*<sup>t</sup>* given *<sup>b</sup><sup>i</sup>* , *θ* � � and *Y<sup>i</sup> t* <sup>0</sup> ð Þ¼ *yt*0, *t* <sup>0</sup> <*t*, has a strictly positive density with regard to the Lebesgue measure on E:

$$\forall y \to P\_Y(y, t - t'|\mathcal{Y}\_{t'}, b^i, \theta) > 0, \ y \in E. \tag{3}$$

*W<sup>i</sup>* ð Þ*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 Σð Þ� : *E* � Θ � *B* ! <sup>þ</sup> represent, respectively, the drift and the diffusion term of the model and are assumed to have some properties sufficiently regular to ensure a unique solution to the model [33].

In this work, we propose an approximate method to estimate the transition density for a time-inhomogeneous NLME model with SDEs in a closed approximate form. So, we suggest to derive the transition density by solving the motion equations of the process using the Fokker-Planck equation, and we deal with the use of its solution given in [34]. Then, using the expression obtained, we get a closed form approximation of the likelihood function that we maximize to obtain the approximated MLEs ^*θ* and Ψ^. The approximated transition density obtained from the

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

� � that we substitute in Eq. (4)

<sup>j</sup><sup>Ψ</sup> � �*db<sup>i</sup>*

, the likelihood function Eq. (4)

<sup>2</sup> *log* j�*<sup>H</sup>* <sup>~</sup>

*b i* <sup>j</sup>*θ*, <sup>Ψ</sup>Þj � i,

*b i* jΨ � � (8)

<sup>j</sup>*θ*, <sup>Ψ</sup> � �:

*iT* (9)

*b i* jΨ (7)

<sup>2</sup> *log* ð Þ� <sup>2</sup>*<sup>π</sup>* <sup>1</sup>

j ~ *b i* , *θ* � � <sup>þ</sup> *logPB* <sup>~</sup>

*:* (6)

*<sup>Y</sup> Y<sup>i</sup> j* , Δ*<sup>i</sup> j* j*Yi <sup>j</sup>*�1, *bi* , *θ*

> *P*ð Þ *<sup>a</sup> <sup>Y</sup> y<sup>i</sup> j* , Δ*<sup>i</sup> j* j*yi <sup>j</sup>*�1, *bi* , *θ*

For a multidimensional vector of random parameters, if the exact transition density or its closed form approximation is available, we can use the Laplace approximation method [16, 20, 35] to obtain an explicit expression of the likelihood function to maximize, despite the fact that the integral in Eq. (4) has no closed

> *b i* <sup>j</sup>ΨÞ þ *<sup>q</sup>*

<sup>j</sup>*θ*, <sup>Ψ</sup> � � <sup>¼</sup> *logPY <sup>y</sup><sup>i</sup>*

j ~ *b i* , *θ* � � <sup>þ</sup> *logPB* <sup>~</sup>

**2.2 Approximate transition density for multidimensional and nonlinear SDME**

Usually, a formal general solution of the stochastic differential equations as in Eq. (1) cannot be given, which makes the calculation of the transition density for this process more complicated. Moreover, this process has a lot of fluctuations that its exact position cannot be determined but can be known given a region by its probability density; with the FP equation, such a probability density can be determined. The FP equation is a differential equation for the distribution function describing a Brownian motion by which the probability density of the stochastic process can be calculated in a much simpler way by solving this equation. This motion equation is usually used for variables describing a macroscopic but small

*∂*~ *b i ∂*~ *b*

h i � �

to obtain the following approximated likelihood function for the SDME model

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

proposed method is denoted by *P*ð Þ *<sup>a</sup>*

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

*2.1.1 Laplace approximation*

can be approximated as:

*M*

*i*¼1

*logL*ð Þ *<sup>θ</sup>*, <sup>Ψ</sup> <sup>≃</sup> <sup>X</sup>

<sup>¼</sup> *argmaxbi f b<sup>i</sup>*

where

**model**

~ *b i*

**95**

*<sup>L</sup>*ð Þ *<sup>a</sup>* ð Þ¼ *<sup>θ</sup>*, <sup>Ψ</sup> <sup>Y</sup>

*M*

<sup>ð</sup> <sup>Y</sup>*ni*

*j*¼1

*i*¼1

solution. So, for a q-dimensional random vector *b<sup>i</sup>*

<sup>j</sup>*θ*, <sup>Ψ</sup> � � � � *and f bi*

and <sup>∣</sup> � <sup>∣</sup> denotes the determinant of the Hessian matrix *H b<sup>i</sup>*

*<sup>∂</sup>*<sup>2</sup> *logPY yi*

� �

*logPY yi* j ~ *b i* , *θ* � � <sup>þ</sup> *logPB* <sup>~</sup>

*H* ~ *b i* j*θ*, Ψ � � <sup>¼</sup>

Eq. (1):

According to the model Eq. (1), the process Y is the same and follows the same form for each individual of the population, and the SDME model of type Itô expresses the dynamic of the individual i perturbed by the Brownian motion. So, the differences between the subjects are modeled on the one hand by the different realizations of the Brownian motion paths *W<sup>i</sup> t* � � *t*≥*t i* 0 and, on the other hand, by the

incorporation of the random parameters *b<sup>i</sup>* in the model. Therefore, the introduction of parameters varying randomly between subjects allows the model Eq. (1) to explain the variability between individuals.

The goal is to estimate the vector of fixed parameters *θ* and the parameter vector Ψ characterizing the distribution of random parameters *b<sup>i</sup>* s, but the statistical inference for such models is a difficult issue, and the level of difficulty is not the same whether the transition density is explicit or not and whether the process is observed directly or with measurement noise; in this work we assume that the process was exactly observed and no observation noise was considered.

#### **2.1 Maximum likelihood estimation in SDME model**

The likelihood function of an SDME model is expressed as follows:

$$L(\theta, \Psi) = \prod\_{i=1}^{M} P\left(\underline{\boldsymbol{\nu}}^{i} | \theta, \Psi\right) = \prod\_{i=1}^{M} \int P\_{\underline{\mathbf{\varnu}}}\left(\underline{\boldsymbol{\nu}}^{i} | b^{i}, \theta\right) P\_{B}\left(b^{i} | \Psi\right) db^{i} \tag{4}$$

with

$$P\_{\underline{Y}}\left(\underline{y}^{i}|b^{i},\theta\right) = \prod\_{j=1}^{n\_{i}} P\_{Y}\left(y\_{j}^{i},\Delta\_{j}^{i}|y\_{j-1}^{i},b^{i},\theta\right),\tag{5}$$

where *ni* is the number of observations for the subject i at discrete points of time *t i* 0, *t i* 1, … , *t i ni* n o, *<sup>i</sup>* <sup>¼</sup> 1, … , *<sup>M</sup>* and <sup>Δ</sup>*<sup>i</sup> <sup>j</sup>* ¼ *t i <sup>j</sup>* � *t i <sup>j</sup>*�1, *<sup>j</sup>* <sup>¼</sup> 1, … , *ni*. The conditional density *PY y<sup>i</sup>* j� � � is equal to the product of the transition densities Eq. (5) for given random effects and *θ*, but the availability of the explicit transition density is the second constraint for the statistical issue of model Eq. (1) to obtain an exact likelihood function and exact estimators, since computing the transition density is not always obvious and requires approximation methods. However, there are some cases where the exact likelihood function is known, and the exact MLEs of *θ* are obtained (see references in the introduction). In fact, to compute the likelihood function in a closed form of an SDME model, we can encounter two types of problems that require approximate methods to overcome them: First, when the transition density *PY yi j* , Δ*<sup>i</sup> j* j*yi <sup>j</sup>*�<sup>1</sup>, *bi* , *θ* � � is known but the integral in Eq. (4) has no solution, in this case, the numerical methods of approximation of the integral are required. Or, second, when *PY y<sup>i</sup> j* , Δ*<sup>i</sup> j* j*yi <sup>j</sup>*�<sup>1</sup>, *bi* , *θ* � � cannot even be expressed explicitly and must also be approximated, see next paragraphs. Usually, in realistic examples, we have both an unknown transition density and an integral that is difficult to solve analytically. In theory, several methods for approximating transition densities and integrals have been proposed (see references cited in the introduction).

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

In this work, we propose an approximate method to estimate the transition density for a time-inhomogeneous NLME model with SDEs in a closed approximate form. So, we suggest to derive the transition density by solving the motion equations of the process using the Fokker-Planck equation, and we deal with the use of its solution given in [34]. Then, using the expression obtained, we get a closed form approximation of the likelihood function that we maximize to obtain the approximated MLEs ^*θ* and Ψ^. The approximated transition density obtained from the proposed method is denoted by *P*ð Þ *<sup>a</sup> <sup>Y</sup> Y<sup>i</sup> j* , Δ*<sup>i</sup> j* j*Yi <sup>j</sup>*�1, *bi* , *θ* � � that we substitute in Eq. (4) to obtain the following approximated likelihood function for the SDME model Eq. (1):

$$L^{(a)}(\theta,\Psi) = \prod\_{i=1}^{M} \left[ \left( \prod\_{j=1}^{n\_i} P\_Y^{(a)} \left( y\_j^i, \Delta\_j^i | y\_{j-1}^i, b^i, \theta \right) \right) P\_B \left( b^i | \Psi \right) db^i. \tag{6} \right]$$

## *2.1.1 Laplace approximation*

Σð Þ� : *E* � Θ � *B* ! <sup>þ</sup> represent, respectively, the drift and the diffusion term of the model and are assumed to have some properties sufficiently regular to ensure a

incorporation of the random parameters *b<sup>i</sup>* in the model. Therefore, the introduction of parameters varying randomly between subjects allows the model Eq. (1) to

inference for such models is a difficult issue, and the level of difficulty is not the same whether the transition density is explicit or not and whether the process is observed directly or with measurement noise; in this work we assume that the process was exactly observed and no observation noise was considered.

The likelihood function of an SDME model is expressed as follows:

<sup>¼</sup> <sup>Y</sup>*ni j*¼1

*<sup>j</sup>* ¼ *t i <sup>j</sup>* � *t i*

random effects and *θ*, but the availability of the explicit transition density is the second constraint for the statistical issue of model Eq. (1) to obtain an exact likelihood function and exact estimators, since computing the transition density is not always obvious and requires approximation methods. However, there are some cases where the exact likelihood function is known, and the exact MLEs of *θ* are obtained (see references in the introduction). In fact, to compute the likelihood function in a closed form of an SDME model, we can encounter two types of problems that require approximate methods to overcome them: First, when the

solution, in this case, the numerical methods of approximation of the integral are

and must also be approximated, see next paragraphs. Usually, in realistic examples, we have both an unknown transition density and an integral that is difficult to solve analytically. In theory, several methods for approximating transition densities and

� �

<sup>¼</sup> <sup>Y</sup> *M*

*i*¼1

*PY y<sup>i</sup> j* , Δ*<sup>i</sup> j* j*yi <sup>j</sup>*�<sup>1</sup>, *bi* , *θ*

where *ni* is the number of observations for the subject i at discrete points of time

is equal to the product of the transition densities Eq. (5) for given

� �

ð *PY y<sup>i</sup>* j*bi* , *θ* � �

The goal is to estimate the vector of fixed parameters *θ* and the parameter vector

form for each individual of the population, and the SDME model of type Itô expresses the dynamic of the individual i perturbed by the Brownian motion. So, the differences between the subjects are modeled on the one hand by the different

According to the model Eq. (1), the process Y is the same and follows the same

*t* � � *t*≥*t i* 0

and, on the other hand, by the

s, but the statistical

<sup>j</sup><sup>Ψ</sup> � �*db<sup>i</sup>* (4)

, (5)

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

*<sup>j</sup>*�1, *<sup>j</sup>* <sup>¼</sup> 1, … , *ni*. The conditional den-

is known but the integral in Eq. (4) has no

cannot even be expressed explicitly

unique solution to the model [33].

*Numerical Modeling and Computer Simulation*

realizations of the Brownian motion paths *W<sup>i</sup>*

explain the variability between individuals.

*<sup>L</sup>*ð Þ¼ *<sup>θ</sup>*, <sup>Ψ</sup> <sup>Y</sup>

with

n o

j� � �

transition density *PY yi*

required. Or, second, when *PY y<sup>i</sup>*

sity *PY y<sup>i</sup>*

*t i* 0, *t i* 1, … , *t i ni*

**94**

Ψ characterizing the distribution of random parameters *b<sup>i</sup>*

**2.1 Maximum likelihood estimation in SDME model**

*M*

*P yi* j*θ*, Ψ � �

*i*¼1

*PY yi* j*bi* , *θ* � �

, *<sup>i</sup>* <sup>¼</sup> 1, … , *<sup>M</sup>* and <sup>Δ</sup>*<sup>i</sup>*

*j* , Δ*<sup>i</sup> j* j*yi <sup>j</sup>*�<sup>1</sup>, *bi* , *θ*

� �

*j* , Δ*<sup>i</sup> j* j*yi <sup>j</sup>*�<sup>1</sup>, *bi* , *θ*

integrals have been proposed (see references cited in the introduction).

For a multidimensional vector of random parameters, if the exact transition density or its closed form approximation is available, we can use the Laplace approximation method [16, 20, 35] to obtain an explicit expression of the likelihood function to maximize, despite the fact that the integral in Eq. (4) has no closed solution. So, for a q-dimensional random vector *b<sup>i</sup>* , the likelihood function Eq. (4) can be approximated as:

$$\log \mathcal{L}(\boldsymbol{\theta}, \boldsymbol{\Psi}) \simeq \sum\_{i=1}^{M} \left[ \log \boldsymbol{\mathcal{P}}\_{\underline{Y}} \left( \underline{\boldsymbol{y}}^{i} | \boldsymbol{\bar{b}}^{i}, \boldsymbol{\theta} \right) + \log \boldsymbol{\mathcal{P}}\_{\boldsymbol{\theta}} \left( \boldsymbol{\bar{b}}^{i} | \boldsymbol{\Psi} \right) + \frac{\boldsymbol{q}}{2} \log \left( 2\pi \right) - \frac{1}{2} \log \left| -H \left( \boldsymbol{\bar{b}}^{i} | \boldsymbol{\theta}, \boldsymbol{\Psi} \right) \right| \right], \tag{7}$$

where

$$\bar{b}^i = \arg\max\_{b^i} \left( f\left( b^i | \theta, \Psi \right) \right) \quad \text{and} \quad f\left( b^i | \theta, \Psi \right) = \log \underline{P}\_{\underline{Y}} \left( \underline{y}^i | \bar{b}^i, \theta \right) + \log P\_B \left( \bar{b}^i | \Psi \right) \tag{8}$$

and <sup>∣</sup> � <sup>∣</sup> denotes the determinant of the Hessian matrix *H b<sup>i</sup>* <sup>j</sup>*θ*, <sup>Ψ</sup> � �:

$$H\left(\tilde{b}^i|\theta,\Psi\right) = \frac{\partial^2 \left[\log P\_{\underline{Y}}\left(\underline{\boldsymbol{y}}^i|\bar{\boldsymbol{b}}^i,\theta\right) + \log P\_B\left(\bar{\boldsymbol{b}}^i|\Psi\right)\right]}{\partial \bar{\boldsymbol{b}}^i \partial \bar{\boldsymbol{b}}^{iT}}\tag{9}$$
