¿Cómo utiliza el algoritmo EM para calcular los MLE para una formulación variable latente de un modelo de Poisson inflado a cero?

El modelo de regresión de Poisson inflado a cero se define para una muestra por y además supone que los parámetros \ mathbf {\ lambda} = (\ lambda_1, \ dots, \ lambda_n) y \ textbf {p} = (p_1, \ dots, p_n) satisfacen$(y_1,\ldots,y_n)$λ = ( λ 1 , ... , λ n

$$Y_i = \begin{cases} 0 & \text{with probability} \ p_i+(1-p_i)e^{-\lambda_i}\\ k & \text{with probability} \ (1-p_i)e^{-\lambda_i} \lambda_{i}^{k}/k! \end{cases}$$ $\mathbf{\lambda} = (\lambda_1, \dots, \lambda_n)$$\textbf{p} = (p_1, \dots, p_n)$

$$\eqalign{ \log(\mathbf{\lambda}) &= \textbf{B} \beta \\ \text{logit}(\textbf{p}) &= \log(\textbf{p}/(1-\textbf{p})) = \textbf{G} \mathbf{\gamma}. }$$

La probabilidad de registro correspondiente del modelo de regresión de Poisson inflado a cero es

$$\eqalign{ L(\gamma,\mathbf{\beta}; \mathbf{y}) &= \sum_{y_i=0} \log(e^{G_i \gamma}+\exp(-e^{\textbf{B}_i \mathbf{\beta}})) +\sum_{y_i >0} (y_i \textbf{B}_i \mathbf{\beta}-e^{\textbf{B}_i \mathbf{\beta}})\\ &\quad -\sum_{i=1}^{n} \log(1+e^{G_{i} \gamma})-\sum_{y_i >0} \log(y_{i}!)}$$

Aquí, $\mathrm{B}$ y $\mathrm{G}$ son las matrices de diseño. Estas matrices podrían ser las mismas, dependiendo de las características que uno desee usar para los dos procesos de generación. Sin embargo, tienen el mismo número de filas.

Suponiendo que pudiéramos observar $Z_i = 1$ cuando $Y_i$ es del estado perfecto, cero y $Z_i = 0$ cuando $Y_i$ es del estado de Poisson, la probabilidad logarítmica sería

$$L(\gamma,\mathbf{\beta}; \mathbf{y}, \mathbf{z}) = \sum_{i=1}^{n} \log(f(z_i|\mathbf{\gamma}))+\sum_{i=1}^{n} \log(f(y_i|z_i, \mathbf{\beta}))$$

$$= \sum_{i=1}^{n} z_{i} (\textbf{G}_i \gamma-\log(1+e^{G_{i} \gamma}))+ -\sum_{i=1}^{n} (1-z_{i})\log(1+e^{G_{i} \gamma})+ \sum_{i=1}^{n} (1-z_i)[y_{i} \textbf{B}_i \beta-e^{\textbf{B}_i \beta} - \log(y_{i}!)]$$ Los dos primeros términos son la pérdida en una regresión logística para separar $z_i=0$ de $z_i=1$ . El segundo término es una regresión a los puntos generados por el proceso de Poisson.

¿Pero no son las variables latentes inobservables? El propósito es maximizar la primera probabilidad de registro. Pero tenemos que introducir variables latentes y derivar un nuevo log-verosimilitud. Luego, utilizando el algoritmo EM, podemos maximizar la segunda probabilidad de registro. ¿Pero esto supone que sabemos que o ?$Z_i = 0$$Z_i = 1$

La raíz de la dificultad que tienes está en la oración:

Luego, utilizando el algoritmo EM, podemos maximizar la segunda probabilidad de registro.

Como has observado, no puedes. En cambio, lo que maximiza es el valor esperado de la segunda probabilidad de registro (conocida como la "probabilidad de registro de datos completa"), donde el valor esperado se toma sobre . $z_i$

Esto lleva a un procedimiento iterativo, donde en la iteración se calculan los valores esperados de dadas las estimaciones de los parámetros de la iteración ( (esto se conoce como el "paso E ",) luego sustitúyalos en la probabilidad de registro de datos completa (consulte EDITAR a continuación para ver por qué podemos hacer esto en este caso), y maximícela con respecto a los parámetros para obtener las estimaciones para la iteración actual (el" paso M " .)$k^{th}$$z_i$$(k-1)^{th}$

La probabilidad log-datos completa para la Poisson inflado de cero en el caso más simple - dos parámetros, por ejemplo y - permite simplificar de manera importante cuando se trata de la M-paso, y esto se traslada en cierta medida a su formulario. Le mostraré cómo funciona eso en el caso simple a través de un código R, para que pueda ver la esencia del mismo. No simplificaré tanto como sea posible, ya que eso podría causar una pérdida de claridad cuando piense en su problema:$\lambda$$p$

# Generate data
# Lambda = 1,  p(zero) = 0.1
x <- rpois(10000,1)
x[1:1000] <- 0

# Sufficient statistic for the ZIP
sum.x <- sum(x)

# (Poor) starting values for parameter estimates
phat <- 0.5
lhat <- 2.0

zhat <- rep(0,length(x))
for (i in 1:100) {
  # zhat[x>0] <- 0 always, so no need to make the assignment at every iteration
  zhat[x==0] <- phat/(phat +  (1-phat)*exp(-lhat))

  lhat <- sum.x/sum(1-zhat) # in effect, removing E(# zeroes due to z=1)
  phat <- mean(zhat)   

  cat("Iteration: ",i, "  lhat: ",lhat, "  phat: ", phat,"\n")
}

Iteration:  1   lhat:  1.443948   phat:  0.3792712 
Iteration:  2   lhat:  1.300164   phat:  0.3106252 
Iteration:  3   lhat:  1.225007   phat:  0.268331 
...
Iteration:  99   lhat:  0.9883329   phat:  0.09311933 
Iteration:  100   lhat:  0.9883194   phat:  0.09310694 

En su caso, en cada paso hará una regresión ponderada de Poisson donde los pesos son 1-zhatpara obtener las estimaciones de y, por lo tanto, , y luego maximizar:$\beta$$\lambda_i$

$\sum (\mathbb{E}z_i\log{p_i} + (1-\mathbb{E}z_i)\log{(1-p_i)})$

con respecto al vector coeficiente de su matriz para obtener las estimaciones de . Los valores esperados , nuevamente calculados en cada iteración.$\mathbf{G}$$p_i$$\mathbb{E}z_i = p_i/(p_i+(1-p_i)\exp{(-\lambda_i)})$

Si desea hacer esto para datos reales, en lugar de simplemente comprender el algoritmo, los paquetes R ya existen; Aquí hay un ejemplo http://www.ats.ucla.edu/stat/r/dae/zipoisson.htm usando la psclbiblioteca.

EDITAR: Debo enfatizar que lo que estamos haciendo es maximizar el valor esperado de la probabilidad de registro de datos completos, NO maximizar la probabilidad de registro de datos completos con los valores esperados de los datos faltantes / variables latentes conectadas. Como sucede, si la probabilidad de registro de datos completos es lineal en los datos faltantes, ya que está aquí, los dos enfoques son iguales, pero de lo contrario, no lo son.