# Finite mixtures in capture-recapture surveys for modelling residency
patterns in marine wildlife populations

Gianmarco Caruso
Dpt. of Statistical Sciences
Sapienza University of Rome
<EMAIL_ADDRESS>&Pierfrancesco Alaimo Di Loro
Dpt. GEPLI
LUMSA
<EMAIL_ADDRESS>&Marco Mingione
Dpt. of Political Sciences
Roma Tre University
<EMAIL_ADDRESS>&Luca Tardella
Dpt. of Statistical Sciences
Sapienza University of Rome
<EMAIL_ADDRESS>&Daniela Silvia Pace
Dpt. of Environmental Biology
Sapienza University of Rome
<EMAIL_ADDRESS>&Giovanna Jona Lasinio
Dpt. of Statistical Sciences
Sapienza University of Rome
<EMAIL_ADDRESS>
###### Abstract

This work aims to show how prior knowledge about the structure of a
heterogeneous animal population can be leveraged to improve on its abundance
estimation from capture-recapture survey data. We combine the Open Jolly-Seber
(JS) model with finite mixtures and propose a parsimonious specification
tailored to the residency patterns of the common bottlenose dolphin. We employ
a Bayesian framework for our inference, discussing the appropriate choice of
priors to mitigate label-switching and non-identifiability issues, commonly
associated with finite mixture models. We conduct a series of simulation
experiments to illustrate the competitive advantage of our proposal over less
specific alternatives. The proposed approach is applied to data collected on
the common bottlenose dolphin population inhabiting the Tiber River estuary
(Mediterranean Sea). Our results provide novel insights into this population’s
size and structure, shedding light on some of the ecological processes
governing its dynamics.

_K_ eywords Abundance estimation; Bayesian modelling; Capture-recapture;
Finite mixtures.

##  1 Introduction

Capture-recapture (CR) methods are statistical techniques widely employed to
estimate the size of an elusive population for which it is impossible to get a
complete enumeration. This task, applied initially to ecology for the study of
fish and wildlife populations (Otis et al., 1978; Wu and Holan, 2017; Matechou
and Argiento, 2022), is now common to many other application fields such as
epidemiology (Chao et al., 2001; Böhning et al., 2020; Maruotti et al., 2023)
and social sciences (Böhning and van der Heijden, 2009; Brittain and Böhning,
2009; Böhning et al., 2018; Silverman, 2020; Di Cecco et al., 2020; Farcomeni,
2022). The term capture is inherited from the traditional way wild animals
have been identified for decades - namely through capture, marking and release
- but it is not necessarily intended for its physical sense anymore.
Researchers increasingly adopt non-invasive methods for monitoring wild
populations to minimise the costs and impact on the population of interest.
Among those, photo-identification (Royle et al., 2009; Pace et al., 2021) and
DNA sampling (Bravington et al., 2016; Morin et al., 2016) are becoming
increasingly popular as they minimise behavioural responses that may bias the
final estimates (see Alunni Fegatelli and Tardella, 2013, and references
therein).

Original applications of such methods date back to the beginning of the 20th
century and were based on standard homogeneity assumptions on the population
structure and the identification process (Le Cren, 1965; Amstrup et al.,
2010). The literature is now rich in alternatives that can address a large
variety of deviations from such basic model assumptions and suit situations
where, for example, individuals exhibit heterogeneous behaviours (Pledger,
2000), sampling occurs in continuous time (Altieri et al., 2022), stop-over
sites are present (Matechou et al., 2013; Worthington et al., 2019; Wu et al.,
2021), temporary emigration is allowed (Zhou et al., 2019), and so on. For an
exhaustive review, see King and McCrea (2019) and references therein.

Our work concentrates on the abundance estimation of a common bottlenose
dolphin population inhabiting a delimited area over multiple years. Such
population is known to be open and calls for using the Open Jolly-Seber model,
a standard CR framework for open populations (Amstrup et al., 2010).
Furthermore, the established ecological literature affirms that it comprises
individuals with different residency patterns and results in a population
clustered into groups with different levels of site-fidelity. Hence, the
homogeneity assumption of the standard CR framework cannot hold, and the
behaviour of individuals belonging to different groups (i.e. entrance, capture
and survival) must be described by different parameters. Accurate estimation
of the clustering structure and parameters is of utmost interest to ecologists
who want to describe the population’s dynamics and inform conservation
policies. A widespread practice in such a heterogeneous setting is to consider
the inclusion of individual covariates that can help explain the differences
among the population’s members. However, informative covariates are often
unavailable and the heterogeneity is entirely latent, as it is in the case
study under consideration. Finite Mixture Models (FMM) represent the natural
solution to this impasse. Each individual is assigned to a different mixture
component with its own set of common and distinct parameter values. FMM
approaches to CR models have been successfully employed in a likelihood-based
framework both in closed population (Pledger, 2000; Dorazio and Royle, 2003;
Pledger, 2005) and open population (Pledger et al., 2010; Guéry et al., 2017)
settings. From the Bayesian perspective, attempts have been made to model
heterogeneity in detection and behavioural effects in closed populations by
Ghosh and Norris (2005). More recently, Turek et al. (2021) proposed a non-
parametric finite mixture model with an unknown number of components and
different capture probabilities.
We build on the Pledger et al. (2003, 2010)’s FMM extension to the Open JS
model to account for latent heterogeneity. We embed this finite mixture
approach in the Royle and Dorazio (2008, 2012)’s parameter-expanded data-
augmentation (PX-DA) framework, which turns out to be particularly convenient
to fit Bayesian CR models via standard Markov Chain Monte Carlo (MCMC)
algorithms. We discuss its implementation challenges and introduce suitable
prior specifications that mitigate the label-switching and non-identifiability
issues of FMM. The model is tested through an extended simulation study and
applied to photo-identification CR survey data of the common bottlenose
dolphins (_tursiops truncatus_) population inhabiting the area of the Tiber
River estuary in the Mediterranean Sea. In particular, we show how the prior
scientific knowledge on the population of interest can be leveraged to specify
a parsimonious FMM tailored to its supposed structure and check its validity.
When this is the case, we show it can sensibly improve the performances of
more comprehensive specifications.

The remainder of the paper is organised as follows: Section 2 describes the
motivation and intuitions behind our work, with a brief description of the
data that will be considered later on; Section 3 illustrates a Bayesian
hierarchical formulation of the Pledger et al. (2003, 2010)’s Jolly-Seber (JS)
class of mixture models within the Royle and Dorazio (2008, 2012)’s data-
augmentation framework and introduces our modelling proposal as a parsimonious
alternative; Section 4 reports a simulation experiment to investigate
different aspects of the model and its estimation; Section 5 provides the
model choice and results on the set of data originally introduced as the
motivating example.

##  2 Motivating example

This work is motivated by the need to get an abundance estimate of the
population of the common bottlenose dolphins inhabiting the area of the Tiber
River estuary in the Mediterranean Sea. Boat-based daily surveys have been
conducted between 2018 and 2020 under favourable weather conditions to collect
photographic and acoustic data of the specimens encountered in the study area
during the search (Papale et al., 2021; Pace et al., 2022c, a, b). The photo-
identification technique was used to identify unique individuals over multiple
sampling occasions and to build single capture histories. For this analysis,
we focus on the so-called well-marked individuals since the probability of
their misidentification can be assumed to be negligible111Photo-identification
is based on matching the same marks across different pictures; the level of
the mark (poor, fair, well) is an index of how identifiable a dolphin is.. As
a consequence, the final estimates are related to the subset of well-marked
individuals only, representing a portion of the population visiting the study
area. Figure 1(a) shows the cumulative number of identified individuals across
the different sampling occasions, where the size of each point is proportional
to the number of newly identified individuals. Its trend is known as the
discovery rate, which is maximum during the first year (notice that 50% of new
identifications occurred in 2018, with a maximum of 25 newly identified
individuals registered in August) and slowly decreases over time. Further
details about the study area, the data collection process and the analysis are
available in Pace et al. (2021).
Most of the recent literature about common bottlenose dolphins converges
towards the identification of three groups characterised by different levels
of site-fidelity to a specific area: from the most to the least frequently
present (Dinis et al., 2016; Hunt et al., 2017; Haughey et al., 2020; La Manna
et al., 2022). This feature is of utter interest to biologists interested in
disentangling the permanent or semi-permanent population from the transient
one. Generally speaking, one group is composed of individuals that (almost)
never leave the study area; these are usually referred to as resident
individuals, observable on many occasions and for long periods of time, and
expected to have the largest number of captures (see Figure 1(b)). The other
group includes individuals who are not continuously present in the study area
but regularly visit it; these are called part-time individuals, observable
throughout a wide time window but usually encountered at occasions far apart
in time. The last group comprises individuals that enter the study area only
once in their lifetime for a short time window and whose captures are rare;
these are transient individuals that are observable only on occasions
occurring on close dates.

(a)

(b)

Figure 1: (a) Cumulative number of individual identifications with size
proportional to the number of newly identified individuals; (b) Total number
of captures by individual.

This inherent heterogeneity cannot be neglected without hindering the
reliability of the final abundance estimates (Gimenez et al., 2018). To this
end, Pace et al. (2021) applies a hierarchical clustering and includes the
group labels in the POPAN-JS framework (Schwarz and Arnason, 1996) to model
heterogeneity in the entrance, capture and survival probabilities. They find
signs of the expected structure and estimate significant differences in most
components of the JS model. While practical for understanding the underlying
structure of the population of interest, this two-step approach has some
relevant issues. First of all, there is neither quantification nor propagation
of the uncertainty of the classification step onto the modelling step; this
can bias the final estimates and yield over-confident conclusions. Second, the
same information set is used twice in two different statistical procedures,
where the latter is performed after conditioning on the first; this can lead
to a confirmation bias in favour of the original hypothesis.

This paper proposes to unify the two steps into a joint statistical procedure.
We embed FMM into the Open-JS framework in a Bayesian hierarchical setting,
allowing for the estimation of cluster labels and population size altogether
while properly propagating the uncertainty at all levels (Clark and Gelfand,
2006).

##  3 State of the art and proposed model

Let $D$ be the number of distinct individuals observed at least once during
$T$ distinct sampling occasions. Data are collected in a ($D\times T$) matrix,
say $\mathbf{Y}$, with the generic element recording whether individual
$i=1,\dots,D$ individual has been observed ($y_{it}=1$) or not ($y_{it}=0$) at
occasion $t=1,\dots,T$. Hence, the rows of $\mathbf{Y}$ contain the capture
histories of all the encountered individuals. The open-population JS model
(Schwarz and Arnason, 1996) envisions individuals entering (i.e., via birth or
immigration) and exiting (i.e., via death or emigration) the population during
the sampling occasions. The emigration is assumed to be permanent for
identifiability purposes, i.e. once an individual has left the population, it
cannot return to it. Furthermore, the JS models assume that all captures are
independent across individuals and over time. The latter is achieved by
considering the existence of a super-population of unknown size, say
$N_{super}$, accounting for all the individuals potentially available
(encountered or not) in the study area between the first and the last sampling
period. The super-population size is the main parameter of interest and
determines the dimension of the model parameter space. Its typical estimation
in the Bayesian framework involves jumping between spaces of different
dimensions, hence requiring the implementation of Reversible Jump MCMC
algorithms (Brooks et al., 2000). Here, we adopt the data-augmentation
approach proposed by Royle and Dorazio (2008) to keep the dimension of the
parameter space fixed throughout the iterations, thus bypassing the need for
the Reversible Jump MCMC (Arnold et al., 2010).

###  3.1 Parameter Expanded Data-Augmentation (PX-DA) formalisation

The Parameter Expanded Data-Augmentation (PX-DA) approach by Royle and Dorazio
(2008) sets the dimension of the parameter space equal to a fixed value $M\gg
N_{super}\geq D$. In this way, solving the varying dimension issue by
converting it into a more manageable missing data problem in the multi-state
processes context (see Section A of the Appendix). The rows of the observed
data matrix $\mathbf{Y}$ are augmented to $M$, hence defining
$\mathbf{Y}_{\text{aug}}=\\{\mathbf{Y},\boldsymbol{0}_{M-D}\\}$, where
$\boldsymbol{0}_{M-D}$ is an $(M-D)\times T$ matrix of zeroes. $M$ must be set
such that $N_{super}\in\left\\{D,D+1,\dots,M\right\\}$ and, consequently,
$N_{super}-D$ among the $M-D$ rows of zeroes correspond to individuals who
belong to the super-population but have never been encountered. The remaining
$M-N_{super}$ correspond to pseudo-individuals that have never been part of
the population during the observation window and hence do not belong to the
super-population. The data generation process assumes that an individual can
be recruited into the population at the beginning of each sampling period only
if it has never been recruited on previous occasions. On the other hand,
individuals who have already been recruited can leave the population between
two subsequent sampling periods. All the $D$ observed individuals will
eventually be recruited in the population between the first and the last
capture occasion, as they have been captured at least once. This dynamic is
controlled through two time-varying latent binary variables: the first one,
$r_{it}$, is equal to $1$ if and only if individual $i$ is recruitable at time
$t$ ($0$, otherwise); the second one, $z_{it}$, is equal to $1$ if and only if
individual $i$ belongs to the population at time $t$ ($0$, otherwise).

All $M$ individuals are recruitable at the first occasion (i.e. $r_{i1}=1$,
$\,i=1,\dots,M$), while they become permanently non-recruitable once they have
entered the population. Let $\rho_{t},\,t=1,\dots,T$ be the recruitment
probabilities, i.e. the probability that an available (not yet entered)
individual in the augmented dataset is recruited in the population at time
$t$. Let $\phi_{t},\,t=2,\dots,T$, be the apparent survival probability222The
term apparent is included as permanent emigration and mortality are
indistinguishable and, thus, treated as the same phenomenon in JS-type
models., i.e. the probability that a recruited individual is in the population
at the following sampling occasion. The following rules govern the latent
label process: $z_{i1}\sim Bern(\rho_{1}),\,i=1,\dots,M$,
$r_{it}=\text{min}\\{r_{i,t-1},1-z_{i,t-1}\\}$ and
$z_{it}\,|\,z_{i,t-1},r_{it}\sim Bern(\phi_{t}\cdot z_{i,t-1}+\rho_{t}\cdot
r_{it})$, $t=2,\dots,T$. Thus, the distribution of the generic element of the
augmented data matrix can be expressed conditionally on $z_{it}$ as:

$y_{it}|z_{it}\sim Bern(p_{t}\cdot z_{it})\,,$ (1)

where $p_{t}$ is the capture probability at occasion $t$. Since $y_{it}=0$
almost surely when $z_{it}=0$, (1) is a zero-inflated binomial model. Notice
that, given $M$, the marginal likelihood of the capture histories can be
expressed as:

$\mathcal{L}_{M}(\boldsymbol{y}\,|\,\boldsymbol{\theta})=\int_{\mathcal{Z}}\prod_{i=1}^{M}\prod_{t=1}^{T}p(y_{it}\,|\,z_{it})\cdot
p(z_{it}\,|\,z_{i,t-1},r_{it})\,p(dz_{it}),$ (2)

where $\mathcal{Z}\subset\\{0,1\\}^{M\times T}$ contains all possible $z_{it}$
configurations and $\boldsymbol{\theta}$ is the full set of model parameters
$\boldsymbol{\theta}=\left\\{p_{t},\phi_{t},\rho_{t}\right\\}_{t=1}^{T}$. The
integral in the likelihood expression has no general closed-form solution and
numerically marginalising out the latent components can be computationally
intensive. This problem is common to the Hidden Markov Models (HMMs)
literature, in which some efficient techniques based on the forward algorithm
have been proposed to efficiently solve the integration problem Worthington et
al. (2019). Alternatively, MCMC algorithms provide a viable solution to work,
iteration after iteration, with the conditional likelihood:

$\mathcal{L}^{c}_{M}(\boldsymbol{y}\,|\,\boldsymbol{\theta},\boldsymbol{z})=\prod_{i=1}^{M}\prod_{t=1}^{T}p(y_{it}\,|\,z_{it}),$

and approximate the posterior distribution of all the quantities of interests
in a Bayesian setting that provides full uncertainty quantification.
Furthermore, standard MCMC methods are now relatively easy to implement by
practitioners due to the availability of software like JAGS (Plummer, 2003) or
NIMBLE (de Valpine et al., 2017).

Marginally, the hierarchical model specification resulting from Equation (1)
implies that $N_{super}\sim Binom(M,\,\psi)$, where $\psi$ is the overall
inclusion probability (throughout all occasions) of any individual in the
super-population. Royle and Dorazio (2008) show that $\psi$ is linked to the
recruitment probabilities through the following equation:

$\psi=1-\prod_{t=1}^{T}\,(1-\rho_{t}).$ (3)

In particular, (3) implies that
$\mathbb{E}[N_{super}\,|\,M,\rho_{1},\dots,\rho_{T}]=M\,\Bigl{[}1-\prod_{t=1}^{T}\,(1-\rho_{t})\Bigr{]}=M\,\psi$.
Hence, choosing the prior distribution for $\rho_{1},\dots,\rho_{T}$ is
crucial to determining the prior on $N_{super}$. Dorazio (2020) demonstrates
that the prior:

$\rho_{t}\sim Beta\biggl{(}\frac{1}{T},\,2-\frac{t}{T}\biggr{)},\,\quad
t=1,\dots,T,$ (4)

induces an objective prior on $N_{super}$. In terms of practical inference,
the estimated population size at each time $t$ and the overall super-
population size can be derived through the latent variables $z$’s, namely
$N_{t}=\sum_{i=1}^{M}\,z_{it}$ and
$N_{{super}}=\sum_{i=1}^{M}\,\mathbb{1}_{\\{\sum_{t=1}^{T}\,z_{it}>0\\}}$.

###  3.2 Jolly-Seber (JS) finite mixture modelling for open populations

Most populations are composed of individuals with heterogeneous behaviours. In
some cases, the heterogeneity can be assumed to be well-described by a finite
number (say $G$) of different patterns. These depend on the individual’s
latent traits, and their consideration requires complex modelling tools, such
as FMM. The underlying assumption of FMM is that each unit can belong to only
one group $g=1,\dots,G$, with unknown prior probabilities $w_{g}$
($\sum_{g}w_{g}=1$). The individuals belonging to different groups in open
population CR studies may have different capture, recruitment, or survival
parameters. In the most general specification, the relative order among the
parameters of different groups can change at each time $t$ (e.g. group one
could have the highest detection rate at the first sampling occasion but the
lowest at the second one). This JS-type mixture model is known as the
Interactive Heterogeneous Model (IHM) and its very rich specification depends
on too many parameters for successful model fitting (see Pledger et al., 2010,
for further discussion). Pledger et al. (2003, 2010) explore simpler
specifications that could adequately represent the population structure and
introduce a convenient notation to navigate through all possible sub-models.
Let $t$ and $h$ be the time and group heterogeneity effects, respectively.
Different expressions correspond to different modelling structures: constant
in time and homogeneous across groups $(\cdot)$; time-varying but homogeneous
across groups $(t)$; constant in time but heterogeneous across groups $(h)$;
time-varying and heterogeneous across groups with separable interaction
$(t+h)$; time-varying and heterogeneous across group with non-separable
interaction $(t\times h)$. For example, the IHM corresponds to
$\left\\{\left[\rho_{t\times h},\,\phi_{t\times h},\,p_{t\times
h}\right]_{G}\right\\}$, where the underlying population is supposed to be
composed by $G$ classes. If we want to specify a model whose heterogeneous
group effect lies in the capture probabilities only, we write
$\\{\left[{\rho}_{t},\,{\phi}_{t},\,{p}_{t\times h_{G}}\right]\\}$. The
subscript is moved to highlight that the mixture of $G$ components is related
only to detection. We will take advantage of this notation in the sequel of
the paper.

###  3.3 Modelling class heterogeneity using finite mixtures within the PX-DA
approach

We embed the PX-DA formalisation of the open JS model into FMM by adding one
layer of hierarchy in the original hierarchical specification.

Now, let $c_{i}\in\\{1,\dots,G\\}$ be the latent membership label of each
individual $i=1,\dots,M$ in $\mathbf{Y_{aug}}$. The full hierarchical
specification is as follows:

$\displaystyle y_{it}\,|\,z_{it},c_{i}=g\sim Bern(p_{gt}\cdot z_{it}),$ (5)
$\displaystyle z_{it}\,|\,z_{i,t-1},r_{it},c_{i}=g\sim Bern(\phi_{gt}\cdot
z_{i,t-1}+\rho_{gt}\cdot r_{it}),\quad
r_{it}=\text{min}\\{r_{i,t-1},1-z_{i,t-1}\\},$ $\displaystyle
p_{gt}\sim\pi_{p_{g}}(\cdot),\quad\phi_{gt}\sim\pi_{\phi_{g}}(\cdot),\quad\rho_{gt}\sim\pi_{\rho_{g}}(\cdot),\quad
c_{i}\sim Multinom\left(1,(w_{1},\dots,w_{G})\right),$
$\displaystyle(w_{1},\dots,w_{G})\sim\pi_{\mathbf{w}}(\cdot)\,$

where $\pi_{.}(\cdot)$ refers to a generic prior distribution for the
parameter. The hierarchical formulation of (5) further complicates the
likelihood expression of Equation (2) to:

$\mathcal{L}_{M}(\boldsymbol{y}\,|\,\widetilde{\boldsymbol{\theta}})=\sum_{g=1}^{G}w_{g}\cdot\int_{\mathcal{Z}}\prod_{i=1}^{M}\prod_{t=1}^{T}p_{g}(y_{it}\,|\,z_{it})\cdot
p_{g}(z_{it}\,|\,z_{i,t-1},r_{it})dz_{it},,$

where
$\widetilde{\boldsymbol{\theta}}=\left\\{\textbf{w},\left\\{\boldsymbol{\theta}_{g}\right\\}_{g=1}^{G}\right\\}$
includes the prior weights of the cluster components and the component-
specific sets of parameters. Its evaluation needs further marginalisation with
respect to the latent group labels, and, once again, MCMC methods are a viable
solution to achieve Bayesian estimation of all quantities of interest.

This modelling framework includes many possible specifications, according to
what varies with time and across groups. We further generalise the model
specification by adapting the survival mechanism to describe not equally
spaced capture occasions. Indeed, the assumption of constant survival across
identical time scales does not transfer to the not equally spaced scenario.
When this is the case - i.e. $l_{t}=(\tau_{t}-\tau_{t-1}),t=2,\dots,T$ are the
time differences between subsequent occasions - the survival probabilities
should be appropriately compounded. Once the time scale is set (e.g., days,
weeks, months, years, etc.), we have $\phi_{gt}=\phi_{g}^{l_{t}}$, where
$\phi_{g}$ represents the survival probability across a single time-unit on
the chosen scale. In addition, following Pledger et al. (2003), a convenient
and parsimonious way to express the time-varying capture probabilities is
through the logit link, i.e.
$\text{logit}(p_{gt})=\mu_{g}+\tau_{t}\,,\,t=1,\dots,T\,,$ where $\mu_{g}$
determines the overall average capture probability of each group and
$\tau_{t}$ is an occasion-specific differential effect.

The last step of the Bayesian model specification involves the choice of prior
distributions for all parameters and their hyper-parameters. The natural prior
for the mixture weights in finite mixture models is the Dirichlet distribution
$Dir_{G}(\alpha_{1},\dots,\alpha_{G})$. It corresponds to a uniform
distribution over the $G$-dimensional simplex when $\alpha_{g}=1,\;\forall g$,
x is usually adopted as weakly informative prior. For $\rho_{gt},t=1,\dots,T$,
we follow Dorazio (2020) and use the prior in (4). General-purpose and weakly
informative priors can be ascribed to $\mu_{g}$. At the same time, we suggest
considering a $N(0,\,\sigma^{2})$ on each parameter $\tau_{t}$, with
$\sigma^{2}$ small, to induce only small time variations (in the logit-scale)
on the capture probabilities. We further impose that $\sum_{t}\tau_{t}=0$ to
favour the otherwise weak identifiability of the $\mu_{g}$ and the $\tau$’s
(Pledger et al., 2003). Nevertheless, such choices do not defend the model
from the label-switching problem of the class-specific parameters typical of
the FMM framework because of the likelihood invariance under permutations of
the components’ labels.

#### Choosing priors on class-specific parameters

The choice of the prior distribution on component-specific parameters of FMMs
can be seen as a challenge and an opportunity. This class of models suffers
from several sources of non-identifiability. In our setting, the most affected
are recruitment and survival probabilities. However, the final inference on
the population size is usually robust with respect to their proper
identification (see Mena and Walker, 2015, and references therein). The
recruitment probabilities $\rho_{gt}$ are model devices that allow the
entrance of new individuals in the population, accounting for the openness of
the population (Royle and Dorazio, 2008). They can be treated as nuisance
parameters without a solid biological meaning and their identifiability is not
of great concern. On the other hand, the capture and survival probabilities
are of biological interest. They can be seen as indicators of how long and how
often units of different types visits the sampling area.

A weakly informative or regularising choice of their prior distributions can
favour identifying these components. In turn, this can ease the identification
of the membership labels already affected by the well-known label-switching
problem. Among the different solutions that have been proposed in the
literature, one relies on imposing an ordering constraint among the component-
specific parameters, i.e. $u_{1}<\dots<u_{G}$ (Richardson and Green, 1997;
Diebolt and Robert, 1994; Chung et al., 2004). This task is trivial in the
context of parameters belonging to the whole real domain, where truncating or
shifting Gaussian distributions allows full control of the position and scale
of the priors. The latter is a reasonable choice for the
$\mu_{g},\,g=1,\dots,G$ parameters. However, when
$u_{g}\in\left(0,1\right),\;\forall\,g$ as in the case of the survival
probabilities, the conditional prior specification must account for their
bounded domains. The standard solution would consider conditionally specified
Uniform distributions. However, we exploit the alternative proposal of Alaimo
Di Loro et al. (2022) that allows for more flexibility. It is based on the
Beta distribution and its generalisations (truncated or restricted), and it is
both effective in imposing the constraint and controlling for the shape and
first moments of the induced marginal priors. In other words, one may set
$u_{1}\sim Beta(\alpha_{1},\,\beta_{1})$ and

$\begin{aligned} &u_{g}|u_{g-1}\sim
tBeta(\alpha_{g},\,\beta_{g};u_{g-1},1),\,\text{ or }&u_{g}|u_{g-1}\sim
rBeta(\alpha_{g},\,\beta_{g};u_{g-1},1)\end{aligned}\qquad g=2,\dots,G,$ (6)

where tBeta and rBeta represent the Truncated Beta and Restricted Beta,
respectively, and it is possible to derive closed-form expressions for the
marginal prior distribution of $u_{g},g\geq 2$, or at least its first and
second moment. Further computational details are given in Section B of the
Appendix, along with formal proofs of the original results.

###  3.4 The RPT model for the common bottlenose dolphin population

The population structure illustrated in Section 2 can be translated to reflect
variations in the parameters of the JS-type PX-DA framework. We have that the
resident (R) individuals, showing high site-fidelity, should be characterised
by high survival and capture probabilities; the part-time resident (P)
individuals, with average site-fidelity, should be characterised by high
survival probabilities but eventually undetectable on some occasions; the
transient individuals (T) only shortly visit the study area and hence should
have very low survival. We want to build a suitably flexible encompassing
model such that it is possible to establish the amount of evidence in favour
of this assumption. In other words, we shall not enforce this exact structure
in our model, but we can make it so that it is recognisable if present. For
instance, we could allow the three groups to have group-specific survival and
capture probabilities as in the full specification of (5). That might be a too
flexible model with a huge number of parameters, therefore yielding highly
uncertain estimates because of weak identifiability. Alternatively, we can
propose a more parsimonious specification tailored to the supposed behavioural
differences among individuals from different groups. In this way, some
components can be assumed to be common to different groups, reducing the
model’s complexity and favouring its identifiability.

First of all, we model the population recruitment dynamic assuming that each
group is characterised by its own set of time-varying recruitability
parameters, i.e. $\rho_{R,\,t},\,\rho_{P,\,t},\,\rho_{T,\,t}$. This setting
naturally induces a clustering in the relationship between the recruitment
process and the super-population size $N_{super}$, yielding group-specific
inflation parameters that express the total population as the sum of three
sub-populations. This modifies the analytical expression of the expected
super-population size as follows:

$\mathbb{E}\left[N_{super}|\psi_{R},\psi_{P},\psi_{T},w_{R},w_{P},w_{T}\right]=M\cdot\sum_{g=R,P,T}\,w_{g}\psi_{g}\,,$

where $w_{g}$ and $\psi_{g}$ ($g=R,P,T$) are, respectively, the $g$-th
component-specific mixture weight and inflation parameter.

Second, we separate the three sub-populations into two distinct classes,
characterised by different tendencies to stay in the study area. We can
distinguish between short-term survivors, i.e. transient individuals who visit
the area for narrow time windows, and long-term survivors, i.e. non-transient
individuals (resident or part-time) who visit the area for wide time windows.
To model such behavioural heterogeneity, we assume that one group (T,
transient) has a smaller survival probability ($\phi_{T}$) than the other two
groups ($\phi_{NT}$): $\phi_{T}<\phi_{NT}$. Notice that these two parameters
represent the survival probabilities across a single time unit on the chosen
scale. When the capture occasions are not equally spaced, the survival
probabilities should be appropriately compounded along different lengths (cfr.
Section 3.3). In practice, given the two base survival probabilities
$\phi_{T}$ and $\phi_{NT}$, the actual survival probabilities at time $t$ are:

$\phi_{T,\,t}=\phi_{T}^{l_{t}},\quad\phi_{NT,\,t}=\phi_{NT}^{l_{t}},$

where $l_{t}=|t-t^{\prime}|,\,t=1,\dots,T$ are the time lags between the $T$
subsequent occasions. Finally, the proposed model retains the time-varying
structure of the detectability, as already introduced in Section 3.3, but
discriminates between part-time (P) and non-part-time (NP) individuals
introducing a partial undetectability component in the latter group. Loosely
speaking, part-time individuals are allowed to be undetectable while alive
(because temporarily not present in the study area) on some occasions chosen
at random with probability $\delta\in(0,1)$. This corresponds to modelling the
capture probability on occasion $t$ as:

$\begin{aligned}
&p_{NP,\,t}=\text{logit}^{-1}\left(\mu+\tau_{t}\right)\,\text{ and
}&p_{P,\,t}=(1-\delta)\,\cdot\,p_{NP,\,t},\end{aligned}\quad\,t=1,\dots,T,$

for the NP and P individuals, respectively; of course,
$p_{P,\,t}<p_{NP,\,t},\,\forall t$. Notice that the larger $\delta$ is, the
more the part-time group is separated from the resident group. This parameter
plays a similar role to the completely random emigration parameter of Kendall
et al. (1997) and it is needed to control the temporary emigration pattern of
part-time individuals discussed in Section 2. Appendix C shows that such
parametrisation is indeed equivalent to temporary emigration under the
simplifying assumption of emigration occurring at random. We name this model
RPT as it encompasses three specific types of behaviour. However, we would
like to point out how the model does not enforce this interpretation. For
instance, both survivals could be estimated to be high, or the undetectability
parameter could be estimated as approximately equal to $1$, etc.

##  4 Simulation experiments

We conduct a simulation experiment to assess the performances of the RPT model
whenever it is well-specified, i.e. the data are generated according to the
structure described in Section 3.4. We generate multiple sets of artificial
data under alternative scenarios from the RPT with fixed parameters and then
estimate a pool of models on them. Our main objective is twofold: i) to
evaluate the ability to recover the true values of the parameters, with a
particular focus on $N_{super}$; ii) to verify whether the RPT is chosen as
the best among other alternatives based on some model selection criterion.

We consider four scenarios characterised by an increasing number of sampling
occasions, i.e. $T\in\left\\{10,20,30,40\right\\}$, to verify the model
performances for different time horizons. We suppose that in the first
scenario (i.e. $T=10$), all the captures are recorded within a relatively
short period (e.g. within a year). Longer time horizons are included in the
other scenarios, where a larger time gap (year gap) is assumed to occur every
$10$ occasions. Further details about the time lags are available in Section D
of the Appendix. The month (and portion of months) is taken as the basic time
unit to avoid the possible numerical instability related to the large values
of the lags in terms of days. Note that this affects the interpretation of the
survival probability parameter as it must be interpreted as the probability of
surviving one month. We adopt the following parameters’ values in all the
scenarios. The survival probabilities are set to $\phi_{T}=0.01$ and
$\phi_{NT}=0.997$. These two values may seem quite extreme at first glance.
However, they guarantee that the short-time survivors (transient individuals)
almost surely stay in the population for less than a year and that the long-
term survivors (non-transient individuals) stay in the population for more
than three months with a very high probability ($>0.99$). Furthermore, notice
that a monthly survival probability equal to $0.01$ corresponds to a survival
probability equal to $0.86$ on a daily scale and equal to $0.34$ on a weekly
scale. On the other hand, a monthly survival probability equal to $0.997$
corresponds to a probability of $0.87$ on a four-years scale. Therefore, the
monthly scale appears as a good compromise to avoid a value that is too low
for $\phi_{T}$ and a value that is too high for $\phi_{NT}$. The capture
probabilities are obtained by setting $\mu=0$ and $\delta=0.7$, and generating
$\tau_{t}\sim N(0,0.25)$ in each scenario. The recruitment parameter for
transient individuals is fixed to $\rho_{T,\,t}=0.02$ for all capture
occasions (assuming a negligible temporal variation for this group). At the
same time, we consider the following structure for resident and part-time
individuals:

$\rho_{R,\,t}=\begin{cases}0.4\,,\quad t=1\\\ 0.02\,,\quad
t=10k+1\,,\,k=1,2,3\\\
0.0025\,,\quad\text{otherwise}\end{cases}\qquad\text{and}\qquad\rho_{P,\,t}=\begin{cases}0.4\,,\quad
t=1\\\ 0.04\,,\quad t=10k+1\,,\,k=1,2,3\\\
0.005\,,\quad\text{otherwise}\end{cases}.$

The mixture weights for the three groups are set to $w_{R}=0.2$, $w_{P}=0.45$
and $w_{T}=0.35$. We envision an augmented super-population of $M^{*}=500$,
that yields an expected super-population size
$\mathbb{E}\left[N_{super}\right]\in\\{170,209,243,271\\}$ for
$T=10,20,30,40$, respectively. Notice that $N_{super}$ increases with $T$ as
more individuals can visit the study area during a longer time horizon.

We simulate independent encounter histories for $K=50$ pseudo-populations.
Along with the RPT model, we consider 10 different alternatives in the class
of JS-type models described in Section 3.2 and 3.3, all having time-varying
recruitment parameters. In the simplest case, we consider a model with
homogeneous capture and survival probability (M1:
$\\{\rho_{t},\,\phi,\,p_{t}\\}$). When the population is supposed to be
structured in $G$ groups, we suppose that each group has its own time-varying
recruitment probability and that mixture components may vary by capture
probability (M2-M4: $\\{\rho_{t\times h_{G}},\,\phi,\,p_{t+h_{G}}\\}$, for
$G=2,3,4$), by survival probability (M5-M7: $\\{\rho_{t\times
h_{G}},\,\phi_{h_{G}},\,p_{t}\\}$, for $G=2,3,4$) or by both (M8-M10:
$\\{[\rho_{t\times h_{G}},\,\phi_{h},\,p_{t+h}]_{G}\\}$, for $G=2,3,4$). Each
simulated dataset was augmented by $500$ all-zero capture histories to
implement the PX-DA approach so that $D_{k}+500=M_{k}\neq M^{*}$ for all
simulated sets. We consider the prior setting described in Section 3. Notably,
we specify a $N(0,\,10)$ for the intercept $\mu$. When the survival
probability is the same for all the individuals, a standard Uniform prior is
placed on that parameter. When two survival probabilities are considered
(generically, $\phi_{1}$ and $\phi_{2}$), we use $\phi_{1}\sim Beta(1,\,2)$
and $\phi_{2}|\phi_{T}\sim tBeta(1,\,1,\,\phi_{T},\,1)$, with the latter
marginally yielding $\phi_{2}\sim Beta(2,\,1)$. Then enforcing
$\phi_{1}<\phi_{2}$ induces a slight repulsion between the two parameters.
When one considers more than two survival parameters, constrained Uniform
priors are instead chosen. These different prior choices for the models do not
substantially affect the model selection criterion associated with each model
but - conversely - are useful to induce a better separation between couples of
survival parameters. Estimation is carried out using JAGS, in which we run $2$
chains with $20,000$ iterations each, discarding $5,000$ as burn-in and
thinning by $2$ the remainder to save storage space (Brooks et al., 2004).

We chose median posterior estimates instead of simple averages, mitigating the
effect of anomalies that can result in occasionally low-informative datasets.
In the same spirit, we rely on the Mean Absolute Error (MAE) as an accuracy
measure instead of the widely employed Root Mean Squared Error. Interval
estimation is assessed through the percentage of times the $95\%$ credible
intervals contain the true value of the parameter (i.e. the coverage) and the
average $95\%$ credible interval width (CIW). The overall goodness-of-fit is
measured via the Watanabe-Akaike Information Criterion (WAIC, Watanabe and
Opper (2010)), following the good practice of Gelman et al. (2014) whenever
finite mixtures are fitted. We also report the overlapping index ($OV$,
Pastore (2018); Pastore and Calcagnì (2019)) between the posterior
distributions of the two survival parameters, averaged over the $K=50$
replicas. If $OV=0$, the two distributions are completely separated, while if
$OV=1$, the two distributions perfectly overlap. This metric is particularly
appealing to understand whether the posterior distributions of the group-
specific parameters are well-separated or not, justifying the related model
parameterisation. Finally, we investigate the fuzzy classification ability of
the RPT model using multi-class AUC (mAUC) (Hand and Till, 2001). Let us
remark that the classification performances can only be evaluated in the
simulation setting and for the RPT model, as the true group labels are known
and are consistent with the estimated ones. Notice that individuals who have
not been observed are not provided with a capture history, thus, it is
impossible to infer the group they belong to. Therefore, the estimated mixture
weights (referred to as the full $M_{k}$-sized pseudo-population) are not
expected to align with the mixture weights employed in the data generation
process. Indeed, the latter was used to generate the cluster labels of the
whole pseudo-population, of which some individuals (i.e. the pseudo-
individuals) never become part of the real population.

Figure 2 shows the differences between the estimated ($\hat{N}_{{super}}$) and
true ($N_{{super}}$) super-population size for each of the $K=50$ replicas.
The error is divided by the expected value of $N_{\text{super}}$ in the
corresponding scenario to allow for a meaningful comparison between scenarios
having expected super-population sizes of different magnitudes.

Figure 2: Relative estimation error of the super-population ($N_{super}$)
abundance for increasing number of sampling occasions ($T$), calculated for
each of the $K=50$ independent replicas, by the RPT and the ten alternative
Pledger et al. (2010)’s mixture models.

We notice that the RPT uniformly provides the best results overall. On the
contrary, its competing models consistently underestimate or overestimate the
super-population size. Notably, as $T$ increases, the underestimation is more
evident for those models that do not account for heterogeneity in survival
probabilities (i.e. M1-M4) - cfr. Pledger et al. (2003, 2010) -, while the
overestimation is substantial for those that consider an excessive number of
parameters (i.e. M8-M10). Intuitively, a greater number of parameters
controlling the capture probabilities tends to infer an excessive number of
uncaught individuals from the zero-histories.

Table 1 reports useful summaries to assess the models’ performances in
estimating $N_{\text{super}}$. The MAE and CIW are scale-dependent measures
and do not allow for scenario comparisons. Therefore, we consider relative
versions of these measures by dividing both by the corresponding expected
value of $N_{\text{super}}$.

Table 1: Estimates of relative MAE (MAE${}_{\text{rel}}$), coverage (Cov.) and
relative average width of the $95\%$ credible intervals (CIW${}_{\text{rel}}$)
for $N_{super}$, median WAIC and percentage of times each competing model has
achieved best WAIC (%waic). All these summaries have been obtained when data
are simulated from the RPT model.

| $T=10$ | $T=20$ | $T=30$ | $T=40$
---|---|---|---|---
Model | MAE${}_{\text{rel}}$ | Cov. | CIW${}_{\text{rel}}$ | WAIC | %waic | MAE${}_{\text{rel}}$ | Cov. | CIW${}_{\text{rel}}$ | WAIC | %waic | MAE${}_{\text{rel}}$ | Cov. | CIW${}_{\text{rel}}$ | WAIC | %waic | MAE${}_{\text{rel}}$ | Cov. | CIW${}_{\text{rel}}$ | WAIC | %waic
RPT | 0.07 | 0.96 | 0.28 | 1433.96 | 6 | 0.04 | 0.96 | 0.17 | 2953.01 | 50 | 0.04 | 0.86 | 0.15 | 4353.83 | 86 | 0.04 | 0.92 | 0.15 | 5840.58 | 88
$\text{M}_{1}:\,\left\\{\rho_{t},\,\phi,\,p_{t}\right\\}$ | 0.06 | 0.70 | 0.15 | 1458.18 | 2 | 0.08 | 0.28 | 0.10 | 3175.72 | – | 0.13 | 0 | 0.06 | 4875.18 | – | 0.14 | 0 | 0.06 | 6668.57 | –
$\text{M}_{2}:\,\left\\{\rho_{t\times h_{2}},\,\phi,\,p_{t+h_{2}}\right\\}$ | 0.08 | 0.94 | 0.33 | 1448.27 | 8 | 0.05 | 0.92 | 0.16 | 2991.11 | – | 0.10 | 0.10 | 0.09 | 4495.79 | – | 0.11 | 0.24 | 0.10 | 6136.03 | –
$\text{M}_{3}:\,\left\\{\rho_{t\times h_{3}},\,\phi,\,p_{t+h_{3}}\right\\}$ | 0.15 | 0.88 | 0.69 | 1435.05 | 10 | 0.05 | 0.94 | 0.25 | 3012.28 | – | 0.08 | 0.60 | 0.15 | 4532.3 | – | 0.09 | 0.42 | 0.13 | 6133.42 | –
$\text{M}_{4}:\,\left\\{\rho_{t\times h_{4}},\,\phi,\,p_{t+h_{4}}\right\\}$ | 0.16 | 0.84 | 0.65 | 1439.58 | 32 | 0.06 | 0.92 | 0.28 | 2985.86 | – | 0.07 | 0.70 | 0.18 | 4510.22 | – | 0.08 | 0.64 | 0.17 | 6169.65 | –
$\text{M}_{5}:\,\left\\{\rho_{t\times h_{2}},\,\phi_{h_{2}},\,p_{t}\right\\}$ | 0.07 | 0.90 | 0.26 | 1448.64 | 14 | 0.08 | 0.90 | 0.27 | 3089.37 | 10 | 0.08 | 0.96 | 0.27 | 4710.57 | – | 0.10 | 0.76 | 0.27 | 6436.28 | –
$\text{M}_{6}:\,\left\\{\rho_{t\times h_{3}},\,\phi_{h_{3}},\,p_{t}\right\\}$ | 0.06 | 0.90 | 0.25 | 1438.10 | 4 | 0.08 | 0.90 | 0.26 | 3071.13 | 18 | 0.08 | 0.96 | 0.26 | 4706.34 | 6 | 0.10 | 0.74 | 0.27 | 6426.04 | 8
$\text{M}_{7}:\,\left\\{\rho_{t\times h_{4}},\,\phi_{h_{4}},\,p_{t}\right\\}$ | 0.06 | 0.92 | 0.26 | 1435.73 | 12 | 0.08 | 0.84 | 0.26 | 3094.53 | 22 | 0.09 | 0.92 | 0.27 | 4725.57 | 8 | 0.11 | 0.72 | 0.27 | 6452.76 | 4
$\text{M}_{8}:\,\left\\{\left[\rho_{t\times h},\,\phi_{t\times h},\,p_{t\times h}\right]_{2}\right\\}$ | 0.12 | 0.90 | 0.56 | 1440.04 | – | 0.05 | 0.96 | 0.21 | 3041.38 | – | 0.10 | 0.68 | 0.47 | 4668.08 | – | 0.12 | 0.64 | 0.64 | 6453.97 | –
$\text{M}_{9}:\,\left\\{\left[\rho_{t\times h},\,\phi_{t\times h},\,p_{t\times h}\right]_{3}\right\\}$ | 0.20 | 0.84 | 0.83 | 1430.92 | 4 | 0.28 | 0.62 | 0.76 | 2967.91 | – | 0.36 | 0.40 | 0.74 | 4444.86 | – | 0.41 | 0.20 | 0.66 | 6071.77 | –
$\text{M}_{10}:\,\left\\{\left[\rho_{t\times h},\,\phi_{t\times h},\,p_{t\times h}\right]_{4}\right\\}$ | 0.25 | 0.78 | 0.84 | 1437.25 | 8 | 0.32 | 0.52 | 0.85 | 2951.84 | – | 0.38 | 0.40 | 0.75 | 4424.14 | – | 0.42 | 0.14 | 0.59 | 6055.66 | –

The RPT model returns the lowest relative MAEs in all the scenarios that
involve more than 1 year of observation, while its competitors are associated
with errors, which mostly increase with $T$. The WAIC seems to fail in
selecting the RPT model when the number of sampling occasions is rather small
(i.e. $T=10$) by attaining the lowest values in correspondence of models that
do not account for survival heterogeneity (i.e. M1-M4 in $52\%$ of the
replicas). This is a reasonable outcome since it is extremely complicated (if
not impossible) to distinguish individuals with low and high survivals when
the occasions span only 1 year. Notably, in such a short-term scenario, the
low survival of a transient individual is likely to be confounded with an
exceptionally low capture probability. This would also explain why model
$\,\left\\{\rho_{t\times h_{4}},\,\phi,\,p_{t+h_{4}}\right\\}$ with four
component-specific capture probabilities is selected in almost $1$ out of $3$
replicas according to the WAIC. However, as $T$ increases, the WAIC tends to
favour the true model, yielding the lowest median score and selecting it
(i.e., returning the lowest WAIC) in most replicas. The classification
performances of the RPT model are quite good in all the scenarios, with the
mAUC improving as $T$ increases. Notably, the resulting median mAUC is always
$\geq 0.82$ and above $0.95$ when a year change occurs. Furthermore, by
assigning the labels to each encountered individual according to the Maximum a
Posteriori (MAP) rule, the median accuracy (across replicas) lies between
$75\%$ and $95\%$ in all scenarios (once again improving as $T$ increases).

Table 2 shows the RPT model performance in estimating some time-constant
parameters. The estimates of the undetectability and the survival parameters
have a very small mean absolute error, although they fail to attain the
nominal credible interval coverage of $95\%$ in most scenarios. This indicates
a good accuracy of the point estimates associated with over-confidence
(visible in the low average CIWs), thus reducing the nominal coverage.
Nonetheless, it settles to a fair and acceptable level. The result is not
particularly surprising since these component-specific parameters govern
latent ecological processes and are potentially mutually confounded.

Table 2: Estimates of MAE, coverage (Cov.) and average width of the $95\%$
credible intervals (CIW) for some time-constant parameters of model RPT. All
these summaries have been obtained when data are simulated from the RPT model.

| $\phi_{T}$ | $\phi_{NT}$ | $\delta$ | $\mu$
---|---|---|---|---
$T$ | MAE | Cov. | CIW | MAE | Cov. | CIW | MAE | Cov. | CIW | MAE | Cov. | CIW
10 | 0.30 | 0.82 | 0.88 | 0.04 | 0.62 | 0.09 | 0.06 | 0.94 | 0.25 | 0.27 | 0.94 | 0.98
20 | 0.08 | 0.88 | 0.35 | 0.005 | 0.68 | 0.01 | 0.03 | 0.82 | 0.10 | 0.13 | 0.82 | 0.43
30 | 0.07 | 0.80 | 0.25 | 0.002 | 0.76 | 0.007 | 0.03 | 0.62 | 0.07 | 0.08 | 0.90 | 0.32
40 | 0.06 | 0.82 | 0.20 | 0.001 | 0.80 | 0.05 | 0.02 | 0.74 | 0.06 | 0.07 | 0.96 | 0.27

Finally, we observe that the two survival probabilities (i.e. $\phi_{T}$ and
$\phi_{NT}$) are well separated in all scenarios, with an $OV$ equal to
$0.078$ when $T=10$ and equal to $0$ for larger $T$.

##  5 Real data analysis

We apply the RPT model to estimate the total population size of the common
bottlenose dolphins inhabiting the Tiber River estuary, as introduced in
Section 2. The data are the detection histories of $D=195$ well-marked
dolphins that have been sighted in the area between June 2018 and November
2020, for a total of $T=87$ occasions. After some preliminary runs with
different values of $M$, we finally set $M-D=500$ rows of pseudo-individuals
(i.e. with null capture histories), thus yielding $M=695$. Other choices of
$M$ led to similar results, with larger values of $M$ only straining the
computational burden, in terms of runtime and storage.

We run $2$ parallel chains, each with $20,000$ iterations with a burn-in of
$5,000$ iterations and no thinning. We compare the performances of the RPT
model with the wide range of alternative models illustrated in Section 4,
using the same prior setting specified in the simulation study for all the
considered models. Table 3 reports the results on abundance estimation along
with the WAIC associated with each competing model.

Table 3: WAIC, estimated total super-population abundance and by group ($95\%$
credible intervals).

Model | WAIC | $\hat{N}_{super}$
---|---|---
RPT | 4725.1 | 311 (266, 373)
$\left\\{\left[\rho_{t\times h},\,\phi_{t\times h},\,p_{t\times h}\right]_{3}\right\\}$ | 4801.7 | 600 (467, 669)
$\left\\{\left[\rho_{t\times h},\,\phi_{t\times h},\,p_{t\times h}\right]_{4}\right\\}$ | 4811.9 | 594 (409, 664)
$\left\\{\rho_{t\times h_{2}},\,\phi_{h_{2}},\,p_{t}\right\\}$ | 4844.6 | 341 (287, 418)
$\left\\{\rho_{t\times h_{3}},\,\phi,\,p_{t+h_{3}}\right\\}$ | 4870.3 | 274 (243, 299)
$\left\\{\rho_{t\times h_{2}},\,\phi,\,p_{t+h_{2}}\right\\}$ | 4877.9 | 276 (240, 307)
$\left\\{\rho_{t\times h_{4}},\,\phi_{h_{4}},\,p_{t}\right\\}$ | 4893.1 | 352 (281, 458)
$\left\\{\rho_{t\times h_{4}},\,\phi,\,p_{t+h_{4}}\right\\}$ | 4903.8 | 271 (244, 306)
$\left\\{\rho_{t\times h_{3}},\,\phi_{h_{3}},\,p_{t}\right\\}$ | 4932.2 | 334 (269, 441)
$\left\\{\left[\rho_{t\times h},\,\phi_{t\times h},\,p_{t\times h}\right]_{2}\right\\}$ | 5088.5 | 356 (267, 628)
$\left\\{\rho_{t},\,\phi,\,p_{t}\right\\}$ | 5108.7 | 230 (216, 242)

We notice that the RPT model does yield the lowest WAIC score. Interestingly,
the second best choice according to the WAIC is model
$\left\\{\left[\rho_{t\times h},\,\phi_{t\times h},\,p_{t\times
h}\right]_{3}\right\\}$, which indeed resembles the same structure as the RPT
model but without allowing for mixture components sharing common parameters.
This lack of parsimony results in an overestimation of the individuals in the
super-population, similar to what has been observed in the simulation study
(cfr. Section 4). Thus, these results suggest the presence of unobserved
heterogeneity in the population, which seems to be better described by the
more parsimonious structure of the RPT model than by a generic three-group
specification.

The two annual survival probabilities are well separated with an $OV=0$ and
posterior estimates are $\hat{\phi}_{T}=1.06\times 10^{-8}$
($CI_{0.95}=[0,2.8\times 10^{-6}]$) for the group of transient individuals and
$\hat{\phi}_{NT}=0.71$ ($CI_{0.95}=[0.62,0.80]$) for the resident and part-
time individuals. Notice that the estimated survival parameter
$\hat{\phi}_{T}$ is of little interpretability on the annual scale. However,
it corresponds to a probability of $0.73$ on the weekly scale and $0.26$ on
the monthly scale. The average capture probability of the resident and
transient individuals is $\hat{p}_{NP}=\text{logit}^{-1}(\hat{\mu})\approx
0.19$; the corresponding temporal variations, captured by $\tau_{t}$, result
in the time-dependent posterior distributions of $p_{NP,\,t}$ reported in
Figure 8 in the Appendix.

The parameter $\delta$ regulating the undetectability of the part-time
individuals is estimated to be $\hat{\delta}\approx 0.74$. This means that
individuals in that group are present in the area approximately for the $26\%$
of their lifetime. Although the estimates of the recruitment parameters have
little interpretation in the Royle and Dorazio (2008)’s considered framework,
it is worth noticing that, on average, the recruitment probabilities are
higher during the first year for the resident individuals, while approximately
constant for part-time and transient individuals (cfr. Figure 9 in the
Appendix). This is a model artefact motivated by all the individuals already
present in the population before the start of the survey (i.e. mostly
residents) and that are virtually recruited on the first occasion. If we
consider the second year only (most reliable in terms of recruitment
probability estimation and individual classification), the average recruitment
probabilities are approximately $7\times 10^{-4},1\times 10^{-3},6\times
10^{-3}$ for the three groups, respectively. This aligns with our
expectations, as recruiting more stable individuals is a slower process than
recruiting less stable ones.

The final abundance estimate of the super-population in the whole observation
window is of $\hat{N}_{super}=311$ ($CI_{0.95}=[266,373]$), with yearly
variations $\hat{N}_{y}$ ($y=2018,2019,2020$) that show a peak in 2019 and a
decrease in the last year of observation (see Figure 3(a)). It is, however,
interesting to look at Figure 3(b), where we report the yearly abundance
estimates by group to better understand the behaviour of the aggregated yearly
pattern in Figure 3(a)).

(a)

(b)

Figure 3: Point estimates and posterior $95\%$ credible intervals of the
yearly super-population size (a) and by group (b).

Indeed, it seems that transients’ abundance ($\hat{N}_{T}=174$
($CI_{0.95}=[123,236]$) is the main factor affecting the overall yearly
counts, especially in 2019 and 2020. They show a substantial decrease in 2020,
while the abundances of residents ($\hat{N}_{R}=57$ ($CI_{0.95}=[48,65]$) and
part-times ($\hat{N}_{P}=81$ ($CI_{0.95}=[55,101]$) group tend to remain
stable across years. Recall that these estimates only refer to well-marked
individuals, who represent only a proportion of the whole population ($\approx
60\%$ of all sighted individuals across the three years)333Inferring the
overall size from the well-marked abundance estimate is not the focus of this
paper. For details on this passage, see Wilson et al. (1999). The posterior
distributions of the yearly sizes of the three groups were obtained by
counting all the individuals present in the population at each iteration by
year and group.

By assigning the $195$ well-marked individuals observed between June 2018 and
November 2020 to a single group according to the Maximum a Posteriori (MAP)
allocation, we have that $51$ are assigned to the group of residents, $54$ to
the group of part-time and $90$ to the group of transients. Notice that the
MAP is a straightforward and well-established method to attain classification
in finite mixture modelling. However, it is known to be sub-optimal in some
contexts (Stephens, 2000; McLachlan et al., 2019) and other appealing methods
have been proposed in the recent literature (e.g. Wade and Ghahramani (2018)).
The problem of summarising membership probabilities into a crisp
classification is a challenging and long-debated issue on which there is no
general agreement. Nevertheless, results from the simulation study (see
Section 4) showed that the MAP is sufficiently reliable in the proposed
framework. On the other hand, we are aware that any procedure deriving a crisp
classification from membership probabilities over-simplifies the complexity of
the inferred results as it does not give any information about the strength
each individual is assigned to a specific group. The latter is one of the main
advantages of the soft clustering returned by finite mixture models, and we
wish to exploit its fuzziness to quantify how decisively each individual is
assigned to one group or the other. For instance, $90\%$ of the individuals
classified into the Resident group have been assigned to it with a probability
greater than $0.9$. This probability is greater than $0.54$ for the Part-Times
and $0.53$ for the Transient.

In this regard, it is possible to visualise the classification results in
Figure 4 using a ternary diagram, which allows quantifying the probability
that each individual belongs to each group. The figure highlights three well-
distinguishable capture history patterns that again comply with the typical
RPT behaviour. Residents are available in the area for the whole study period
and are spotted very frequently in subsequent sampling occasions; their
capture histories are very informative. Hence, their identification is clear-
cut (intense colour) in most cases. Part-time individuals are available in the
area for most of the study period and are frequently spotted, but not as often
as residents. Their classification is crisp if they have been encountered for
the first time toward the beginning of the study period, while it is dicey
when they have been encountered toward the end for the first time. Finally,
the transients show short capture histories with few captures that never cross
two years; this reflects that they are spotted only a few times and do not
visit (survive) the area for long. Given the little information provided by
their capture histories, their classification is slightly vaguer compared to
the other groups, especially for individuals only observed in the last
sampling occasions.

Figure 4: Individual cumulative frequencies of capture for all the encountered
individuals divided into the three groups defined by the RPT model. Posterior
allocation was based on MAP.

Finally, it is important to remark that the proportions of the three groups do
not align with the estimated weights $\hat{w}=(0.09,0.39,0.52)$. This is
because the prior weights are not a good indicator of the true impact of each
component in the super-population, but only of their impact on the augmented
dataset (as mentioned in Section 4). However, we can infer on the composition
of the $N-D$ uncaught individuals, namely those individuals that were
recruited (i.e. not pseudo-individuals) but were never captured. The model
allocates (on average) the uncaught individuals to the resident group for the
$2\%$, the part-time group for the $15\%$ and the transient group for the
$83\%$. This result complies with the short-term survival of the transient
individuals and the more elusive nature of the part-time ones. On the other
hand, the residents are always present in the area and, therefore, more
susceptible to capture in a broad sense. The results suggest that most of the
resident individuals ($98\%$) have already been observed and documented, a
result which is in line with the progressive decrease of the discovery rate
highlighted in Section 2 (cfr. Figure 1(a)).

##  6 Discussion

Estimating the abundance of marine wildlife species is a challenging but
critical activity that can tell much about the undergoing ecological
processes. Thus, combining high-quality data with solid analytical approaches
is essential to improve our knowledge of these dynamics and increase the
potential for management actions (Vella et al., 2021; Lin et al., 2022). In
this work, we pursued a Bayesian estimation of the size of a common bottlenose
dolphin population, which is organised into three groups with different
residency and site-fidelity patterns. Accounting for such unobserved
heterogeneity is a very common problem in the environmental literature.
However, few papers approach the problem from the Bayesian perspective and
develop ad-hoc solutions based on prior scientific knowledge of the population
of interest.

We proposed a parsimonious specification of a finite mixture model within the
PX-DA setting for CR analysis, which we named RPT. This specification reflects
the typical bottlenose dolphin residency pattern, with individuals showing
high, partial or low site-fidelity (Dinis et al., 2016; Hunt et al., 2017;
Haughey et al., 2020; La Manna et al., 2022). Its adoption, characterised by
fewer free parameters, simplifies the identification of all the model
components compared to more generic and flexible alternatives. Furthermore,
while the application of finite mixture models to bring evidence about the
true existence of different population groups has sometimes been discouraged
(Pledger, 2000; Pledger et al., 2003), here we have shown how they can be
exploited to identify classes of individuals sharing similar profiles whenever
strong scientific evidence of population groups’ existence is available.

We devoted part of the model description to the discussion of a suitable prior
elicitation preventing the label-switching issue of FMM. Formal derivations
have enriched the results on constrained Beta priors from Alaimo Di Loro et
al. (2022). The performances of our proposal have been assessed through a
simulation study that considered scenarios with the increasing length of the
observation window and, hence, of the capture histories (i.e. mimicking
continuous monitoring of a population of interest over time). When data were
generated from the RPT model, we evaluated the influence of the number of
sampling occasions on the estimates’ quality. As expected, yielded better
accuracy and coverage thanks to the larger sample size. Comparison with
alternative model versions over many replicas showed that the model estimation
exhibited satisfying and robust performances when the observation window was
long enough (more than $1$ year of monitoring).

We estimated the RPT model and other alternatives on the data that motivated
our work. The results showed that, in terms of the WAIC score, our model
outplays several well-established competitors in the class of JS-type models.
The results correspond to biologically meaningful findings and align with
previous research work. The estimated abundance is of
$\hat{N}^{*}_{super}=311$, of which $51$ are estimated to be resident members
of the population. This quantity is particularly important as it is a proxy of
the breeding population in the area. The yearly trend shows an increase
between 2018 and 2019 and a decrease between 2019 and 2020. While the
estimated overall abundances are close to the ones obtained by Pace et al.
(2021), the trend is slightly in contrast with their results (non-significant
increasing trend between 2019 and 2020). This difference must be due to the
better quantification of the overall uncertainty in our analysis. This is
particularly important in this context as the major interest lies in
estimating a reasonable range that accounts for both the worst and best case
scenario. Indeed, we estimate that the driving factor of the variations in the
population size across the years is due to the abundance variation of the
transient individuals, the proportion of which is dynamic and not static.
Changes in the presence of such dolphins in an area may reflect changes in
driving features like habitat quality, prey distribution, or anthropogenic
disturbance. The year 2020 was the COVID-19 pandemic, with changes in
ecological conditions and human pressures on coastal marine waters (Carome et
al., 2022). Although difficult to assess, the influence of these factors on
the number of transient individuals cannot be ruled out, as several dolphins
may have been induced to reduce their mobility toward the study area by
improved conditions in their native habitat.

One basic assumption of CR experiments that we do not drop is that captures of
different individuals are independent. However, it is well known that
bottlenose dolphin populations can form structured societies with complex
social networks (Pace et al., 2022a). In future studies, we would like to drop
the independence assumption and include information on the population’s social
structure and the consequent statistical dependence in the capture histories.
Including the effect of external or individual covariates would also be
interesting. For instance, it is possible to recognise the gender and age
class in high-quality pictures. This partial information may be incorporated
in the Bayesian framework and could help for a better assessment, for example,
of the membership of different individuals in different groups or of their
marking probability (Wu et al., 2021). The capture-recapture model to estimate
the abundance of the bottlenose dolphins could be extended by adopting a
stopover model (Pledger et al., 2009; Worthington et al., 2019), which allows
capture and survival probabilities to depend both on time and time since
arrival in the population. In that case, the model could be formulated as a
multi-state Hidden Markov model (Worthington et al., 2019), where the
different states refer to the different groups of individuals in the
population. Last but not least, it would be extremely interesting to conduct
the survey on a larger spatial scale and include external information to
account for the spatial heterogeneity (Wu and Holan, 2017).

## Conflict of interest statement

The authors declare no conflict of interest.

## Data availability statement

The data that supports the findings of this study are available in the
supplementary material of this article.

## References

  * Alaimo Di Loro et al. (2022) P. Alaimo Di Loro, G. Caruso, M. Mingione, G. Jona Lasinio, and L. Tardella.  Specification of informative priors for capture-recapture finite mixture models.  In _Book of short paper – SIS 2022_. Pearson, 2022.
  * Altieri et al. (2022) L. Altieri, A. Farcomeni, and D. Alunni Fegatelli.  Continuous time-interaction processes for population size estimation, with an application to drug dealing in Italy.  _Biometrics_ , 2022.
  * Alunni Fegatelli and Tardella (2013) D. Alunni Fegatelli and L. Tardella.  Improved inference on capture recapture models with behavioural effects.  _Statistical Methods & Applications_, 22(1):45–66, 2013.
  * Amstrup et al. (2010) Steven C Amstrup, Trent L McDonald, and Bryan FJ Manly.  _Handbook of capture-recapture analysis_.  Princeton University Press, 2010.
  * Arnold et al. (2010) R. Arnold, Y. Hayakawa, and P. Yip.  Capture–recapture estimation using finite mixtures of arbitrary dimension.  _Biometrics_ , 66(2):644–655, 2010.
  * Böhning and van der Heijden (2009) Dankmar Böhning and Peter G. M. van der Heijden.  Recent developments in life and social science applications of capture-recapture methods.  _AStA Advances in Statistical Analysis_ , 93(1):1, 2009.
  * Böhning et al. (2018) Dankmar Böhning, John Bunge, and Peter G. M. van der Heijden.  _Capture-recapture methods for the social and medical sciences_.  CRC Press Boca Raton, 2018.
  * Böhning et al. (2020) Dankmar Böhning, Irene Rocchetti, Antonello Maruotti, and Heinz Holling.  Estimating the undetected infections in the covid-19 outbreak by harnessing capture–recapture methods.  _International Journal of Infectious Diseases_ , 97:197–201, 2020.  ISSN 1201-9712.  doi:https://doi.org/10.1016/j.ijid.2020.06.009.  URL https://www.sciencedirect.com/science/article/pii/S1201971220304446.
  * Bravington et al. (2016) M. V. Bravington, H. J. Skaug, and E. C. Anderson.  Close-kin mark-recapture.  _Statistical Science_ , 31(2):259–274, 2016.
  * Brittain and Böhning (2009) Sarah Brittain and Dankmar Böhning.  Estimators in capture–recapture studies with two sources.  _AStA Advances in Statistical Analysis_ , 93(1):23–47, 2009.
  * Brooks et al. (2000) S. P. Brooks, E. A. Catchpole, Byron J. T. Morgan, and S. C. Barry.  On the Bayesian analysis of ring-recovery data.  _Biometrics_ , 56(3):951–956, 2000.
  * Brooks et al. (2004) SP Brooks, Ruth King, and BJT Morgan.  A bayesian approach to combining animal abundance and demographic data.  _Animal biodiversity and conservation_ , 27(1):515–529, 2004.
  * Carome et al. (2022) William Carome, Elisabeth Slooten, William Rayment, Trudi Webster, Lindsay Wickman, Tom Brough, and Stephen M Dawson.  A long-term shift in the summer distribution of Hector’s dolphins is correlated with an increase in cruise ship tourism.  _Aquatic Conservation: Marine and Freshwater Ecosystems_ , 32(10):1660–1674, 2022.
  * Chao et al. (2001) A. Chao, P. K. Tsay, S. H. Lin, W. Y. Shau, and D. Y. Chao.  The applications of capture-recapture models to epidemiological data.  _Statistics in medicine_ , 20(20):3123–3157, 2001\.
  * Chung et al. (2004) H. Chung, E. Loken, and J. L. Schafer.  Difficulties in drawing inferences with finite-mixture models: a simple example with a simple solution.  _The American Statistician_ , 58(2):152–158, 2004\.
  * Clark and Gelfand (2006) James S Clark and Alan E Gelfand.  _Hierarchical modelling for the environmental sciences: statistical methods and applications_.  OUP Oxford, 2006.
  * de Valpine et al. (2017) Perry de Valpine, Daniel Turek, Christopher J Paciorek, Clifford Anderson-Bergman, Duncan Temple Lang, and Rastislav Bodik.  Programming with models: writing statistical algorithms for general model structures with nimble.  _Journal of Computational and Graphical Statistics_ , 26(2):403–413, 2017.
  * Di Cecco et al. (2020) Davide Di Cecco, Marco Di Zio, and Brunero Liseo.  Bayesian latent class models for capture–recapture in the presence of missing data.  _Biometrical Journal_ , 62(4):957–969, 2020.
  * Diebolt and Robert (1994) Jean Diebolt and Christian P Robert.  Estimation of finite mixture distributions through bayesian sampling.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 56(2):363–375, 1994.
  * Dinis et al. (2016) A. Dinis, F. Alves, C. Nicolau, C. Ribeiro, M. Kaufmann, A. Cañadas, and L. Freitas.  Bottlenose dolphin Tursiops truncatus group dynamics, site fidelity, residency and movement patterns in the Madeira Archipelago (North-East Atlantic).  _African Journal of Marine Science_ , 38(2):151–160, 2016.
  * Dorazio (2020) R. M. Dorazio.  Objective prior distributions for Jolly-Seber models of zero-augmented data.  _Biometrics_ , 76(4):1285–1296, 2020.
  * Dorazio and Royle (2003) R. M. Dorazio and J. A. Royle.  Mixture models for estimating the size of a closed population when capture rates vary among individuals.  _Biometrics_ , 59(2):351–364, 2003.
  * Farcomeni (2022) Alessio Farcomeni.  How many refugees and migrants died trying to reach europe? joint population size and total estimation.  _The Annals of Applied Statistics_ , 16(4):2339–2351, 2022.
  * Gelman et al. (2014) A. Gelman, J. Hwang, and A. Vehtari.  Understanding predictive information criteria for Bayesian models.  _Statistics and computing_ , 24(6):997–1016, 2014\.
  * Ghosh and Norris (2005) S. K. Ghosh and J. L. Norris.  Bayesian capture-recapture analysis and model selection allowing for heterogeneity and behavioral effects.  _Journal of Agricultural, Biological, and Environmental Statistics_ , 10(1):35–49, 2005.
  * Gimenez et al. (2018) Olivier Gimenez, Emmanuelle Cam, and Jean Michel Gaillard.  Individual heterogeneity and capture–recapture models: what, why and how?  _Oikos_ , 127(5):664–686, 2018.
  * Guéry et al. (2017) L. Guéry, S. Descamps, R. Pradel, S. A. Hanssen, K. E. Erikstad, G. W. Gabrielsen, H. G. Gilchrist, and J. Bêty.  Hidden survival heterogeneity of three Common eider populations in response to climate fluctuations.  _Journal of Animal Ecology_ , 86(3):683–693, 2017\.
  * Hand and Till (2001) D. J. Hand and R. J. Till.  A simple generalisation of the area under the ROC curve for multiple class classification problems.  _Machine learning_ , 45(2):171–186, 2001.
  * Haughey et al. (2020) Rebecca Haughey, Tim Hunt, Daniella Hanf, Robert William Rankin, and Guido J Parra.  Photographic capture-recapture analysis reveals a large population of Indo-Pacific bottlenose dolphins (Tursiops aduncus) with low site fidelity off the North West Cape, Western Australia.  _Frontiers in Marine Science_ , 6:781, 2020.
  * Hunt et al. (2017) T. N. Hunt, L. Bejder, S. J. Allen, R. W. Rankin, D. Hanf, and G. J. Parra.  Demographic characteristics of Australian humpback dolphins reveal important habitat toward the southwestern limit of their range.  _Endangered Species Research_ , 32:71–88, 2017.
  * Kendall et al. (1997) William L Kendall, James D Nichols, and James E Hines.  Estimating temporary emigration using capture–recapture data with pollock’s robust design.  _Ecology_ , 78(2):563–578, 1997.
  * King and McCrea (2019) R. King and R. McCrea.  Chapter 2 - Capture–Recapture Methods and Models: Estimating Population Size.  In Arni S.R. Srinivasa Rao and C.R. Rao, editors, _Integrated Population Biology and Modeling, Part B_ , volume 40 of _Handbook of Statistics_ , pages 33–83. Elsevier, 2019.  doi:https://doi.org/10.1016/bs.host.2018.09.006.
  * La Manna et al. (2022) G. La Manna, N. Rako-Gospić, D. S. Pace, S. Bonizzoni, L. Di Iorio, L. Polimeno, F. Perretti, F. Ronchetti, G. Giacomini, G. Pavan, G. Pedrazzi, H. Labach, and G. Ceccherelli.  Determinants of variability in signature whistles of the Mediterranean common bottlenose dolphin.  _Scientific Reports_ , 12(1):1–16, 2022.
  * Le Cren (1965) E David Le Cren.  A note on the history of mark-recapture population estimates.  _Journal of animal ecology_ , 34:453–454, 1965.
  * Lin et al. (2022) W. Lin, R. Zheng, B. Liu, S. Chen, M. Lin, M. Liu, W. Liu, and S. Li.  Low Survivals and Rapid Demographic Decline of a Threatened Estuarine Delphinid.  _Frontiers in Marine Science_ , 9, 2022.  ISSN 2296-7745.  doi:10.3389/fmars.2022.782680.
  * Maruotti et al. (2023) Antonello Maruotti, Dankmar Böhning, Irene Rocchetti, and Massimo Ciccozzi.  Estimating the undetected infections in the monkeypox outbreak.  _Journal of Medical Virology_ , 95(1):e28099, 2023\.
  * Matechou and Argiento (2022) Eleni Matechou and Raffaele Argiento.  Capture-recapture models with heterogeneous temporary emigration.  _Journal of the American Statistical Association_ , pages 1–14, 2022\.
  * Matechou et al. (2013) Eleni Matechou, Byron JT Morgan, Shirley Pledger, JA Collazo, and JE Lyons.  Integrated analysis of capture–recapture–resighting data and counts of unmarked birds at stop-over sites.  _Journal of Agricultural, Biological, and Environmental Statistics_ , 18:120–135, 2013.
  * McLachlan et al. (2019) G. J. McLachlan, Sharon X. Lee, and S. I. Rathnayake.  Finite mixture models.  _Annual review of statistics and its application_ , 6:355–378, 2019.
  * Mena and Walker (2015) R. H. Mena and S. G. Walker.  On the Bayesian mixture model and identifiability.  _Journal of Computational and Graphical Statistics_ , 24(4):1155–1169, 2015.
  * Morin et al. (2016) D. J. Morin, M. J. Kelly, and L. P. Waits.  Monitoring coyote population dynamics with fecal DNA and spatial capture–recapture.  _The Journal of Wildlife Management_ , 80(5):824–836, 2016.
  * Otis et al. (1978) D. L. Otis, K. P. Burnham, G. C. White, and D. R. Anderson.  _Statistical inference from capture data on closed animal populations_ , 1978.
  * Pace et al. (2021) D. S. Pace, C. Di Marco, G. Giacomini, S. Ferri, M. Silvestri, E. Papale, E. Casoli, D. Ventura, M. Mingione, P. Alaimo Di Loro, G. Jona Lasinio, and G. Ardizzone.  Capitoline Dolphins: Residency Patterns and Abundance Estimate of Tursiops truncatus at the Tiber River Estuary (Mediterranean Sea).  _Biology_ , 10(4):275, 2021.
  * Pace et al. (2022a) D. S. Pace, S. Ferri, G. Giacomini, C Di Marco, E. Papale, M. Silvestri, G. Pedrazzi, D. Ventura, E. Casoli, and G. Ardizzone.  Resources and population traits modulate the association patterns in the common bottlenose dolphin living nearby the Tiber river estuary (Mediterranean Sea).  _Frontiers in Marine Science_ , page 1506, 2022a.
  * Pace et al. (2022b) D. S. Pace, G. Panunzi, A. Arcangeli, S. Moro, G. Jona-Lasinio, and S. Martino.  Seasonal distribution of an opportunistic apex predator (Tursiops truncatus) in marine coastal habitats of the Western Mediterranean Sea.  _Frontiers in Marine Science_ , page 1785, 2022b.
  * Pace et al. (2022c) D. S. Pace, C. Tumino, M. Silvestri, G. Giacomini, G. Pedrazzi, G. Pavan, E. Papale, M. Ceraulo, G. Buscaino, and G. Ardizzone.  Bray-Call Sequences in the Mediterranean Common Bottlenose Dolphin Tursiops truncatus Acoustic Repertoire.  _Biology_ , 11(3):367, 2022c.
  * Papale et al. (2021) E. Papale, G. Alonge, F. Caruso, R. Grammauta, S. Mazzola, B. Mussi, D. S. Pace, and G. Buscaino.  The higher, the closer, the better? Influence of sampling frequency and distance on the acoustic properties of short-beaked common dolphins burst pulses in the Mediterranean Sea.  _Aquatic Conservation: Marine and Freshwater Ecosystems_ , 31:51–60, 2021.
  * Pastore (2018) M. Pastore.  Overlapping: a R package for estimating overlapping in empirical distributions.  _Journal of Open Source Software_ , 3(32):1023, 2018.
  * Pastore and Calcagnì (2019) M. Pastore and A. Calcagnì.  Measuring distribution similarities between samples: a distribution-free overlapping index.  _Frontiers in psychology_ , 10:1089, 2019.
  * Pledger (2000) S. Pledger.  Unified maximum likelihood estimates for closed capture–recapture models using mixtures.  _Biometrics_ , 56(2):434–442, 2000.
  * Pledger (2005) S. Pledger.  The performance of mixture models in heterogeneous closed population capture–recapture.  _Biometrics_ , 61(3):868–873, 2005.
  * Pledger et al. (2003) S. Pledger, K. H. Pollock, and J. L. Norris.  Open capture-recapture models with heterogeneity: I. Cormack-Jolly-Seber model.  _Biometrics_ , 59(4):786–794, 2003.
  * Pledger et al. (2010) S. Pledger, K. H. Pollock, and J. L. Norris.  Open capture–recapture models with heterogeneity: II. Jolly–Seber model.  _Biometrics_ , 66(3):883–890, 2010.
  * Pledger et al. (2009) Shirley Pledger, Murray Efford, Kenneth Pollock, Jaime Collazo, and James Lyons.  Stopover duration analysis with departure probability dependent on unknown time since arrival.  _Modeling demographic processes in marked populations_ , pages 349–363, 2009.
  * Plummer (2003) M. Plummer.  JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling.  In _Proceedings of the 3rd international workshop on distributed statistical computing_ , volume 124, pages 1–10. Vienna, Austria., 2003.
  * Richardson and Green (1997) Sylvia Richardson and Peter J Green.  On bayesian analysis of mixtures with an unknown number of components (with discussion).  _Journal of the Royal Statistical Society: series B (statistical methodology)_ , 59(4):731–792, 1997.
  * Royle and Dorazio (2008) J. A. Royle and R. M. Dorazio.  _Hierarchical modeling and inference in ecology: the analysis of data from populations, metapopulations and communities_.  Elsevier, 2008.
  * Royle and Dorazio (2012) J. A. Royle and R. M. Dorazio.  Parameter-expanded data augmentation for Bayesian analysis of capture–recapture models.  _Journal of Ornithology_ , 152(2):521–537, 2012\.
  * Royle et al. (2009) J. A. Royle, K. U. Karanth, A. M. Gopalaswamy, and N S. Kumar.  Bayesian inference in camera trapping studies for a class of spatial capture–recapture models.  _Ecology_ , 90(11):3233–3244, 2009.
  * Schwarz and Arnason (1996) C. J. Schwarz and A. N. Arnason.  A general methodology for the analysis of capture-recapture experiments in open populations.  _Biometrics_ , pages 860–873, 1996.
  * Silverman (2020) B. W. Silverman.  Multiple-systems analysis for the quantification of modern slavery: classical and Bayesian approaches.  _Journal of the Royal Statistical Society: Series A (Statistics in Society)_ , 183(3):691–736, 2020.
  * Stephens (2000) M. Stephens.  Dealing with label switching in mixture models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 62(4):795–809, 2000.
  * Turek et al. (2021) D. Turek, C. Wehrhahn, and O. Gimenez.  Bayesian non-parametric detection heterogeneity in ecological models.  _Environmental and Ecological Statistics_ , 28(2):355–381, 2021.
  * Vella et al. (2021) A. Vella, S. Murphy, J. Giménez, R. de Stephanis, B. Mussi, J. G. Vella, K. Larbi Doukara, and D. S. Pace.  The conservation of the endangered Mediterranean common dolphin (Delphinus delphis): Current knowledge and research priorities.  _Aquatic Conservation: Marine and Freshwater Ecosystems_ , 31:110–136, 2021.
  * Wade and Ghahramani (2018) Sara Wade and Zoubin Ghahramani.  Bayesian cluster analysis: Point estimation and credible balls (with discussion).  _Bayesian Analysis_ , 13(2):559–626, 2018.
  * Watanabe and Opper (2010) S. Watanabe and M. Opper.  Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory.  _Journal of machine learning research_ , 11(12), 2010.
  * Wilson et al. (1999) B. Wilson, P. S. Hammond, and P. M. Thompson.  Estimating size and assessing trends in a coastal bottlenose dolphin population.  _Ecological applications_ , 9(1):288–300, 1999\.
  * Worthington et al. (2019) Hannah Worthington, Rachel McCrea, Ruth King, and Richard Griffiths.  Estimating abundance from multiple sampling capture-recapture data via a multi-state multi-period stopover model.  _Annals of Applied Statistics_ , 13(4), 2019.
  * Wu and Holan (2017) Guohui Wu and Scott H Holan.  Bayesian Hierarchical Multi-Population Multistate Jolly–Seber Models With Covariates: Application to the Pallid Sturgeon Population Assessment Program.  _Journal of the American Statistical Association_ , 112(518):471–483, 2017.
  * Wu et al. (2021) Guohui Wu, Scott H Holan, Alexis Avril, and Jonas Waldenström.  A Bayesian semiparametric Jolly–Seber model with individual heterogeneity: An application to migratory mallards at stopover.  _The Annals of Applied Statistics_ , 15(2):813–830, 2021.
  * Zhou et al. (2019) Ming Zhou, Rachel S McCrea, Eleni Matechou, Diana J Cole, and Richard A Griffiths.  Removal models accounting for temporary emigration.  _Biometrics_ , 75(1):24–35, 2019.

## Appendix

##  Appendix A Multi-state interpretation of the ecological process

The ecological model described by Royle and Dorazio (2008) can be also seen as
a multi-state process. At each occasion $t$ each individual in the augmented
matrix can be in one and only one of the following three states:

  1. 1.

it has never been part of the population

  2. 2.

it is part of the population

  3. 3.

it was part of the population, but now it is not

We first notice that when an individual becomes part of the population, it
cannot be recruited any more: following the notation introduced in Section
3.1, this implies that for $t>1$, $r_{it}$ and $z_{it}$ cannot be
simultaneously equal to $1$.

In the JS modelling framework individuals that leave the population cannot
return to it. Hence, state $3$ is an absorbing state. Let us momentarily
ignore the population heterogeneity (clustering structure) for the sake of
clarity. If we allow temporal heterogeneity, then the transition probability
matrix associated with the three states at times $t=2,\dots,T$ is:

$\blockarray{c@{\hskip 1.0pt}ccc@{\hskip 3.0pt}}&123\\\ \block{c@{\hskip
3.0pt}|@{\hskip 2.0pt}ccc@{\hskip 2.0pt}|}11-\rho_{t}\rho_{t}0\\\
20\phi_{t}1-\phi_{t}\\\ 3001\\\ $

where rows and columns represent the states at time $t$ and $t+1$,
respectively. At period $t=1$, all individuals can be recruited in the
population. As $t$ increases, more and more individuals enter the population
or, equivalently, are removed from state 1. All the observed individuals will
eventually be recruited in the population before time $T$, but not all of them
will leave it (may have survived to future, unobserved, periods). Remember
that individuals are exposed to capture, with probability $p_{t}$, only during
their transitory stay in state 2. Hence, also a portion of never observed
individuals may have been recruited into the population at some point without
ever being captured. They represent the unknown part we aim to estimate.

##  Appendix B More details on the prior specification

The standard solution is that of concatenating conditionally specified Uniform
distributions as follows:

$u_{1}\sim Unif(0,1),\quad u_{j}|u_{j-1}\sim Unif(u_{j-1},1)\qquad
j=2,\dots,G,$ (7)

where the $u_{j}$’s are, for instance, the survival probabilities (Turek et
al., 2021). While effective in imposing the constraint, the priors in (7) do
not allow for the inclusion of previous information that can ease parameters’
identification. Alaimo Di Loro et al. (2022) explore alternative conditional
prior specification that, while implementing the ordering constraint, allows
to control for the shape and first moments of the induced marginal prior
distributions:

$\pi_{u_{g}}(u_{g})=\int_{\mathcal{S}_{g-1}}\pi_{u_{1}}(u_{1})\prod_{j=2}^{g}\pi_{u_{j}|u_{j-1}}(u_{j}\,|\,u_{j-1})\,du_{1},\dots,u_{g-1}\,,$

where $\mathcal{S}_{g-1}$ is a simplex of order $g-1$. Possible choices are
the Beta and Truncated Beta distributions. The latter corresponds to the
following set of prior distributions:

$\begin{aligned} &u_{1}\sim
Beta(\alpha_{1},\,\beta_{1}),\,\quad&u_{g}|u_{g-1}\sim
tBeta(\alpha_{g},\,\beta_{g};u_{g-1},1)\end{aligned}\qquad g=2,\dots,G,$ (8)

where $tBeta(\alpha_{g},\,\beta_{g};l,1)$ denotes the Truncated Beta
distribution in $(l,1)$. Note that for $\alpha_{g}=\beta_{g}=1$ we obtain (7).
When $G=2$ and $\alpha_{2}=1$, the prior specification in (8) induces the
following marginal prior density on $u_{2}$:

$\pi_{u_{2}}(u_{2})=\frac{B(\alpha_{1},\beta_{1}-\beta_{2})}{B(\alpha_{1},\beta_{1})}\,\beta_{2}(1-u_{2})^{\beta_{2}-1}\,F_{Beta(\alpha_{1},\beta_{1}-\beta_{2})}(u_{2})\,,\quad\text{with
}\beta_{1}>\beta_{2}\,\,,$ (9)

where $F_{Beta(\alpha_{1},\beta_{1}-\beta_{2})}(\cdot)$ is the cdf of a
$Beta(\alpha_{1},\beta_{1}-\beta_{2})$. See Section B.1 for the formal proof.
When $\alpha_{1}=\beta_{2}=k$ and $\beta_{1}=k+1$, we have that
$\pi_{u_{1}}(u_{1})=Beta(u_{1}\,|\,k,\,k+1)$ and the density in (9) is a
$Beta(k+1,\,k)$. The two distributions are mirrored with respect to the
vertical line $v=0.5$ (see Figure 5(a) of the Appendix as an example). Hence,
we have $\mathbb{E}[u_{2}]=1-\mathbb{E}[u_{1}]$, for all $k>0$; if $k>1$ we
also have $Mode(u_{2})=1-Mode(u_{1})$. The low-parametrised structure induces
well-separated prior means or prior modes for $u_{1},u_{2}$ marginal
distributions, favoring the mixture components’ separation. Alternative
settings inducing well separated modes in beta-type priors are reported in
Section B.2 of the Appendix.

###  B.1 The Beta and Truncated Beta distribution

The density of a Beta random variable with shape $\alpha$ and rate $\beta$ is:

$Beta\left(z\,|\,\alpha,\beta\right)=\frac{\Gamma\left(\alpha+\beta\right)}{\Gamma(\alpha)\Gamma(\beta)}\cdot
z^{\alpha-1}\cdot(1-z)^{\beta-1},\quad z\in(0,1),$

where $\Gamma(\cdot)$ is the Euler Gamma function. The truncation requires
normalising the same distribution over the truncated domain. We say
$u_{2}|u_{1}$ is a truncated Beta in $(u_{1},1)$ when:

$\pi(u_{2}|u_{1})=tBeta(\alpha_{2},\beta_{2},\,u_{1},1)=\frac{1}{B(\alpha_{2},\beta_{2})}\,\frac{u_{2}^{\alpha_{2}-1}\,(1-u_{2})^{\beta_{2}-1}}{1-F_{Beta(\alpha_{2},\beta_{2})}(u_{1})}\,\mathbb{1}_{\\{u_{1},1\\}}(u_{2})\,,$

with $F_{Beta(\alpha_{2},\beta_{2})}(\cdot)$ being the cdf of a
$Beta(\alpha_{2},\beta_{2})$. Now, suppose that $u_{1}\sim
Beta(\alpha_{1},\beta_{1})$ and $u_{2}|u_{1}\sim
tBeta(\alpha_{2},\beta_{2},\,u_{1},1)$, then if $\alpha_{2}=1$ the marginal
prior distribution induced on $u_{2}$ is given by:

$\displaystyle\pi(u_{2})$
$\displaystyle=\int_{0}^{1}\,\pi(u_{2}|u_{1})\pi(u_{1})\,du_{1}$
$\displaystyle=\int_{0}^{1}\,\frac{\beta_{2}\,(1-u_{2})^{\beta_{2}-1}}{1-F_{Beta(1,\beta_{2})}(u_{1})}\,\mathbb{1}_{\\{u_{1},1\\}}(u_{2})\frac{1}{B(\alpha_{1},\beta_{1})}\times
u_{1}^{\alpha_{1}-1}(1-u_{1})^{\beta_{1}-1}\,du_{1}$
$\displaystyle=\frac{\beta_{2}\,(1-u_{2})^{\beta_{2}-1}}{B(\alpha_{1},\beta_{1})}\,\int_{0}^{u_{2}}\,\frac{u_{1}^{\alpha_{1}-1}(1-u_{1})^{\beta_{1}-1}}{1-F_{Beta(1,\beta_{2})}(u_{1})}\,du_{1}$
$\displaystyle=\frac{\beta_{2}\,(1-u_{2})^{\beta_{2}-1}}{B(\alpha_{1},\beta_{1})}\,\int_{0}^{u_{2}}\,u_{1}^{\alpha_{1}-1}(1-u_{1})^{(\beta_{1}-\beta_{2})-1}\,du_{1}$
$\displaystyle=\frac{B(\alpha_{1},\beta_{1}-\beta_{2})}{B(\alpha_{1},\beta_{1})}\,\beta_{2}\,(1-u_{2})^{\beta_{2}-1}\,\int_{0}^{u_{2}}\,\frac{u_{1}^{\alpha_{1}-1}(1-u_{1})^{(\beta_{1}-\beta_{2})-1}}{B(\alpha_{1},\beta_{1}-\beta_{2})}\,du_{1}$
$\displaystyle=\frac{B(\alpha_{1},\beta_{1}-\beta_{2})}{B(\alpha_{1},\beta_{1})}\,\beta_{2}(1-u_{2})^{\beta_{2}-1}\,F_{Beta(\alpha_{1},\beta_{1}-\beta_{2})}(u_{2})\,,\quad{with}\,\,\beta_{1}>\beta_{2},$

where we exploited the fact that
$F_{Beta(1,b)}(z)=\int_{0}^{z}\,(1-t)^{b-1}\,dt=1-(1-z)^{b}$. Observe that the
constraint $\beta_{1}>\beta_{2}$ is essential to avoid the divergence of the
beta function.

In particular, it is straightforward to show that if $\alpha_{1}=\beta_{2}=k$
and $\beta_{1}=k+1$, that is $u_{1}\sim Beta(k,k+1)$ and $u_{2}|u_{1}\sim
tBeta(1,k,\,u_{1},1)$, then $u_{2}\sim Beta(k+1,k)$. Notice that, in this
particular case, the marginal distribution induced on $u_{2}$ is symmetrical
with respect to the distribution of $u_{1}$ around the vertical line $v=0.5$;
equivalently, $u_{2}\overset{d}{=}1-u_{1}$.

###  B.2 The Beta and Restricted Beta

An alternative conditional specification that allows for a properly informed
marginal prior can be obtained using Restricted Beta distribution. Also known
as 4-parameters Beta, it is a Beta r.v. that has been shifted and scaled to
reside on the domain $(l,u)$:

$rBeta\left(z\,|\,\alpha,\beta,l,u\right)=\frac{\Gamma\left(\alpha+\beta\right)}{\Gamma(\alpha)\Gamma(\beta)}\cdot\frac{(z-l)^{\alpha-1}\cdot(u-z)^{\beta-1}}{(u-l)^{\alpha+\beta-1}}\,\mathbb{1}_{\\{l,u\\}}(z).$

We can use it to specify recursively a set of $G$ conditional priors as:

$\begin{aligned} &u_{1}\sim
Beta(\alpha_{1},\,\beta_{1}),\quad&u_{g}|u_{g-1}\sim
rBeta(\alpha_{g},\,\beta_{g};u_{g-1},1)\end{aligned}\qquad g=2,\dots,G,$

where $rBeta(\alpha_{g},\,\beta_{g};u_{g-1},1)$ denotes the Beta restricted to
$(u_{g-1},1)$. The corresponding joint prior is:

$\pi\left(u_{1},u_{2}\right)=\frac{\Gamma(\alpha_{1}+\beta_{1})\Gamma(\alpha_{2}+\beta_{2})}{\Gamma(\alpha_{1})\Gamma(\beta_{1})\Gamma(\alpha_{2})\Gamma(\beta_{2})}\cdot
u_{1}^{\alpha_{1}-1}(1-u_{1})^{\beta_{1}-\beta_{2}-\alpha_{2}}(u_{2}-u_{1})^{\alpha_{2}-1}(1-u_{2})^{\beta_{2}-1}\,\mathbb{1}_{\\{u_{1},1\\}}(u_{2}),$

which does not allow for an analytical marginalisation to get $\pi(u_{2})$ in
the general scenario. However, the $rBeta$ expected value and variance are
available in closed form and hence we can use the law of total expectation to
derive the marginal expected value and variance of all the components. For
$g=2,\dots,G$:

$\displaystyle\mathbb{E}\left[u_{g}\right]=\frac{\alpha_{g}}{\alpha_{g}+\beta_{g}}+\mu_{g-1}\frac{\beta_{g}}{\alpha_{g}+\beta_{g}}$
$\displaystyle\mathbb{V}\left[u_{g}\right]=\sigma^{2}_{g-1}\cdot\frac{\beta_{g}^{2}}{(\alpha_{g}+\beta_{g})^{2}}\left(1+\frac{\alpha_{g}}{\beta_{g}(\alpha_{g}+\beta_{g}+1)}\right)+\left(1-\mu_{g-1}\right)^{2}\cdot\frac{\alpha_{g}\beta_{g}}{(\alpha_{g}+\beta_{g})^{2}(\alpha_{g}+\beta_{g}+1)},$

where $\mu_{g-1}=\mathbb{E}[u_{g-1}]$ and
$\sigma^{2}_{g-1}=\mathbb{V}[u_{g-1}]$. Therefore, one can define a system of
equations to find the combination of $\alpha_{g},\beta_{g}$ that complies with
a prior knowledge on the moments of the parameters (see Figure 5(b) as an
example).

(a)

(b)

Figure 5: (a) Marginal priors induced on $u_{1}$ (solid red) and $u_{2}$
(dashed orange) by the conditional specification $u_{1}\sim
Beta(.5,1.5),\;u_{2}|u_{1}\sim tBeta(1,.5,u_{1},1)$. (b) Marginal priors
induced on $u_{1}$ (solid red) and $u_{2}$ (dashed orange) with
$\mathbb{E}[u_{1}]=0.2$, $\mathbb{E}[u_{2}]=0.4$ with $\mathbb{V}[u_{1}]=0.02$
and $\mathbb{V}[u_{2}]=0.04$ by the conditional prior specification $u_{1}\sim
Beta(1.4,5.6),\;u_{2}|u_{1}\sim rBeta(0.826,2.478,u_{1},1)$ .

#### A convenient parameter setting

A convenient parameter setting arises when $G=2$ and if we fix
$\alpha_{1}=\alpha$, $\alpha_{2}=1$, $\beta_{1}=\beta+1$ and
$\beta_{2}=\beta$. Indeed, if $u_{1}\sim Beta(\alpha,\beta+1)$ and
$u_{2}\,|\,u_{1}\sim rBeta\left(1,\beta,u_{1},1\right)$, then the joint
distribution of $(u_{1},u_{2})$ is:

$\pi\left(u_{1},u_{1}\right)=\frac{(\alpha+\beta)\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\cdot
u_{1}^{\alpha-1}(1-u_{2})^{\beta}.$

Marginalising with respect to $u_{1}$ we get:

$\displaystyle\pi\left(u_{2}\right)$
$\displaystyle=\int_{0}^{1}\pi(u_{1},u_{2})\,\mathbb{1}_{\\{u_{1},1\\}}(u_{2})\,du_{1}=\int_{0}^{1}\pi(u_{1},u_{2})\,\mathbb{1}_{\\{0,u_{2}\\}}(u_{1})\,du_{1}=\int_{0}^{u_{2}}\pi(u_{1},u_{2})\,du_{1}=$
$\displaystyle=\frac{(\alpha+\beta)\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\cdot(1-u_{2})^{\beta-1}\cdot\int_{0}^{u_{2}}u_{1}^{\alpha-1}\,du_{1}=$
$\displaystyle=\frac{\Gamma(\alpha+\beta+1)}{\Gamma(\alpha)\Gamma(\beta)}\cdot(1-u_{2})^{\beta-1}\cdot\frac{u_{2}^{\alpha}}{\alpha}=$
$\displaystyle=\frac{\Gamma(\alpha+\beta+1)}{\Gamma(\alpha+1)\Gamma(\beta)}\cdot
u_{2}^{(\alpha+1)-1}(1-u_{2})^{\beta-1},$

which is a standard Beta density
$\pi\left(u_{2}\right)=Beta(u_{2}\,|\,\alpha+1,\beta)$. Expected value and
variance can then be derived from basic properties of the Beta distribution.

##  Appendix C Modelling the temporary emigration

Here, we show that the introduction of the undetectability parameter $\delta$
on the part-time individuals is equivalent to allowing for random temporary
emigration. We focus on the higher hierarchy level of the part-time model
specification, as all other components are not affected. For the sake of
clarity, we drop the $g$-subscript and let the reference to the part-time
group be implied.

The proposed model specification for the part-time group is as follows:

$\displaystyle y_{it}\,|\,z_{it}\,\sim\,Bern(p_{t}\cdot(1-\delta)\cdot
z_{it}),$ (10) $\displaystyle
z_{it}\,|\,z_{i,t-1},r_{it}\,\sim\,Bern(\phi_{t}\cdot z_{i,t-1}+\rho_{t}\cdot
r_{it}),\quad r_{it}=\text{min}\\{r_{i,t-1},1-z_{i,t-1}\\},$

which is different from the other groups only through the introduction of the
parameter $\delta\in(0,1)$ in the detection process.

Let us recall that $z_{it}$ is a latent variable indicating whether individual
$i$ is “alive” at time $t$. This is generally confounded with permanent
emigration, but it cannot account for temporary emigration as exited
individuals cannot ever re-enter the study area and return susceptible to
captures. Therefore, this first latent variable is only able to model the time
at which individual $i$ starts visiting the area (is born) and the time at
which it stops visiting it for good (dies). The explicit modelling of
temporary emigration within this time window requires the introduction of an
additional latent variable $v_{it}$ denoting whether individual $i$ is present
given that $z_{it}=1$. The specification of Equation (10) arises if we assume
that temporary emigration occurs randomly and with equal probability $\delta$
while individual $i$ is alive. This corresponds to the following hierarchical
specification:

$\displaystyle y_{it}\,|\,v_{it}\,\sim\,Bern(p_{t}\cdot v_{it}),$ (11)
$\displaystyle v_{it}\,|\,z_{i,t}\,\sim\,Bern((1-\delta)\cdot z_{i,t}),$
$\displaystyle z_{it}\,|\,z_{i,t-1},r_{it}\,\sim\,Bern(\phi_{t}\cdot
z_{i,t-1}+\rho_{t}\cdot r_{it}),\quad
r_{it}=\text{min}\\{r_{i,t-1},1-z_{i,t-1}\\},$

where $y_{it}\perp z_{it}$ if $v_{it}$ is known. We can easily marginalise
$v_{it}$ out of Equation (11) by noting that each $y_{it}|z_{it}$ is a
Bernoulli random variable with probability of success:

$\displaystyle\widetilde{p}_{t}$
$\displaystyle=\mathbb{E}\left[y_{it}\,|\,z_{it}\right]\,=\,\mathbb{E}\left[\mathbb{E}\left[y_{it}\,|\,v_{it}\right]\,|\,z_{it}\right]=$
$\displaystyle=\mathbb{E}\left[p_{t}\cdot
v_{it}\,|\,z_{it}\right]=p_{t}\cdot\mathbb{E}\left[v_{it}\,|\,z_{it}\right]=p_{t}\cdot(1-\delta)\cdot
z_{it},$

from which:

$y_{it}\,|\,z_{it}\,\sim\,Bern(p_{t}\cdot(1-\delta)\cdot z_{it})$

This equivalence and the corresponding interpretation are what motivates the
use of a multiplicative parametrisation on the capture probability at the
invlogit scale in place of a more straightforward group-specific intercept
within the logit specification.

##  Appendix D Details about time lags used in the simulation study

The number of days between two consecutive capture occasions (daily time lags,
henceforth) within a single year has been simulated from a shifted geometric
distribution with probability $0.05$, which has an expected value equal to
$20$ and a standard deviation equal to $19.5$. The resulting random sequence
of time lags is

$(20,1,12,15,56,9,9,12,10)$

and, for scenarios that contemplate more than one year of study, the same
sequence is repeated during each new year. The shift of year occurring each
$10$ occasions is achieved by using a higher constant time lag (i.e. $240$
days) between the $(10k)$th occasion and the $(10k+1)$th occasion, with
$k=1,2,3$. This results in a scenario $k$ composed by $k$ years of study, for
$k=1,2,3,4$. For example, scenario 2 ($T=20$) is composed of the following
sequence of time lags, resulting in $2$ years of capture occasions:

$(20,1,12,15,56,9,9,12,10,240,20,1,12,15,56,9,9,12,10)$

##  Appendix E Convergence of relevant parameters estimated on real data

We checked the convergence of the parameter chains in the real data
application through the general-purpose Gelman diagnostic. All potential scale
reduction factors $\hat{R}$ are below $1.01$, which suggests a good mixing of
all parameter chains.

We show the behaviour of the traceplots and density plots of the most relevant
time-static parameters for the sake of saving space. We can observe how the
two chains explore the same parameters space in all cases and produce a well-
shaped posterior distribution, with no bad behaviour.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 6: Traceplots of posterior samples for the main parameters of interest
of model RPT.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 7: Densities of posterior samples for the main parameters of interest
of model RPT.

##  Appendix F Estimated time-varying parameters of the RPT model on the real
data application

Figure 8: Posterior estimates and $95\%$ credible intervals for capture
probabilities of resident and transient individuals at each sampling occasion.

(a)

(b)

(c)

Figure 9: Estimated recruitment probabilities (a) $\rho_{R,t}$, (b)
$\rho_{P,t}$ and (c) $\rho_{T,t}$.