Como você usa o algoritmo EM para calcular MLEs para uma formulação de variável latente de um modelo de Poisson inflado com zero?

O modelo de regressão de Poisson inflado com zero é definido para uma amostra por e assume ainda que os parâmetros e atendem$(y_1,\ldots,y_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}. }$$

A probabilidade de log correspondente do modelo de regressão de Poisson inflado a zero é

$$\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}!)}$$

Aqui, e são as matrizes de design. Essas matrizes podem ser as mesmas, dependendo dos recursos que se deseja usar para os dois processos geradores. Eles têm o mesmo número de linhas, no entanto.$\mathrm{B}$$\mathrm{G}$

Supondo que pudéssemos observar quando é do estado perfeito, zero e quando é do estado de Poisson, a probabilidade de log seria$Z_i = 1$$Y_i$$Z_i = 0$$Y_i$

$$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}!)]$$ Os dois primeiros termos são a perda em uma regressão logística para separar de . O segundo termo é uma regressão aos pontos gerados pelo processo de Poisson.$z_i=0$$z_i=1$

Mas as variáveis ​​latentes não são observáveis? O objetivo é maximizar a primeira probabilidade de log. Mas temos que introduzir variáveis ​​latentes e derivar uma nova probabilidade de log. Em seguida, usando o algoritmo EM, podemos maximizar a segunda probabilidade de log. Mas isso pressupõe que sabemos que ou ?$Z_i = 0$$Z_i = 1$

A raiz da dificuldade que você está tendo está na frase:

Em seguida, usando o algoritmo EM, podemos maximizar a segunda probabilidade de log.

Como você observou, você não pode. Em vez disso, o que você maximiza é o valor esperado da segunda probabilidade de log (conhecida como "probabilidade completa do log de dados"), em que o valor esperado é assumido sobre o . $z_i$

Isso leva a um procedimento iterativo, onde na iteração você calcula os valores esperados de considerando as estimativas de parâmetro da iteração ( (conhecida como "E-step ",) substitua-os na probabilidade completa do log de dados (consulte EDIT abaixo para saber por que podemos fazer isso neste caso) e maximize isso em relação aos parâmetros para obter as estimativas para a iteração atual (a" etapa M " .)$k^{th}$$z_i$$(k-1)^{th}$

O log probabilidade-de dados completa para o Poisson zero inflado no caso mais simples - dois parâmetros, dizer e - permite simplificar de forma significativa quando se trata de o M-passo, e isso leva por cima de alguma forma ao seu formulário. Vou mostrar como isso funciona no caso simples, através de um código R, para que você possa ver a essência dele. Não simplificarei o máximo possível, pois isso pode causar uma perda de clareza quando você pensa no seu 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 

No seu caso, em cada etapa, você fará uma regressão de Poisson ponderada em que os pesos devem 1-zhatobter as estimativas de e, portanto, e, em seguida, maximizar:$\beta$$\lambda_i$

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

com relação ao vetor de coeficiente da sua matriz para obter as estimativas de . Os valores esperados , calculados novamente a cada iteração.$\mathbf{G}$$p_i$$\mathbb{E}z_i = p_i/(p_i+(1-p_i)\exp{(-\lambda_i)})$

Se você deseja fazer isso com dados reais, em vez de apenas entender o algoritmo, já existem pacotes R; aqui está um exemplo http://www.ats.ucla.edu/stat/r/dae/zipoisson.htm usando a psclbiblioteca.

EDIT: Devo enfatizar que o que estamos fazendo é maximizar o valor esperado da probabilidade do log de dados completos, NÃO maximizar a probabilidade do log de dados completos com os valores esperados dos dados ausentes / variáveis ​​latentes plugados. Como acontece, se a probabilidade do log de dados completos é linear nos dados ausentes, como aqui, as duas abordagens são iguais, mas, caso contrário, não são.