El modelo de infra-reportado

El fenómeno del infra-reporte en el contexto de la investigación biomédica ha sido ampliamente estudiado, al menos desde los trabajos de Bernard et. al. en 2011 [2]. Una referencia un poco más reciente es Fernández-Fontelo et. al. en 2016 [4], donde se propone un mecanismo muy simple para modelar el fenómeno. Consideremos un proceso oculto (el número real diario de casos positivos, que no es observable puesto que no se hacen tests diarios a toda la población de interés) X_n con una estructura autoregresiva de orden 1, para datos enteros y con innovaciones Poisson, Po-INAR(1):

     $$X_n=\alpha \circ X_{n-1}+W_n(\lambda)$$

donde  0<\alpha<1 es un parámetro fijo,  W_n \sim Poisson( \lambda ), i.i.d., independientes de  X_n y  \circ es el operador de thinning binomial:

     \begin{equation*}\alpha \circ X_{n-1}=\sum_{i=1}^{X_{n-1}}Z_i \end{equation*}

con  Z_i variables aleatorias i.i.d. con distribución de Bernoulli( \alpha).
El proceso INAR(1) es una cadena de Markov homogénea con probabilidades de transición

     $$\pr (X_n = i | X_{n-1} = j) = \sum_{k=0}^{i\wedge j} {j \choose k}  \alpha^k (1-\alpha)^{j-k} \pr(W_n=i-m)$$

La esperanza y varianza del operador de thinning binomial son

    \[\textrm{{\bf E}}\left(\alpha \circ X_{n-1}|X_{n-1}=x_{n-1}\right)=\alpha x_{n-1}\]


y

    \[\textrm{{\bf Var}}\left(\alpha \circ X_{n-1}|X_{n-1}=x_{n-1}\right)=\alpha (1-\alpha)x_{n-1}\]

Un esquema simple de infra-reporte

 

El fenómeno de infra-reporte se modela asumiendo que los recuentos observados tienen la forma siguiente: 

(1)   \begin{align*} Y_n = \left\{ \begin{array}{ll} X_n & \text{con probabilidad } 1-\omega\\ q \circ X_n & \text{con probabilidad } \omega, \end{array} \right.\nonumber \end{align*}

donde \omega y q representan la frecuencia e intensidad del proceso infra-reportado. Esto es, para cada n, observamos X_n con probabilidad 1-\omega, y un q-thinning de X_n con probabilidad \omega, independientemente del pasado \{X_j: j \leq n\}. Por tanto, lo que se observa (los recuentos reportados) es

    \[Y_n = (1- \uno_n)X_n + \uno_n \sum_{j=1}^{X_n} \xi_j \quad\quad \uno_n\sim\mbox{Bern}(\omega), \quad \xi_j\sim\mbox{Bern}(q)\]

Seguramente habrá que modificar este modelo para hacer que tanto q como \lambda dependan del tiempo.

Propiedades del modelo

La esperanza y la variancia del proceso INAR(1) X_n con innovaciones Poisson(\lambda) es \mu_X=\sigma_X^2=\frac{\lambda}{1-\alpha}. Las auto-covariancias y autocorrelaciones son

    \[\gamma_X(k)=\alpha^{|k|}\lambda \textrm{ y }\rho_{X}(k)=\alpha^{|k|}\]

respectivamente. Por otro lado, \E Y_n =\mu_{Y}=\mu_X(1-\omega(1-q)).  La autocovariancia del proceso observado Y_n es

    \[\gamma_Y(k)=(1-\omega(1-q))^2\alpha^{|k|}\mu_X\]

Por lo tanto, la función de auto-correlación de Y_n es un múltiplo de \rho_{X}(k):

    \[\rho_Y(k)=\frac{(1-\alpha)(1-\omega(1-q))^2}{(1-\alpha)(1-\omega(1-q))+\lambda(\omega(1-\omega)(1-q)^2)}\alpha^{|k|}=c(\alpha,\lambda,\omega,q)\alpha^{|k|}.\]

La forma de la esperanza de Y_n permite además obtener predicciones para Y_{n+k} habiendo observado hast Y_n.

Estimación de parámetros

La distribución marginal de  Y_n es una mezcla de dos distribuciones de Poisson

(2)   \begin{align*} Y_n &\sim \begin{cases} \textrm{Poisson}\left(\frac{\lambda}{1-\alpha}\right) &\quad \mbox{con probabilida } 1-\omega,\\ \textrm{Poisson}\left(\frac{q\lambda}{1-\alpha}\right) &\quad \mbox{con probabilidad } \omega.\nonumber \end{cases} \end{align*}

Cuando q=0 la distribución de \{Y_n\} Poisson cero-inflada. A partir de la mezcla podemos tener estimaciones preliminares para los parámetros que luego se usarán en un cálculo de los estimadores por máxima verosimilitud. La verosimilitud de Y es muy complicada,

    \[P(Y)=\sum_XP(X,Y)=\sum_xP(Y|X=x)P(X=x)\]

y para calcularla se usa el forward algorithm ~\cite{lystig02} muy usado en el contexto de modelos de Markov escondidos. Las  probabilidades de avance (forward)  son

(3)   \begin{align*} \alpha_k(X_k)=P(Y_k|X_k)\sum_{X_{k-1}}P(X_k|X_{k-1})\alpha_{k-1}(X_{k-1}) \nonumber, \end{align*}

con \alpha_1(X_1)=P(X_1)P(Y_1|X_1). Así que la verosimilitud es

    \[P(Y)=\sum_{n}\alpha_n(X_n).\]

P(Y_k|X_k) y P(X_k|X_{k-1}) son las llamadas probabilidades de emisión y de transición. Las probabilidades de transción se calculan así:

(4)   \begin{align*} P(X_n=x_n \mid X_{n-1}=x_{n-1})=e^{-\lambda}\sum_{j=0}^{x_n\wedge x_{n-1}}\binom{x_{n-1}}{j}\alpha^j(1-\alpha)^{x_{n-1}-j}\frac{\lambda^{x_n-j}}{(x_n-j)!} \nonumber \end{align*}

y las de emisión

(5)   \begin{align*} P(Y_i=j \mid X_i=k) & = \left\{ \begin{array}{lcc} 0 & if & k < j \\ (1-\omega) + \omega q^k& if & k=j\\ \omega \binom{k}{j} q^j (1-q)^{k-j}& if & k > j, \end{array} \right. \nonumber \end{align*}

Reconstrucción de la cadena oculta X_n

Para reconstruir la serie oculta X_n, se usa el algoritmo de \textbf{Viterbi} [10] descrito anteriormente.

La idea es reconstruir la cadena latente X_1^*=x_1^*,\dots, X_N^*=x_N^* que maximice la probabilidad del proceso latente dadas las observaciones Y_n, suponiendo que los parámetros son conocidos.

Bibliografía

[2] Bernard H,Werber D, Hoehle M. Estimating the under-reporting of norovirus illness in Germany utilizing enhanced awareness of diarrhoea during a large outbreak of Shiga toxin-producing E. coli O104: H4 in 2011-a time series analysis. BMC
Infectious Diseases 2014; 14:116. DOI:10.1186/1471-2334-14-116.

[3] Fan, C. et al. Prediction of Epidemic Spread of the 2019 Novel Coronavirus Driven by Spring Festival Transportation in China: A Population-Based Study. Int. J. Environ. Res. Public Health 17, (2020).

[4] Fernández-Fontelo A., Cabaña A., Puig P., Moriña D. Under-reported data analysis with INAR-hidden Markov chains
Statistics in Medicine (2016), vol. 35, Issue 26, 4875-4890.

[5] Fernández‐Fontelo, A., Cabaña, A., Joe, H., Puig, P. & Moriña, D. Untangling serially dependent underreported count data for gender‐based violence. Stat. Med. 38, 4404–4422 (2019).

[6] Moriña, D., Fernández-Fontelo, A., Cabaña, A.  & Puig, P. New statistical model for misreported data with application to current public health challenges. Submitted to Statistical Methods in Medical Research, (2020).

[7] Moriña, D., Fernández-Fontelo, A., Cabaña, A., Puig, P., Monfil, L., Brotons, M. & Diaz, M. Quantifying the underreporting of genital warts cases. Submitted to the European Journal of Epidemiology (2020).

[8] Kermack W. O., McKendrick A. G. (1927). A Contribution to the Mathematical Theory of Epidemics. Proceedings of the Royal Society A. 115 (772): 700–721.

[9] T.C. Lystig, J.P. Hughes (2002),  Exact computation of the observed information matrix for hidden Markov models, J. of Comp.and Graph. Stat.

[10] Viterbi, A.J. (1967), Error bounds for convolutional codes and an asymptotically optimum decoding algorithm. IEEE Transactions on Information Theory,  13, 260D269.

Comments are closed.

  • Relevant Posts

    • How Coronavirus Spreads
      Person-to-person spread The virus is thought to spread mainly from person-to-person. Between people who are in close contact with one another (within about 6 feet). Through respiratory droplets […]