ZIP Fit Indices Calculated from an EM Algorithm I am working through @ben-bolker's owls example available here:https://groups.nceas.ucsb.edu/non-linear-modeling/projects/owls/WRITEUP/owls.pdf
In particular, I am making use of the R zipme function which is defined separately as follows: 
zipme <- function(cformula, zformula, cfamily=poisson,
                  data, maxitr=20, tol=1e-6, verbose=TRUE) {
  #####################################
  ## EM algorithm for fitting ZIP mixed-effects model
  ##
  ##   y is the observation from the distribution:
  ##           P(Y=0)=p+(1-p)F(0,lambda)
  ##           P(Y=k)=(1-p)F(k,lambda).
  ##
  ##   data : the owl data frame with covariates; assumes data
  ## already pre-processed per pg 333 of Zuur et al 2009;
  ##        column order/names: Nest, FoodTreatment, SexParent, ArrivalTime,## NCalls, BroodSize, NegPerChick,
  ##        logbrdsze; logbrdsze is log(BroodSize).
  ##
  ##   formlog : formula for logistic regression. left side should be: z~
  ##   formpoi : formula for Poisson or NB regression. left side should be: y~
  ##
  ##   maxitr  : maximum number of iterations
  ##
  ## 2011.3.14 modified from Mihoko's GAMZINB to run ZIP mixed-effect model
  ##
  #############
  # number of observations
  m<-nrow(data)
  rname <- as.character(cformula)[2]

  ## initialize z and probz (z=1 -> perfect state; probz is probability of 0 in imperfect state for poisson)

  z<-numeric(m)
  probz<-numeric(m)
  z[data[[rname]]==0]<- 1/(1+exp(-1))  ## starting value

  ## n.b. we are looking for [3] since zformula has a LHS
  randz <- length(grep("\\(.*\\|.*\\)",as.character(zformula)[3]))>0
  ## delta is used to gauge convergence. after initialization, it is the abs. difference between current z and new z.    
  itr <- 1
  delta <- 2
  deltainfo <- numeric(maxitr)
  while(delta>tol & itr <= maxitr){
    if (verbose) cat("itr:",itr,"\n")
    ## make (update) working data frame
    bydataw <- data.frame(z=z,data)
    ##
    ## Maximization 1: logistic
    old.z<-z
    if (randz) {
      uu <- glmer(zformula, family=binomial, data=bydataw)
    } else {
      ## suppress warnings 
      uu <- suppressWarnings(glm(zformula, family=binomial, data=bydataw))
    }
    ## save current logistic model output
    u <- fitted(uu)
    ##
    ## Maximization 2: poisson loglinear with weights
    vv <- glmer(cformula, family=cfamily, weights=(1-z), data=bydataw)   
    ## save Poisson model output
    v <- fitted(vv)
    ##
    ## Expectation: used to update z with conditional expectation;only need to update at y=0.
    zdat <- data[[rname]]==0
    z[zdat] <- u[zdat]/( u[zdat]+(1-u[zdat])*exp(-v[zdat]))
    new.z<-z
    ## updated convergence indicator
    delta<-max(abs(old.z-new.z))
    ## save delta for this iteration; to be output
    deltainfo[itr] <- delta
    itr <- itr+1
  }            
  L <- list("zfit"=uu, "cfit"=vv, itr=itr, deltainfo=deltainfo, z=z)
  ##    uu.binom : output object of logistic regression; 
  ##    vv.flm   : output object of poisson regression
  class(L) <- "zipme"
  L
}

In essence, the function is an EM wrapper for fitting a ZIP mixed-effects model. The output of the function is, among other elements related to the maximization, two model objects - a mixed or fixed effect binomial GLM (for the zero component of the ZIP model) and a mixed or random effect Poisson GLM (for the Poisson component). 
I would like to compare the output of zipme (m2 below) to other models of the same data - a regular Poisson in particular (m1 below).
For example, using Bolker's Owls data...
library(lme4)
download.file("https://groups.nceas.ucsb.edu/non-linear-modeling/projects/owls/DATA/Owls.rda", destfile = "Owls.rda")
load("Owls.rda")

Owls$NCalls <- Owls$SiblingNegotiation


Owls <- transform(Owls,ArrivalTime=scale(ArrivalTime,center=TRUE,scale=FALSE))


m1 <- glmer(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
              offset(logBroodSize)+(1|Nest), family = "poisson", data=Owls)

m2 <- zipme(cformula=NCalls~(FoodTreatment+ArrivalTime)*SexParent+
        offset(logBroodSize)+(1|Nest),
      zformula=z ~ 1,
      data=Owls,maxitr=20,tol=1e-6,
      verbose=FALSE)

BIC(m1)
##[1] 5040.163

BIC(m2$cfit)
##[1] 3373.52

Question: Is it ok for me to compare the output of BIC(m1) with the output of BIC(m2$cfit) or should I be calculating a summary fit index which includes both the binomial and Poisson components to compare back to m1?
 A: Thanks to @gregor and @ben-bolker for getting me thinking along the following lines (@gregor via email and @ben-bolker above).
@ben-bolker noted that I should start from scratch with my own likelihood and @gregor noted 

My instinct is that adding the BIC's from each model component together sounds pretty good - and pretty easy. ZI model equations (and I think likelihoods) are two models multiplied together, so added together when logged... with correlations it might not be exactly that but seems like an acceptable estimate. 

Thinking about this more, I think that this is not only an acceptable estimate, but indeed the correct method for calculating BIC. 
Recall that the ZIP predicts $Y_i$ as 
$$
Y_i=\left\{
                \begin{array}{ll}
 0 \text{ with probability } \pi_i + (1 - \pi_i) e^{-\lambda_i} \\
 h \text{ with probability } (1 - \pi_i) \frac{\lambda^h_i e^{-\lambda_i}} {h_i!}.
                 \end{array}
              \right.
$$
The likelihood for this model is thus set up in two separate components: 
$$
\begin{align}
 \mathcal{L_1} &= \prod_{y_i=0} \pi_i + (1 - \pi_i) e^{-\lambda_i},\qquad y_i = 0 \\
 \mathcal{L_2} &= \prod_{y_i>0} (1 - \pi) \frac{\lambda_i^{y_i} e^{-\lambda_i}} {y_i!},\qquad y_i \ge 1. 
\end{align}
$$
The total log-likelihood, $\ln\mathcal{L_T}$, is then simply given as 
$$
\ln\mathcal{L_T} = \ln\mathcal{L_1} + \ln\mathcal{L_2}.
$$
Returning to BIC, we have the formula  
$$
 \mathrm{BIC} = {-2 \cdot \ln{\mathcal{L_T}} + k \cdot \ln(n)}
$$
where $k$ is the number of estimated parameters and $n$ is the number of observations in our model. 
We can thus calculate the BIC for a ZIP, as 
$$
\begin{align}
 \mathrm{BIC_{ZIP}} &= {-2 \cdot \ln{\mathcal{L_T}} + (k_1 + k_2) \cdot \ln(n)} \\
&= (-2 \cdot \ln\mathcal{L_1} + (k_1) \cdot \ln(n)) + (-2 \cdot \ln\mathcal{L_2} + (k_2) \cdot \ln(n)) \\
&= BIC_1 + BIC_2
\end{align}
$$
where the subscripts $1$ and $2$ continue to indicate information from the Binomial and Poisson components of the ZIP respectively. 
Unless I'm missing something, it appears to me that $gregor's intuition is correct. 
