A hierarchical approach for estimating state‐specific mortality and state transition in dispersing animals with incomplete death records

Unbiased mortality estimates are fundamental for testing ecological and evolutionary theory as well as for developing effective conservation actions. However, mortality estimates are often confounded by dispersal, especially in studies where dead‐recovery is not possible. In such instances, missing individuals (i.e. individuals with unobserved time of death) may have died or permanently emigrated from a study area, making inferences about their fate difficult. Mortality before and during dispersal, as well as the decision to disperse, usually depend on a suite of individual, social and environmental covariates, which in turn can be used to draw conclusions about the fate of missing individuals. Here, we propose a Bayesian hierarchical model that takes into account time‐varying covariates to estimate transitions between life‐history states and mortality in each state using mark‐resighting data with missing individuals. Specifically, our framework estimates mortality rates in two states (resident and dispersing state) by treating the fate of missing individuals as a latent (i.e. unobserved) variable that is statistically inferred based on information from individuals with a known fate and given the individual, social and environmental conditions at the time of disappearance. Our model also estimates rates of state transition (i.e. emigration) to assess whether a missing individual was more likely to have died or survived due to unobserved emigration from the study area. We used simulations to check the validity of our model and assessed its performance with data of varying degrees of uncertainty. Our modelling framework provided accurate mortality and emigration estimates for simulated data of different sample sizes, proportions of missing individuals, and resighting intervals. Variation in sample size appeared to affect the precision of estimated parameters the most. Our approach offers a solution to estimating unbiased mortality of both resident and dispersing individuals as well as the probability of emigration using mark‐resighting data with incomplete death records. Conditional on the availability of data on known‐fate individuals and relevant time‐varying covariates, our model can reconstruct the fate (death or emigration) of missing individuals. The modularity of our framework allows mortality analyses to be tailored to a variety of species‐specific life histories.


| INTRODUC TI ON
Mortality is a key demographic process and its unbiased estimation is paramount to understand population dynamics (Hanski, 1998;Saether & Bakke, 2000). Unbiased information about the magnitude of mortality, how it varies with age, geographical location and behaviour, and how this variation is influenced by environmental factors, is fundamental to demographic studies. In wild animal populations, however, empirical information on mortality is often incomplete as individuals frequently disappear under unexplained circumstances (Colchero & Clark, 2012;Schaub & Royle, 2014).
Assuming death or survival for these missing individuals (terms in bold at first occurrence are defined in the Glossary) can become problematic, especially for species that tend to disperse beyond well-monitored study areas . The resulting bias in mortality estimates can have severe implications for our ability to test ecological and evolutionary theory and to predict population dynamics for wildlife management purposes (Griffith et al., 2016).
A variety of methods have been developed to address the challenges of unbiased mortality estimation under state uncertainty or incomplete records. Multistate mark-resighting models (Arnason, 1973;Schwarz et al., 1993) were first introduced to deal with variation in mortality and detection probabilities arising from individuals in different states (Lebreton et al., 2009). Subsequent extensions of these models allowed estimating survival under state uncertainty (Dupuis, 1995;Pradel, 2005), particularly in cases where it is not possible to distinguish between death and permanent emigration from the study area (Dupont et al., 2022;Ergon & Gardner, 2014;Schaub & Royle, 2014). In field studies that collect resighting data on known individuals, additional uncertainties in the estimation of mortality can arise due to incomplete records of age and sex. For such datasets, approaches have been developed to infer mortality by estimating latent variables within a Bayesian framework (Clark et al., 2005;Colchero & Clark, 2012). For studies involving missing individuals, this Bayesian framework has been expanded to estimate mortality by treating the fate (i.e. death or emigration) of missing individuals as a latent variable, which was inferred based on the posterior probabilities of death and emigration Colchero et al., 2016;Colchero et al., 2021). Emigration, in this context, was derived from the observed age distributions of immigrants because empirical information on emigrants was unavailable. In most animal species, however, the decision to emigrate is usually not only driven by age but also related to a combination of intrinsic and extrinsic factors (Clobert et al., 2012). Therefore, if additional information is available on how and under what circumstances individuals disperse, a more reliable emigration probability can be assigned to missing individuals . Consequently, bias in mortality estimates can be further reduced.
Mortality typically changes with age but may also depend on an individual's 'state' (e.g. geographical location, behaviour, physiology or reproductive status) because such states are usually associated with different threats (Lebreton et al., 2009). Dispersing individuals face different mortality risks than resident conspecifics; for instance, unfamiliarity with a new area may result in reduced feeding efficiency (Pinter-Wollman et al., 2009), higher susceptibility to predation (Garrett & Franklin, 1988) and increased intraspecific aggression (Packer et al., 1988). Males of polygynous species generally have higher mortality during the mating period due to decreased feeding time and increased male-male competition for gaining access to oestrous females (Toïgo & Gaillard, 2003). State-specific differences in mortality are evident also in bird species. For instance, breeders usually suffer higher mortality than non-breeders due to increased exposure to predators during incubation (Low et al., 2010). An understanding and quantification of state-specific mortality is thus crucial as the contribution of individuals to the overall population viability varies substantially depending on their life-history state (e.g. dispersing vs. resident, breeding vs. non-breeding; Catlin et al., 2019).
State-specific mortality usually depends on a range of individual, social and environmental covariates. In yellow-bellied marmots Marmota flaviventer for example, state-specific mortality can be influenced by covariates such as body mass (Ozgul et al., 2010), the strength of social relationships (Montero et al., 2020) and drought severity (Cordes et al., 2020). Such covariates are, however, usually not constant over time and an omission of this temporal variability can seriously hamper predictions of population dynamics (Paniw et al., 2019). Considering time-varying covariates in mortality analyses is particularly important for species that live in seasonal (Cordes et al., 2020) or social environments (Angulo et al., 2018), or for highly mobile species that are found in spatially heterogeneous landscapes (Basille et al., 2013).
Monitoring free-ranging individuals under natural conditions is often difficult, leading to long and irregular time intervals between data on known-fate individuals and relevant time-varying covariates, our model can reconstruct the fate (death or emigration) of missing individuals. The modularity of our framework allows mortality analyses to be tailored to a variety of species-specific life histories.

K E Y W O R D S
Bayesian inference, capture-recapture, capture-recovery, interval censoring, multistate markresighting, Siler model, state-based mortality, true survival successive observations. In such situations, the exact time of an event (i.e. death, disappearance or transition between states) is frequently unobserved but known to lie within a certain time interval. The estimation of mortality and state transition can thus be biased if interval censoring is not explicitly modelled (Kalbfleisch & Prentice, 2011).
Here, we present a Bayesian hierarchical model to estimate state-specific mortality and state transition in relation to timevarying covariates for mark-resighting data with incomplete death records. We extend previous methods by considering (1)

| MATERIAL S AND ME THODS
We present a hierarchical model to estimate mortality of both resident and dispersing individuals as well as the probability of emigration. While doing so, we use latent variables to account for the fate of missing individuals, and time-varying covariates to account for the individual heterogeneity. In the following, we provide a general overview of all model components and assumptions, present a simulation reflecting a hypothetical life history of a social species that allows demonstrating the full potential of the model, describe the mathematical details of the model, and finally, validate and assess the performance of the model using the simulated data.

| Model overview
Our model is designed for mark-resighting studies characterized by the absence of dead-recovery and species that disperse beyond focal study areas, resulting in the unobserved disappearance of individuals. Consequently, our model is applicable to any study system to which the following properties and assumptions apply: (1) individuals can be in either a resident or dispersing state, (2) at least one of the two sexes may disperse, (3) individuals are identifiable upon resight, (4) individual identity and individual state are always recorded without error, and (5) individuals are sighted with certainty if they are alive and present in the focal study area.
In this study, we extended previous frameworks Colchero et al., 2016;Colchero et al., 2021) in several ways. Our model accommodates interval-censored data where the exact time of an event (i.e. death, emigration or disappearance from the study area) is unobserved but known to lie within a certain time interval. In addition to age and sex, the effects of other time-varying covariates on mortality are included in our framework.
Our model estimates emigration probabilities based on a set of proximate causes (reflected by using time-varying covariates) associated with emigration events instead of solely using age-based information from immigration events.
Based on the life history of the focal species, the properties of the empirical data, and the research question in mind, our model can be adapted to accommodate analyses of varying detail. A fully parameterised model allows: (1) estimating mortality probability of individuals that are in either a resident or dispersing state, (2) quantifying emigration probability, (3) accounting for the effects of timevarying covariates on these life-history events, (4) inferring the fate of individuals that disappeared with unknown fate, (5) estimating all previously mentioned parameters separately by sex, and (6) inferring the sex of individuals that died before their sex could be determined.

| Simulated life-history data
To demonstrate the potential of our model to estimate mortality while accounting for the fate of missing individuals, we simulated individual-based data of a social species that loosely follows the life history of lions, African wild dogs or dwarf mongooses. We did so because of a peculiarity of these species that allows inferring the fate not only of individuals that disappeared before dispersal but also of individuals that disappeared during dispersal. To understand how our model infers the fate of missing individuals in either state, we thus first qualitatively describe the simulated life history, while the computational details of the simulation are provided later (Section 2.4).
Our simulation consists of individuals of known age who can either be in a resident or a dispersing state at a given age. After birth, individuals belong to the resident state unless they emigrate and move to the dispersing state. In our simulation, individuals can emigrate as part of a multiple-member dispersing coalition, as is common in lions (Packer & Pusey, 1987), wild dogs  or dwarf mongooses (Waser et al., 1994). For such species, it is often observed that individuals who emigrate together can separate during dispersal, which results in a breakup of the dispersing coalition (Waser, 1996). Our simulation thus generates individual life histories (i.e. time of death, emigration and dispersing coalition breakup), while taking into account the effect of a single time-varying covariate on these life-history events. Although our model can handle cases in which individuals could disperse more than once and consequently be in the resident state more than twice (see Section 2.3.1), in our simulation we have only considered life histories with at most one dispersal event. A schematic representation of such a life history is shown in Figure 1a.
In our simulation, possible fates for missing individuals are as follows: A resident individual may have disappeared because it had died in the resident state or because it had permanently emigrated from the study area ( Figure 1b). A dispersing individual may have disappeared because it had died during dispersal or because it had permanently separated from its dispersing coalition as a result of a coalition breakup (Figure 1c). Therefore, the explicit modelling of emigration and dispersing coalition breakup allows inferring the fate of individuals that disappeared before or during dispersal. While our framework is designed as a two-sex model, only one sex is considered in our simulation for simplicity. Nonetheless, a schematic F I G U R E 1 Schematic representation of exemplary life histories (LHs): (a) LH of an individual that dispersed sometime between age y left and y right as part of a multiple-member dispersing coalition, separated from the original coalition between age z left and z right , settled between age t R left and age t R right , and died following its last observation as resident at age x left ; (b) LH of an individual that was a potential disperser when it disappeared under unexplained circumstances following its last observation as resident; (c) LH of an individual that was in the dispersing state at the time of last observation and part of a multiple-member dispersing coalition when it disappeared under unexplained circumstances; (d) LH of an individual that died prior to sex determination. Times τ refer to ages at which the focal individual was observed alive and relevant covariates were updated. A description of all remaining variables and parameters is provided in Table 1.  representation of how our model would infer the sex of an individual that died prior to sex determination is shown in Figure 1d.

| Model details
Below, we provide a description of the mathematical details un-

State-specific mortality
To describe how mortality changes as an animal ages, our model estimates the mortality hazard rate μ (Ergon et al., 2018) dependent on an individual's age t using the Siler model (Siler, 1979). Although other functions could be considered for modelling animal mortality (Colchero & Clark, 2012), the Siler model was chosen because it adequately reflects the bathtub-shaped mortality typical of most species (Gage & Dyke, 1986;Siler, 1979). The Siler model, with parameters a 0 , b 0 ∈ ℝ and a 1 , b 1 , c ∈ ℝ >0 , is composed of three additive components where the first component captures decreasing mortality of young individuals, the second component reflects age-independent mortality and the third component captures increasing mortality of old individuals (Siler, 1979). To estimate separate mortality parame-  (Ergon et al., 2018), which describes the probability to survive from birth to age x, as where X is a random variable for age at death. To model state-specific survival, the integral in Equation 2 (i.e. the cumulative hazard) is decomposed into periods during which the focal individual was either in the resident state or the dispersing state. Since individuals have different life histories, the survival function varies case by case. For example, if an individual dispersed sometime between age y left and age y right and settled between age t R left and age t R right (Figure 1a), the probability to survive from birth to the age of last observation as resident x left is given by

GLOSSARY
State: The behavioural state to which an individual belongs at a given time. Here, we distinguish between a resident state and a dispersing state.

State-specific mortality: Mortality hazard rates estimated separately for resident and dispersing individuals.
Missing individuals: Individuals lost to follow-up whose time of death was not observed.
State transition: An individual's transition from the resident state to the dispersing state (i.e. emigration or dispersal) or vice versa (i.e. settlement).

Emigration:
The movement of an individual away from its site of birth or natal group, corresponding to the transition from resident to dispersing state.

Dispersing coalition breakup:
For individuals that disperse in multiple-member dispersing coalitions, this corresponds to the separation of an individual from the original coalition.
Interval censoring: This refers to the case where the exact timing of a life-history event (e.g. death, emigration, disappearance from the study area) is unobserved but known to be within a certain time interval. For comparison, right censoring refers to the case where the exact time of death is also unknown but, in contrast to interval censoring, no information is available about the time interval in which death occurred. Thus, right censoring truncates (i.e. censors) the true survival time on the right side of the observed survival time, which corresponds to the age at last observation.

Indicators:
We consider three binary indicators: The sex (s i ) of an individual i, whether an individual died or dispersed after the last observation (d i ), and finally, whether or not an individual experienced a dispersing coalition breakup (b i ) after the last observation. If the value of such an indicator is known, we declare this with a subscript k (e.g. s k, i ); if it is unknown and hence estimated as a latent variable, we add a subscript u.

Latent variable:
Unobserved indicator of an individual that is inferred retrospectively. Based on the three indicators used in our model, we consider the following latent variables: The fate (d u, i ) of a missing resident individual, the fate (b u, i ) of a missing dispersing individual, and the sex (s u, i ) of a resident individual that died too young to be sexed.

Observed variable:
A variable that is recorded when individuals are sighted, such as an individual's age at first and last detection, as well as time-varying covariates.

Random variable:
A variable (e.g. age at death) that has a set of possible values (e.g. different ages), each with a different probability, that depend on a random event (e.g. mortality).
where y = 1 2 y left + y right and t R = 1 2 t R left + t R right . The terms of Equation 3 that are highlighted in blue and orange represent resident and dispersing states, respectively. Equation 3 can be expanded analogously if the individual disperses even more times than stated above. While our framework allows for an arbitrary number of state changes, we present it in this paper with at most one dispersal event for simplicity. Nonetheless, a description of how multiple dispersal events are implemented in our model is provided in SI Appendix S1.
Finally, the survival function can be used to estimate the probability that death occurs before age x, that is, the cumulative distribution function (CDF) for age at death For interval-censored data (Sparling et al., 2006), this CDF can be employed to calculate the probability of death occurring during the

Effects of time-varying covariates
Our model allows accounting for effects of time-varying covariates on mortality, emigration and dispersing coalition breakup. Two options for considering covariate effects on mortality are implemented: The user can choose between (1) assuming that covariates have a constant effect on mortality across all possible ages, or (2) assuming that covariates have an age-dependent effect on mortality. In the first option, mortality is described with a proportional hazard model (Kleinbaum & Klein, 2010) where v(t) is a vector function of time-varying covariates evaluated at age t, and β is a coefficient vector describing the effect size (i.e. the log hazard ratio) of these covariates. Figure  For individuals in the dispersing state, our model assumes that covariates have an age-independent effect on mortality because for most species the dispersal duration (i.e. the age range during which individuals disperse) is usually short (Clobert et al., 2012); hence, option 1 is implemented and the vector of mortality parameters associated with the dispersing For individuals in the resident state, our model offers the choice between both options, with option 2 implemented as the default, thus allowing the user to test whether covariates affect mortality differently for young and old individuals. Choosing option 2, the vector of mortality parameters associated with the resident state is To include the effects of time-varying covariates in the survival function, the integrals of the state-specific mortality hazards (Equation 3) are decomposed based on discrete times when covariates were recorded. Considering covariates that were recorded at update times = 1 , 2 , … , n (i.e. discrete ages at which the focal individuals were observed alive), the survival function for the life-history example in Equation 3 can be rewritten as Depending on the lower and upper limit of the separate integrals in Equation 8, we approximated them differently with regard to the covariate values: For integrals where the upper and lower limits are update times (e.g. τ j − 1 and τ j ), we used the mean of the covariate values recorded at those times, e.g. 1 2 v a j−1 + v a j . For integrals containing only one covariate update time (either in the lower or upper limit), we used covariate values recorded at the corresponding update time for the evaluation of the entire integral (e.g. v a (τ j + 1 ) for the integral with limits τ j + 1 and y). The exact forms of the integrals of mortality hazards for both age-dependent (Equation 7) and ageindependent (Equation 6) effects of time-varying covariates are given in SI Appendix S2.

Dispersal and dispersing coalition breakup
Of the two possible state transitions emigration (i.e. dispersal) and settlement, only the former is modelled within our framework. The reason for this is that information on emigration is needed to infer whether a missing individual that was last observed in the resident state has died or emigrated. Settlement is not explicitly modelled because it is assumed that individuals that disappeared in the dispersing state would be observed if they settled; thus, settlement is not considered as a reason for disappearance. However, our model assumes that for species with a similar life history as in our simulation and application example (Behr et al., in review), a missing dispersing individual could have disappeared due to death or due to a separation from the dispersing coalition. Therefore, dispersing coalition breakup is explicitly modelled. In the following, at most one emigration event per individual is considered, while multiple breakup events in the dispersing state are possible. Hence, the indicator d obs i , which describes if emigration was observed, can be 0 or 1, whereas b obs i , which describes the number of observed breakup events, can take (in principle) any non-negative integer value. A description of how multiple dispersal events are implemented in the full model is provided in SI Appendix S1.
Emigration and dispersing coalition breakup are modelled as fully parametric proportional hazard (FFPH) models. Hence, the emigration rate g(y) is estimated dependent on an individual's age y as where g 0 (y − y min ) represents the baseline emigration rate (i.e. the emigration rate that results when all time-varying covariates are at mean values, or no such covariates are considered) shifted by the minimum age y min individuals must reach to be able to disperse, v em (y) is a vector function of time-varying covariates and β em is a coefficient vector describing the effect size of these covariates. Following the findings on the species considered in the application example , our model calculates the baseline emigration rate using a log-logistic function, although other functions could be considered as outlined in Munda et al. (2012). The exact form of the baseline hazard can be found in SI Appendix S3. Hereinafter, we will omit the explicit dependence on the time-varying covariates for the sake of notational simplicity.
Based on the emigration rate, the probability that dispersal occurs before age y, or the CDF for age at emigration, can be calculated as  1). In a proportional hazards framework (Kleinbaum & Klein, 2010), the covariate effects are estimated as hazard ratios (HRs), which have a multiplicative effect on the baseline mortality rate (solid line): HR > 1 shifts the baseline mortality curve upwards (dotted line), which results in higher mortality compared to the 'baseline' (for a given age and keeping all other covariates unchanged). The opposite applies to HR < 1 (dashed line). For covariates with an age-independent effect on mortality (Equation 6), panel (a) shows how such covariates would alter the mortality curve proportionally with age. For covariates with an age-dependent effect (Equation 7), panels (b) and (c) exemplify how the mortality curve would change at 'young' and 'old' ages, respectively.

Mortality rate
Baseline mortality Hazard ratio > 1 Hazard ratio < 1 where Y is a random variable for age at emigration, and γ is a vector of dispersal parameters (i.e. β em and the parameters describing the baseline emigration rate). For interval-censored data (Sparling et al., 2006), this CDF can be employed to estimate the probability that dispersal occurs during the age interval (y left , y right ] where the individual is known to have been resident at age y left and to have emigrated sometime before age y right ≥ y left . For dispersing coalition breakups, a similar approach was followed. However, the difference is that our model assumes the breakup hazard to depend on the time an individual has been dispersing, rather than on the age of the individual. To reflect the dispersal duration, the estimation of the dispersing coalition breakup rate r(z| , y) and the corresponding CDF R(z| , y), where ρ is a vector of breakup parameters, are carried out similarly to the ones for emigration and are described in detail in SI Appendix S3. Further details on the calculations of the integrals in Equation 11 can be found in SI Appendix S3 as well. The sex of resident individuals that died prior to sex determination, s u (Figure 1d): This latent variable applies to individuals that died too young to be sexed.

Sex-specific estimation
With our model, all parameters (i.e. θ, γ and ρ) can be estimated separately by sex. The corresponding parameter vectors are indicated with a subscript f for female and m for male (e.g. θ f and θ m ).

| Full Bayesian model
The  (Clark, 2020;Gelfand & Smith, 1990) is implemented in our framework. This involves using the Metropolis-Hastings algorithm to sample from the posterior distribution of one of the unknown parameters, conditional on the current values of the others, and repeating this for each parameter (see Section 2.3.3). In the following, we will therefore introduce all six conditional posterior distributions of θ, γ, ρ, s u , d u and b u , and we will again omit the explicit dependence on time-varying covariates for notational simplicity.

Posterior for mortality parameters
Mortality influences the distribution of X (age at death), which is unobserved but in turn determines the observations x left and x right . To construct the likelihood function of the mortality parameters s i of an individual i with known or assumed sex s i that was alive at the age of the first observation t F i , a different probability is assigned depending on whether that individual died during the age interval the second line of the likelihood in Equation 12, otherwise the first one. Divided by the prior p(θ | θ p ), where θ p is a vector of prior hyper-parameters, the posterior of mortality parameters is

Posterior for dispersal and breakup parameters
The dispersal parameters determine the distribution of Y (age at dispersal), which is unobserved but in turn determines the observations x left , x right , y left , and y right . For an individual i with (assumed or known) sex s i that is a potential disperser (i.e. it disappeared later than the minimum dispersal age y min and did not enter the study as an immigrant), our model defines the dispersal likelihood by distinguishing between two cases: (1) the focal individual dispersed during the age interval y left

Posterior for unknown sexes
The conditional posterior of the unknown sex indicator s u, i can be calculated for an individual i whose sex was unknown at the age of the last observation x left i . Note that in our model, individuals can only possibly be unsexed if they disappeared before the minimum dispersal age y min , which means they are neither potential dispersers nor candidates for a dispersing coalition breakup.
Therefore, the posterior of s u, i can be estimated using only the mortality likelihood as where the prior s p for the probability of being a female can be specified using empirical data such as species-specific sex ratios of the focal population  or information on litter sex ratio of the focal individual (Behr et al., in review). In the latter case, a different prior s p, i would be used for each individual instead of s p .

| Model fitting
To sample the conditional posteriors, a Markov Chain Monte Carlo (MCMC) algorithm is used (Figure 3). Specifically, each iteration of the MCMC algorithm consists of the following: for a set of parameters, e.g. the dispersal parameters γ = (α s = 0 , α s = 1 , κ s = 0 , ⋯), the algorithm proposes a new value for the first 'sub-parameter', α s = 0 , giving γ new .  all parameters θ, γ, ρ, s u , d u and b u . This is a Metropolis-within-Gibbs scheme (Clark, 2020;Gelfand & Smith, 1990), and it can be shown that the sequence of samples obtained in this way constitute a Markov chain, the stationary distribution of which is the joint posterior distribution (Gelman et al., 2013). The widths of the distributions from which proposals for new parameter values are drawn, are dynamically updated every 100 steps depending on the acceptance ratio. In case time-varying covariates are considered and assumed to have an age-dependent effect on mortality, the mortality parameters can be updated in blocks to achieve faster convergence of MCMC chains. The modelling algorithm was implemented in R, version 4.0.2 (R Core Team, 2020), and C + + (ISO/ IEC, 2020).

| Model validation
To validate our modelling framework, we used a simulation-based validation technique (Cook et al., 2006). Considering the hypothetical life history of the social species described in Section 2.2, we simulated life histories of 5000 individuals 72 times with a timestep of one day and a single covariate value generated at each timestep from a standard normal distribution. We also repeated the entire procedure without covariates. We chose a normal distribution for each parameter (e.g. In each simulation, 5000 individual life histories (i.e. time of death, dispersal, and dispersing coalition breakup) were generated in the following way. First, we sampled the age at which an individual dispersed from the probability density function (PDF) of the dispersal age using rejection sampling. This PDF is the product of the dispersal 'hazard' (Equation 9) and the probability that dispersal occurs before age y (Equation 10), g y | , y min 1 − G y | , y min .
Next, we sampled the age at which an individual died from the PDF of the age at death. This PDF is the product of the mortality hazard (Equations 6 and 7) and survival function (Equation 2), which both depend on the previously sampled age at dispersal (see Equation   3). If the drawn age at dispersal was smaller than the drawn age at death, we considered the individual to have dispersed, otherwise, we considered it to have died without dispersal. Lastly, if the individual had dispersed, we sampled the age at which a dispersing coalition breakup occurred in the same fashion as described above. If the age at breakup was smaller than the drawn age at death, we considered the individual to have separated from the dispersing coalition, otherwise we considered it to have died without a coalition breakup.
For each of the 72 simulations, we then applied our modelling framework with the corresponding life histories and the distributions from which we had drawn the simulation parameters as priors. We generated 48,000 samples from the posterior distribution of each parameter in each simulation. Next, for each simulation and for each parameter, we calculated the empirical posterior quantile of the true parameter, that is, the proportion of samples from the posterior distribution that are larger than the true parameter value. If our modelling framework was implemented correctly, the 72 quantile values that were obtained for each parameter should follow a uniform distribution (Cook et al., 2006). For each parameter, we checked this visually and by assessing goodness-of-fit using a Cramér-von Mises test.

| Model performance
We assessed the performance of our modelling framework using different degrees of uncertainty in the data, namely different sample sizes (large sample n = 5000 and small sample n = 1000), different resighting intervals (1 day and 30 days), different proportions of individuals that disappeared under unexplained circumstances (all fates known and partially unknown fates) resulting in a total of eight combinations. We chose a value for each parameter and ran one simulation (as described above) for every combination. To generate unknown fates, we selected 50% of the individuals that were simulated F I G U R E 3 Flow-chart representation of the procedure used to sample the conditional posteriors (see Section 2.3.3). A Markov chain Monte Carlo algorithm is used with parallel chains in which each chain is initiated with different randomly drawn starting values for the unknown parameters. Following the initialisation of all chains, the algorithm involves sampling conditional posteriors of the unknown parameters in successive steps. At each step, the algorithm updates the unknown parameters (analogous to Equation 18) successively. Each parameter is updated by looping over all individuals of the focal sample (illustrated by the blue boxes). to be potential dispersers (i.e. to have died after the minimum dispersal age). We considered their fate unknown, and if they had dispersed according to our simulation, we set their last observation to their actual time of dispersal. If they had not dispersed according to our simulation, we left their age at death as their last observation. Next, for those individuals known to have dispersed, we repeated the procedure to generate unknown fates with regard to dispersing coalition breakup. We used vaguely informed, normally distributed priors (σ = 10) that were centered on the chosen parameter values and fitted the eight simulated datasets using 16 parallel MCMC-chains with 8000 steps in each chain. MCMC convergence was visually assessed by inspecting trace plots.

| Sensitivity to censoring
To assess the sensitivity of mortality estimates to censoring of individuals with a known fate, we fitted our model to four datasets subject to substantial censoring. Specifically, we used the previously simulated dataset with 5000 individuals, 1 day resighting interval, and all fates known, and applied right-censoring to a random subset of 50% of individuals in each state. We did this four times, with varying proportions of individuals that were censored at the time of their death (0%, 25%, 50%, and 75%). All remaining right-censored individuals were censored at a random point in time (i.e. prior to death). To fit the model, we employed the same MCMC settings as described above.

| RE SULTS
For both sets of simulated data (one with and one without covariates), there was one out of fourteen parameters for which the Cramér-von Mises test gave a p-value of less than 0.05, indicating that the quantiles of the corresponding parameter did not follow a uniform distribution ( Figure S1). It should be noted, however, that the two 'problematic' parameters were not the same for the two sets of simulated data (a 0, R in case of no covariate and κ disp with one covariate). Furthermore, there were two issues potentially compromising the accuracy of our approach. First, due to computational restraints, we were only able to generate each of the (in total 2 × 72) quantile values using 16 chains of 3000 steps each, which led to some chains not being fully converged. This was especially the case for the resident mortality parameters, where we used a blocking scheme (see Section 4). Secondly, especially when using covariates, the number of data points per parameter becomes low for the dispersal and coalition breakup parameters, since only a fraction of all individuals undergo dispersal and only a fraction of those experience coalition breakup. Nonetheless, we checked and tested the corresponding sections of our code, but were unable to find any mistakes, leading us to conclude that our modelling framework performed satisfactorily.
For all eight simulated datasets, the mortality rates used to simulate the data lay within the 95% credible intervals of the estimated mortality posteriors for both resident and dispersing states, suggesting accurate model estimates (Figure 4). While the estimated posterior means were close to the true mortality means in almost all simulations, the variation in sample size appeared to affect the precision and accuracy of our modelling framework the most (Figure 4a-d versus Figure 4e-h). As expected, smaller sample sizes led to wider credible intervals, particularly for older ages in resident individuals and younger ages in dispersing individuals (Figure 4e-h). Neither higher proportions of unknown fates (Figure 4a,b,f,h), nor longer resighting intervals ( Figure 4c,d,g,h) alone significantly affected model performance.
Only in combination, these two variations resulted in lower precision (i.e. wider credible intervals; Figure 4d of missing individuals can be regarded as excellent (Hosmer & Lemeshow, 2000).
For datasets subject to substantial censoring, fitted models underestimated the mortality rates used to simulate the data. As expected, the magnitude of this underestimation increased as the proportion of individuals that were censored at the time of death increased ( Figure S2).

| DISCUSS ION
We presented a hierarchical model for estimating mortality of both resident and dispersing individuals as well as rates of state transition from mark-resighting data containing individuals that disap- Finally, although there is no universal guideline as to which of the simulated scenarios shown in Figure 4 is the most realistic (due to strong differences in the characteristics of each study), based on our application example (Behr et al., in review) with a sample of roughly 2000 individuals, resighting intervals in the range of 1 month, and a proportion of missing individuals of 16%, we can at least assume that for this example a simulation scenario located somewhere between the ones shown in Figure 4c,d,g,h is most likely to apply.
When using our modelling framework in studies with limited sample size (e.g. due to an imbalanced age distribution of marked individuals) and insufficient data (e.g. due to a lack of confirmed mortality and emigration events), we advise employing a data collection

F I G U R E 4
Comparison of estimated mortality rates (shaded areas) with mortality rates used to simulate data (solid lines) for individuals in a resident state (blue) or dispersing state (orange). Shaded areas depict 95% credible intervals (CI) and dashed lines indicate the mean of estimated posteriors. Results are shown for simulated data of different sample sizes (a-d vs. e-h), proportions of missing individuals (a,c,e,g vs. b,d,f,h), and resighting intervals (a,b,e,f vs. c,d,g,h). Note that for computational reasons, we deliberately chose parameters that made simulated mortality of dispersers smaller than mortality of residents (see main text for details). While a minimum dispersal age of 1 year was used in the simulation, the youngest dispersal age of all simulated life histories was 1.3 years; thus, the mortality rates shown for the dispersing state were truncated at that age. For datasets with unknown fates, the ability of the model to correctly infer the fates of missing individuals is indicated by the receiver operator characteristic (ROC) curve and the corresponding area under the curve (AUC), where true positive rates (TPR) are plotted against false positive rates (FPR). When AUC ≥ 0.8, the predictive performance of the model can be considered as excellent (Hosmer & Lemeshow, 2000).   .
Although studies with long resighting intervals could alternatively censor missing individuals at their last observation, our simulations advise against this and instead suggest inferring their fate to avoid an underestimation of mortality rates, which occurs even when individuals are right-censored at a random time ( Figure S2).
The widths of the credible intervals of the estimated mortality rates showed two noteworthy patterns: (1) for dispersing individuals, the parameters that model mortality in young individuals (a 0, D and a 1, D ) are estimated with less precision than the parameters describing mortality in old individuals, resulting in wide credible intervals for mortality rates at young dispersal ages ( Figure 4). This was due to the fact that most individuals do not disperse immediately after reaching the minimum dispersal age and the pool of dispersed individuals increases only slowly with age. Hence, the majority of dispersing individuals do not contribute to the estimation of the dispersal mortality parameters at young ages. This issue becomes even more apparent with small sample sizes (Figure 4e-h), thus reinforcing the conclusion drawn in the previous paragraph.
(2) For resident individuals, fitted mortality rates at old ages underestimated true mortality rates in simulations with long resighting intervals and unknown fates (Figure 4d,h) Colchero & Clark, 2012). If information on confirmed emigration events is limited, but at the same time information on the age of immigrants is available, immigration could be integrated in our model to parameterize emigration as described in SI Appendix S6. If individuals that are alive and present in the focal study area are not resighted with certainty, imperfect detection probabilities could be included (Kéry & Royle, 2016). If the state of an individual cannot be determined with certainty at any resighting occasion, state uncertainty could be explicitly modelled (Kendall, 2009). While our model includes two states (resident state and dispersing state), it could be extended by additional states (e.g. reproductive or behavioural states), making our framework suitable for investigating various hypotheses related to mortality. However, implementing such extensions would introduce additional unknown parameters into the model, potentially exacerbating some of the limitations noted above and requiring more computational resources.
Although our modelling framework is based on a robust statistical approach, it is not without limitations. First, the complexity of our model requires high computational power to achieve convergence, possibly restricting the applicability to other studies. The computational time for one MCMC step scales linearly with the number of individuals and increases far slower than linearly with the number of covariates. However, a larger number of covariates means that the time to convergence will be longer, and more or longer chains have to be run. For example, on a cluster of AMD EPYC 7702 processors, one of the 72 simulations to validate our modelling framework consumed approximately 640 CPU hours. Secondly, the effectiveness of fitting our model depends on the availability of information on known deaths and emigration events for at least a proportion of all marked individuals-a requirement that may not be met in datasets spanning only a few years. Finally, the reliability of model estimates for right-censored data rests on the assumption of independent and non-informative censoring, a limitation inherent in any mortality analysis (Kleinbaum & Klein, 2010). Thus, our approach (and likely any other mortality model) would underestimate mortality if censoring was more likely in the case of deaths.
In summary, we introduced a hierarchical model that offers a framework for estimating unbiased mortality and emigration rates for mark-resighting data with missing individuals while accounting for interval-censored observations, effects of time-varying covariates, and different behavioural states. The modularity of our framework allows the model to be tailored to a variety of speciesspecific life histories by adapting the provided R code. Following the workflow presented in Figure

CO N FLI C T O F I NTE R E S T S TATE M E NT
None.

PEER R E V I E W
The peer review history for this article is available at https://   (Cook et al., 2006). This null hypothesis was tested with a Cramér-von Mises goodness-of-fit test and we report here the resulting p-values.  the predictive performance of the model can be considered as good (Hosmer and Lemeshow, 2000).