# Nc-mixture occupancy model

Huu-Dinh Huynh  Institute of Statistics, National Chung Hsing University,
Taichung, Taiwan  Wen-Han Hwang  Institute of Statistics, National Tsing Hua
University, Hsinchu, Taiwan
<EMAIL_ADDRESS>
###### Abstract

A class of occupancy models for detection/non-detection data is proposed to
relax the closure assumption of N$-$mixture models. We introduce a community
parameter $c$, ranging from $0$ to $1$, which characterizes a certain portion
of individuals being fixed across multiple visits. As a result, when $c$
equals $1$, the model reduces to the N$-$mixture model; this reduced model is
shown to overestimate abundance when the closure assumption is not fully
satisfied. Additionally, by including a zero-inflated component, the proposed
model can bridge the standard occupancy model ($c=0$) and the zero-inflated
N$-$mixture model ($c=1$). We then study the behavior of the estimators for
the two extreme models as $c$ varies from $0$ to $1$. An interesting finding
is that the zero-inflated N$-$mixture model can consistently estimate the
zero-inflated probability (occupancy) as $c$ approaches $0$, but the bias can
be positive, negative, or unbiased when $c>0$ depending on other parameters.
We also demonstrate these results through simulation studies and data
analysis.

Keywords: Community parameter; N$-$mixture models; Zero-inflated model

##  1 Introduction

Estimating the occupancy and abundance of species in a specific region,
despite imperfect detection, is a crucial problem in ecology conservation and
management. In reality, a species may exist at a survey site but go undetected
due to limitations in the survey method or timing. Zero counts in a site
sampling survey can be caused by the species not being present or by detection
errors. If this issue is not addressed, the occupancy rate will be
underestimated. Site occupancy models (MacKenzie et al., 2002), based on
multiple-visit detection/non-detection (occurrence, presence/absence) data,
can estimate the occurrence rate by accounting for detection errors. These
models, which can be temporal or spatial replication surveys, are cost-
effective as they do not require individual identification or marking. They
are widely used in species distribution modeling, with various extensions such
as multi-season open models, multi-species models, dynamic models, and
spatial-temporal models (MacKenzie et al., 2017; Hogg et al., 2021; MacKenzie
et al., 2009; Johnson et al., 2013) having been developed. This study focuses
on single-species, single-season occupancy models.

In the standard occupancy model (MacKenzie et al., 2002), detection
probability and species abundance are confounded and cannot be distinguished
from one another. To overcome this limitation, Royle and Nichols (2003); Royle
(2004) introduced N$-$mixture occupancy models that allow for the separation
of detection probability and species abundance. These models enable the
estimation of species abundance through multiple-visit occurrence or count
data. N$-$mixture models have received a lot of attention in the literature
(Haines, 2016a, b; Joseph et al., 2009; Gomez et al., 2018) as they can
estimate population size like capture-recapture models, but without the need
for capturing and marking individuals. However, the performance of these
models depends on the assumptions made (Dennis et al., 2015; Link et al.,
2018; Barker et al., 2018). Therefore, our goal is to improve these models by
relaxing the closure model assumption, which is often violated even in single-
season surveys (Kendall et al., 2013; Otto et al., 2013).

The closure assumption states that the number of individuals at a site remains
constant during multiple visits. Surveys are often conducted through temporal
or spatial replication, and sometimes with the use of multiple detectors or a
combination of these methods (MacKenzie et al., 2017; Kendall and White,
2009). Temporal replication involves surveying the same sites at different
times, while spatial replication involves selecting random sampling units
within a larger area at a single site. While this closure assumption is
reasonable for surveys conducted in the same locations over a short period, it
may not be appropriate for highly mobile species (Hayes and Monfils, 2015). We
also note that the study in Royle (2004) is an example of temporal
replication, despite the paper’s title referencing spatial replication to
indicate the distribution of sites.

The inference of the N$-$mixture occupancy model has been shown to produce
biased point estimates and incorrect interval estimates when the closure
assumption is violated. While a few studies have highlighted these effects
(Dénes et al., 2015; Duarte et al., 2018; Dail and Madsen, 2011; Ke et al.,
2022), most evidence is derived from simulation studies. To the best of our
knowledge, there is a lack of theoretical results that explain the behavior of
estimation bias under the N$-$mixture model. However, we have recently
provided theoretical evidence under the proposed N${}_{c}-$mixture model,
which is an extension of the N$-$mixture model. It is important to note that
Dail and Madsen (2011) also proposes a class of generalized N$-$mixture
models, which allows for the immigration and emigration of species populations
and estimates year-to-year immigration and emigration rates, providing
valuable insights for conservation managers. However, this model is a multi-
season open population model that goes in a different direction than our
extension.

The N${}_{c}-$mixture model is designed to create a framework that can unify
both temporal and spatial replicating surveys. For example, consider a triple-
visit survey conducted at a site where $N_{1}$, $N_{2}$, and $N_{3}$ represent
the number of observable individuals during each visit. These variables are
assumed to be identically distributed random variables. In the case of
temporal replication, where the surveys are conducted in a short period, the
three $N_{j}$ variables can be considered equal and meet the closure
assumption of the N$-$mixture model. In contrast, for spatial replication, the
$N_{j}$ variables are treated as independent. To account for both scenarios,
we decompose the $N_{j}$ variables into two components: $N_{j}=K+M_{j}$ for
$j=1,2,3$, where $K$ represents the number of common individuals during the
triple-visit survey, and $M_{j}$ represents the number of non-common
individuals. We assume that $K$ and $M_{j}$ are independent, with
$E(K)=cE(N_{j})$ for some $0\leq c\leq 1$. This parameter $c$, referred to as
the community parameter, indicates the proportion of individuals who are
residents and remain fixed during the triple-visit. It also allows us to
easily let $K$ degenerate to $0$ if $c=0$, or let $M_{j}$ degenerate to $0$ if
$c=1$. Figure 1 provides an illustration of this decomposition.

Figure 1:  An illustration of the N${}_{c}-$mixture model. In a triple-visit
survey of a site, some individuals are present in the site for each of the
three visits (denoted as $v_{1}$, $v_{2}$, and $v_{3}$). The figure
illustrates the situations corresponding to the community parameter $c=0$
(left), $c=0.5$ (middle), and $c=1$ (right). The colored circles represent
those individuals who are residents (or, equivalently, are fixed) in the site
during the three visits. For $c=0$, the number of individuals can differ from
visit to visit, mimicking the scenario of spatial replication. For $c=1$, the
number of individuals is constant, and the closure assumption of the
N$-$mixture model is satisfied. For $c=0.5$, around half of the individuals
were residents during the survey.

Under the N${}_{c}-$mixture model, we are able to demonstrate that if the
community parameter $c$ is incorrectly specified as $1$ (i.e., the N$-$mixture
model is used instead), the mean abundance will be overestimated.
Additionally, our results indicate that the bias increases as the community
parameter $c$ decreases, reaching infinity as $c$ approaches $0$.

We propose an extension to the N${}_{c}-$mixture model by incorporating a
zero-inflated component. This allows us to bridge the standard occupancy model
when $c=0$ (MacKenzie et al., 2002) and the zero-inflated N$-$mixture model
when $c=1$ (Haines, 2016a). We then investigate the behavior of estimators for
these two extreme models as the community parameter $c$ ranges from $0$ to
$1$. Our findings reveal that the standard occupancy model underestimates the
zero-inflated probability (occupancy), and the bias increases as the community
parameter $c$ increases. However, an interesting finding is that the zero-
inflated N$-$mixture model can estimate the occupancy consistently as $c$
approaches $0$, but the bias can be positive or negative depending on the
other parameters.

The paper is organized as follows: In Section 2, we develop the
N${}_{c}-$mixture model and present estimation methods. In Section 3, we
consider the zero-inflated N${}_{c}-$mixture model. Section 4 includes
simulation studies to evaluate estimator performance. Two real data examples
are presented in Section 5, and a discussion is in Section 6. Web Appendix A
includes proofs for all propositions.

##  2 N${}_{c}-$mixture model

###  2.1 Notation and estimation

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)

###  2.2 Occupancy rate

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.

###  2.3 Behavior of the N$-$mixture model estimators

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.

Under the N${}_{c}-$mixture model with a double-visit survey, as the number of
sites increases to infinity, the estimators $\tilde{\mu}_{1}$ and
$\tilde{r}_{1}$ converge to $\mu/c$ and $cr$ respectively with probability
one, for all $0<c\leq 1$.

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.

######  Proposition 2.

The method of moment estimators of the N$-$mixture model,
$\tilde{\mu}_{1\rm{M}}$ and $\tilde{r}_{1\rm{M}}$, are obtained by solving the
equations

$\bar{Y}=Tp~{}~{}\mbox{and~{}}~{}\bar{Y}^{2}=Tp+T(T-1)\\{2p-1+(1-p)^{2-r}\\},$

where $\bar{Y}$ is the sample average and
$\bar{Y}^{2}=\sum_{i=1}^{n}Y_{i}^{2}/n$. Furthermore, if the N${}_{c}-$mixture
model is true, as the number of sites increases to infinity, the estimators
$\tilde{\mu}_{1\rm{M}}$ and $\tilde{r}_{1\rm{M}}$ converge to $\mu/c$ and $cr$
respectively with probability one for all $0<c\leq 1$.

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$.

######  Proposition 3.

Under the N${}_{c}-$mixture model with $c$ approaching zero, as the number of
sites $n$ increases to infinity, the estimators $\tilde{\mu}_{1}$ and
$\tilde{r}_{1}$ converge such that $\tilde{\mu}_{1}\rightarrow\infty$,
$\tilde{r}_{1}\rightarrow 0$, and $\tilde{\mu}_{1}\tilde{r}_{1}\rightarrow\mu
r$ with probability one.

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.

##  3 Zero-inflated N${}_{c}-$mixture model

To account for the species occupancy, we extend the model (1) by incorporating
a zero-inflation component. Following MacKenzie et al. (2002), let $\psi$ be
the site occupancy probability, then the probability of a zero count
($Y_{i}=0$) is $(1-\psi)+\psi f(0;\bm{\theta})$. The likelihood function thus
has an additional parameter $\psi$ so that we may write

$L(\bm{\theta},\psi)=\prod_{i=1}^{n}\left\\{(1-\psi)\mathrm{\bf
I}(y_{i}=0)+\psi f(y_{i};\bm{\theta})\right\\}.$ (5)

We refer the model (5) as a zero-inflated N${}_{c}-$mixture (ZINc) model.

When $c=1$, the model (5) becomes a zero-inflated N$-$mixture (ZIN) model, as
previously described in Haines (2016a). In the case of $c=0$, it simplifies to
a zero-inflated binomial (ZIB) model with the detection probability
$p=1-\exp(-\mu r)$, which is also known as the first occupancy model
(MacKenzie et al., 2002). Thus, the ZINc model unifies both the ZIB and ZIN
models into a single framework, with the special cases of $c=0$ and $c=1$
corresponding to ZIB and ZIN, respectively.

###  3.1 Estimation

When the community parameter $0<c<1$, the ZINc 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
ZINc 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.

###  3.2 Behaviors of $\tilde{\psi}_{0}$ and $\tilde{\psi}_{1}$

Like the N${}_{c}-$mixture model, the ZINc 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}$).

######  Proposition 4.

Under the ZINc model with $c>0$, the ZIB occupancy estimator
$\tilde{\psi}_{0}$ shows an underestimation of $\psi$ with probability one as
$n$ increases to infinity. A linear approximation of this underestimation can
be represented as $\tilde{\psi}_{0}\approx\frac{p}{p+\Delta}\psi$ where

$\Delta=\frac{\\{1-(1-p)^{T}\\}\left\\{f(0;\bm{\theta})-(1-p)^{T}\right\\}}{\left[\frac{1}{p}\\{1-(1-p)^{T}\\}-T(1-p)^{T-1}\right]\\{1-f(0;\bm{\theta})\\}},$
(6)

and $\Delta$ increases as $c$ increases, when $c=0$ then $\Delta=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 ZINc model).

######  Proposition 5.

Under the N${}_{c}-$mixture model with $c$ approaching zero, as the number of
sites $n$ increases, the estimators $\tilde{\mu}_{1}$ and $\tilde{r}_{1}$
converge such that $\tilde{\mu}_{1}\rightarrow\infty$,
$\tilde{r}_{1}\rightarrow 0$, and $\tilde{\mu}_{1}\tilde{r}_{1}\rightarrow\mu
r$ with probability one. Additionally, $\tilde{\psi}_{1}$ is a consistent
estimator for $\psi$ when $c=0$.

The estimator $\tilde{\psi}_{1}$ is also a consistent estimator under the ZIN
(or ZINc 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.

###  3.3 Tests for ZIN and ZIB

The null hypothesis of the ZIN or ZIB model can be tested within the framework
of the ZINc 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 ZINc 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).

##  4 Simulation study

We conducted simulations to evaluate the performance of proposed models and
estimators. We considered two scenarios: the N${}_{c}-$mixture model and the
ZINc model. In the first scenario, we computed the maximum likelihood
estimators for both the N$-$mixture model and N${}_{c}-$mixture model. In the
second scenario, we calculated the maximum likelihood estimators for the ZIB,
ZIN, and ZINc models. All estimates were calculated by using the optim
function in the R software (R Core Team, 2022).

We specified the true parameter values as $\mu=1,2$, $r=0.25,0.5$, and
$c=0,0.05,\ldots,1$ for the simulations. For both scenarios, the number of
sites was set at $n=200,500,1000$, and the number of visits was set at
$T=5,7,10$. In the second scenario, we also set the occupancy probability to
$\psi=0.7$. We generated 1,000 data sets for each parameter setting and
calculated the estimates and associated standard error estimates for each data
set.

To account for outliers in some of the simulated data sets, we present the
median of the parameter estimates (Med), the median of the estimated standard
error (Med.se), and the median absolute deviation (MAD) scaled to align with
the normal distribution. Additionally, we report the coverage percentage (CP)
of the nominal $95\%$ Wald-type confidence intervals. We note that in some
instances, the numerical methods utilized for estimating the parameters for
each model did not converge, with a higher frequency of non-convergence
observed for the ZIN model and a lower frequency for the ZINc model. However,
in most cases, the percentage of failures was minor and not reported.

###  4.1 Simulation study A: N${}_{c}-$mixture model

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:  Median estimates of abundance $\mu$ for N$-$mixture and
N${}_{c}-$mixture models (Simulation study A) as a function of the community
parameter $c$, with the number of visits $T=7$ and sites $n=500$. The sub-
graphs correspond to four combinations of $\mu=1,2$ and individual detection
probabilities $r=0.25,0.5$. Data points with high values were removed for
clarity.

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.

###  4.2 Simulation study B: ZINc model

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.

Figure 3:  Median estimates of abundance $\mu$ for ZIB, ZINc, and ZIN models
(Simulation study B) as a function of the community parameter $c$, with the
number of visits $T=7$ and sites $n=500$. The sub-graphs correspond to four
combinations of $\mu=1,2$ and individual detection probabilities $r=0.25,0.5$.
Data points with high values were removed for clarity.

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 ZINc 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.

Figure 4: Median estimates of occupancy rate $\psi$ for ZIB, ZINc and ZIN
models (Simulation study B) plotted against the community parameter $c$, with
the number of visits $T=7$ and sites $n=500$. The sub-graphs correspond to
four combinations of $\mu=1,2$ and individual detection probabilities
$r=0.25,0.5$.

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 ZINc
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 ZINc 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 ZINc 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 ZINc
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.

##  5 Examples

###  5.1 Example 1. Fisher data

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 ZINc 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 ZINc had the smallest
AIC value. Furthermore, the ZINc 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 ZINc 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.

Model | $\mu$ | $r$ | $c$ | $\psi$ | AIC
---|---|---|---|---|---
| Fisher
ZINc | 0.93 (0.21) | 0.49 (0.10) | 0.51 (0.13) | 0.154 (0.020) | 626.9
ZIN | 1.79 (0.68) | 0.22 (0.05) | $1^{*}$ | 0.175 (0.030) | 633.2
ZIB | 0.51 (0.04) | $1^{*}$ | $0^{*}$ | 0.140 (0.016) | 663.1
N${}_{c}-$mix. | 0.13 (0.02) | 0.46 (0.03) | 0.93 (0.02) |  | 636.8
N$-$mix. | 0.16 (0.02) | 0.37 (0.02) | $1^{*}$ | 0.150 (0.017) | 650.6
| Blue Jay
ZIN | 48.96 (427.97) | 0 (0.04) | $1^{*}$ | 0.729 (0.091) | 170.1
ZIB | 0.22 (0.03) | $1^{*}$ | $0^{*}$ | 0.723 (0.077) | 168.1
N$-$mix. | 1.64 (0.56) | 0.09 (0.03) | $1^{*}$ | 0.806 (0.109) | 169.5
| Catbird
ZINc | 0.60 (0.32) | 0.42 (0.21) | 0.19 (0.14) | 0.421 (0.079) | 139.1
ZIN | 0.85 (1.31) | 0.16 (0.09) | $1^{*}$ | 0.740 (0.691) | 137.8
ZIB | 0.26 (0.04) | $1^{*}$ | $0^{*}$ | 0.403 (0.074) | 138.5
N${}_{c}-$mix. | 0.52 (0.17) | 0.19 (0.06) | 0.97 (0.12) |  | 137.9
N$-$mix. | 0.55 (0.14) | 0.18 (0.04) | $1^{*}$ | 0.421 (0.079) | 135.9
| Common Yellow-Throat
ZINc | 0.93 (0.23) | 0.47 (0.09) | 0.56 (0.10) | 0.775 (0.075) | 222.9
ZIN | 2.06 (0.81) | 0.20 (0.05) | $1^{*}$ | 0.851 (0.123) | 227.8
ZIB | 0.49 (0.04) | $1^{*}$ | $0^{*}$ | 0.723 (0.064) | 254.5
N${}_{c}-$mix. | 1.09 (0.38) | 0.31 (0.09) | 0.88 (0.15) |  | 226.7
N$-$mix. | 1.45 (0.25) | 0.23 (0.03) | $1^{*}$ | 0.765 (0.059) | 226.8
| Tree Swallow
ZINc | 0.75 (0.22) | 0.45 (0.08) | 0.61 (0.12) | 0.682 (0.106) | 196.4
ZIN | 1.81 (0.75) | 0.18 (0.05) | $1^{*}$ | 0.726 (0.135) | 204.0
ZIB | 0.40 (0.04) | $1^{*}$ | $0^{*}$ | 0.587 (0.071) | 228.0
N${}_{c}-$mix. | 0.48 (0.11) | 0.46 (0.07) | 0.72 (0.08) |  | 198.4
N$-$mix. | 1.02 (0.19) | 0.22 (0.03) | $1^{*}$ | 0.641 (0.069) | 203.9
| Song Sparrow
ZINc | 1.01 (0.59) | 0.37 (0.08) | 0.90 (0.09) | 0.678 (0.240) | 184.7
ZIN | 2.23 (0.98) | 0.21 (0.06) | $1^{*}$ | 0.573 (0.105) | 187.9
ZIB | 0.54 (0.05) | $1^{*}$ | $0^{*}$ | 0.501 (0.072) | 210.5
N${}_{c}-$mix. | 0.58 (0.13) | 0.42 (0.05) | 0.93 (0.04) |  | 183.4
N$-$mix. | 0.81 (0.16) | 0.31 (0.04) | $1^{*}$ | 0.555 (0.071) | 190.5
Table 1:  Parameter estimates (standard errors in parenthesis) and AIC values
for each model fitted to the fisher data and the five bird species data of the
BBS survey. Note that ∗ indicates that the parameter value is a fixed constant
without estimation. The occupancy rate $\psi$ of the N$-$mixture model is a
deriving parameter defined as $\psi=1-\exp(-\mu)$; there is no occupancy rate
for the N${}_{c}-$mixture model (Section 2.2). For the blue jay data, the
results for the ZINc and N${}_{c}-$mixture models are omitted because they are
reduced to ZIB and N$-$mixture models, respectively. The five bird species are
sorted according to the estimate of $c$ from the ZINc model.

###  5.2 Example 2: Breeding Bird Survey (BBS) data

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 ZINc model. The first
two species show small estimates of $c$, with the ZINc 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 ZINc 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 ZINc 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.

##  6 Discussion

The N$-$mixture model has been extended to allow for the number of individuals
at a sample site to vary with each visit. In this extension, the number of
individuals at a site per visit is treated as a latent variable that is
decomposed into two components, with only one of these components being
considered in the original N$-$mixture model. It should be noted that the
extended N${}_{c}-$mixture model is only applicable to closed populations
(single season) due to the assumption of constant $\mu$. While it is possible
to relax this assumption using regression models with the site and/or visit
covariates, it is unclear whether this approach is effective for some
multiseasonal data.

Unlike the N$-$mixture model, the occupancy of the N${}_{c}-$mixture model may
not be well defined by the abundance parameter. To address this, we propose a
zero-inflated N${}_{c}-$mixture model that uses a zero-inflated probability to
define occupancy explicitly. This new model unifies the commonly used standard
occupancy and zero-inflated N$-$mixture models into one framework. As a
result, we discover interesting properties of both models, such as the use of
standard occupancy estimators under the zero-inflated N$-$mixture model or
vice versa.

Our extension models offer greater flexibility in fitting data, but they also
increase the complexity of the model. These models may also experience issues
such as numerical instability, parameter identifiability, and model
identification problems, particularly when data is sparse, such as in cases of
low detection probability, low abundance, few replicates of visits, or a small
sample of sites. These issues have been observed in both simulation studies
and data analysis. Indeed, these problems were initially raised in N$-$mixture
models (Dennis et al., 2015; Barker et al., 2018; Link et al., 2018) as a
caution against their use. Our extension models require an additional
community parameter, which may exacerbate these problems. However, these
issues may be mitigated by incorporating additional capture-recapture data
(Barker et al., 2018) or including covariates in the models (Link et al.,
2018). Further research is necessary to improve the estimation in these
contexts.

There are several potential avenues for extending this work, each of which
would require further implementation and investigation.

  * •

Using the negative binomial distribution to model abundance instead of the
Poisson distribution. This distribution is often more appropriate for
biological and ecological data analysis, but it also includes an additional
aggregation parameter, increasing model complexity and estimation difficulty.

  * •

Examining count data instead of occurrence data. This would require a
likelihood function without an explicit form, making estimation
computationally intensive, particularly for large datasets.

  * •

Incorporating covariates into the models. Including covariates is important
for real-world applications. This can be done by using a log link for
abundance and a logit link for detection and occupancy probability.
Incorporating species covariates for the community parameter in a joint
species model may also be meaningful.

  * •

Extending to time-to-detection data. Time-to-detection occupancy models have
been developed recently, which also show some advantages over occurrence data
models (Priyadarshani et al., 2022). Strebel et al. (2021) propose an
N$-$mixture time-to-detection occupancy model that enables estimating the
abundance without marking individuals. An analogous N${}_{c}-$mixture
occupancy model for time-to-detection data is currently under development.

## Acknowledgements

This work was supported by the National Science and Technology Council of
Taiwan.

## References

  * Barker et al. (2018) Barker, R.J., Schofield, M.R., Link, W.A., and Sauer, J.R. (2018). On the reliability of $\text{N}-$mixture models for count data. Biometrics, 74, 369--377.
  * Dail and Madsen (2011) Dail, D. and Madsen, L. (2011). Models for estimating abundance from repeated counts of an open metapopulation. Biometrics, 67, 577--587.
  * Dénes et al. (2015) Dénes, F.V., Silveira, L.F., and Beissinger, S.R. (2015). Estimating abundance of unmarked animal populations: accounting for imperfect detection and other sources of zero inflation. Methods in Ecology and Evolution, 6, 543--556.
  * Dennis et al. (2015) Dennis, E.B., Morgan, B.J., and Ridout, M.S. (2015). Computational aspects of N$-$mixture models. Biometrics, 71, 237--246.
  * Duarte et al. (2018) Duarte, A., Adams, M.J., and Peterson, J.T. (2018). Fitting N$-$mixture models to count data with unmodeled heterogeneity: Bias, diagnostics, and alternative approaches. Ecological Modelling, 374, 51--59.
  * Gomez et al. (2018) Gomez, J.P., Robinson, S.K., Blackburn, J.K., and Ponciano, J.M. (2018). An efficient extension of N$-$mixture models for multi-species abundance estimation. Methods in Ecology and Evolution, 9, 340--353.
  * Haines (2016a) Haines, L.M. (2016a). A note on the Royle--Nichols model for repeated detection--nondetection data. Journal of Agricultural, Biological, and Environmental Statistics, 21, 588--598.
  * Haines (2016b) Haines, L.M. (2016b). Maximum likelihood estimation for N$-$mixture models. Biometrics, 72, 1235--1245.
  * Hayes and Monfils (2015) Hayes, D.B. and Monfils, M.J. (2015). Occupancy modeling of bird point counts: Implications of mobile animals. The Journal of Wildlife Management, 79, 1361--1368.
  * Hogg et al. (2021) Hogg, S.E., Wang, Y., and Stone, L. (2021). Effectiveness of joint species distribution models in the presence of imperfect detection. Methods in Ecology and Evolution, 12, 1458--1474.
  * Johnson et al. (2013) Johnson, D.S., Conn, P.B., Hooten, M.B., Ray, J.C., and Pond, B.A. (2013). Spatial occupancy models for large data sets. Ecology, 94, 801--808.
  * Joseph et al. (2009) Joseph, L.N., Elkin, C., Martin, T.G, and Possingham, H.P (2009). Modeling abundance using N$-$mixture models: the importance of considering ecological mechanisms. Ecological Applications, 19, 631--642.
  * Karavarsamis and Huggins (2020) Karavarsamis, N. and Huggins, R. (2020). Two-stage approaches to the analysis of occupancy data I: The homogeneous case. Communications in Statistics - Theory and Methods, 49, 4751--4761.
  * Ke et al. (2022) Ke, A., Sollmann, R., Frishkoff, L.O., and Karp, D.S. (2022). A hierarchical N$-$mixture model to estimate behavioral variation and a case study of neotropical birds. Ecological Applications, 32, e2632.
  * Kendall and White (2009) Kendall, W.L. and While, G.C. (2009). A cautionary note on substituting spatial subunits for repeated temporal sampling in studies of site occupancy. Journal of Applied Ecology, 46, 1182--1188.
  * Kendall et al. (2013) Kendall, W.L., Hines, J.E., Nichols, J.D., and Grant, E.H.C. (2013). Relaxing the closure assumption in occupancy models: staggered arrival and departure times. Ecology, 94, 610--617.
  * Link (2003) Link, W.A. (2003). Nonidentifiability of population size from capture-recapture data with heterogeneous detection probabilities. Biometrics, 59, 1123--1130
  * Link et al. (2018) Link, W.A., Schofield, M.R., Barker, R.J., and Sauer, J.R. (2018). On the robustness of N$-$mixture models. Ecology, 99, 1547--1551.
  * MacKenzie et al. (2002) MacKenzie, D.I., Nichols, J.D., Lachman, G.B., Droege, S., Royle, J.A., and Langtimm, C.A. (2002). Estimating site occupancy rates when detection probabilities are less than one. Ecology, 83, 2248--2255.
  * MacKenzie et al. (2009) MacKenzie, D.I., Nichols, J.D., Seamans, M.E., and Gutiérrez, R.J. (2009). Modeling species occurrence dynamics with multiple states and imperfect detection. Ecology, 90, 823--835.
  * MacKenzie et al. (2017) MacKenzie, D.I., Nichols, J.D., Royle, J.A., Pollock, K.H., Bailey, L.L., and Hines, J. E. (2017). Occupancy estimation and modeling: inferring patterns and dynamics of species occurrence, 2nd Edition. London: Elsevier/Academic Press.
  * Otto et al. (2013) Otto, C.R.V., Bailey, L.L., and Roloff, G.J. (2013). Improving species occupancy estimation when sampling violates the closure assumption. Ecography, 36, 1299--1309.
  * Priyadarshani et al. (2022) Priyadarshani, D., Altwegg, R., Lee, A.T.K., and Hwang, W.H. (2022). What can occupancy models gain from time-to-detection data? Ecology, 103, e3832.
  * R Core Team (2022) R Core Team (2022). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.
  * Royle and Nichols (2003) Royle, J.A. and Nichols, J.D. (2003). Estimating abundance from repeated presence-absence data or point counts. Ecology, 84, 777--790.
  * Royle (2004) Royle, J.A. (2004). N$-$mixture models for estimating population size from spatially replicated counts. Biometrics, 60, 108--115.
  * Royle (2006) Royle, J.A. (2006). Site occupancy models with heterogeneous detection probabilities. Biometrics, 62, 97--102.
  * Royle and Dorazio (2008) Royle, J.A. and Dorazio, R.M. (2008). Hierarchical modeling and inference in ecology: The analysis of data from populations, metapopulations and communities. Elsevier/Academic Press.
  * Self and Liang (1987) Self, S.G., and Liang, K.L. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association, 82, 605--610.
  * Strebel et al. (2021) Strebel, N., Fiss, C.J., Kellner, K.F., Larkin, J.L., Kéry, M., and Cohen, J. (2021). Estimating abundance based on time-to-detection data. Methods in Ecology and Evolution, 12, 909--920.
  * Zielinski et al. (2005) Zielinski, W.J., Truex, R.L., Schlexer, F.V., Campbell, L.A., and Carroll, C. (2005). Historical and contemporary distributions of carnivores in forests of the Sierra Nevada, California, USA. Journal of Biogeography, 32, 1385--1407.