\begin{matrix}c \\ ◯ \end{matrix} Revista de Matemática: Teoría y Aplicaciones 2019 26(1) : 83–98 cimpa – ucr issn: 1409-2433 (Print), 2215-3373 (Online)

doi: https:// doi. org/ 10.15517/ rmta. v26i1. 35527

estimación de un modelo de efectos mixtos usando un proceso de difusión

parcialmente observado

estimation of a mixed effects model using a partially observed diffusion

process

José Soto^{∗}  Saba Infante^{†}  Franklin Camacho^{‡}

Isidro R. Amaro^{§}

Received: 2/Mar/2018; Revised: 2/Nov/2018; Accepted: 4/Nov/2018

∗Escuela de Ciencias Matemáticas y Computacionales Universidad de Investigación de Tecnología Experimental Yachay Tech, Ecuador y Universidad de Los Andes Venezuela. E-Mail: jsoto@yachaytech.edu.ec , soto@ula.ve.

†Escuela de Ciencias Matemáticas y Computacionales Universidad de Investigación de Tecnología Experimental Yachay Tech, Ecuador y Universidad de Carabobo, Venezuela . E-Mail: sinfante@yachaytech.edu.ec

‡Misma dirección que/Same address as: J.Soto in Ecuador. E-Mail: fcamacho@yachaytech.edu.ec §Misma dirección que/Same address as: J.Soto in Ecuador. E-Mail: iamaro@yachaytech.edu.ec

84 j. soto – s. infante – f. camacho – i. amaro

Resumen

Consideramos un modelo de efectos mixtos general, donde la variabilidad de los efectos aleatorios de los individuos o unidades experimentales son incorporados a través de una ecuación diferencial estocástica. Estos modelos son útiles para analizar simultáneamente datos de medidas repetidas tomadas en tiempo discreto y con errores. Se implementó un algoritmo Monte Carlo por cadenas de Markov para hacer la inferencia a posteriori. Se realizó un análisis de diagnóstico sobre los parámetros estimados para detectar si el modelo es adecuado y mostrar su convergencia, además se muestran las trazas y las densidades estimadas a posteriori. La metodología se ilustró empleando datos sintéticos.

Palabras clave: modelos de efectos mixtos; ecuaciones diferenciales estocásticas; algoritmos Monte Carlo por cadenas de Markov.

Abstract

We consider a general mixed-effects model, where the variability of random effects of experimental individuals or units is incorporated through a stochastic differential equation. These models are useful for simultaneously analysing data from repeated measurements taken in discrete time and with errors. A Markov chain Monte Carlo algorithm was implemented to make the statistical inference a posteriori. A diagnostic analysis was carried out on the estimated parameters to detect if the model is suitable and show its convergence, in addition to the traces and posterior densities are shown. The methodology is illustrated using simulated data.

Keywords: mixed effects models; stochastic differential equations; Markov chains; Monte Carlo algorithms.

Mathematics Subject Classification: 62MXX.

1 Introducción

Una ecuación diferencial estocástica (SDE), por sus siglas en inglés, es una clase de procesos estocásticos de tiempo continuo que puede ser utilizada para modelar una amplia gama de fenómenos del mundo real. Cuando el proceso se observa en forma independiente en múltiples individuos o unidades experimentales en un tiempo t, la variabilidad entre las trayectorias pueden ser modeladas mediante la introducción de un efecto aleatorio a través del drift de la SDE, estos modelos se conocen en la literatura estadística como el modelo de efectos mixtos impulsado por una SDE [15, 4]. El modelo de difusión de efectos mixtos puede ser considerado como una extensión de los modelos de regresión para datos longitudinales, en el cual se introduce un componente de efectos aleatorios para explicar las variaciones “entre” y “dentro” de individuos. Estos modelos son una herramienta estándar para describir simultáneamente medidas repetidas tomadas en tiempos discretos y con errores sobre una serie de individuos. La estimación de parámetros en estos modelos no es clara, excepto en casos muy simples. Un método natural es utilizar la estimación por máxima verosimilitud, pero la densidad de transición pocas veces es conocida; la verosimilitud no tiene una forma explícita, lo que hace intratable la estimación de parámetros por este procedimiento. Estimar la verosimilitud implica tratar con densidades de transición complicadas, y manipular con integrales sobre el espacio de los parámetros no observados que pocas veces se pueden calcular en forma cerrada.

Recientemente, en la literatura estadística se desarrollan metodologías basadas en los modelos de difusión de efectos mixtos donde el proceso de difusión satisface la ecuación diferencial estocástica tipo Itô; es decir, se tiene una clase de procesos estocásticos de tiempo continuo que modelan la solución de la SDE como un proceso de Markov de primer orden. Esta metodología ha sido aplicada en dinámica de sistemas térmicos [1], experimentos en la industria farmacéutica ([4], [13], [9], entre otros), en neurociencia [8, 6], en pronóstico de energía solar y eólica [12], en dinámica neuronal [5] y en crecimiento de grietas [11]. Otras aplicaciones incluyen análisis de enfermedades de transmisión, tales como epidemias, series financieras, dinámica de crecimiento de población, y procesos intracelulares [16]

En partícular, en [7] llevaron acabo un análisis de un proceso continuo de curvas de crecimientos usando modelos mixtos definidos en términos de una SDE, ellos propusieron técnicas bayesianas para hacer las inferencias y estimar los parámetros; además, tomaron en consideración los casos cuando la SDE tienen una solución explícita y utilizaron un algoritmo Gibbs; y cuando la distribución condicional del proceso de difusión no tiene forma explícita, utilizan un esquema de Euler-Maruyama para hacer las aproximaciones. Delattre y Lavielle en [3] combinaron un algoritmo SAEM con el filtro de Kalman extendido (FKE) para estimar los parámetros de un modelo de efectos mixtos, adicionalmente, proporcionaron herramientas para la recuperación de las trayectorias de la difusión subyacente del modelo; la metodología es ilustrada mediante simulaciones usando datos provenientes de la industria farmacéutica. En [2] modelaron datos de panel usando un modelo de efectos mixtos con heterogeneidad en los efectos aleatorios individuales. Compararon dos enfoques utilizando mezclas de distribuciones Gaussianas y de Dirichlet; la metodología fue ilustrada utilizando datos sintéticos y datos reales para obtener las estimaciones a posteriori de los parámetros. Con respecto a [16] se propone un esquema Monte Carlo por cadenas de Markov (MCMC) en el escenario de los modelos de difusión de efectos mixtos. La propuesta incluye un algoritmo que permite muestrear eficientemente a partir de las distribuciones marginales, las cuales pueden tener dinámicas lineales o no lineales entre los tiempos de observación, ellos ilustran la propuesta usando datos sintéticos a partir de un modelo SDE de crecimiento de árboles de naranjas y usando datos reales obtenidos del control de plagas en plantas de algodón bajo diferentes regímenes de tratamiento. Picchini y Lyng [14] propusieron un modelo de efectos mixtos de una SDE para modelar el crecimiento de dinámicas de medidas repetidas de volúmenes de tumores en ratones, ellos consideraron dos secciones divididas en una representación a través de fracciones de tumores muertos y los que sobreviven a un tratamiento; ellos estimaron los parámetros del modelo usando una metodología Bayesiana exacta y un método aproximado. En este trabajo se propone implementar un algoritmo Monte Carlo por cadenas de Markov (MCMC), por sus siglas en inglés, para estimar las soluciones y parámetros de un modelo de efectos aleatorios mixto definido a través de una SDE. Esta propuesta tiene la ventaja que no requiere usar la densidad de transición y permite generar muestras de las distribuciones marginales a posteriori de los estados desconocidos del sistema y parámetros del modelo, las cuales se utilizan para aproximar la distribución a posteriori y luego hacer la inferencia. El resto del artículo es como sigue: en la sección 2 se realiza el planteamiento del problema, en la sección 3 se define la metodología a utilizar, esto es, se calculan las distribuciones marginales y se define el algoritmo a utilizar para generar las muestras y llevar a cabo la estimación de estados y parámetros; en la sección 4 se lleva acabo una aplicación, y en la sección 5 se establece una discusión y las conclusiones.

2 Planteamiento del problema

Considérese el caso en que se tiene N unidades experimentales escogidos aleatoriamente de una población teórica y asociada con cada unidad experimental i, se tiene un proceso de Itô de tiempo continuo de dimensión d, dado por:

x_{ij} = x(t_{ij}, ϕ_{i}) ; ϕ_{i}|μ, ν ∼ N (μ, \frac{1}{ν} )
dx(t, ϕ_{i}) = μ (ϕ_{i}, t, x_{t}) dt + σ (κ^{2}, t, x_{t}) dB_{i}(t)
y_{ij} = g (x_{ij}, t_{ij}) + ϵ_{ij} , ϵ_{ij} ∼ N (0, σ_{ϵ}^{2}) ,  (1)

donde: x = x_{ij} := x_{i}(t_{j}), x_{i} = {x_{ij}}_{j=0,...,n_{i}} denotan los estados (latentes) desconocidos sobre el individuo i, y = y_{ij} := y_{i}(t_{j}), y y_{i} = {y_{ij}}_{j=0,...,n_{i}} denota las observaciones sobre el individuo i, i = 1, …, I.

Cada proceso x(t) ∈ X tiene dimensión d_{x} y está definido sobre un espacio X, y_{ij} ∈ Y tiene dimensión d_{y}. El modelo tiene la siguiente interpretación: para cada sujeto i, {x_{i}(t)} representa el estado oculto no observable de un proceso biológico de interés en el tiempo t. Los estados verdaderos {x_{i}(t)} son imposible de obtenerlos y en su lugar se observan{y_{i}(t)} en un tiempo discreto t, las cuales son medidas con errores. Las dinámicas del proceso latente son modeladas a través de funciones de drift μ(.) de dimensión d y el coeficiente de difusión σ(.) representa una matriz definida positiva de orden d × d, ambas expresiones se suponen conocidas, salvo los parámetros del modelo y el efecto aleatorio ϕ_{i} = (ϕ_{1}^{i}, …, ϕ_{q}^{i}). Cada efecto aleatorio individual ϕ_{i} sigue un modelo de probabilidad p(.|ψ) dependiendo sobre el parámetro de la población ψ y dB_{i}(t) es un vector de dimensión d de un movimiento Browniano estándar no correlacionados. La solución de la ecuación diferencial {x_{i}(t)} para cada individuo i está sujeta a las mismas condiciones de regularidad (condición de crecimiento global y de Lipschitz en la variable x) para garantizar la existencia y unicidad de la solución.

Cuando las observaciones y_{ij} son condicionalmente independientes dado los estados x_{ij}, entonces se tiene un modelo en la forma espacio estado (o un modelo de Markov oculto). El objetivo de la inferencia es estimar los estados latentes y los parámetros θ = (μ, ν, κ^{2}, σ_{ϵ}^{2}). Una de las características importantes de este modelo propuesto es la capacidad de discriminar la variabilidad dentro los individuos (parámetro κ), entre individuos (parámetro ϕ_{i}) y la variabilidad residual σ_{ϵ}. Los autores de [14] señalan que cada observación y_{ij} depende de tres fuentes de variabilidad: la variabilidad de las dinámicas estocásticas individuales x_{i}(t_{j}), los efectos aleatorios ϕ_{i} y el error residual ϵ_{ij}.

Los métodos de estimación clásicos y Bayesianos requieren de la evaluación de la función de verosimilitud:

p(y|θ) = ∏_{i=1}^{I} p(y_{i}|θ),  (2)

donde

p(y_{i}|θ) = ∫ p(y_{i}|ϕ_{i}, θ)p(ϕ_{i}|θ)dϕ_{i} = ∫ (∫ p(y_{i}|x_{i}, θ)p(x_{i}|ϕ_{i}, θ)dx_{i}) p(ϕ_{i}|θ)dϕ_{i}
= ∫ (∫ ∏_{j=0}^{n_{i}} p(y_{ij}|x_{ij}, ϕ_{i}, θ)p(x_{ij}|x_{i,j−1}, ϕ_{i}, θ)p(x_{i,0}|ϕ_{i}, θ)dx_{ij}) p(ϕ_{i}|θ)dϕ_{i}. (3)

El término p(x_{i,0}|ϕ_{i}, θ) puede eliminarse en los casos cuando x_{i,0} se incluye entre los efectos aleatorios.

En general, la función de verosimilitud p(y_{i}|θ) no es analíticamente tratable. En problemas de interés, la forma de la SDE asociada con cada unidad experimental i no permite una solución explícita, entonces una manera de abordar el problema es discretizando los estados usando una aproximación Euler-Maruyama. Supóngase que se tienen observaciones x_{i0}, …, x_{in_{i}} en los puntos t_{i0}, …, t_{in_{i}}, tal que

t_{i0} = t_{0} < t_{1} < … < t_{m_{i}} < … < t_{M_{i}} = t_{in_{i}}  (4)

denotan discretizaciones igualmente espaciadas del intervalo [t_{i0}, t_{in_{i}}]. Sea Δ el tamaño del paso de la discretización. Entonces el proceso latente denotado por x_{ij} de la SDE se puede discretizar mediante la siguiente estructura jerárquica.

x_{io} = x_{0}(ϕ_{i})
x_{ij} = x_{ij−1} + μ (ϕ_{i}, t_{ij}, x_{ij}) Δ_{ij} + σ (κ^{2}, t_{ij}, x_{ij}) z_{i,j},

z_{i,j} ∼ N(0, 1); ϕ_{i}|μ, ν ∼ N (μ, \frac{1}{ν} ); μ ∼ N(a, b)

σ_{ϵ}^{2} ∼ IG(g, h); ν ∼ G(c, d); κ^{2} ∼ IG(e, f); Δ_{ij} := t_{ij} − t_{ij−1}

y_{ij} = x_{ij} + ϵ_{ij} , ϵ_{ij} ∼ N(0, σ_{ϵ}^{2}),  (5)

donde G y IG denotan distribuciones Gamma e inversa Gamma respectivamente.

Sean β es una constante y

μ (ϕ_{i}, t_{ij}, x_{ij}) = ϕ_{i} − βx_{ij−1}; σ (κ^{2}, t_{ij}, x_{ij}) = κ\sqrt{Δ_{ij}},  (6)
𝔼 (x_{ij}) = x_{ij−1} + (ϕ_{i} − βx_{ij−1}) Δ_{ij}; 𝕍ar (x_{ij}) = κ^{2}Δ_{ij}. (7)

3 Estimación de estados y parámetros

En el modelo (1) se observan tres fuentes de errores: la variabilidad entre individuos que se mide a través de \frac{1}{ν}, el ruido del proceso dinámico κ^{2} que refleja las fluctuaciones aleatorias del proceso dinámico teórico y el ruido de medición de las observaciones σ_{ϵ}^{2}. El problema consiste en determinar la distribución a posteriori de los estados x_{1:t} = (x_{1}, …, x_{t}) y los parámetros θ = (μ, ν, κ^{2}, σ_{ϵ}^{2}), dado los datos observados y_{1:t} = (y_{1}, …, y_{t}); es decir, p(x_{1:t}, θ|y_{1:t}).

Como los modelos definidos en (1), y (5) son no lineales, las distribuciones a posteriori no son explícitas, entonces se requiere de un procedimiento de estimación recursivo, en este trabajo se propone usar algoritmo Monte Carlo por cadenas de Markov (muestreador de Gibbs).

3.1 Distribuciones full condicionales

La distribución a priori conjunta es dada por

p(μ, ν, ϕ, κ^{2}, σ_{ϵ}^{2}) = p (ϕ|μ, ν, κ^{2}, σ_{ϵ}^{2}) p (μ, ν, κ^{2}, σ_{ϵ}^{2})
= p(ϕ|μ, ν)p(μ)p(ν)p(κ^{2})p(σ_{ϵ}^{2}),  (8)

donde

p(ϕ|μ, ν) = p(ϕ_{1}|μ, ν), …, p(ϕ_{I}|μ, ν) ∝ ν^{\frac{I}{2}} exp (−\frac{ν}{2} ∑_{i=1}^{I}(ϕ_{i} − μ)^{2}) . (9)

La distribución a posteriori:

p(μ, ν|ϕ) = p(ϕ|μ, ν)p(μ)p(ν)
∝ ν^{\frac{I}{2} +c−1} exp\{−(\frac{ν}{2} ∑_{i=1}^{I}(ϕ_{i} − μ)^{2}+ \frac{b}{2}(μ − a)^{2}+dν)\}. (10)

Las distribuciones marginales son:

• La marginal p (μ|ϕ, ν) es

p (μ|ϕ, ν) ∝ exp \{−\frac{1}{2} (ν ∑_{i=1}^{I}(ϕ_{i} − μ)^{2} + b(μ − a)^{2})\}

= exp \{−\frac{b + Iν}{2} (μ − \frac{ν ∑_{i=1}^{I} ϕ_{i} + ab}{b + Iν})^{2}\}

∼ N (\frac{ab + ν ∑_{i=1}^{I} ϕ_{i}}{b + Iν} , \frac{1}{b + Iν} ) .  (11)

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 26(1): 83–98, Jan–Jun 2019

90 j. soto – s. infante – f. camacho – i. amaro

• La marginal p (ν|ϕ, μ) es

p (ν|ϕ, μ) = ν^{\frac{I}{2}} exp \{−\frac{ν}{2} ∑_{i=1}^{I}(ϕ_{i} − μ)^{2}\} ν^{c−1} exp \{−dν\}

∝ ν^{c+ \frac{I}{2} −1} exp \{−ν (d + \frac{1}{2} ∑_{i=1}^{I}(ϕ_{i} − μ)^{2})\}

∼ G (c + \frac{I}{2}, d + \frac{1}{2} ∑_{i=1}^{I}(ϕ_{i} − μ)^{2}) .  (12)

Por otra parte, se tiene

x_{ij}|x_{ij−1}, ϕ, κ^{2} ∼ N (x_{ij−1} + (ϕ_{i} − βx_{ij−1})Δ_{ij}, κ^{2}Δ_{ij}) .  (13)

La verosimilitud de los datos x_{ij}, es dada por

L (x_{ij}; ϕ, κ^{2}) = ∏_{i=1}^{I} ∏_{j=0}^{n_{i}} p (x_{ij}|x_{ij−1}, ϕ, κ^{2})

= ∏_{i=1}^{I} ∏_{j=0}^{n_{i}} \frac{1}{\sqrt{2πκ^{2}Δ_{ij}}} exp \{− \frac{1}{2κ^{2}Δ_{ij}} [x_{ij} − x_{ij−1} − (ϕ_{i} − βx_{ij−1})Δ_{ij}]^{2}\}

∝ \frac{1}{(κ^{2}Δ_{ij}) \frac{∑_{i=1}^{I}(n_{i}+1)}{2}} exp\{\frac{−1}{2κ^{2}} ∑_{i=1}^{I}∑_{j=0}^{n_{i}} \frac{1}{Δ_{ij}} [x_{ij} −x_{ij−1}−(ϕ_{i}−βx_{ij−1})Δ_{ij}]^{2}\} . (14)

Además como κ^{2} ∼ IG(e, f), entonces:

• La marginal p (κ^{2}|x_{ij}, ϕ) es dada por

p (κ^{2}|x_{ij}, ϕ) ∝ L (x_{ij}; ϕ, κ^{2}) p(κ^{2})

∝ \frac{1}{(κ^{2}Δ_{ij})^{^{\frac{1}{2}} ∑{I \atop i=1}(n_{i}+1)}} exp\{− \frac{1}{2κ^{2}} ∑_{i=1}^{I}∑_{j=0}^{n_{i}} \frac{1}{Δ_{ij}} [x_{ij} −x_{ij−1}−(ϕ_{i}−βx_{ij−1})Δ_{ij}]^{2}\}

× ( \frac{1}{κ^{2}} )^{e+1} exp (− \frac{f}{κ^{2}} )

∼IG(e+2+ \frac{1}{2} ∑_{i=1}^{I} (n_{i}+1) , f + \frac{1}{2} ∑_{i=1}^{I}∑_{j=0}^{n_{i}} \frac{1}{Δ_{ij}} [x_{ij} −x_{ij−1}−(ϕ_{i}−βx_{ij−1})Δ_{ij}]^{2}) . (15)

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 26(1): 83–98, Jan–Jun 2019

estimación de un modelo de efectos mixtos usando … 91

Por otra parte, se tiene que

y_{ij}|x_{ij}, ϕ, κ^{2}, σ_{ε}^{2} ∼ N (x_{ij}, σ_{ε}^{2}) .

La verosimilitud de los datos y_{ij} es dada por

L (y_{ij}; ϕ, κ^{2}, σ_{ε}^{2}) = ∏_{i=1}^{I} ∏_{j=0}^{n_{i}} p (y_{ij}|x_{ij}, ϕ, κ^{2}, σ_{ε}^{2})

= ∏_{i=1}^{I} ∏_{j=0}^{n_{i}} \frac{1}{\sqrt{2πσ_{ϵ}^{2}}} exp [−\frac{1}{2σ_{ϵ}^{2}} (y_{ij} − x_{ij})^{2}]

∝ \frac{1}{(σ_{ε}^{2})^{\frac{1}{2} ∑_{i=0}^{I}(n_{i}+1)}} exp \{−\frac{1}{2σ_{ϵ}^{2}} ∑_{i=1}^{I} ∑_{j=0}^{n_{i}} (y_{ij} − x_{ij})^{2}\} . (16)

Además, se supone que σ_{ε}^{2} ∼ IG(g, h), entonces:

• La marginal de p (σ_{ε}^{2}|x_{ij}, y_{ij}, ϕ, κ^{2}) es dada por

p (σ_{ε}^{2}|x_{ij}, y_{ij}, ϕ, κ^{2}) = L (y_{ij}; ϕ, κ^{2}, σ_{ε}^{2}) p(σ_{ε}^{2})

∝ \frac{1}{(σ_{ε}^{2})^{\frac{1}{2} ∑_{i=0}^{I}(n_{i}+1)}} exp [− \frac{1}{2σ_{ϵ}^{2}} ∑_{i=1}^{I} ∑_{j=0}^{n_{i}} (y_{ij} − x_{ij})^{2}] \frac{1}{σ_{ε}^{2}} exp (− \frac{h}{σ_{ε}^{2}} )

= ( \frac{1}{σ_{ε}^{2}} )^{g+ \frac{1}{2} ∑_{i=0}^{I}(n_{i}+1)+2−1}exp \{− \frac{1}{σ_{ε}^{2}} [h + \frac{1}{2} ∑_{i=1}^{I} ∑_{j=0}^{n_{i}} (y_{ij} − x_{ij})^{2}]\}

∼ IG (g + \frac{1}{2} ∑_{i=1}^{I}(n_{i} + 1), h + \frac{1}{2} ∑_{i=1}^{I} ∑_{j=0}^{n_{i}} (y_{ij} − x_{ij})^{2}) .  (17)

• La marginal p (ϕ_{i}|x_{i}, y_{i}) se obtiene como sigue, nótese que la densidad

conjunta p(x, y, ϕ) se puede escribir como sigue

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 26(1): 83–98, Jan–Jun 2019

p(x, y, ϕ) = p(y|x, ϕ)p(x, ϕ) = p(y|x, ϕ)p(x|ϕ)p(ϕ)

= ∏_{i=1}^{I} p(ϕ_{i}|μ, ν) ∏_{j=0}^{n_{i}} p(y_{ij}|x_{ij})p(x_{ij}|x_{ij−1}, ϕ_{i})

= ∏_{i=1}^{I} ν^{\frac{1}{2}} exp \{−\frac{ν}{2} (ϕ_{i} − μ)^{2}\} ∏_{j=0}^{n_{i}} \frac{1}{2πσ_{ϵ}^{2}} exp \{−\frac{1}{2σ_{ϵ}^{2}} (y_{ij} − x_{ij})^{2}\}

× \frac{1}{\sqrt{2πκ^{2}Δ_{ij}}} exp [−\frac{1}{2κ^{2}Δ_{ij}} (x_{ij} − (x_{ij−1} + (ϕ_{i} − βx_{ij})Δ_{ij})]^{2} . (18)

Así que

p(x, y, ϕ) = p(ϕ|x, y)p(x, y) ∝ p(ϕ|x, y).  (19)

En consecuencia

p(ϕ_{i}|x_{i}, y_{i}) = \frac{ν^{\frac{1}{2}}}{(κ^{2}Δ_{ij})^{\frac{I}{2}}} exp \{− \frac{v}{2} (ϕ_{i} − μ)^{2}+∑_{j=0}^{n_{i}} \frac{1}{κ^{2}Δ_{ij}} [x_{ij} −x_{ij−1}(ϕ_{i}−βx_{ij−1})Δ_{ij}]^{2}\} , (20)

desarrollando el exponente y completando cuadrado en los ϕ, se obtiene:

− \frac{v}{2} (ϕ_{i} − μ)^{2} + ∑_{j=0}^{n_{i}} \frac{1}{κ^{2}Δ_{ij}} [x_{ij} − x_{ij−1}(ϕ_{i} − βx_{ij−1})Δ_{ij}]^{2}

= \frac{1}{\frac{κ^{2}}{νκ^{2}+∑_{j=0}^{n_{i}} Δ_{ij}}} [ϕ_{i} − (\frac{μνκ^{2} + ∑_{j=0}^{n_{i}} βx_{ij−1}Δ_{ij}^{2} + (x_{ij} − x_{ij−1})Δ_{ij}}{νκ^{2} + ∑_{j=0}^{n_{i}} Δ_{ij}} )]^{2} ,

de donde se deduce que

p(ϕ_{i}|x_{i}, y_{i})∝exp\{\frac{1}{\frac{κ^{2}}{νκ^{2}+∑_{j=0}^{n_{i}} Δ_{ij}}} [ϕ_{i}−(\frac{μνκ^{2}+∑_{j=0}^{n_{i}} βx_{ij−1}Δ_{ij}^{2} +(x_{ij} −x_{ij−1})Δ_{ij}}{νκ^{2}+∑_{j=0}^{n_{i}} Δ_{ij}} )]^{2}\}
∼N( \frac{μνκ^{2} + ∑_{j=0}^{n_{i}} βx_{ij−1}Δ_{ij}^{2} + (x_{ij} − x_{ij−1})Δ_{ij}}{νκ^{2} + ∑_{j=0}^{n_{i}} Δ_{ij}} , \frac{1}{\frac{κ^{2}}{νκ^{2}+∑_{j=0}^{n_{i}} Δ_{ij}}} ) . (21)

• La marginal p (x|θ, ϕ, y), se obtiene como sigue

p (x|θ, ϕ, y) = p (y|θ, ϕ) p (y|x, σ_{ϵ}^{2}) ,  (22)

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 26(1): 83–98, Jan–Jun 2019

donde

p (x|θ, ϕ, y) = p(x|θ, ϕ)p(y|x, σ_{ϵ}^{2})

= ∏_{i=1}^{I} ∏_{j=0}^{n_{i}} p(y_{i}|x_{i}, θ)p(x_{i}|ϕ_{i}, θ)

= ∏_{i=1}^{I} ∏_{j=0}^{n_{i}} p(y_{ij}|x_{ij}, θ)p(x_{ij}|x_{ij−1}, ϕ_{i}, θ)

= ∏_{i=1}^{I} ∏_{j=1}^{n_{i}} \frac{1}{\sqrt{2πσ_{ϵ}^{2}}} exp \{− \frac{1}{2σ_{ϵ}^{2}} (y_{ij} − x_{ij})^{2}\}

× \frac{1}{\sqrt{2πκ^{2}Δ_{ij}}} exp \{− \frac{1}{2κ^{2}Δ_{ij}} [x_{ij} − (x_{ij−1} − (ϕ_{i} − βx_{ij−1}))]^{2}\}

∝exp\{−∑_{i=1}^{I} ∑_{j=0}^{n_{i}} [\frac{1}{2σ_{ϵ}^{2}} (y_{ij} −x_{ij})^{2}+ \frac{1}{2κ^{2}Δ_{ij}} [x_{ij} −(x_{ij−1}−(ϕ_{i}−βx_{ij−1}))]^{2}]\} , (23)

desarrollando el exponente y completando cuadrado en los x_{ij}, se obtiene:

\frac{1}{2σ_{ϵ}^{2}} (y_{ij} − x_{ij})^{2} + \frac{1}{2κ^{2}Δ_{ij}} [x_{ij} − (x_{ij−1} − (ϕ_{i} − βx_{ij−1}))]^{2}
∝ \frac{1}{2} \frac{1}{\frac{σ_{ϵ}^{2}κ^{2}Δ_{ij}}{κ^{2}Δ_{ij}}} [x_{ij} − (\frac{y_{ij}κ^{2}Δ_{ij}}{κ^{2}Δ_{ij} + σ_{ϵ}^{2}} + \frac{σ_{ϵ}^{2} (x_{ij−1} + (ϕ_{i} − βx_{ij−1})) Δ_{ij}}{κ^{2}Δ_{ij} + σ_{ϵ}^{2}} )]^{2} ,

así que la marginal de p (x|θ, ϕ, y) es

p (x|θ, ϕ, y) ∝ exp \{∑_{i=1}^{I} ∑_{j=0}^{n_{i}} −\frac{1}{2} \frac{1}{\frac{σ_{ϵ}^{2}κ^{2}Δ_{ij}}{κ^{2}Δ_{ij}}} [x_{ij} − (\frac{y_{ij}κ^{2}Δ_{ij}}{κ^{2}Δ_{ij} + σ_{ϵ}^{2}}
+\frac{σ_{ϵ}^{2} (x_{ij−1} + (ϕ_{i} − βx_{ij−1})) Δ_{ij}}{κ^{2}Δ_{ij} + σ_{ϵ}^{2}} )]^{2}\}

∼N (\frac{y_{ij}κ^{2}Δ_{ij} +σ_{ϵ}^{2} (x_{ij−1}+(ϕ_{i}−βx_{ij−1})) Δ_{ij}}{κ^{2}Δ_{ij} + σ_{ϵ}^{2}} , \frac{σ_{ϵ}^{2}κ^{2}Δ_{ij}}{κ^{2}Δ_{ij} + σ_{ϵ}^{2}} ) . (24)

−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

Algoritmo Gibbs: para el modelo de efectos mixtos usando una SDE. −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− Paso 1. Generar: p (μ|ϕ, ν) como en la ecuación dada en (11).

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 26(1): 83–98, Jan–Jun 2019

Paso 2. Generar: p (ν|ϕ, μ) como en la ecuación dada en (12).

Paso 3. Generar: p (κ^{2}|x_{ij}, ϕ, κ^{2}) como en la ecuación dada en (15).

Paso 4. Generar: p (σ_{ε}^{2}|x_{ij}, y_{ij}, ϕ, κ^{2}) como en la ecuación dada en (17).

Paso 5. Generar: p (ϕ_{i}|x_{i}, y_{i}) usando la ecuación (21).

Paso 6. Generar: p (x|θ, ϕ, y) usando la ecuación (24).

−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

4 Resultados

Para realizar el proceso de inferencia; esto es, llevar a cabo la estimación de las soluciones y parámetros del modelo de efectos mixtos a través de una SDE, el procedimiento es el siguiente, se simularon N = 4000 observaciones sobre i = 10 individuos. En este sentido, se implementó el muestreador de Gibbs con las distribuciones a posteriori marginales obtenidas para este modelo (ver el algoritmo Gibbs), La programación del algoritmo fue implementado usando el software R. Se obtuvieron las muestras simuladas para cada uno de los parámetros y estas salidas son analizadas usando el paquete BOA disponible en R para chequear la convergencia de los parámetros.

En la Tabla 1, se muestra un resumen estadístico de los estimados a posteriori de los parámetros que contiene: la media (mean), la desviación estándar (SD), un estimador Naive (que representa la desviación típica de la muestra dividida por la raíz cuadrada del tamaño de la muestra), un error Monte Carlo (que se obtiene usando una estimación de una serie de tiempo, calculando la raíz cuadrada de la estimación de la varianza de la densidad espectral, dividida por el tamaño de la muestra, esto permite obtener error estándar asintótico), una estimación Batch (calculada usando la desviación típica de la muestra de dos promedios de lotes consecutivos de tamaño 50, dividida por la raíz cuadrada del número de lotes) y la función autocorrelación estimada llamada Batch ACF (se obtiene usando dos lotes de muestras consecutivas, se espera que la función Batch ACF debe estar cerca de cero).

En la Figura 1, se muestra las trazas a posteriori para los valores simulados de μ, ν, ϕ, y k^{2} respectivamente.

La Figura 2, muestra las trazas estimadas para cuatro estados soluciones de las simuladas para los 10 individuos considerados, donde se observa convergencia.

Mean  SD  Naive SE  MC Error  Batch SE  Batch ACF
μ 0.11480131 0.0727719731 1.627231e-03 8.666184e-03 1.137195e-02 0.9087463
ν 0.02483368 0.0006701769 1.498561e-05 2.566849e-05 2.698157e-05 0.6452355
ϕ -0.18328172 6.3800266347 1.426617e-01 1.227308e-01 1.282178e-01 -0.1134824
k^{2} 0.16059294 0.0028727866 6.423746e-05 3.073844e-04 4.079278e-04 0.9164281
σ^{2} 0.93764967 0.1824914933 4.080634e-03 2.246261e-02 2.918704e-02 0.9253488

Tabla 1: Resumen estadístico de los parámetros estimados.

Figura 1: Trazas simuladas para los parámetros μ, ν, ϕ, y k^{2}.

En la Figura 3, se muestran las densidades estimadas a posteriori para los parámetros μ, ν, ϕ, y k^{2}. En cada uno de estos gráficos se observa que los parámetros convergen y se estabilizan casi simétricamente alrededor de un valor que corresponde a su media a posteriori. En la Figura 3 se muestra un gráfico que representa la densidad a posteriori estimada a partir de las cadenas de cada parámetro usando muestras de tamaño 2000 y 4000, respectivamente.

Figura 2: Trazas simuladas para X.

Figura 3: Densidades empíricas a posteriori de μ, ν, ϕ y k^{2}.

En la Tabla 2, se muestran los resultados de la prueba de diagnóstico de convergencia de [10] que consiste en una prueba de estacionaridad basada en la teoría del puente Browniano y el estadístico de Cramer-Von-Mises. Si hay evidencia de no estacionariedad, la prueba se repite después de descartar el primer 10% de las iteraciones. Este proceso continúa hasta que la cadena resultante pasa la prueba o más del 50% de las iteraciones hayan sido descartadas. La falla de la cadena para pasar esta prueba indica que se necesita una ejecución más larga de la muestra de MCMC para lograr la convergencia. Los resultados obtenidos indican la convergencia de las trazas simuladas de los parámetros ν, k^{2}, ϕ y σ^{2} (es decir, pasan la prueba).

Stationarity Keep Discard C-Von-M Halfwidth  Mean  Halfwidth
Test  Test
μ  failed  400  600 1.1377662 failed 0.08136377 0.0146534384
ν  failed  400  600 2.4806535 passed 0.02483932 0.0001961795
ϕ  passed  1000  0  0.1157022 failed 0.14841167 0.3634571976
k^{2}  failed  400  600 0.9448437 passed 0.15681077 0.0007604466
σ^{2}  failed  400  600 5.7027727 passed 0.92951749 0.0460030231

Tabla 2: Prueba de diagnóstico de convergencia para las muestras simuladas.

5 Conclusiones y discusiones

El modelaje estocástico propuesto en este trabajo incorpora estructuras más generales para determinar el error de medición y los errores de estimación del modelo. Los modelos mixtos definidos por una ecuación diferencial estocástica pueden ser vistos como una extensión de una cadena de Markov en un espacio estado finito, de tiempo continuo y cubren una amplia clase de problemas de procesos biológicos, en donde los tiempos de observación son irregulares y a diferentes escalas. El algoritmo propuesto permite estimar a posteriori las soluciones y parámetros del modelo, incluyendo ruido de observación. Se demuestra la convergencia de los parámetros estimados cuando se tiene una densidad de transición explícita. La metodología propuesta nos proporciona una solución flexible, robusta bajo una implementación secuencial automática.

Agradecimientos

Agradecemos al Editor de la Revista y a los árbitros por sus valiosas sugerencias y contribuciones para mejorar este manuscrito. Esté trabajo fue parcialmente apoyado por Instituto de Estadística Aplicada y Computación de la Universidad de los Andes, Venezuela, y la Escuela de Ciencias Matemáticas y Computacionales, Universidad Yachay Tech, Ecuador.

Referencias

[1] Bacher, P.; Madsen, H. (2011) “Identifying suitable models for the heat dynamics of buildings”, Energy and Buildings 43(7): 1511–1522.

[2] Dahlin, J.; Kohn, R.; Schön, T. (2016) “Bayesian inference for mixed effects models with heterogeneity”, Technical Report from Automatic Control at Linköpings Universitet, (LiTH-ISY-R-3091). Linköping, Sweeden.

[3] Delattre, M.; Lavielle, M. (2013) “Coupling the SAEM algorithm and the extended Kalman filter for maximum likelihood estimation in mixedeffects diffusion models”, Statistics and Its Interface 6(4): 519–532.

[4] Ditlevsen, S.; De Gaetano, A. (2005) “Mixed effects in stochastic differential equation models”, Statistical Journal 3(2): 137–153.

[5] Ditlevsen, S.; Samson, A. (2013) “Introduction to stochastic models in biology”, in: M. Bachar & J. Batzel & S. Ditlevsen (Eds.) Stochastic Biomathematical Models, Springer, Berlín, Heidelberg: 3–35.

[6] Doucet, A.; De Freitas, N.; Gordon, N. (2001) Sequential Monte Carlo Methods in Practice. Springer, New York.

[7] Donnet, S.; Foulley, J.; Samson, A. (2010) “Bayesian analysis of growth curves using mixed models defined by stochastic differential equations”, Biometrics 66(3): 733–741.

[8] Faugeras, O.; Touboul, J.; Cessac, B. (2009) “A constructive meanfield analysis of multi-population neural networks with random synaptic weights and stochastic inputs”, Frontiers in Computational Neuroscience 3(1): 1–28.

[9] Hansen, A.; Duun-Henriksen, A.; Juhl, R.; Schmidt, S.; Norgaard, K.; Jorgensen, J.; Madsen, H. (2014) “Predicting plasma glucose from interstitial glucose observations using bayesian methods”, Journal of Diabetes Science and Technology 8(2): 321–330.

[10] Heidelberger, P.; Welch, P. D. (1983) “Simulation run length control in the presence of an initial transient”, Operations Research 31(6): 1109–1144.

[11] Hermann, S.; Ickstadt, K.; Müller, C. (2016) “Bayesian prediction of crack growth based on a hierarchical diffusion model”, Applied Stochastic Models in Business and Industry 32(4): 494–510.

[12] Iversen, J.; Morales, J.; Moller, J.; Madsen, H. (2014) “Probabilistic forecasts of solar irradiance by stochastic differential equations”, Environmetrics 25(3): 152–164.

[13] Overgaard, R.; Jonsson, N.; Tornoe, C.; Madsen, H. (2005) “Non-linear mixed-effects models with stochastic differential equations: implementation of an estimation algorithm”, Journal of Pharmacokinetics and Pharmacodynamics 32(1): 85–107.

[14] Picchini, U.; Lyng, J. (2016) “Stochastic differential equation mixed effects models for tumor growth and response to treatment”, Preprint, Universidad de Lund, Suecia.

[15] Pinheiro, J.; Bates, D. (2000) Mixed-Effects Models in S and S-PLUS. Springer-Verlag, New York.

[16] Whitaker, G.A.; Golightly, A.; Boys, R.; Sherlock, C. (2017) “Bayesian inference for diffusion-driven mixed-effects models”, Bayesian Analysis 12(2): 435–463.