Joint estimation of growth and survival from mark‐recapture data to improve estimates of senescence in wild populations

Understanding age-dependent patterns of survival is fundamental to predicting population dynamics, understanding selective pressures, and estimating rates of senescence. However, quantifying age-specific survival in wild populations poses significant logistical and statistical challenges. Recent work has helped to alleviate these constraints by demonstrating that age-specific survival can be estimated using mark-recapture data even when age is unknown for all or some individuals. However, previous approaches do not incorporate auxiliary information that can improve age estimates of individuals. We introduce a survival estimator that combines a von Bertalanffy growth model, age-specific hazard functions, and a Cormack-Jolly-Seber mark-recapture model into a single hierarchical framework. This approach allows us to obtain information about age and its uncertainty based on size and growth for individuals of unknown age when estimating age-specific survival. Using both simulated and real-world data for two painted turtle (Chrysemys picta) populations, we demonstrate that this additional information substantially reduces the bias of age-specific hazard rates, which allows for the testing of hypotheses related to aging. Estimating patterns of senescence is just one practical application of jointly estimating survival and growth; other applications include obtaining better estimates of the timing of recruitment and improved understanding of life-history trade-offs between growth and survival.


Accepted Article
This article is protected by copyright. All rights reserved. conflict) as they age (Medawar, 1952;Williams, Day, Fletcher, & Rowe, 2006). Moreover, this susceptibility itself can also change with age (Williams & Day, 2003). The ubiquity of demographic senescence (i.e., changes in age-specific survival and fecundity due to agespecific deterioration) motivates the need for statistical models that estimate age-specific demography (Siler 1979, Pletcher 1999, Colchero and Clark 2012. We focus here on characterizing age-specific mortality. The ability to measure the relationship between age and mortality is essential for biodemographic modeling (Miller et al. 2014), informs evolutionary theory (Ricklefs 2010), and aids conservation and population management efforts (Lynch & Fagan 2009).
Laboratory experiments and theoretical work on senescence and aging often rely on assumptions and constraints that are unlikely to be true for free-living wild populations (Williams et al. 2006, Nussey et al. 2008. Thus, it is important to measure and understand patterns of senescence in the wild where natural selection occurs. Mark-recapture studies provide data on animals in the wild and these types of longitudinal data sets are necessary for understanding inter-and intra-specific variation in senescence (Nussey et al. 2008). For example, despite it long being thought that turtles senesce slowly or not at all, Warner et al.
(2016) detected significant declines in age-specific reproductive success and survival by measuring reproductive output, embryo survival rates, and post-hatching mortality in painted turtles (Chrysemys picta) using a 20-year mark-recapture data set.
Mark-recapture methods are a common approach for measuring survival in wild populations. However, estimating age-specific patterns of survival has traditionally required following known-age individuals for long periods of time to track their mortality. This process requires knowledge of the age of individuals, which is often only available for

Accepted Article
This article is protected by copyright. All rights reserved.
We demonstrate that estimating growth and survival simultaneously improves estimates of age-specific survival. We first describe a hierarchical model that incorporates growth, age-specific hazards, and mark-recapture survival estimates and use simulated data to validate the model. Then, using mark-recapture data sets for males and females in two populations of painted turtles, Chrysemys picta, we show how the model can reduce bias in the estimates of the rate of senescence. In doing so, we shed light on an ongoing debate on the generalizability of female advantage in lifespan and survival (e.g., Austad and Fischer 2016). By focusing our case study on this important model taxon, we broaden the comparative landscape of wild aging biodemography and life-history studies.

Model:
We develop a general model for estimating age-specific survival that combines three distinct components: a model to estimate age and associated uncertainty using a mark-recapture version of a von Bertalanffy growth curve model (Fabens 1965, Schofield et al. 2013, an age-specific hazard function to characterize the functional relationship between age and mortality that estimates patterns of demographic senescence (Siler 1979, Pletcher 1999, Colchero and Clark 2012, and a mark-recapture survival model to estimate hazard function parameters dependent on known or estimated ages (Cormack 1964, Seber 1965, Lebreton et al. 1992. We combine the three components into a single hierarchical model to share information and distribute uncertainty among the components. The Colchero and Clark (2012) model estimates age-specific survival when age is unknown for a subset of the population. Our model incorporates additional information available from sizes and growth rates for the unknown-age component of the population. Our goal is to improve estimates by reducing uncertainty in the age-based component of the model.

Accepted Article
This article is protected by copyright. All rights reserved.
Growth model: We use a von Bertalanffy growth model to predict size at capture Li t for individual i during survey occasion t. We use an age-at-first-capture formulation of the mark-recapture model (Schofield et al. 2013). We denote the initial size at age 0 for individual i as L0i . Similarly, asymptotic size of an individual is denoted as LIi and represents the maximum size attained by an individual at age infinity. Finally, each individual is assumed to grow at a growth rate Ki, which represents the proportion of growth from L0i to LIi that remains after a year of time. The size of individual i at time t is calculated as: where AFCi is the age of the individual at first capture (may be known or unknown) and Δit is the number of years since first capture. At the first capture, Δi t = 0 and (AFCi + Δit) is the age of the individual i at time t.
We allow for individual variation in each of the estimated parameters. For L0 and LI, we assume normal variation where L0i ~ Normal(µL0, ) and LIi ~ Normal(µLI, ). For K, we assume logit normal variation where logit(Ki) ~ Normal(µK, ). For AFC, we assume a negative binomial distribution, which constrains AFC to be a positive integer. It is also possible to incorporate AFC as a continuous variable, using a distribution such as a lognormal (see Schofield et al. 2013). We chose to model AFC as discrete since most animals have limited birthing seasons that are amenable to assignment to a year. We also let each parameter vary by sex, though it is additionally possible to vary them by time in the case of K, or by other useful covariates that explain within-population growth patterns (e.g., site).
The final component of the growth model relates Lit , the expected size of individual i at time t, to the observed size value, obsLit . This accounts for sampling error as well as lackof-fit to the vB function. We assume a normal error structure where obsLit ~ Normal (Lit,, ).

Accepted Article
This article is protected by copyright. All rights reserved.
When fitting models in subsequent sections, we use uniform priors for µL0 and µLI across reasonable values for each. For µK, we use a uniform prior on the probability scale where logit(µK) ~ Uniform(0,1). For the variance parameters , , , and , we use a t-distributed prior where σ ~ t(0,0.0004,3) but is constrained to be greater than 0, as in Schofield et al. 2013. Using the t-distribution allows for robust estimates and occasional extreme values (Gelman and Hill 2007). Note that Eaton and Link (2011) provide another formulation to incorporate individual variation into a vB growth function, which has been successfully applied to estimate growth for other mark-recapture data sets (Fellers et al. 2013, Rose et al. 2018a).

Age-specific hazard function:
In studies that explore rates of demographic senescence, it is typical to estimate the relationship between age and survival using hazard functions (Bronikowski et al. 2011, Miller et al. 2014. Unlike most estimates of survival using markrecapture data, mortality is treated as a continuous rather than as a discrete process, such that the hazard is an estimate of the instantaneous ability to an animal to survive at a given time, rather than over a determined time period. In this case, if the mortality rate μ is constant, then survival for one year is given by e -μ and survival from birth to age T is given by e -μT .
In the case of senescence, mortality is expected to depend on the age, a, of the individual, and is denoted by m(μ|a). This relationship may be described by Gompertz, Weibull, Siler, and Logistic equations, with or without a constant Makeham mortality term, all of which describe how mortality changes as an animal ages (Siler 1979, Pletcher 1999, Colchero and Clark 2012. If the continuous mortality rate varies as a function of age, then estimating the survival over a given interval requires determining the average mortality during the interval multiplied by the length of that interval. To do so, one calculates the

Accepted Article
This article is protected by copyright. All rights reserved.
integral of m(μ|a) over the interval of interest and as before, takes the negative exponent of the natural logarithm for this value. For example, the probability of surviving, φ, from age a to a + 1 is: The form used for m(μ|a) is flexible. We focus on exponential (i.e., constant) mortality using a Gompertz function (Gompertz 1825), where:

= *
In most cases, we expect mortality to be high for young individuals and then to fit a Gompertz-type function once animals reproductively mature. In this case, we can separate the early ages and use an exponential mortality function for these younger animals and then use a Gompertz function for older animals; we subsequently call the age at which this separation is made the truncation age. Importantly, truncation age is not necessarily equivalent to the age at maturation, but rather is the age at which a population transitions from exponential to Gompertz mortality. Other functions can be used, as described above, and we include code

Accepted Article
This article is protected by copyright. All rights reserved.
The CJS formulation conditions on the initial capture, where z i0 = 1, which denotes that the individual is alive (i.e., z = 1). The probability zit +1 = 1 (i.e., the animal is still alive in time t + 1 given that it is alive in time t) for all subsequent years is given by zit +1 = Bernoulli(zit * φi[ai t]).
The value for φi[ai t] comes from the age-specific hazard function. The probability of observing individual i at time t, yit, is the probability of detecting an animal given that it is alive at time t. We express this probability as yit ~ Bernoulli(zit * θt). We use θt to denote the probability of detecting an individual, given that it is alive at time t, and we allowed θt to vary as a function of time, age, sex, or other factors thought to affect capture probabilities in the system.

Simulated data set:
We use Monte Carlo simulations to test model performance and to assess whether including information about individual growth increases the accuracy of our approach. We simulate the capture of 25 new individuals (i.e., not previously captured) each year for 20 years where animals still alive in subsequent years were recaptured with a moderate 50% detection probability. Age at first capture (AFC) is assumed to come from a negative binomial distribution. Length at birth and asymptotic size for each individual come from normal distributions with means of 30 and 150 and standard deviations of 1 and 10, respectively. Additionally, we ran simulations with varied numbers of animals and a lower detection probability and found that the results are decently robust (Appendix S1).
We use a truncated mortality function for early ages where survival is allowed to vary until age 6, and thereafter follow a Gompertz model where ß 0 = -4, and ß 1 = 0.2. We generate 100 data sets and then analyze each data set with and without growth information included in the estimator to determine whether our growth model improved estimates. In addition, we varied the proportion of known ages to estimate the impact of unknown age percentages on

Accepted Article
This article is protected by copyright. All rights reserved.
universally age faster than females. Specifically, we found that female and male painted turtles demonstrate similar adult mortality acceleration (Appendix S2: Figure S3), although the amount of uncertainty merits further investigation. Mammals appear to have a general phenomenon of enhanced female life expectancy relative to males, often accompanied by reduced mortality acceleration, whereas birds have a reported male advantage (Bronikowski et al. 2011, reviewed in Austad andFischer 2016). Interestingly, one of the hypotheses for sex-specific variation in lifespan centers on an advantage to the homogametic sex, with some available data supporting this contention (Tower 2006, Maklakov andLummaa 2013).
Painted turtles, like many species of reptiles, do not have sex chromosomes and sex is determined by the temperature experienced during egg incubation (Janzen and Krenz 2004).
Furthermore, across the reptile clade, both forms of genotypic sex determination (i.e., heterogamety and homogamety) occur in addition to environmental sex determination (Janzen and Phillips 2006). Thus, studying sex-specific aging across non-avian reptiles has the potential to disentangle causes of aging from the presence of sex chromosomes.
We find that incorporating growth increases the accuracy of models for estimating senescence in populations of animals of unknown ages. The method is likely to be most effective in species where size is indicative of age; species that reach asymptotic size early in life may not be appropriately fit with the vB function. The general approach used here--simultaneously estimating growth and survival---has many additional applications not highlighted here. The approach allows for improved estimates of age and uncertainty about age, life-span of individuals, average life-span of populations, and size-specific survival rates.
In addition, the model could be extended to improve estimates of recruitment as well as covariation between growth and survival. Improved methods for estimating these parameters are necessary not only for understanding the evolutionary ecology of long-lived species, but also for aiding their management and conservation.

Accepted Article
This article is protected by copyright. All rights reserved.

Accepted Article
This article is protected by copyright. All rights reserved.