N_c-mixture occupancy model

Consider a multiple-visit sampling survey consisting of $n$ sites and $T$ visits. Let $N_{ij}$ be the number of individuals at site $i$ during the $j$-th visit. Following the motivation in the Introduction (see Figure 1), we decompose $N_{ij}$ into two independent Poisson distributed random variables: $K_{i}$ and $M_{ij}$. The expected value of $K_{i}$ is $c\mu$ and the expected value of $M_{ij}$ is $(1-c)\mu$, where $0\leq c\leq 1$ and $\mu>0$ are constant parameters. Therefore, the number of individuals at each site, $N_{ij}$, follows an identical Poisson distribution with mean $\mu$, representing the abundance parameter over sites (per each visit). However, the numbers of individuals from multiple visits at a site, $N_{ij},j=1,\cdots,T$, are not independent as they share a common variable $K_{i}$. The community parameter $c$ characterizes the degree of dependence between visits, as it also represents the correlation between each visit of $N_{ij}$. Note that $M_{ij}$ and $K_{i}$ are degenerate to 0 when $c=1$ and $c=0$, respectively.

To model the data observation process, we assume that each individual is independently detectable with a probability of detection $r$. If the species is detected at site $i$ during visit $j$, let $Y_{ij}=1$, otherwise $Y_{ij}=0$. We denote the probability of detection at site $i$ during visit $j$ as $P(Y_{ij}=1)=p_{ij}$, then $p_{ij}=1-(1-r)^{N_{ij}}$ when $N_{ij}$ is known. This forms an N${}_{c}-$mixture model.

It is clear that $Y_{ij},j=1,\cdots,T$, are exchangeable variables, though not independent unless $c=0$. Let $Y_{i}=\sum_{j=1}^{T}Y_{ij}$ be the frequency of occurrence at site $i$. In Web Appendix A, we show that the probability function of $Y_{i}$ under the N${}_{c}-$mixture model is given by

$$f(y_{i};\bm{\theta})={T\choose y_{i}}\exp(-d\mu rT-c\mu)\sum_{k=0}^{y_{i}}{y_{i}\choose k}(-1)^{k}\exp\left\{c\mu{(1-r)}^{T-y_{i}+k}+d\mu r(y_{i}-k)\right\},$$

(1)

where $\bm{\theta}=(\mu,r,c)$ is the vector of model parameters and $d=1-c$. Equation (1) is also referred to as the N${}_{c}-$mixture model.

When the community parameter $c=1$, the model in equation (1) simplifies to

$$f(y_{i};\bm{\theta}_{1})={T\choose y_{i}}\exp(-\mu)\sum_{k=0}^{y_{i}}{y_{i}\choose k}(-1)^{k}\exp\left\{\mu{(1-r)}^{T-y_{i}+k}\right\},$$

(2)

where $\bm{\theta}_{1}=\bm{\theta}$ with $c=1$. This is the probability function of the N$-$mixture model (Royle and Nichols, 2003), but the explicit form (2) is firstly given in Haines (2016a).

Similarly, when the community parameter $c=0$, the model (1) reduces to

$$f(y_{i};\bm{\theta}_{0})={T\choose y_{i}}\left\{\exp(-\mu r)\right\}^{T-y_{i}}\left\{1-\exp(-\mu r)\right\}^{y_{i}},$$

(3)

where $\bm{\theta}_{0}=\bm{\theta}$ with $c=0$. The reduced model (3) is a binomial distribution denoted as $Y_{i}\sim\text{Bino}(T,p)$, where $p=1-\exp(-\mu r)$. We also note that the binary variables $Y_{ij}$ are independently and identically distributed with a mean of $E\{1-(1-r)^{N_{ij}}\}$ when $c=0$. A calculation using the Poisson moment generating function shows that $E\{1-(1-r)^{N_{ij}}\}=1-\exp(-\mu r)=p$, which leads to the same result as (3). However, the model parameters $\mu$ and $r$ in (3) cannot be separated, making the model unidentifiable, and only their product $\mu\times r$ is identifiable.

When the community parameter $0<c<1$, the N${}_{c}-$mixture model (1) is identifiable for $T\geq 3$. The log-likelihood function of the parameter vector $\bm{\theta}$ is given by $\ell(\bm{\theta})=\sum\nolimits_{i=1}^{n}\log\{f(y_{i};\bm{\theta})\}$, and the maximum likelihood estimation is straightforward. For later use, the likelihood function can be represented in the framework of the multinomial model. Let $m_{j}=\sum_{i=1}^{n}\mathrm{\bf I}(y_{i}=j)$ for $j=0,\cdots,T$, where $\mathrm{\bf I}(\cdot)$ is the indicator function. These statistics $m_{j}$ can be viewed as a result of the multinomial model with cell probabilities $f(j;\bm{\theta})$. Therefore, we can write $\ell(\bm{\theta})=\sum_{j=0}^{T}m_{j}\log\left\{f(j;\bm{\theta})\right\}$, and the score function for $\bm{\theta}$ is

$$S(\bm{\theta})=\sum\limits_{j=0}^{T}\frac{\partial f(j;\bm{\theta})}{\partial\bm{\theta}}\left\{\frac{m_{j}-nf(j;\bm{\theta})}{f(j;\bm{\theta})}\right\}.$$

(4)

Royle and Nichols (2003) define the occupancy rate $\psi$ as a derived parameter under the N$-$mixture model. Specifically, $\psi=P(K_{i}>0)=1-\exp(-\mu)$ in terms of our notations. Under the N${}_{c}-$mixture framework, the occupancy rate per each visit can be defined as $P(N_{ij}>0)=1-\exp(-\mu)$. Alternatively, the rate can be defined as $1-P(N_{ij}=0,~{}\forall~{}1\leq j\leq T)$, where if some $N_{ij}>0$, the site $i$ is considered occupied by the species. In this way, the occupancy rate is $1-\exp\left\{-\mu-(T-1)d\mu\right\}$, which depends on the number of visits $T$ and converges to $1$ when $T$ increases infinitely. Both definitions, therefore, differ from the current concept of site occupancy. The problem is that the number of individuals at site $i$ may vary from visit to visit because $M_{ij},j=1,\ldots,T,$ are random in the N${}_{c}-$mixture model. Therefore, rather than directly defining occupancy, we would like to include a zero-inflated parameter to determine site occupancy under the N${}_{c}-$mixture model. This extension will be addressed in Section 3.

In this subsection, we examine the behavior of N$-$mixture model estimators for $\mu$ and $r$ when the number of individuals at a site during multiple visits is not a fixed constant (i.e., $c<1$).

First, we examine the scenario where the community parameter $c$ is known and $T=2$ (double-visit). In this scenario, we use the notation $\tilde{\mu}_{c}$ and $\tilde{r}_{c}$ to represent the (restricted) maximum likelihood estimators (MLEs) of $\mu$ and $r$, respectively. As shown in Web Lemma 1 of Web Appendix A, the MLEs have closed forms when $T=2$, which are given by $\tilde{\mu}_{c}=cz_{1}^{2}/(2z_{1}+z_{2})$ and $\tilde{r}_{c}=(2z_{1}+z_{2})/(cz_{1})$, where $z_{1}=\log\left\{2n/(2m_{0}+m_{1})\right\}$ and $z_{2}=\log(m_{0}/n)$. These closed forms allow us to easily determine the limits of the MLEs of the N$-$mixture model, $\tilde{\mu}_{1}$ and $\tilde{r}_{1}$, as given in the following proposition.

Proposition 1 delivers insights into the behavior of N$-$mixture model estimators when the community parameter is less than $1$. The results show that the MLEs are consistent when $c=1$, which is a common property of the maximum likelihood approach. Additionally, the results indicate that the N$-$mixture model overestimates the abundance parameter $\mu$, and the bias increases as the community parameter $c$ decreases. In contrast, the N$-$mixture model exhibits the opposite bias behavior for the detection probability $r$. Despite these biases, the estimator $\tilde{\mu}_{1}\times\tilde{r}_{1}$ can consistently estimate $\mu\times r$ under the framework of N${}_{c}-$mixture model, which is noteworthy.

We initially believed that the results of Proposition 1 would hold for surveys with more than two visits ($T>2$), but further investigation revealed that this is only partially true. The correct part is that there are moment estimators of the N$-$mixture model that exhibit similar behaviors to those described in Proposition 1. To demonstrate this, we derived the moments of $Y_{i}$ under the N$-$mixture model and defined the resulting estimators as $\tilde{\mu}_{1\rm{M}}$ and $\tilde{r}_{1\rm{M}}$. The following proposition summarizes the results of this analysis.

In the case of multiple-visit surveys, our simulation study (Section 4) found that the limit of the MLE of the N$-$mixture model, $\tilde{\mu}_{1}$, exhibits a similar pattern to $\mu/c$ when $0<c<1$. However, there is a discrepancy between them, and we can only determine the behavior when the community parameter $c$ is close to $0$.

Based on the results presented, we suggest that when $c<1$, the MLE of the N$-$mixture model tends to overestimate the parameter $\mu$ and underestimate the parameter $r$, and the bias increases as $c$ moves further away from $1$. However, the estimator $\tilde{\mu}_{1}\times\tilde{r}_{1}$ may still be able to consistently estimate $\mu\times r$ at certain ranges of $c$, such as when $c$ is close to 0, within the framework of the N${}_{c}-$mixture model.

When the community parameter $0<c<1$, the ZIN_c model is identifiable for $T\geq 4$. By defining $p_{0}({\mbox{\boldmath$\theta$}},\psi)=(1-\psi)+\psi f(0;\bm{\theta})$ and $f(+;\bm{\theta})={1-f(0;\bm{\theta})}$, the likelihood function can be written as

$$L(\bm{\theta},\psi)=L_{0}(\bm{\theta},\psi)\times L_{1}(\bm{\theta}),$$

where $L_{0}(\bm{\theta},\psi)=\{p_{0}({\mbox{\boldmath$\theta$}},\psi)\}^{m_{0}}\{1-p_{0}({\mbox{\boldmath$\theta$}},\psi)\}^{n-m_{0}}$ and $L_{1}(\bm{\theta})=\prod\limits_{j=1}^{T}\left\{f(j;\bm{\theta})/f(+;\bm{\theta})\right\}^{m_{j}}.$ Note that $L_{0}$ reflects the probability function of $\mathrm{\bf I}(Y_{i}>0)$ while $L_{1}$ is the conditional likelihood based on $Y_{i}>0$. As the conditional likelihood function is independent of the occupancy probability $\psi$, it allows us to estimate $\theta$ without the confounding of $\psi$ (Karavarsamis and Huggins, 2020). Specifically, we can find the conditional MLE of $\theta$, denoted as $\hat{\mbox{\boldmath$\theta$}}$, by solving the score function of $L_{1}({\mbox{\boldmath$\theta$}})$. Additionally, by taking the MLE of the profile likelihood $L_{0}(\hat{\mbox{\boldmath$\theta$}},\psi)$ we can find $\hat{\psi}=(n-m_{0})/\{nf(+;\hat{\mbox{\boldmath$\theta$}})\}$. As shown in Web Appendix A, Web Lemma 2 confirms that the estimators $\hat{\mbox{\boldmath$\theta$}}$ and $\hat{\psi}$ resulting from this method are also the usual MLEs based on (5). The asymptotic variances of $\hat{\mbox{\boldmath$\theta$}}$ and $\hat{\psi}$ can be derived in the usual way.

In the zero-inflated type models, the main focus is on the estimation of the occupancy probability. Let $\tilde{\psi}_{c}$ be the MLE of $\psi$ for the ZIN_c model, given that the community parameter $c$ is known. We next examine the behavior of $\tilde{\psi}_{0}$ and $\tilde{\psi}_{1}$, which correspond to the occupancy estimators for the ZIB and ZIN models, respectively.

Like the N${}_{c}-$mixture model, the ZIN_c model also has an identifiability issue for $\mu$ and $r$ when $c=0$. As a result, under the ZIB model, only the parameter $p$ (i.e., $1-\exp(-\mu r)$) and the occupancy probability $\psi$ can be estimated. In practice, to fit ZIB models, it is common to set $r=1$ to estimate $\mu\times r$ with the resulting abundance estimator ($\tilde{\mu}_{0}$).

We notice that as either $T$ or $\mu\times r$ increases, the value of $\Delta$ decreases to zero. This result is expected, as when there are more visits or when the species abundance is high, the observed occupancy approaches $\psi$. The linear approximation bias provides a reasonable representation of the trend of $\tilde{\psi}_{0}$ in various aspects, depicting that as $c$ moves away from $0$ (or as the correlation between visits increases), the underestimation becomes more significant.

Remark 1. As a direct consequence of Proposition 4, we can also see that the corresponding abundance estimator $\tilde{\mu}_{0}$ increases as $c$ increases, and that $\tilde{\mu}_{0}$ at $c=0$ is a consistent estimator of the product of species abundance and detection probability ($\mu\times r$).

Estimators of the ZIN model tend to have bias when the community parameter $c$ is less than $1$. However, the behavior of these biases is complex. For example, the bias of the occupancy estimator $\tilde{\psi}_{1}$ does not vary monotonically with decreasing $c$. Despite this, when $c\approx 0$, the estimators of the ZIN model (except $\tilde{\psi}_{1}$) behave similarly to those of the N$-$mixture model, as shown in Proposition 3. Interestingly, $\tilde{\psi}_{1}$ can consistently estimate $\psi$ at $c=0$, as shown in the next proposition. For clarity, if there is no confusion, $\tilde{\mu}_{1}$ and $\tilde{r}_{1}$ in this proposition also refer to the MLE of the ZIN model (or the restricted MLE of the ZIN_c model).

The estimator $\tilde{\psi}_{1}$ is also a consistent estimator under the ZIN (or ZIN_c with $c=1$) model. Therefore, Proposition 5 shows that $\tilde{\psi}_{1}$ behaves like a bridge with both ends (at $c=0$ and 1) at the same level; however, it is not clear whether the bridge deck is always above or below this level, or if it varies in different sections. We will further investigate this behavior through a simulation study.

The null hypothesis of the ZIN or ZIB model can be tested within the framework of the ZIN_c model, as it is equivalent to testing the values $\{c=1\}$ or $\{c=0\}$. A simple method to justify the hypothesis is to use a Wald-type confidence interval based on the estimate of $c$. A more formal approach is to use the likelihood ratio test, as the null hypothesis is a submodel of the full ZIN_c model. However, it is important to note that the asymptotic distribution of the likelihood ratio test under the null hypothesis is a mixture of $0$ and chi-square distributions, rather than the usual chi-square distribution (Self and Liang, 1987). In practice, we also suggest generating bootstrap samples under the null hypothesis to find the $p-$value. This is similar to the conclusion of Dail and Madsen (2011).

The median estimates of $\mu$ for $T=7$ and $n=500$ are displayed in Figure 2. A comprehensive examination of the simulation results for $c=(0.25,0.5,0.75)$ can be found in Web Tables 1-6 in Web Appendix B.

Figure 2 illustrates the behavior of the N$-$mixture abundance estimator $\tilde{\mu}_{1}$ as a function of the community parameter $c$. The figure shows that the estimator exhibits a positive bias for all values of $c$ less than $1$. This bias follows a monotonically decreasing pattern and only approaches the true value when $c$ equals $1$. When $c$ is close to $0$, the bias can be substantial but has been removed from the figure for ease of visualization. Additionally, the figure shows that $\tilde{\mu}_{1}$ is greater than $\mu/c$ for all values of $c$ less than 1, although the difference has not been explicitly explored. In a separate simulation study (not reported), we calculated the moment estimator for the N$-$mixture model and found that its pattern closely resembled the reference curve $\mu/c$. These simulation results agree with the findings outlined in Propositions 1-3. Corresponding conclusions can also be drawn from the results of the estimator $\tilde{r}_{1}$ shown in Web Figure 1 of Web Appendix B. The MLE of the N${}_{c}-$mixture model can consistently estimate $\mu$, but it also frequently exhibits bias around $c=0$. This is mainly because the parameters of $\mu$ and $r$ are almost unidentifiable when $c\approx 0$. The bias is more pronounced when $\mu$ is large and $r$ is small, but it becomes smaller when $n$ or $T$ increases; see Web Table 6 (for $n=1000$ and $T=10$ cases).

In Web Tables 1-6, it can be seen that the N$-$mixture abundance estimator has a relative bias ranging from $40\%$ to $400\%$ when $c$ varies from $0.75$ to $0.25$. The bias is more pronounced when $T$ or $r$ increases but less severe when $\mu$ increases. On the other hand, the N${}_{c}-$mixture model shows nearly unbiased estimates in all cases, except when $n=200,~{}T=5$, and $r=0.25$, where the relative bias can reach up to $7.5\%$-$16\%$ when $\mu$ varies from $2$ to $1$; see Web Tables 1 and 4.

In terms of mean absolute deviation (MAD), both models show a decrease in variation as $c$ increases, with the N$-$mixture model showing a much steeper decrease compared to the N${}_{c}-$mixture model. Specifically, when $c=0.25$, the MAD of $\tilde{\mu}_{1}$ can be $2$ to $4$ times that of $\hat{\mu}$, but the former is usually smaller than the latter when $c=0.75$. The asymptotic standard error estimates, as measured by Med.se, generally match the results of the corresponding MAD, making them reasonably reliable for the scenarios considered. The Wald-type confidence interval estimator of the N${}_{c}-$mixture model performs well when data information is sufficient, with a close match to the nominal $95\%$ confidence level at $n\geq 500$ and $T\geq 7$. However, in some cases with small $r$ and $c$, the coverage probability (CP) can be lower than $80\%$ (as seen in Web Tables 1, 4, and 5), indicating an unsatisfactory performance. The N$-$mixture model estimator often reaches $0$ for the coverage probability due to the severe bias problem of $\tilde{\mu}_{1}$ when $c\leq 0.75$.

The results from Web Tables 1-6 on the detection probability parameter $r$ are similar to $\mu$, with the only difference being that the N$-$mixture model estimator’s bias is in the opposite direction. Finally, it is found that the median of the product $\tilde{\mu}_{1}\times\tilde{r}_{1}$ presents nearly unbiased results for estimating $\mu\times r$ in all cases. Note that this property is only proven for the range $c\approx 0$ (Proposition 5), but the simulation results suggest that the range of $c$ can be extended somewhere.

In Figures 3 and 4, we present the median estimate results for the parameters $\mu$ and $\psi$ when $T=7$ and $n=500$. The detailed simulation results for $c=(0.25,0.5,0.75)$ can be found in Web Tables 7-12 of the Web Appendix B. In the ZIB model, we have fixed the value of $r=1$. This is because the parameter $\mu\times r$ is non-separable in this case.

In Figure 3, we can see that the ZIB estimator $\tilde{\mu}_{0}$ consistently underestimates $\mu$. In fact, $\tilde{\mu}_{0}$ is approximately equal to $\mu\times r$ when $c$ is close to 0 and increases as $c$ increases (as stated in Remark 1). The ZIN estimator $\tilde{\mu}_{1}$ has a similar trend as the N$-$mixture model when $c$ is close to $0$, but it falls off quickly, rises slowly, and becomes consistent with $\mu$ at $c=1$. In general, $\tilde{\mu}_{1}$ mainly exhibits positive bias but occasionally also shows negative bias at some $0<c<1$. The ZIN_c estimator $\hat{\mu}$ typically exhibits some bias at both ends of the range of $c$, which can be very large at $c\approx 0$ when $r$ is small. As expected, increasing $n$ and $T$ can mitigate bias issues of the maximum likelihood estimator.

In Figure 4, the ZIB estimator $\tilde{\psi}_{0}$ was found to underestimate $\psi$ when $c>0$, with a greater bias observed when $\mu=1$ compared to $\mu=2$. The relative bias reached a maximum of $-35\%$ when $c=1$ and $\mu=1$. The approximate formula of Proposition 4 was found to be consistent with the behavior trend of $\tilde{\psi}_{0}$, although it is not shown in the figure. The ZIN occupancy estimator $\tilde{\psi}_{1}$ was found to fit well with $\psi$ at both ends of the $c$ range, as stated in Proposition 5. A positive bias was observed around $c\approx 0$ and a negative bias around $c\approx 1$. The partial derivative of $\frac{\partial}{\partial c}\tilde{\psi}_{1}$ was proven to be positive when $c\approx 0$, indicating overestimation of $\psi$ in that region. However, at $c\approx 1$, the partial derivative was found to be positive in most situations, with some negative signs observed in some instances, such as when $\mu=2$ and $r=0.15$. In this case, the plot of $\tilde{\psi}_{1}$ against $c$ exhibited behavior similar to an arch bridge, with only positive biases observed for all $0<c<1$. The ZIN_c estimator $\hat{\psi}$ was found to be unbiased in general, except when $c\approx 1$. The bias was found to be smaller when $\mu$ or $r$ is increased.

Next, we present a more detailed summary of the performance of the three occupancy estimators based on the results in Web Tables 7-12. The upper half of Web Table 7 ($\mu\times r=0.25$ and $T=5$) shows that the ZIB occupancy estimator $\tilde{\psi}_{0}$ has a relative bias of $-11\%$ at $c=0.25$, which becomes more pronounced as $c$ increases, reaching a maximum of $-31\%$ at $c=0.75$. However, increasing $T$ or $\mu\times r$ can reduce the bias of $\tilde{\psi}_{0}$ as seen in Web Tables 7-12, which is consistent with Proposition 4. In Web Tables 7-9 ($\mu=1$), the ZIN occupancy estimator $\tilde{\psi}_{1}$ has a relative bias of around $28\%$-$40\%$ at $c=0.25$. Specifically, the ZIN model often estimated $\psi$ as one at $c=0.25$ when $T\geq 7$ and $r=0.5$ (Web Tables 8-9). The bias of $\tilde{\psi}_{1}$ decreases as $c$ increases, showing a small negative bias at $c=0.75$. When $\mu$ is increased to $\mu=2$, Web Tables 10-12 show that $\tilde{\psi}_{1}$ overestimates $\psi$ in all cases, with the most significant bias around $10\%$-$13\%$ at $c=0.5$, and small or even negligible bias at $c=0.25,0.75$. Web Tables 7-12 also reveal that the bias of the ZIN_c occupancy estimator $\hat{\psi}$ is generally close to $0$, except for a few cases of $r=0.25$, $T=5$, and $c=0.75$ (Web Table 7).

The ZIB occupancy estimator has the lowest MAD, indicating that it is the most stable of the three methods. When $\mu=1$, the MAD of the ZIN estimator decreases with increasing $c$; however, when $\mu=2$, the MAD is more prominent at $c=0.5$ due to the significant bias of the ZIN estimator at this point. In contrast, the MAD of the ZIN_c estimator increases with increasing $c$ values; when $\mu=1$, the MAD at $c=0.75$ can be $1.5$ to $3$ times the MAD at $c=0.25$, but the increase is much smaller when $\mu=2$. Generally, the stabilities of all three estimators improve with increasing values of $n,~{}T$, $r$, and $\mu$.

The Med.se of the ZIB occupancy estimator fits the corresponding MAD reasonably well, except for some negative biases observed at $c=0.75$. However, the resulting interval estimator is only reliable when $\mu=2$ and $c=0.25$; in other cases, the corresponding CP may even reach zero showing very poor performance. For the ZIN occupancy estimator, Med.se often overestimates MAD, particularly with increasing $c$, but this is less of an issue when $r$ or $\mu$ are increased. The resulting CP also shows many abnormal values, even when it appears to be close to the nominal level in some cases. For example, Web Tables 8-9 show that the CP can reach about $99\%$ in the upper panels but may drop to zero in the bottom panels. As a result, the interval estimates for the ZIN model are not considered reliable. The ZIN_c estimator, being a maximum likelihood estimate, can perform well on Med.se and CP measures when the data information is sufficient, but its performance is not guaranteed otherwise.

Lastly, we note that abundance estimators can be similarly affected by issues such as violation of model assumptions or a lack of sufficient data, which can result in substantial bias in some instances, such as when $c=0.25$ and $r=0.25$. In general, compared to abundance estimators, the corresponding occupancy estimators are relatively less affected by these issues and their performance is relatively robust.

The fisher $(Martes~{}pennanti)$ is a carnivorous mammal native to the boreal forests of North America. In 2000, a survey program using noninvasive methods was conducted to collect data on fisher species distribution in northern and central California (Zielinski et al., 2005; Royle and Dorazio, 2008). The data consists of multiple-visit occurrences at $n=464$ sites, each visited $T=8$ times. See Web Table 13 for the data. Out of the eight visits, $400$ sites had zero counts, resulting in a sample occupancy rate of $64/464=13.8\%$. We analyzed the fisher data using the models presented in Sections 2 and 3. The results are summarized in the top panel of Table 1, where we also calculated the Akaike information criterion (AIC) value for each model to evaluate model fit.

The AIC values for the ZIN and ZIN_c models were significantly lower than those of the N$-$mixture and N${}_{c}-$mixture models, respectively, indicating that including zero-inflated probabilities improves the fit of the model to the fisher data. Among the three zero-inflated models, the ZIN_c had the smallest AIC value. Furthermore, the ZIN_c model estimated the community parameter $c$ to be in the middle of $(0,1)$ and the associated Wald-type confidence interval (Web Table 14) provided strong evidence against the null hypotheses $c=0$ and $c=1$. The likelihood ratio statistics for the tests $c=0$ and $c=1$ were $40.21$ and $8.23$ respectively, and the $p-$values based on 10,000 bootstrap samples were $0$ and $0.006$, respectively, leading to the same conclusion for rejecting ZIB and ZIN models. To sum up, the ZIN_c model fits the fisher data better than the other models, however, the occupancy estimates for the three zero-inflated models are similar.

It is worth noting that the ZIN model estimated $\psi$ to be $0.175$, which appears larger than the other models, but the ZIN model has an unoccupied probability of $(1-\psi)+\psi\exp(-\mu)$. Therefore, to compare occupancy rates with other models, we should use $\psi\{1-\exp(-\mu)\}$ with an estimate of $0.175\{1-\exp(-1.79)\}=0.146$ (called the occupied probability of the ZIN model), which is similar to other estimates of $\psi$. In contrast, the resulting abundance estimates are quite different, with ZIN reporting a higher value and ZIB having the lowest value. The standard error estimates for the three models also follow the same order, with the parameter estimates for ZIN showing the greatest estimated variation. This phenomenon is similar to what was observed in Web Table 8 of Simulation Study B.

This example aims to understand how the community parameter $c$ may vary among different species in data from the same survey. We consider the data used in Royle (2006), originally from the North American Bird Breeding Survey (BBS) program, which includes records of occurrence data for five species of birds -- blue jay ($Cyanocitta~{}cristata$), catbird ($Dumetella~{}carolinensis$), common yellow-throat ($Geothlypis~{}trichas$), tree swallow ($Tachycineta~{}bicolor$), and song sparrow ($Melospiza~{}melodia$)-- from $50$ locations with $11$ repeated visits. The sample occupancy rates for the five species (in the above order) are $66\%$, $38\%$, $72\%$, $58\%$, and $52\%$; see Web Table 13 for the data.

Results are reported in Table 1, where the order of the species was actually sorted according to the estimated value of $c$ in the ZIN_c model. The first two species show small estimates of $c$, with the ZIN_c model even degenerating to the ZIB model for the blue jay data due to an estimate of $c=0$. All the considered models showed comparable AIC values for both species, indicating a similar level of model fit. However, the parameter estimates of $\mu,~{}r$, and $c$ for each model produced conflicting results. We believe this is a problem of model identifiability, similar to the one described Royle (2006) and Link (2003). Royle (2006) concluded that this problem is more pronounced in occupancy models when the probability of detection is low, which is the case for these two species.

The next two species, common yellow-throat and tree sparrow, show middle estimates of $c$, where both $95\%$ Wald type confidence interval does not include $0$ or $1$ (Web Table 14 ). The results for model comparisons here can be summarized as similar to Example 1, particularly for zero-inflated models. In terms of AIC, the ZIN_c model performs better than other models.

The last species, song sparrow, has a high estimate of $c$ of $0.90$, which is close to $1$. In fact, the upper limit of its $95\%$ confidence interval exceeds the upper bound of $1$. In this case, the standard errors of $\hat{\mu}$ and $\hat{\psi}$ for the ZIN_c model are relatively high, indicating that estimation uncertainty increases as $c$ approaches $1$. This phenomenon was also observed in Simulation Study B.

In Web Table 14, we also report bootstrap $p-$values for the likelihood ratio tests to the hypotheses of ZIB and ZIN models, which suggest significant evidence against the null hypotheses for the last three species of BBS data.