Open access peer-reviewed chapter

Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential Equations Using the Fokker-Planck Equation Solution

By Bakrim Fadwa, Hamid El Maroufy and Hassan Ait Mousse

Submitted: September 16th 2019Reviewed: December 4th 2019Published: February 19th 2020

DOI: 10.5772/intechopen.90751

Downloaded: 52

Abstract

This chapter is concerned with estimation method for multidimensional and nonlinear dynamical models including stochastic differential equations containing random effects (random parameters). This type of model has proved useful for describing continuous random processes, for distinguishing intra- and interindividual variability as well as for accounting for uncertainty in the dynamic model itself. Pharmacokinetic/pharmacodynamic modeling often involves repeated measurements on a series of experimental units, and random effects are incorporated into the model to simulate the individual behavior in the entire population. Unfortunately, the estimation of this kind of models could involve some difficulties, because in most cases, the transition density of the diffusion process given the random effects is not available. In this work, we focus on the approximation of the transition density of a such process in a closed form in order to obtain parameter estimates in this kind of model, using the Fokker-Planck equation and the Risken approximation. In addition, the chapter discusses a simulation study using Markov Chain Monte Carlo simulation, to provide results of the proposed methodology and to illustrate an application of mixed effects models with SDEs in the epidemiology using the minimal model describing glucose-insulin kinetics.

Keywords

  • stochastic differential equations
  • mixed effects model
  • Fokker-Planck equation
  • transition density
  • maximum likelihood estimators
  • genetic algorithm
  • Markov Chain Monte Carlo simulation

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, 3, 4, 5, 6] or stochastic and linear or nonlinear models (see many applications of stochastic NLME models in biomedical [7, 8, 9] and in pharmacokinetic [10, 11, 12] fields).

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, 25, 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-Carlo-integrate the transition density [27, 28, 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 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.

2. Theoretical tools

Consider an N-dimensional continuous and stochastic process Ytin the state space ERNdescribed by the general first-order nonlinear stochastic differential equations of the Itô type [30]:

dYti=μYtitθbidt+ΣYtiθbidWti,Y0i=y0i,i=1,,M,E1

where Ytiis defined as the solution of the SDME model Eq. (1) that exists under some conditions that we suppose satisfied [31, 32, 33] and represents the observation of individual i from M different experimental units, (i=1,..,M), at the moment tt0i, and Y0i=Yt0iis the initial state of Ytfor each subject. The process Ytit0i=1..Mis assumed to verify the same model structure Eq. (1) according to the individual parameters bi; θΘRpis a p-dimensional fixed effects parameter which represents the same and common characteristics for all subjects, and biBRqare 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 PBbiΨdepending on a population parameter Ψ; in the population approach, this parameter vector allows for a data from several subjects to be considered simultaneously. Each component blimay follow a different distribution, l=1q, and a standard choice for the joint density function PBbiΨof the vector bicould be the Gaussian distribution; however, any other distributions may be considered continuous or discrete:

bii.i.dNϑϕ.E2

The joint density function of the vector biis parameterized by a q-dimensional parameter ϑυRqand a q×q-dimensional matrix ϕΦRq×qrepresenting the covariance matrix of biand specifying the parameters of the marginal distributions of the components bli, 1lq; the components of ϕand ϑrepresent the population parameters Ψ. For Y0i, 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 Yitgiven biθand Yit=yt, t<t, has a strictly positive density with regard to the Lebesgue measure on E:

yPYyttytbiθ>0,yE.E3

Witare the standard Brownian motions, and they are assumed mutually independent with bjfor all 1i,jM. The functions μ:E×R×Θ×BRand Σ:E×Θ×BR+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].

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 Wtitt0iand, on the other hand, by the incorporation of the random parameters biin 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 bis, 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θΨ=i=1MPy¯iθΨ=i=1MPY¯y¯ibiθPBbiΨdbiE4

with

PY¯y¯ibiθ=j=1niPYyjiΔjiyj1ibiθ,E5

where niis the number of observations for the subject i at discrete points of time t0it1itnii,i=1,,Mand Δji=tjitj1i,j=1,,ni. The conditional density PY¯y¯iis 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 PYyjiΔjiyj1ibiθ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 PYyjiΔjiyj1ibiθ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).

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 PYaYjiΔjiYj1ibiθthat we substitute in Eq. (4) to obtain the following approximated likelihood function for the SDME model Eq. (1):

LaθΨ=i=1Mj=1niPYayjiΔjiyj1ibiθPBbiΨdbi.E6

2.1.1 Laplace approximation

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 bi, the likelihood function Eq. (4) can be approximated as:

logLθΨi=1MlogPY¯y¯ib˜iθ+logPB(b˜iΨ)+q2log2π12logH(b˜iθΨ),E7

where

b˜i=argmaxbifbiθΨandfbiθΨ=logPY¯y¯ib˜iθ+logPBb˜iΨE8

and denotes the determinant of the Hessian matrix HbiθΨ:

Hb˜iθΨ=2logPY¯y¯ib˜iθ+logPBb˜iΨb˜ib˜iTE9

2.2 Approximate transition density for multidimensional and nonlinear SDME model

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 system, where the fluctuations are important as for some cases in physics, e.g., the position and the speed of the Brownian motion of a small particle. However, it can be also used for the larger system where, in spite of their small fluctuations, the stochastic description remains necessary when the deterministic equations may not be stable for this type of system.

So, under some assumptions, the transition density PYyjiΔjiyj1ibiθof the process Yin the SDME model Eq. (1) is the solution of the following functional partial differential Equation [23, 34]:

PYyjiΔjiyj1ibiθt=LFPPYyjiΔjiyj1ibiθ,E10

where LFPis the FP operator and LFP=k=1NYkμkYitθbi+12k=1Nl=1N2YkYlΣkland μkis the k-th element of the drift vector and Σklis the kl-th element of the diffusion matrix. We assume that the diffusion matrix is positive definite so that its inverse exists DetΣ0.

In [34], H. Risken deals in this book with the derivation of the FP equation and its solution methods with some of its applications especially for problems of Brownian motion. So, here, we propose to characterize the transition density of the N-dimensional process Ytiin Eq. (1) using the Risken approximation based on the formal solution of Eq. (10) proposed in [34], p: 4.109, denoted by PYayjiΔjiyj1ibiθ:

PYaYjiΔjiYj1ibiθ=2ΠΔjNDetΣ12exp(14Δjl=1Nk=1NΣ1lkYjilYj1ilμlYj1itθbiΔjiYjikYj1ikμkYj1itj1θbiΔji)E11

To test the effectiveness of this approach, we will guide our statistical methodology by simulation studies in order to examine the flexibility of its application to deduce its advantages and disadvantages. So, we substitute PYaYjiΔjiYj1ibiθin Eq. (6), and by solving the integral with respect to the random effects density and maximizing the obtained likelihood function, we get the approximated estimators θ̂and Ψ̂.

2.2.1 Approximated estimators

For a nonlinear SDME model with Gaussian random effects using Eqs. (6), (7), and (11), we obtain the following approximated likelihood function:

logLaθΨi=1M[q2log2π12logdetϕni2logdetΣ+(j=1nilog2ΠΔjN14ΔjΣ1lkYjilYj1ilμlYj1itθb˜iΔjiYjikYj1ikμkYj1itj1θb˜iΔji)12b˜iνϕ1b˜iν+q2log2π12logdetHb˜iθΨ].E12

The MLEs of θΨcan be obtained using one of the optimization tools and numerical computation software, especially when it is complicated to compute the gradients of the likelihood analytically. Here, we propose to use the genetic algorithm as an optimization tool to maximize the approximate likelihood function Eq. (12) using MATLAB software:

θ̂Ψ̂=argminGAlogLaθΨ.E13

2.3 Genetic algorithm

The genetic algorithm is a random search technique to look for an exact or approximated optimum points for optimization problems [36, 37, 38]. It is based on the concepts of natural genetic evolution which contains the following stages: the reproduction, the crossing, and the mutation of a constantly evolving population. It sets up the evolution of a random population of potential solutions of the N cardinal; then, the N simultaneous iterative trajectories interact with each other by following or imitating the biological evolution, for a convergence of certain elements of the population towards an optimal point of the fitness function.

The GA can search in multiple directions to explore all the search space by the possibility of jumping across them, so that the seeds spread uniformly over the whole search space. In this algorithm, we have a diversity of initial populations which gives the global optimum faster than other algorithms, where the initial value is very important and should be enough close to the global optimum. All of these features allows the GA to be regarded as a driving tool of evolution giving good results for optimization processes [37, 39]. In the literature, there were many works about the application of GA in optimizing problems specially for likelihood function [40, 41]. For the use of GA, we must first define some parameters of the algorithm: Population size N, EN, SR, CP, MP, fitness function, and convergence criteria. In the following we present the GA steps:

Steps of GA:

1. Generate initial population β10β20βN0,m=0via an initialization strategy (random generation), in our case β=θΨ.

For m=0:

2. Evaluate the Fitness function logLaθmΨm.

3. While (convergence criteria are not satisfied):

Do:

4. Replacement step (by using SR and EN): At the SR rate, individuals with the worst results in step 2 of fitness function are replaced by new ones randomly generated, and a number EN of individuals is selected and accepted for the next step.

5. Selection operator by using roulette wheel method, based on the fact that the more the individual has a good result of fitness function, the more likely he will be selected.

6. Crossover operator by using CP and mutation operator by using MP: It is a mechanism of perturbation on the candidate individuals (parents) according to CP and MP to generate new groups of individuals, and we obtain a new m+1ndpopulation β1m+1β2m+1βNm+1.

Else:

7. Evolution stops; get GA output.

8. m=m+1.

End For.

In the simulation studies paragraphs, the GA is implemented using the MATLAB software, where the function “ga,” to generate the genetic algorithm, requires inputs that are chosen according to the constraints of each example (see the help window in MATLAB). Moreover, the algorithm parameters are chosen according to what is early used in the literature [39], EN=4,MP=0.2,CP=0.8andSR=1/3, and the search spaces are around the confidence interval of the minimal model parameters (see [42] and references therein).

3. Simulation studies for SDME model

3.1 The two-dimensional Ornstein-Uhlenbeck process

To apply the proposed methodology and evaluate its effectiveness, we consider the two-dimensional OU process that is very useful in pharmacokinetic/pharmacodynamic studies and in biology [43], physics, engineering, finance, and neuroscience applications [14, 44]. Indeed, the choice of this process is due to the fact that it is one of the few known multivariate SDME models with known transition density. For this reason, we choose the OU process to evaluate the methodology presented above, and we perform a comparison study between the results obtained using the proposed transition density in Eq. (11) and those obtained using the exact density. The model is defined as follows:

dY1it=β11b11iY1itα1+β12b12iY2itα2dt+Σ11dW1it,Y0i=y01i,i=1,,MdY2it=β21b21iY1itα1+β22b22iY2itα2dt+Σ22dW2it,Y0i=y02i,i=1,,ME14

With Yit=Y1itY2it; β=β11β12β21β22; α=α1α2; Σ=Σ1100Σ22; Wit=W1itW2it; Yi0=Y1i0Y1i0and bi=b11ib12ib21ib22iwhere bllii.i.dΓrllrll1,l,l=1,2;i=1,,M.

We rewrite the system in matrix notation under the Itô formula; we denote by the elementwise multiplication:

dYit=βbiαYitdt+ΣdWit,Y0i=y0i,i=1,,M.E15

Here, the random effects biare a matrix and not a vector in order to have a uniform dimension writing of the Eq. (15) and are assumed mutually independent and independent of Y0iand Wi. The fixed parameter vector is θ=β11β12β21β22α1α2Σ11Σ22, and the population parameter vector is Ψ=r11r12r21r22. The exact transition density of model Eq. (15) for a given realization of the random effects is a bivariate normal:

PYYtjiΔjiYtj1ibiθ=2π1ς12expYtjiμς1Ytjiμ2,E16

with mean vector μ=α+Ytj1iαexpβ.biΔjiand covariance matrix ς=τexpβ.biΔjiτexpβ.biΔji, where

τ=12trβ.biβ.biβ.biΣΣ+β.bitrβ.bi.IΣΣβ.bitrβ.bi.I.E17

We assume that the matrices β.biand Σhave full rank and the real parts of the eigenvalues of β.biare positive definite in order that a stationary solution to Eq. (15) exists. Under these assumptions, we derive from Eqs. (14) and (11) the following approximated transition density of Y:

PYaYtjiΔjiYtj1ibiθ=2ΠΔj2Σ11Σ2212exp(14ΔjΣ111[Ytj1iYtj11i+(β11b11iYtj11iα1+β12b12iYtj12iα2)Δji]2+Σ221Ytj2iYtj12i+β21b21iYtj11iα1+β22b22iYtj12iα2Δji2).E18

In Figure 1, we illustrate in (a) the simulation of the OU process using the Euler scheme [45] with the following set of parameters: (β11=2.8,β12=2.5,β21=1.8,β22=2,α1=0.8,α2=1.5,Σ11=0.3,Σ22=0.5,r11=45,r12=100,r21=100,r22=125), and a time step of Δj=0.001and represent in (b) the graph of the two transition densities given by Eqs. (16) and (18) for Yjto Yj+1using the same set of parameters and time step.

Figure 1.

A sample path of the OU process in the third graph of (a) for the given parameters set with the initial condition: Y0= (3,3) and time interval [3,10]; and the transition density for a transition from Yj to Yj+1 in (b).

3.1.1 Simulation results

In this simulation study, we generate 1000 artificial datasets of dimension 2×n+1×Mfrom Eq. (15), where M is the number of subjects and n presents the number of repetitions of the experiment for each subject; then we estimate the parameters using the proposed approximated method, and we obtain 1000 sets of parameter estimates. The observations are obtained by linear interpolation from simulated trajectories using the Euler-Maruyama scheme with step size equal to 103.

By plugging Eq. (16) in Eq. (4) and Eq. (18) in Eq. (6), we obtain a huge expression of the likelihood function but without a closed solution of integrals, so the exact estimators of θand Ψare unavailable. Therefore, in both cases, either using the exact or the approximated transition density, the Hessian matrix in Laplace approximation can be obtained analytically after a tedious calculation; then we apply the GA to obtain parameter estimates. But we cannot ignore the time consumed by this algorithm to get the results because of the long and complex expressions.

Table 1 shows that, for the given sample size, the results can be correctly identified using this estimated approach; even if some of parameters are overestimated or underestimated, the results remain acceptable because the results belong to the confidence interval. However, we believe that these results could be further proven by using other sample sizes and by adding alternative assumptions for the model that we did not consider in the methodology proposed in this chapter, which could further complicate the problem and be more time-consuming, as well as the present methodology suffers from some limitations. For the random parameters, the estimates can be provided using the optimization algorithm on Eq. (8) using the obtained estimate results of the parameters vector Ψ. Moreover, we conduct this simulation study using Figure 2, which shows that the empirical distribution of the most approximated estimators seems to be reasonably close to a normal distribution:

True valuesMean and (Std) (M=40,n=10)
β11β12β21β22β̂11β̂12β̂21β̂22
2.82.51.823.10–3.252.48–2.561.72–1.652.11–2.37
(0.164)–(0.283)(0.015)–(0.283)(0.095)–(0.189)(0.153)–(0.255)
α1α2Σ11Σ22α̂1α̂2Σ̂11Σ̂22
0.81.50.30.51.06–1.131.57–1.640.28–0.370.56–0.45
(0.071)–(0.102)(0.109)–(0.158)(0.026)–(0.073)(0.023)–(0.061)
r11r12r21r22r̂11r̂12r̂21r̂22
451001002544.75–52.43100–112.75102.35–89.6424.72–31.02
(0.523)–(9.372)(01.166)–(28.113)(01.04)–(22.05)(2.297)–(11.20)

Table 1.

Ornstein-Uhlenbeck model: maximum likelihood estimates from 1000 simulations of model Eq. (14), using the exact and the approximated transition density.

Figure 2.

Empirical distribution of estimates obtained using the exact and approximated transition density.

3.2 Stochastic minimal model

The minimal model describes the glucose-insulin kinetics and the dynamics of these processes illustrating the diabetes disease mechanisms. The diabetes is one of the most prevalent diseases in individuals and the nature and degree of its assignment changes from an individual to another and depends on certain individual characteristics, which implies that the concepts of stochastic modeling with random effects could be a good approach to modeling this disease. The Diabetes may be due to the insufficient insulin production (type 1 diabetes) or to the fact that the cells do not respond to the secreted insulin (type 2 diabetes) and the T2D patients tend to have substantially lower insulin sensitivity than healthy individuals, so that the T2D can be characterized by the level of insulin sensitivity for each individual. Therefore, to model the T2D, we observe how a person’s body responds to insulin in the process of transporting glucose to tissues by measuring his insulin sensitivity. So, we present in this section the estimation of the minimal model which represents a powerful model describing the glucose-insulin kinetics for an individual’s body in three differential equations, see the mathematical formulation of the model in [46, 47, 48, 49]. So, it is already clear that the model will contain both fixed and random effects, because the study of diabetes disease takes into account the response of each individual according to his own parameters and other common parameters that describe the process of glucose-insulin for the entire population. See the description of the glucose-insulin kinetics in Figure 3.

Figure 3.

At first, glucose and insulin concentrations in the blood are described by two sets of differential equations (see [50]); at a rate p1, glucose leaves and enters the glucose space in proportion to the difference between the plasma glucose concentration G(t) and the basal plasma concentration Gb that represent the known pre-injection glucose level for each individual. Therefore, the parameter p1 represents the glucose’s own ability to be eliminated in muscles, liver, and tissues independently of insulin which is called glucose efficiency and denoted by SG. Then, the glucose disappears from the glucose space at a rate proportional to insulin concentration in the insulin compartment Xt which represents the dynamic of insulin response according to the two rates p2 and p3. These two parameters represent, respectively, the decreased glucose absorption capacity in tissues and its increased insulin dependency. For insulin secretion I(t), it is secreted by the pancreas independently of the glucose concentration, and proportional to a rate n to its own level already in the body and to the glucose level deferred from a threshold h at a rate γ when Gt is above h, the insulin secretion does not only depend on the hyperglycemia level but also to the time spent since glucose injection.

From the mentioned literature and Figure 3, the glucose-insulin disposal can be represented, with respect to time, by the following nonlinear stochastic differential equations, perturbed by the stochastic terms σ1dw1t, σ2dw2t, and σ3dw3t:

dGt=p1+XtGt+p1Gbdt+σ1dw1t,G0=G0dXt=p2Xt+p3ItIbdt+σ2dw2t,X0=0dIt=nItIb+γGthtdt+σ3dw3t,I0=I0,E19

where Gtand Itare, respectively, the concentration of glucose and insulin at time t in the blood. Gband Ibindicate the basal level of glucose and insulin concentration before the glucose injection, this injection will cause a disturbance of the concentrations according to the mechanism described in these equations, and these values are assumed known for each individual. And, G0and I0are the theoretical measure of the concentrations at glucose injection moment at the beginning of the experiment.

The insulin sensitivity is defined by combining the two rates p3and p2as SI=p3/p2, representing insulin’s ability to increase the net glucose utilization [51].

Finally, the stochastic minimal model under the Itô sense re-parameterized by SGand SIcan be rewritten as:

dYit=SGi+XtiGti+SGiGbip2Xti+SIiItiIbinItiIbi+γGtihtdt+ΣdWt;Y0i=G0i0I0i,E20

where Yit=GtiXtiItiand Σis the diagonal diffusion matrix which contains the unknown constants σ1, σ2, and σ3. The parameters SGi,p2,SIi,n,γ,h,G0iand I0iare unknown in the model and will be estimated. The parameters SGi,SIi,I0iand G0iare assumed random, because they represent individual parameters that change from an individual to another. Each person has its own insulin sensitivity SIiwhich allows to know if the cells of his body react correctly or not to the insulin and if the insulin produced by the pancreas is sufficient or not, which can make some people with T2D and others without diabetes. Also, for glucose effectiveness SGi, which represents the glucose’s own ability to be eliminated independently of insulin, it is unique to each individual and changes from a person to another, as well as for the measurement of glucose and insulin concentration. For the rest of the parameters, we consider them fixed since they describe the common side of the glucose-insulin process for the entire population (all individuals). So, we have the following random effect vector bi=SGiSIiI0iG0i, and we assume that:

SGiNμSGσSG,SIiNμSIσSI,I0iNμI0σI0,G0iNμG0σG0.E21

Random effects are assumed to be independent with a multinormal joint density function, with the mean ϑ=μSGμSIμI0μG0and the covariance matrix ϕ=diagσSGσSIσI0σG0, so we have Ψ=μSGμSIμI0μG0σSGσSIσI0σG0and θ=p2nγhσ1σ2σ3.

So, here we deal with a time-inhomogeneous NLME model with SDE describing the glucose-insulin kinetics (see [52] for the implementation of SDE time-inhomogeneous model) and [53] where the maximum likelihood estimation for a time-inhomogeneous stochastic differential model of glucose dynamics was treated. The measurements in the model Eq. (20) are assumed directly observed without measurement errors as well as in the theoretical approach presented above.

So, we wish to estimate θΨgiven the observations y¯=y1¯yM¯from model Eq. (20). By using the approximated transition density (Eq. (11)), we get the following approximated likelihood function for model Eq. (20):

LaθΨ=i=1M2π2σSGσSIσI0σG012j=1ni2ΠΔj3σ1σ2σ312)R4exp(j=1ni[14Δj[1σ1A1ji2+1σ2A2ji2+1σ3A3ji2)]12σSG1SGiμSG2σSI1SIiμSI2σI01I0iμI02σG01G0iμG02)dSGidSIidI0idG0iE22

with

Alji=YjilYj1ilμlYj1itj1θbiΔji.E23

We have no closed form solution to this integral, so exact estimators of θand Ψare unavailable. So, we use the Laplace approximation method described in Section 2 to obtain a closed form approximation to the log-likelihood function logLaθΨfor the model Eq. (20); then by applying the GA, we get the approximate estimates θ̂and Ψ̂.

3.2.1 Simulation results

We start this study with an application on artificial data generated on the intravenous glucose tolerance test principle (see [53] for a mathematical modeling of the test where glucose and insulin concentrations in plasma are subsequently sampled after an intravenous glucose injection). We generate 5000 sets of simulated artificial data of dimension ni+1×Mfrom Eq. (20) using the Euler-Maruyama scheme [45] with a step size of 103and a set of true parameters that are chosen according to [53, 54] representing the normal range of parameters values to simulate healthy subjects (without diabetes), with M being the number of units and nibeing the number of observations or repeated measurements collected on each unit i.

For each data from 5000 generated data sets, we estimate θΨby applying the proposed method. We first assume that the number of repeated measurements collected on each unit is constant ni=mand Δji=Δfor all 1iMand 1jni; then we apply the GA after choosing the good algorithm parameters (N, EN, SR, CP, MP), and we get 5000 sets of parameters estimates. We repeat this for large and small data with different possibilities of repetition of the experiment mM=4060;4010and 1020; for each parameter the sample mean and standard deviation are reported in Table 2. The simulation study on small data is treated in order to see if the size of the sample influences on the results and if the number of measures taken over time has a negligible effect or not, in other words, to see if it is possible to select only the essential measuring moments without repeating the measurements several times to well simulate a subject. We treat this issue in relation to our model and its study context, since in epidemiology the availability of data (measurements) at any point of time is an interesting constraint. We note that quantities of Gband Ibare randomly simulated from the normal range of healthy subjects.

True values(M,m)Mean and (Std)
p2nγhp̂2n̂γ̂ĥ
0.0740.100.000790(40,60)0.0737 (0.0014)0.100 (0.0013)0.00073 (2.02 × 10−4)89.08 (0.031)
(40,10)0.0784 (0.0059)0.209 (0.0295)0.00068 (2.54 × 10−4)91.34 (0.712)
(10,20)0.0794 (0.0402)0.156 (0.046)0.00054 (0.0016)62.11 (1.13)
σ1σ2σ3σSGσ̂1σ̂2σ̂3σ̂SG
0.010.060.030.006(40,60)0.009 (0.0012)0.0616 (0.0011)0.0343 (0.0027)0.0061 (0.00010)
(40,10)0.014 (0.0023)0.0708 (0.0060)0.0340 (0.0035)0.0073 (0.00030)
(10,20)0.015 (0.0031)0.0938 (1.0196)0.0480 (0.0108)0.0067 (0.0006)
σSIσI0σG0μSGσ̂SIσ̂I0σ̂G0μ̂SG
0.00002546500.03(40,60)0.000021 (3 × 10−6)44.98 (1.25)45.96 (2.11)0.0315 (0.0007)
(40,10)0.000029 (4.4 × 10−6)43.94 (1.75)45.16 (2.53)0.0349 (0.0036)
(10,20)0.000016 (0.8 × 10−5)41.09 (2.82)44.64 (3.07)0.0178 (0.0051)
μSIμI0μG0μ̂SIμ̂I0μ̂G0
0.000295320(40,60)0.00021 (1.2 × 10−6)94.20 (1.12)321.2 (1.07)
(40,10)0.00025 (1.6 × 10−6)92.05 (2.25)318.6 (3.11)
(10,20)0.00037 (1.35 × 10−5)122.51 (2.51)281.5 (4.88)

Table 2.

Approximated maximum likelihood estimates and standard deviation from simulations of model Eq. (20), using the approximated transition density Eq. (11) with large and small DATA.

In Table 2, we report the results obtained on large and small data by maximizing Eq. (22) using the GA. We notice that, for the case of the large data, the true values of the parameters are well identified with the exception of some in which, in all the simulations and samples, the true value does not belong to the estimated confidence interval, such as σSI, σI0, and σG0; nevertheless, the results are more satisfactory when the sample size M is large for all parameters. As mentioned before, we cannot ignore the limitations of this method, and the results could be improved by eliminating these limits and improving this approach, but for the given assumptions and tools, we can say that the results are still satisfactory even for a sample of M=10with 20measures taken on each subject. So, we can conclude that we can rely on the results obtained from small samples with a small number of measurement repetitions of at least 20observations. In addition, although we can use a relatively small number of measurements, we need to know how to choose the time points to perform the measurements in the blood, as this could, physiologically, affect selected observations and results. Thus, it is specified here that the essential task is to know how to choose the good moments of measurement after the injection, even in small numbers, chosen according to medical knowledge. Figure 4 shows that, in the case of Mm=4060, the empirical distribution of the most approximated estimates seems to be reasonably close to a normal distribution.

Figure 4.

SMM: empirical distribution of population parameter estimates obtained using (18) for (M, m) = (40,60).

Finally, from this simulation study where we have considered two SDME models, we can conclude that the parameters values of the models appear to be correctly identified using the proposed approach based on the Risken approximation to approximate the transition density of a stochastic process.

4. Conclusions

In this work, we have proposed a procedure to estimate the parameters of a mixed effects model containing stochastic differential equations, known by the SDME models, by proposing an approach to approximate its likelihood function to obtain the MLEs. This method has been evaluated by simulation studies on two SDME models in epidemiology: the two-dimensional Ornstein-Uhlenbeck process and stochastic minimal model. In fact, in models with SDEs instead of ODEs with random effects, the estimation of parameters is still not obvious even for one individual (one trajectory) because of the difficulties in deriving the transition densities, and difficulties become more interesting in using the population approach that treats the entire population simultaneously. The derivation of the exact density is not always possible for a stochastic and continuous process in an SDME model, so the search for an approximation of this density is an important step and requires an expensive calculation. This task is very interesting to give good results with good statistical properties of the estimators obtained by maximizing the likelihood function. In this work, the proposed estimation method has been applied to multidimensional and nonlinear SDME models with many random parameters normally distributed that can be extended to random parameters of any distribution. So, an approximation of the transition density PYais obtained in a closed form using the Risken approximation for the formal solution of the Fokker-Planck equation proposed by Risken [34], and then the approximated likelihood is obtained using the Laplace approximation method and optimized using the genetic algorithm; these calculation procedures can be obtained using any numerical calculation software or with symbolic computing capabilities.

The classical inference of SDME models implies the problem of the numerical evaluation of the integral for the given random effects in the likelihood function, which becomes complicated especially when the model contains more than two random parameters. In the literature, several methods have been proposed and tested for the approximation of the integral (see references in the introduction) and the following examples: [8] which proposes the Gaussian quadrature method to solve the integrals for the case of an SDME models with a single random effect and [20] where the study was revised for a general case with several random parameters using the Laplace approximation to compute the integral in Eq. (4) and Eq. (6) numerically. For the mixed effects framework, see [3, 16, 55, 56]. In the case of using the Laplace method, as in this chapter, the calculation of the Hessian matrix can be done analytically when it is possible, as the examples in Section 2, or with the help of a symbolic calculus software or the automatic differentiation (AD) tools [57].

The results of simulation studies are satisfactory and can be obtained either by using moderate values for the number of experimental units M and of observations n taken for each experimental unit or by using a small sample size but with a number of measurements taken for each subject of 10 at least; this is relevant for applications where large sets of data are not available, such as biomedical applications where the mixed effect theory is widely applied.

The advantages of this approach, compared to those proposed in the literature for multidimensional SDME models with more than one random parameter [26], are that the computation of the approximate density is very easy and does not require a lot of time to calculate it or to program it in a software; the only task that can be time-consuming is in the optimization step to search for the optimum solution of the likelihood function, and also the proposed method is effective even with large data with a MATLAB program on a common PC (Intel Pentium IV 3.0 GHz with 512 MB of RAM). Nevertheless, the method suffers some limitations, for example, when the conditions to use Eq. (11) are not verified when, e.g., the inverse of the diffusion term does not exist and when, in certain cases, it is not obvious to derive the gradients and Hessians terms. Another limitation is that measurement errors are not considered in this work, and for a good stochastic version, it will be better to include noise on process increments and noise on observations that may be significant compared to system noise. These limitations may provide a perspective towards a more elaborate extension of the statistical study for SDME models, particularly in the field of epidemiology.

In conclusion, in this work, we proposed a method of parameter estimation for a mixed effects models with SDEs proposing an approximation method for the transition density in the case when it cannot be obtained in a closed form, with an approximation method of integral computation for the case of an SDME model with several random parameters. We have treated two examples to illustrate the effectiveness of this approach using computer tools. Indeed, we believe that this type of models is very interesting and provides a powerful and flexible modeling approach for repeated measurement studies such as biological and pharmacokinetic/pharmacodynamic and financial studies, as they combine the good characteristics of mixed effects and stochastic increments in intra-subject dynamics for a good modeling of a phenomenon.

Acknowledgments

Z.T., H.E., and F.B conceived the presented idea and contributed to the design and implementation of the research. F.B., H.E., Z.T., and H.A contributed to the analysis of the results. Both F.B. and H.E authors developed the theory and performed the computations and the numerical simulations and verified the analytic methods and calculations.

Conflict of interest

The author(s) declare no competing interests.

Abbreviations

SDEs

stochastic differential equations

MCMC

Markov Chain Monte Carlo

FP

Fokker-Planck equation

MLEs

maximum likelihood estimators

ODE

ordinary differential equations

NLME

nonlinear mixed effects

SDME

stochastic differential mixed effects

OU

Ornstein-Uhlenbeck

GA

genetic algorithm

EN

elite number

SR

selection rate

CP

crossover probability

MP

mutation probability

SMM

stochastic minimal model

T2D

type 2 diabetes

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Bakrim Fadwa, Hamid El Maroufy and Hassan Ait Mousse (February 19th 2020). Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential Equations Using the Fokker-Planck Equation Solution, Numerical Modeling and Computer Simulation, Dragan M. Cvetković and Gunvant A. Birajdar, IntechOpen, DOI: 10.5772/intechopen.90751. Available from:

chapter statistics

52total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

The Possibilities of Modeling Petri Nets and Their Extensions

By Goharik Petrosyan

Related Book

First chapter

Introductory Chapter: Computer Simulation

By Dragan Cvetković

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us