Improving assessments of data- limited populations using life- history theory

1. Predicting how populations may respond to climate change and anthropogenic pressures requires detailed knowledge of demographic traits, such as survival and reproduction. However, the availability of these data varies greatly across space and taxa. Therefore, it is common practice to conduct population assessments by filling in missing values from surrogate species or other populations of the same species. Using these independent surrogate values concurrently with observed data neglects the life- history trade- offs that connect the different aspects of a population's demography. Consequently, this


| INTRODUC TI ON
Predicting how vulnerable populations are likely to respond to future threats requires detailed knowledge of life-history traits, including survival and reproduction (Pearson et al., 2014;Salguero-Gómez et al., 2016;Winemiller, 2005). However, our ability to monitor individual traits differs greatly across space and taxa, such that the availability of data is highly unbalanced (Horswill et al., 2019). This problem is so widespread that even the most optimistic data collection programme is not sufficient as a single solution to fill in all missing data. Therefore, identifying reliable ways to assess population status, even when demographic data are missing, has been highlighted as a key ecological problem facing policy makers (Kindsvater et al., 2018).
A typical approach for assessing populations with missing demographic data involves filling in gaps with values from surrogate species or populations of the same species (Githiru et al., 2007;Hernández-Camacho et al., 2015;Schtickzelle et al., 2005). However, combining independent surrogate estimates with observed demographic data neglects the physiological compromises that connect the different aspects of a population's demography (Stearns, 1992). Covariation among life-history traits is widely recognised across animal and plant groups (Bakewell et al., 2020;Healy et al., 2019;Jeschke & Kokko, 2009;Salguero-Gómez et al., 2016), and reflects the core relationships that connect fitness-related traits, such as body growth, reproduction and survival (Bennett & Owens, 2002;Lande, 1982;Stearns, 1992). At the extremes, species grow slowly, mature late, have low fecundity and high survival, or the reverse pattern. Similar patterns of covariation among life-history traits also exist across populations of the same species (Devenish-Nelson et al., 2013;Frederiksen et al., 2005;Horswill et al., 2019). Consequently, population assessments that fill in missing values with demographic parameters from independent surrogate sources will incorporate biases and imprecisions that may ultimately lead to population mismanagement (Devenish-Nelson et al., 2013;Hernández-Camacho et al., 2015;Johnson et al., 2010).
Predictive approaches that combine available data with established life-history trade-offs between demographic traits offer a new avenue for addressing taxonomic, spatial and temporal imbalances in demographic data (Horswill et al., 2019;Kindsvater et al., 2018;Thorson et al., 2017). This method considers multiple life-history traits, populations and species simultaneously in a hierarchical framework to reconstruct missing values based on the trade-offs that determine an organism's life-history strategy. The approach provides estimates for traits that are poorly monitored or difficult to measure empirically, and supports the estimation of parameter uncertainty based on measurement errors and data availability (Horswill et al., 2019). Consequently, we can use this technique to account for life-history trade-offs when filling in missing demographic data, as well as estimate parameter uncertainty that can be propagated through to population predictions to provide precautionary estimates of population change.
In this study, we use a predictive method based on life-history trade-offs (Horswill et al., 2019) to reconstruct population-specific demographic information across a substantial part of a species breeding range. The reconstruction framework used here can be applied to a wide range of species throughout the plant and animal kingdoms and we apply our analysis to an intensively monitored species of seabird, black-legged kittiwake Rissa tridactyla, where the spatial coverage of demographic monitoring is highly unbalanced between traits ( Figure 1). Black-legged kittiwakes are the most numerous gull in the world (Coulson, 2011); however, long-term population declines mean that this species is currently categorised as globally Vulnerable (Birdlife International, 2019). A lack of population-specific demographic information for this species has also recently been highlighted as a key contributing factor to uncertainty in the consenting process for offshore wind development in the United Kingdom (Ruffino et al., 2020). To address this, we conduct a projection analysis based on age-structured models to examine how reconstructed (i.e. correlated) demographic parameters may improve population assessments in data-limited situations, compared to assessments that fill in missing data with independent surrogate values.

| Data
We considered two demographic parameters for populations of black-legged kittiwakes breeding in the United Kingdom and Ireland: adult survival and fecundity (i.e. reproduction). The number of black-legged kittiwake populations in the United Kingdom and Ireland with empirical estimates for fecundity (n = 68) far exceeds the number of populations with empirical estimates for rates of adult survival (n = 5, Figure 1; Table S1). For the populations that have empirical estimates for both adult survival and fecundity (n = 5), we restricted the datasets to years when both demographic rates were monitored (Table S1). For the survival dataset, we limited our analysis to studies using mark-recapture modelling that had at least 8 years of data. The ability to estimate adult survival rates dramatically increases when mark-recapture time series are between 5 and 10 years in duration . For the wider fecundity dataset, we limited our analysis to populations and study plots in the Seabird Monitoring Programme Database (SMP, 2020) that had more than 50 breeding pairs and at least 5 years of data since 2009. Fecundity can be estimated with considerably shorter time series than are required to estimate rates of survival, and we selected this threshold for population size in order to minimise the influence of ecological processes that differentially affect small populations, such as demographic stochasticity  and depensatory regulation (Horswill, O'Brien, et al., 2017). The field protocols for collecting data at each site followed international guidelines for seabird monitoring (Walsh et al., 1995).

| Potential demographic impact of vital rates
We used elasticity analysis to examine the potential demographic impact of adult survival and fecundity to population growth. We applied this analysis to the five populations with observed data for both demographic rates (Figure 1). The projection matrices (A) followed a pre-breeding census (Caswell, 2001), and were parameterised using the mean observed rates of adult survival (Φ i ) and fecundity (F i ) for each population (Equation 1). The matrices were structured by age and the dimensions matched the mean age of first reproduction for black-legged kittiwake; 4 years (Horswill & Robinson, 2015). We calculated the elasticity of each demographic rate using the statistical package popbio (v. 2.7;Stubben & Milligan, 2007) in Program R (v. 4.0.2; R Core Team, 2020). This analysis estimates the change in population growth given a proportional perturbation in each vital rate (Caswell, 2001).
Like many seabirds, black-legged kittiwakes are largely unobservable during the first years of life. Therefore, population-specific estimates of juvenile (i.e. from fledging to age 1 year) and immature (i.e. before the age of first breeding) survival rates are limited (Horswill & Robinson, 2015). Available estimates for black-legged kittiwakes breeding in France report rates of juvenile survival (J i ) to be 75% of adults, rates of immature survival between age 1 and 2 years (I i ) to be 87.5% of adults, and birds to have survival rates comparable to adults from age 2 onwards (Cam et al., 2005). Therefore, we assumed an additive relationship between age-classes (Cam et al., 2005;Desprez et al., 2011;Horswill et al., 2014Horswill et al., , 2016, and applied age-specific scalars to estimate juvenile and immature survival rates from the imputed population-specific rates of adult survival. We assumed a 1:1 sex ratio in the population and halved fecundity to model female numbers only.

| Reconstructing demographic parameters for populations of black-legged kittiwakes
We constructed a Bayesian hierarchical model (hereafter 'cor- (1) The United Kingdom and Ireland has considerably fewer sites where populations of black-legged kittiwakes are monitored for rates of adult survival and fecundity (black triangles), compared to sites where only fecundity is monitored (red circles). Projection WGS84 The correlated model uses stochastic functions to quantify the life-history trade-off connecting rates of adult survival to fecundity and account for sampling variability in these traits.
Incorporating these sources of variance and covariance into the statistical imputation process means that the reconstructed (i.e. correlated) demographic parameters provide more biologically plausible replacements of missing values, compared to independent surrogate data from similar species or other populations of the same species (e.g. Horswill et al., 2019;Kindsvater et al., 2018).

The demographic component of the correlated model (Equation 2)
contained linear functions for adult survival (ϕ) and fecundity (f): Here, we use lower-case notation to denote reconstructed demographic parameters. Subscript i represents population, and is a multivariate normal (MVN) random effect that hierarchically connects rates of adult survival and fecundity. The covariance structure (Σ) of the random effect allows the relationship (i.e. lifehistory trade-off) between demographic parameters to emerge during model fitting: The mean values of the multivariate normal distribution ( ) represent the national mean values for each vital rate and were resolved on the scale of the linear predictor. Therefore, we assigned these parameters from normal prior distributions (N). We set the expectation of these prior distributions using estimates of black-legged kittiwake demographic rates in France (Cap Sizun); adult survival = logit(0.804), and fecundity = log(0.651; cited by Frederiksen et al., 2005). We also set the variance to allow imputation of the mean value within a realistic range for black-legged kittiwakes in the United Kingdom and Ireland (i.e. 0.03). The observed range of population-specific survival estimates was 0.81 to 0.86, and the 95% confidence interval of the prior distribution for the mean value was 0.69 to 0.89 ( Figure S1a). Likewise, the observed range of population-specific fecundity estimates was 0.24 to 1.25, and the 95% confidence interval of the prior distribution for the mean value was 0.35 to 1.16 ( Figure S1b; see Figure S2 for a comparison between the prior and posterior distributions of these parameters on the scale of the linear predictor).
We parameterised the variance-covariance matrix of the multivariate normal distribution (Σ, Equation 4) using an inverse precision matrix. We set the range of the variance priors ( , f ) to reflect predictions for long-lived species (such as black-legged kittiwakes) based on life-history theory, that is, low variance in rates of adult survival and higher variance in rates of fecundity.
The limited number of data points for adult survival required the shape of the prior distribution for variance in this trait to be more constrained than the equivalent prior distribution for fecundity.
Specifically, we used a truncated (T) normal distribution (N) for survival and uniform distribution (U) for fecundity (Equation 4).
The correlation parameter was drawn from a normal distribution, where the expectation was informed by the life-history trade-off between adult survival and fecundity reported for black-legged kittiwakes by Frederiksen et al. (2005); −0.62. We kept the variance of this prior distribution relatively broad (i.e. 0.05) because we wanted the data to determine the strength of the trade-off (Equation 4, see Figure S3 for a comparison between the prior and posterior distributions of this parameter). We truncated the correlation parameter at −0.9 and 0.9 to prevent unrealistically strong correlations.
We linked the empirical data for survival and fecundity with the reconstructed (i.e. correlated) demographic parameters through an observation component that reflected sampling variability. Survival values were drawn from a beta distribution and fecundity values from a truncated gamma distribution. We selected these distributions to bound inference within a plausible range on the observed scale, 0-1 for survival, and 0-3 for fecundity. We truncated the gamma distribution at three to match the maximum clutch size for black-legged kittiwakes.
Here, upper-case notation denotes empirical data for survival (Φ i ) and fecundity (F i ) on the observed scale. We defined the shape parameters of these distributions ( i , i , r i , i ) so that the population values reconstructed in Equation 2 ( i , f i ) provided the expectation, and published measures of sampling variability for black-legged kittiwakes determined the variance; 6 × 10 −4 for survival and 3 × 10 −3 for fecundity (Frederiksen et al., 2004; Supporting Information Appendix 1). This stochastic node acts on data points with and without observed data to provide populationspecific values with estimates of uncertainty for all survival and fecundity parameters.
To assess the influence of using informative prior distributions for certain parameters in our correlated model, we replaced these with uniform prior distributions that had biologically or theoretically driven bounds (Hobbs & Hooten, 2015). We then re-ran our analysis and compared the output to the model with informative prior distributions (Supporting Information Appendix 3).
Using a bounded uniform prior distribution to assign the expectation for survival, fecundity and the correlation parameter in the multivariate normal distribution, did not prevent model convergence but decreased the minimum effective sample size, and increased the posterior credible intervals of the reconstructed demographic rates (Supporting Information Appendix 3). By contrast, using a bounded uniform prior distribution to assign the variance parameter for survival in the variance-covariance matrix (i.e. replicating the prior distribution used to assign the variance parameter for fecundity, Equation 4) prevented model convergence (Supporting Information Appendix 3).
We implemented all models using JAGS (v. 4.3.0) via the 'jag-sUI' library (v 1.5.1; Kellner, 2019) for program R (v. 4.0.2). We fitted the models by running three Monte Carlo Markov chains (MCMC) for 2 × 10 5 iterations and retaining every 100th step to minimise autocorrelation in the MCMC sampling. We initialised each chain at different points in the parameter space. To confirm convergence of the chains we used the Brooks-Gelman-Rubin diagnostic tool (all values � r < 1.01) and the effective sample size of each parameter (all values >950). Convergence occurred before 5,000 MCMC draws, and therefore we removed the first 5,000 draws as burn-in. We provide a guide on how to assess model convergence in the supplementary material (Supporting Information Appendix 2).

| Population assessments using independent and correlated demographic parameters
We conducted a projection analysis to examine the potential biases and imprecisions generated by assessing populations with independent (i.e. combining observed population-specific data with independent surrogate values) and correlated demographic parameters (i.e. values reconstructed using the correlated model). In long-lived species, such as seabirds, the expectation is that rates of adult survival will experience low levels of variability, such that mean values should be relatively constant between populations. Therefore, it is common practice to use the national mean in place of missing adult survival data to conduct impact assessments on data-limited seabird populations (e.g. . The life-history tradeoff (i.e. covariance) between rates of survival and fecundity means that the difference between the national mean rate of adult survival and the true population-specific value will be greatest for populations with the highest and lowest rates of fecundity. Consequently, we restricted the projection analysis to the 10 populations with the highest and lowest mean values for fecundity, thereby examining the worst-case scenario for introducing bias into population assessments.
We conducted the projection analysis under four scenarios (Table 1). We ran each projection scenario for 25 years and 1,000 iterations. All projection scenarios also included demographic stochasticity using binomial and Poisson distributions for survival and fecundity events respectively. In the first and second projection scenarios, we held the mean demographic rates constant across all iterations. For the first scenario, we projected populations using independent demographic parameters. To do this, we com- To quantitatively compare the four projection scenarios, we examined the distribution of absolute (i.e. final) population sizes generated and measured vulnerability as the percentage of projections that resulted in a population size of less than 50 breeding females. We used this threshold because depensatory regulation, associated with increased predation at low density, differentially affects individual populations of black-legged kittiwakes (Horswill, O'Brien, et al., 2017), and it is not possible to identify which populations experience this regulatory pressure. We used pairwise comparisons between the outputs of our four projection scenarios to examine how independent and correlated demographic parameters and their associated TA B L E 1 The four projection scenarios (PS) used to examine the potential biases and imprecisions generated by assessing populations with independent and correlated demographic parameters uncertainty influence population assessments. For example, by comparing the outcomes of projection scenarios 1 and 2, we examined how independent demographic parameters influence population assessments, compared to correlated values (Table 1).
By comparing the outcomes of projection scenarios 1 and 3, as well as scenarios 2 and 4, we examined how incorporating parameter uncertainty can influence population assessments, and by comparing the outcomes of projection scenarios 2 and 4, we examined the effect of incorporating uncertainty that also accounts for observed levels of covariance between demographic parameters (Table 1).
For each projection scenario (p) we used an age-structured model to predict the size (S) of each population (i) at time t + 1 as a function of population size at time t and constant demographic parameters (Equation 6). We structured the projection matrices for these deterministic models following the same principles that we used to construct the projection matrices in the elasticity analysis

| Potential demographic impact of vital rates
As expected for long-lived species, we found that potential demographic impact (measured as elasticity) was consistently higher in rates of adult survival (mean = 0.61, SD = 0.04), compared to rates of survival from hatching to age 1 year (mean = 0.10, SD = 0.01).
The much greater potential demographic impact of adult survival on population growth highlights the importance of using accurate estimates of this demographic parameter in population assessments.

| Reconstructing demographic parameters for populations of black-legged kittiwakes
We detail the reconstructed (i.e. correlated) rates of survival and fecundity for 68 populations of black-legged kittiwakes in Table S1. This includes 63 populations that are unmonitored for

| Population assessments using independent and correlated demographic parameters
Using independent demographic parameters resulted in faster rates of decline in the low fecundity group, and slower rates of decline in the high fecundity group, compared to assessments using correlated demographic parameters (Table 2;

| D ISCUSS I ON
Life-history trade-offs constrain the range of possible life-history strategies that can evolve across the tree of life (Lande, 1982;Stearns, 1992). Accordingly, the correlated model used in this study can be applied across the plant and animal kingdoms to provide potentially more precise estimates of life-history traits with missing data, compared to values selected from independent surrogate species or populations. These reconstructed (correlated) values can be used to parameterise any analysis of population dynamics based on demographic rates, including population viability analysis, integrated population models, integral projection models and partial differential and ordinary differential equations.  Table 1) that resulted in vulnerability or populations increasing from the initial size F I G U R E 4 (a) The absolute (final) population sizes obtained from the projection analyses strongly differed for populations with low and high rates of fecundity when demographic parameters (adult survival and fecundity) were from independent sources (projection scenario 1, Table 1). (b) This difference was smaller when assessments were based on reconstructed values from our correlated model (projection scenario 2, Table 1). (c) Propagating parameter uncertainty that was independent between rates of adult survival and fecundity (projection scenario 3, Table 1) increased the range of possible outcomes and increased predicted vulnerability in the low fecundity group, compared to projection scenario 1. (d) Propagating parameter uncertainty that maintained observed levels of covariance between traits increased the range of possible outcomes and increased predicted vulnerability in both the low and high fecundity groups, compared to projection scenario 2 (projection scenario 4, West Scotland (Figure 1). Using elasticity analysis, we show that rates of adult survival in this species have a considerably larger demographic impact than rates of fecundity and juvenile survival.
This result agrees with life-history theory and highlights the importance of having robust estimates of adult survival rates when conducting impact assessments on long-lived species, such as seabirds. It may be less critical to use correlated values to fill in missing demographic parameters with low elasticity if assessments are deterministic. However, environmental stochasticity can increase the demographic impact of vital rates with low elasticity (Gaillard et al., 1998), such that reliable estimates of mean values will become more important if assessment models include demographic rates that vary temporally.
In data-limited circumstances, it is common practise to adopt a pre-cautionary approach to population assessment, and favour the method that provides the maximum estimate of potential mortality (Horswill, O'Brien, et al., 2017). We show that filling in missing data with independent surrogate values does not offer a consistently precautionary approach. Using an independent value that underestimates the missing demographic parameter will overestimate rates of population decline, whilst using an independent value that overestimates the missing demographic parameter could seriously underestimate rates of population decline. At present, it is not possible for practitioners to ascertain whether surrogate values taken from independent sources will over or underestimate missing values. Consequently, correlated vital rates reconstructed using life-history trade-offs can also be used as reference estimates to identify potential biases associated with independent surrogate values.
Population-specific data on other life-history traits, such as age of recruitment and body growth rates, are limited in seabirds.
Therefore, we conducted our analysis using two vital rates, adult survival and fecundity, and a one-dimensional parameter space to represent the life-history trade-off (i.e. covariance between traits). As a result, populations with rates of fecundity close to the national average will have reconstructed rates of adult survival that are also close to the national average. In accordance, national averages could be used to fill in missing data for populations that have life-history strategies very close to the average. However, this may change if the modelling framework is extended to incorporate more than two traits, such that inference becomes multidimensional. Multivariate extensions of the correlated model can be easily implemented and are likely to improve the reconstruction of highly fragmented datasets (Horswill et al., 2019). Future reviews and syntheses of how different life-history traits relate to each other by species and groups are likely to be a valuable resource for specifying these models.
In Bayesian statistics, prior distributions can be used to summarise our understanding of how the world works to obtain meaningful inference from small and fragmented datasets (Hobbs & Hooten, 2015). In our correlated model, we used scientific knowledge on the biology of the study species and from lifehistory theory to inform the shape of the prior distributions assigning certain parameters. In the absence of reference studies, it is still possible to incorporate existing scientific knowledge when defining prior distributions by using biologically or theoretically driven bounds. For example, survival can be assigned from a uniform distribution bound by 0 and 1, and the correlation parameter highly limited in our original dataset (n = 5, Figure 1) and evidently lacked the information required to quantify variance in this trait.
Consequently, correlated models incorporating highly fragmented data are likely to benefit from formulating expectations based on existing data and expertise, and using these to assign informative prior distributions on relevant parameters. We recommend defining these prior distributions to reflect predictions from life-history theory. For example, among populations of a fast-lived species, fecundity is expected to exhibit lower variation than adult survival (Lande, 1982;Stearns, 1989).
In large populations, environmental stochasticity is usually a dominant source of demographic variation. However, in populations of <100 individuals, demographic variation is mostly driven by demographic stochasticity . Instead of explicitly modelling demographic stochasticity, deterministic thresholds are often used to predict extinction (Bascompte, 2003). Extinction thresholds, that is, critical densities beyond which a population can no longer persist, are difficult to define because they will vary by taxon, mating system and local environmental pressures Lande et al., 2003;Saether et al., 2013). Consequently, we incorporated demographic stochasticity into our projection analysis. This provides a more mechanistic approach for modelling processes at low density, thus reducing uncertainty in how small or declining populations will respond to environmental pressures.
Other rule-based methods of impact assessment, such as potential biological removal (PBR), have lower data requirements than matrix (e.g. age-structured) models, making them quite appealing in data-limited situations (Dillingham & Fletcher, 2008;Wade, 1998). However, the algorithms underpinning rule-based methods often rely on implicit assumptions that limit their application (Niel & Lebreton, 2005). Consequently, rule-based methods are considered unsuitable for environmental impact assessments, especially on seabirds Miller et al., 2019;. Alternative approaches for overcoming data limitations include using phylogenetic and spatial distance to classify species with similar traits and conservation status (James et al., In press;Jetz & Freckleton, 2015;Penone et al., 2014). However, phylogeny may not be a strong predictor of all species-level demographic parameters (Che-Castaldo et al., 2018). It also becomes less informative when considering population, as opposed to species-level strategies, and when examining closely related species within the same genus. Finally, populations in close proximity may not share similar vital rates if environmental pressures differ, such as density dependence or predation (Table S1).
Including density-dependent regulation in population assessment models remains debated Horswill, O'Brien, et al., 2017;. In our study, including compensatory density dependence would have slowed the rate of population decline, but the relative differences between simulations based on independent and correlated demographic parameters would not have changed. We also used baseline (i.e. mean) demographic parameters in our projection analysis. This approach is common for conducting population and environmental impact assessments (Hernández-Camacho et al., 2015;Johnson et al., 2010;Peck et al., 2008). However, temporal variance and covariance between demographic rates can contribute substantially to population growth (Coulson et al., 2005;Doak et al., 2005) (Horswill et al., 2019). If data are available, separate covariance structures may also be required to account for age-specific differences in carryover effects associated with adverse environmental conditions (Horswill et al., 2014;Horswill, Trathan, et al., 2017).
In addition to supporting internationally important breeding populations of seabirds (Mitchell et al., 2004), North-West Europe also provides the best wind resources in the region (Wind Europe, 2020). Extensive development of offshore wind farms is planned to take place within European waters over the next decade (Corbetta et al., 2015;OGL, 2019), and mortality due to turbine collision has been identified as a serious potential threat to many species of seabird, including black-legged kittiwakes (Bradbury et al., 2014;Dias et al., 2019;Furness et al., 2013).
International legislation requires that most European countries, granting or refusing consent for the construction of wind farms, conduct a full assessment of the possible impacts to local and However, at present, large amounts of missing demographic data are contributing to uncertainty in these impact assessments (Ruffino et al., 2020). While it may be possible to estimate a population mean for fecundity using relatively short time series of data, collecting enough mark-recapture data to robustly estimate rates of adult survival may not be attainable within the consenting process of a wind farm . Therefore, the correlated demographic parameters we detail for black-legged kittiwakes have immediate application for improving environmental management decisions relating to the construction of wind farms in UK and Irish waters. We recommend that impact assessments on data-limited populations evaluate potential future risk by comparing predictions based on independent and correlated demographic parameters, as well as models that propagate correlated parameter uncertainty through to the population predictions. The reconstruction technique applied in this study can easily be adapted to include additional populations. It can also easily be refocused to another species. It is important to highlight that applications to more data-limited species may require informed prior distributions (Appendix 3), and we recommend using the posterior distribution of the correlation parameter reported in this study to set the corresponding prior distribution in correlation models of seabird species with fewer estimates of adult survival. Being able to benefit from several sources of information is a key strength of the reconstruction approach applied in this study. This also highlights the importance of maintaining an extensive network of long-term monitoring datasets to support future applications.

DATA AVA I L A B I L I T Y S TAT E M E N T
Input data and the script for the correlated model available via the Dryad Digital Repository https://doi.org/10.5061/dryad.qnk98 sfg0 . All data and sources are also listed in the Supporting Information Table S1.