Hierarchical model: does leaving out a latent variable (hierarchy level) result in an equivalent model? Say we have a hierarchical model:
$$z_i \sim \mbox{Bernoulli}(\pi_i); \mbox{logit}(\pi_i) =  ... \text{(linear function of covariates for site i)}$$
$$y_{i,j} \sim \mbox{Bernoulli}(z_i \cdot p_{i,j}) ; \mbox{logit}(p_{i,j}) =  ... \text{(linear function of covariates for visit <i,j>)}$$
where $z_i$ is a latent variable, saying whether bird is present at site i ($z_i = 1$) or absent ($z_i = 0$); $\pi_i$ is a probability of bird being present at site $i$, which is a function of some covariates. $y_{i,j}$ is the actual observation of the bird at site $i$ during visit $j$; $p_{i,j}$ being the probability of detection of the bird which is present, i.e. $p_{i,j} = \mbox{P}(y_{i,j} = 1 | z_i = 1)$.
Because the only input data are the $y_{i,j}$ and we actually don't know $z_i$, I thought we could simplify (flatten) the model by removing the latent variable $z$ and having just:
$$y_{i,j} \sim \mbox{Bernoulli}(\pi_i \cdot p_{i,j})$$
The question is: Is this model equivalent to the first one? $y$-prediction-wise, coefficient-wise, everything-wise (of course everything except having $z_i$'s :-) If no, where is the difference?
My thoughts are:


*

*maybe it should be equivalent, because $z$ is latent, so it should be "integrated out" :-)

*maybe it shouldn't, because $z_i$ could actually transfer information between $y_{i,j}$ of the same site $i$. For example, when $y_{i,1} = 1$, the model will immediatelly know $z_i = 1$ without question, so this could somehow influence the inference for nodes $y_{i,2}$, $y_{i,3}$ etc... 


Just thoughts, I'm very likely mistaken somewhere... I'd be grateful for some thoughts on the matter :-)
Thanks and have a beautiful Christmass time! :)
 A: Let's have a thourough look at the likelihood function $\cal{L}(y_{ij}|\pi_i,p_{ij})$ of both models.
The likelihood function of the second model is simple:
$${\cal L}(y_{ij}|\pi_i,p_{ij}) = \mbox{Bernoulli}(y_{ij}|\pi_i \cdot p_{ij})$$
The likelihood function of the first model is the integration (summation) over all possible values of the latent variable $z$:
$$\begin{eqnarray}{\cal L}(y_{ij}|\pi_i,p_{ij}) &=& \sum_{z_i \in \{0,1\}} {\cal L}(y_{ij}|z_i,p_{ij}) \cdot {\cal L}(z_i|\pi_{i}) = \\ &=& {\cal L}(y_{ij}|z_i= 0,p_{ij}) \cdot {\cal L}(z_i=0|\pi_{i}) + {\cal L}(y_{ij}|z_i=1,p_{ij}) \cdot {\cal L}(z_i=1|\pi_{i}) = \\ &=& \mbox{Bernoulli}(y_{ij}|0)\cdot(1-\pi_i) + \mbox{Bernoulli}(y_{ij}|p_{ij})\cdot\pi_i\end{eqnarray}$$
This gives us:
$${\cal L}(y_{ij}|\pi_i,p_{ij}) = \begin{cases} 
y_{ij} = 0: & 1\cdot(1-\pi_i)+(1-p_{ij})\cdot\pi_i =1 - p_{ij}\cdot\pi_i \\ 
y_{ij} = 1: & 0 + p_{ij}\cdot\pi_i
\end{cases}$$
Which is exactly the same likelihood function as in the first model, right? This proves the models are fully equivalent, right? Q.E.D., right?
No.
The models need the joint likelihood to be equivalent (thanks @JimB for comment)!
So, the joint likelihood of the second model is again simple:
$${\cal L}(y_{i*}|\pi_i,p_{i*}) = \prod_{j} {\cal L}(y_{ij}|\pi_i,p_{ij}) = \prod_{j} \mbox{Bernoulli}(y_{ij}|\pi_i \cdot p_{ij})$$
For the first model, again, lets do the summation over all possible values of $z$:
$$\begin{eqnarray}{\cal L}(y_{i*}|\pi_i,p_{i*}) &=& \sum_{z_i \in \{0,1\}} {\cal L}(z_i|\pi_{i}) \prod_{j} {\cal L}(y_{ij}|z_i,p_{ij})  = \\ &=& {\cal L}(z_i=0|\pi_{i}) \prod_{j} {\cal L}(y_{ij}|z_i= 0,p_{ij})  + {\cal L}(z_i=1|\pi_{i}) \prod_{j} {\cal L}(y_{ij}|z_i=1,p_{ij})  = \\ &=& (1-\pi_i) \prod_{j} \mbox{Bernoulli}(y_{ij}|0) + \pi_i \prod_{j} \mbox{Bernoulli}(y_{ij}|p_{ij})\end{eqnarray}$$
It can be easily seen that this is different from the second model; here, we have only singular power of $\pi_i$, whereas in the second model, $\pi_i$ will reach the power of how many $j$ indices we have for the single $i$ index.
So the answer is: the models are not equivalent.
A: If one starts out modeling the joint probability (or equivalently the likelihood) of the status of all visits and mentions the appropriate assumptions, then one can certainly avoid talk of latent variables and construct the likelihood (although the use of latent variables is implied).
Using much of the OP's notation we assume independence of visit status between sites but because the true status of a site is assumed not to change during all visits (the "closure" assumption), a dependence among visits at site $i$ is induced:  either all visits will be zero (at an unoccupied site) or the visit status for the $v_i$ visits will vary independently according to a probability $p_{ij}$ for visits $j=1,2,\ldots,v_i$.
We also assume that there are no false positive visits and that a site with no detections could either be an unoccupied site or an occupied site.
The likelihood at a site for observed visit status $y_{i,1},y_{i,2},\ldots,y_{i,v_i}$ is given by
$$\mathcal{L}_i(y_{i,1},y_{i,2},\ldots,y_{i,v_i})=1-\pi_i+\pi_i \prod_{j=1}^{v_i}(1-p_{ij}) ~~~\text{if}~~\sum_{j=1}^{v_i} y_{ij}=0$$
$$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~=\pi_i \prod_{j=1}^{v_i}p_{ij}^{y_{ij}}(1-p_{ij})^{1-y_{ij}} ~~~~~~\text{if}~~\sum_{j=1}^{v_i} y_{ij}>0$$
A: It does seem that you've answered your own question. In your example with y[i,1]=1, occupancy (z) is known for all replicate visits to site i, so the remaining visits to the site (y[i,2...n.reps]) provide information about detection probability. It is a hierarchical model with site occupancy being estimated first, then detection being estimated within the site conditional on the site being occupied or not). If your model only uses Bernoulli(site occupancy X detection probability), the model might estimate that site i was occupied on visit 1 but not on visits 2:n.reps -- which is not consistent with the closure assumption that makes the model work.
