1 Introduction

Survival and recruitment of juveniles during their first stage of life can drive fish and wildlife population trajectories (Altwegg et al. 2005; Benson et al. 2018; Eberhardt and Siniff 1977; Gaillard et al. 1998; Hatter and Janz 1994; Middleton et al. 2013; Pojar and Bowden 2004; Unsworth et al. 1999). Methods for estimating juvenile survival often rely on intense data collection efforts to capture and uniquely mark individuals during an early stage of life and tracking them through time (Brownie and Robson 1983; Pradel 1996; White 1983). These data can provide high-quality information on juvenile survival. However, marking juveniles can affect survival or behavior (Côté et al. 1998; White et al. 1972), and capturing juveniles in a timeframe sufficient to estimate early-stage survival can be challenging due to narrow and stochastic windows of opportunity.

Juvenile survival of many species relies on postnatal parental care (e.g., lactating mammals). The window for capturing parents is wider than that of juveniles, where adult animals can be targeted at more convenient spatial, temporal, or behavioral domains. When a marked adult provides care for juvenile animals, investigators can use the adult mark to locate and count unmarked juveniles (Hastings et al. 2009). Technological devices (e.g., vaginal implant transmitters) can also aid in identifying birth time and location for several species (Bishop et al. 2007; Bowman and Jacobson 1998; Monteith et al. 2014).

Due to the importance of reliable estimates of juvenile survival for understanding population dynamics, and the challenges associated with marking juveniles, investigators have been developing innovative ways to obtain survival estimates from unmarked juveniles and marked parents for at least 50 years (e.g., Mayfield 1961). Current methods assume knowing the initial number of juveniles born (Lukacs et al. 2004), assume perfect detection of marked adults (Schmidt et al. 2015), perfect detection of unmarked juveniles (Flint et al. 1995; Manly and Schmutz 2001), do not distinguish between adult detection probability and juvenile detection probability (Hastings et al. 2009), or assume equal survival across groups of a given age (Lukacs et al. 2004). Innovation that helps relax these assumptions will offer investigators more flexibility and provide more realistic models of juvenile survival, especially if we can leverage this innovation with advances occurring at the forefront of aerial wildlife survey technology, including unmanned aerial systems, and advances in photographic methodology.

Aerial surveys are an important tool for estimating abundance and distribution of vertebrate populations. They permit observations of wildlife over landscape-scale areas that are often remote or inaccessible to ground-based observation and can be used in rugged or dangerous terrain (e.g., areas affected by recent wildfires). Available technology now permits investigators to study animals using high-resolution video and thermal imagery systems (Cilulko et al. 2013; Conn et al. 2014; Linchant et al. 2015; Spaan et al. 2019), permitting increased detection of inconspicuous animals or animals with elusive life history traits (Fig. 1). Platforms for using high-resolution video and thermal imagery systems include fixed wing aircrafts, helicopters, and unmanned aerial systems (UAS; Chabot and Bird 2015; Hodgson et al. 2016; Ritchie et al. 2012; Spaan et al. 2019). While aerial surveys have been used extensively for estimating abundance and distribution (e.g., Williams et al. 2017), and to a lesser extent, estimating survival of uniquely marked individuals (e.g., Bradford et al. 2006; Langtimm et al. 1998), lack of an appropriate statistical framework has often precluded using aerial surveys for estimating survival of unmarked animals or cases where parents are marked and juveniles are unmarked. Schmidt et al. (2015) developed novel tools to combine known-fate data of marked individuals with counts of unmarked individuals in an open-population N-mixture modeling framework (sensu Dail and Madsen 2011) and demonstrated its application with data consisting of individually marked wolves and counts of unmarked wolves in a pack with the marked individual. Here, we extend this framework by applying it to obtain juvenile survival, information that is critically important to wildlife population dynamics and often difficult or expensive to obtain. We also extend the methods of Schmidt et al. (2015) to permit both juvenile and adult detection probability, an important consideration when designing wildlife field studies.

Fig. 1
figure 1

Aerial photographs (top) and infrared images (bottom) of a mule deer taken during our field study in a remote portion of central Nevada. The images on the left-hand side contain a radio-marked adult doe, and two unmarked juvenile fawn mule deer that were obtained by homing in on the radio-marked adult

Our objective was to develop a cohesive framework to estimate reproduction rates, juvenile survival and detection probability, and adult detection probability using repeated observations of radio-marked adults with known fates and unmarked juveniles with unknown fates. We first developed a statistical framework to estimate our parameters of interest. Second, we evaluated the consistency of our statistical methods using simulated data with known parameter values. Third, we conducted a field study to examine the feasibility of using a fixed-wing aircraft mounted camera system to monitor marked adult female mule deer (Odocoileus hemionus), and their juvenile offspring, in a remote section of central Nevada, and fit our statistical model to the data. While we were interested in developing methods that support using piloted or unmanned aircrafts for efficient data collection in hard-to-access environments, the methods we developed can be generally applied to many systems and taxa, including camera traps or other remote image devices.

2 Methods

2.1 Model Development

Following the terminology of Berliner (1996), we present a Bayesian hierarchical model consisting of three levels (Appendix 1 of ESM). At the top level, we developed a data model linking the observed data and associated variation to latent ecological processes. Next we discuss a process model describing the underlying ecological processes (i.e., birth and death rates). Finally, we developed parameter models representing our prior knowledge about the parameter inputs in the ecological process model and data model. We use the square-bracket notation [a|b] to represent the probability density or mass function of variable a given variable b (Gelfand and Smith 1990).

2.1.1 Data Model

Investigators commonly radio-mark individuals and subsequently attempt to relocate them by homing in on their signal (e.g., Williams et al. 2011, 2014). Relocating marked individuals typically occurs to obtain visual confirmation that the radio is still functioning, to verify the individual is still alive, or to collect other relevant information (e.g., breeding state; Nichols et al. 2007). During relocation, observations of marked adults and their unmarked juvenile offspring fall into one of three categories: (1) positive identification of the adult and an accurate count of all living juveniles, (2) positive identification of the adult and an inaccurate count of all juveniles (typically consisting of false negative errors because some juveniles are missed), or (3) failure to accurately identify the adult, and therefore, failure to obtain any information on juveniles. We require a data model that appropriately quantifies each potential outcome.

We let \(y_{i,t}\) represent binary data for the ith marked individual during the tth (\(t=1,\ldots ,T\)) sampling period in which \(y_{i,t}=1\) if the marked individual was positively identified, and \(y_{i,t}=0\) otherwise. After detecting a marked adult, investigators have the opportunity to count the number of juveniles accompanying the adult. We let \(c_{i,t}\) represent counts of juveniles. If an adult is not detected, \(c_{i,t}=0\) with probability one. If an adult is detected, the data \(c_{i,t}\) often under-represent true number of juveniles present \(N_{i,t}\), due to imperfect detection. Therefore, it is necessary to incorporate detection probability (\(p_{\text {juv}}\)) to estimate the true number of juveniles, and subsequently, juvenile survival. We represent the relationship among juvenile count data, true latent number of juveniles, and observed adults as

$$\begin{aligned} \begin{aligned} c_{i,t} \sim&{\left\{ \begin{array}{ll} [c_{i,t}|N_{i,t},p_{juv,i}], &{} \quad y_{i,t}=1 \\ 0, &{} \quad y_{i,t}=0 \\ \end{array}\right. }, \end{aligned} \end{aligned}$$
(1)

where \([c_{i,t}|N_{i,t},p_{juv,i}]\) is a probability mass function relating the counts \(c_{i,t}\) of juveniles to the latent, true number of juveniles \(N_{i,t}\), and potentially, a detection probability \(p_{juv,i}\). The binomial distribution is the most commonly used model in ecology to address false negative errors. However, if juvenile counts contain false positive errors (e.g., the potential for double counting individuals), other distributions (e.g., Poisson) may be considered (standard model checking guidelines apply; c.f., Conn et al. 2018).

It is not always possible to obtain a positive identification of marked animals. For example, vegetation or rugged terrain can obstruct views of animals, even when they are in close proximity (e.g., Hoffman et al. 2010). Further, marked individuals may not be identifiable from other adult con-specifics in a group, herd, raft, or pod. Therefore, their juvenile offspring cannot be uniquely identified. Assuming the marked individual is alive, we model the observation process of the marked adult using a Bernoulli distribution

$$\begin{aligned} y_{i,t}&\sim \text {Bernoulli}(p_{\text {adult}}), \end{aligned}$$
(2)

where \(p_{\text {adult}}\) is the probability of positive identification of marked parents. We assumed marked parents remained alive from periods \(t=1\) to \(t=T\) (which can be confirmed with GPS technology for many wildlife species). If a parent dies during the course of a study, it can either be removed from the data prior to analysis or Eq. 2 can be modified to address adult survival, given appropriate data exist.

2.1.2 Process Model

To model survival of juveniles from time \(t-1\) to t, we used a binomial distribution with survival probability \(\phi \) that may depend on a vector (denoted with a bold lower-case letter) of exogenous covariates x (hypothesized to influence survival), estimated parameters \(\varvec{\beta }\), and a link function f (e.g., the logit or probit link function)

$$\begin{aligned} \begin{aligned} N_{i,t}&\sim \text {binomial}(N_{i,t-1},\phi _i),\\ f(\phi _{i})&=\mathbf{x }_{i}'\varvec{\beta }. \end{aligned} \end{aligned}$$
(3)

An advantage of the process model in Eq. 3 is that if we do not know the number of juveniles born to each parent initially, we can estimate it probabilistically. For example,

$$\begin{aligned} N_{i,1}&\sim [N_{i,1}|\mu ,\nu ], \end{aligned}$$
(4)

where \(N_{i,1}\) is the number of juveniles born to the ith parent in the first time period, \(\mu \) represent the mean fecundity rate, and \(\nu \) describes the uncertainty of mean fecundity. Common ecological models for count data include the Poisson and negative binomial distributions. For our field study, described in Mule Deer Application below, and for modeling mammal reproduction in general, counts of juveniles are under-dispersed relative to these distributions. One model for under-dispersed data is the Conway–Maxwell Poisson distribution, with mean fecundity rate \(\mu \), and under-dispersion parameter \(\nu \)

$$\begin{aligned} N_{i,1}&\sim \text {Conway-Maxwell Poisson}(\mu ,\nu ) \end{aligned}$$

(Wu et al. 2013; Lynch et al. 2014). The Conway–Maxwell Poisson distribution permits realistic representation of ungulate reproduction, which we assessed in our mule deer application using prior and posterior predictive distributions (Appendix 1 of ESM). For other taxa (e.g., birds), the Poisson or negative binomial distribution may be appropriate.

2.1.3 Parameter Models

To complete the Bayesian specification of our model, we describe probability models for the parameters discussed in the data and process components, above. Parameters requiring prior distributions include \(p_{\text {juv}}, p_{\text {adult}}\), \(\varvec{\beta }\), \(\mu \), and \(\nu \). The support of the detection parameters should be restricted to real numbers between 0 and 1. In our simulations and field study, we used conjugate beta(q, r) distributions for detection probability. This results in a beta-binomial model for the observations \(c_{i,t}\) (Gelman et al. 2013). The beta-binomial model can be used to model heterogeneity in detection probabilities between groups or correlation within groups (Hisakado et al. 2006; Martin et al. 2011). The correlation parameter \(\rho \) is

$$\begin{aligned} \rho&=\frac{1}{q+r+1}, \end{aligned}$$

and can be fixed, by selecting values of q and r (e.g., beta(\(q=1, r=1\)); \(\rho =1/3\)), or can be estimated as a derived quantity by placing prior distributions on q and r and estimating their posterior distributions (Gelman et al. 2013). For \(\varvec{\beta }\), we used a normal distribution with mean equal to zero and variance equal to 1.5\(^2\) (Northrup and Gerber 2018). For our simulated data, we used a uniform(a, b) distribution for \(\mu \) and \(\nu \). Full conditional distributions and associated R code for each parameter are provided in Appendix 1 of ESM. The joint Bayesian posterior distribution for our data, process, and parameter models is

$$\begin{aligned}{}[\mathbf{N} , \mathbf{p} _{\text {juv}},q,r,p_{\text {adult}},\varvec{\beta },\mu ,\nu |\mathbf{C} ,\mathbf{Y} ] \propto&\nonumber \\ \prod _{i=1}^n\prod _{t=1}^{T}&[c_{i,t}|N_{i,t},p_{\text {juv},i}]^{y_{i,t}}\varvec{1}_{\{c_{i,t}=0\}}^{1-y_{i,t}}[y_{i,t}|p_{\text {adult}}][p_{\text {juv},i}|q,r] \nonumber \\&\times [N_{i,t}|N_{i,t-1},\varvec{\beta }][N_{i,1}|\mu ,\nu ][q,r,p_{\text {adult}},\varvec{\beta },\mu ,\nu ]. \end{aligned}$$
(5)

2.2 Simulated Data

To validate and assess model performance, we simulated data using the model described in Eq. 5, with known parameter values. The simulated data allowed us to assess consistency of posterior distributions with true parameter values. R script for conducting the simulation is provided in Appendix 1 of ESM. We simulated 100 data sets, each containing 100 adult deer, and a random number of offspring. Each data set was simulated using different parameter values obtained randomly from pre-specified distributions (Table 1). We used distributions that we felt would result in data that were realistic for many wildlife applications. For example, we let adult and juvenile detection probability range between 0.5 and 1 (Table 1). We also let \(\mu \) range from 1–3.5, and \(\nu \) range from 1–5. Sampling from the resulting uniform distributions for \(\mu \) and \(\nu \), and then using the sampled values of \(\mu \) and \(\nu \) to sample observations from the Conway–Maxwell Poisson distribution resulted in a distribution for \(N_{i,1}\) with median equal to 2, 0.025 \(\text {quantile} = 0\), and 0.975 \(\text {quantile} = 5\). This parameterization permits \(N_{i,1}\) to reasonably cover all possible outcomes for the number of mule deer offspring. For our simulation, we let \(\varvec{\beta }\equiv (\beta _0,\beta _1)'\) to represent two survival rates corresponding to two groups of interest [e.g., a control, \(f(\phi _i)=\beta _0\), and a treatment, \(f(\phi _i)=\beta _0+\beta _1\)]. For our simulations, we also set \(r=q=1\). After simulating data, we fit our model to the data using a Markov Chain Monte Carlo (MCMC) algorithm. An R function to fit the model to simulated and real data is provided in Appendix 1 of ESM. After we fit the model to each data set, we calculated a score (using a 0–1 score function; Williams and Hooten 2016) to describe the proportion of times the 95% credible intervals of each parameter overlapped the true parameter value used to generate the data. A summary of simulation information and the resulting scores are provided in Table 1.

Table 1 Parameters in Eq. 5, and the distribution used to generate true parameter values to simulate 100 different data sets

2.3 Mule Deer application

We fit our model to data obtained from high-resolution infrared (IR) camera surveys (Fig. 1) conducted on a marked sample of mule deer in the Ruby Mountains and Diamond Mountains in central Nevada (Elko, Eureka, and White Pine County, Nevada, USA). Marked adult female mule deer were selected for this study based on known distributions of mule deer from a concurrent study on migration ecology and the effects of predators on mule deer survival in Nevada. We used a sample of 22 mule deer. Fifteen mule deer were selected from a control area where no predator management was taking place, and 7 mule deer were selected from a treatment area where active predator removal via trapping and aerial gunning was occurring. We were interested in estimating whether there was a difference in survival rates of juveniles between the control and treatment areas. Based on our model and data, we defined survival as the probability of surviving between two successive surveys. The duration between surveys were different lengths; 22 days between the first two surveys, and 42 days between the last two surveys. Thus, our estimated survival rates were a weighted average of the first two durations. We initially allowed survival to vary between the two sampling periods, but this did not seem to have a large effect in our data. Sampling designs can be tailored to estimate necessary survival periods optimal for the study system. The number of adult deer marked and the survey days for our study were chosen based on criteria of the study design for the concurrent study and therefore were not necessarily optimal for our investigation. In Sect. 3, below, we discuss and implement a power analyses to estimate the optimal sample size required to detect a difference in survival rates for our study, if a difference existed.

We conducted 3 aerial surveys during summer 2018 (June 22–23, July 18–19, August 29–30) using a gimbal-mounted IR camera system (Schoenecker et al. 2018) attached to a fixed-wing aircraft (P68 Observer 2; Vulcanair, Naples, ID, USA). We chose our initial survey period to follow the peak birthing time in northern Nevada, so we could minimize the probability fawns would be born after our initial survey, but also have the opportunity to monitor fawns during their early stage of life when survival is lowest. Although it is possible that some births may have occurred after the initiation of the first survey on June 22nd, we chose June 22nd because less than 5% of mule deer births occurred after this date across study areas with similar latitudes in Utah and Colorado (Freeman et al. 2014). We estimated June 7–15th to be the peak birthing time in northern NV and therefore estimated that fawns would be 1–2 weeks old during the first survey. Time between the first and second survey was 22 days, time between the second and third survey was 42 days, and therefore, average time between surveys was 32 days.

We used recent telemetry fixes from global position system (GPS) radio-collars (Lotek Wireless Iridium Litetrack 420) to narrow the search area for adult marked animals and used very high frequency (VHF) telemetry to locate marked animals. All survey flights were conducted during the first 2 h of daylight to minimize thermal noise from solar radiation (Schoenecker et al. 2018). We attempted to locate all marked individuals (\(n=22\)) using the IR camera system during each survey. We spent approximately 5–10 min searching for each marked animal after we were near it. Marked adult mule deer locations and the number of juveniles with the associated adult were recorded in real-time during the survey and a video of each encounter was recorded for post-flight analysis and verification (e.g., Fig. 1). The resulting data from our surveys are shown in Table 2.

We used the same prior distributions as those used for our simulation study with one exception to tailor our analysis to our data. Instead of fixing \(q=r=1\) in our beta prior distribution for \(p_\mathrm{juv}\), we placed gamma (\(\text {shape}=0.01\), \(\text {rate}=0.01\)) priors on q and r and allowed them to be estimated (Gelman et al. 2013). We did this because observations of siblings might not be independent and estimating q and r, as opposed to fixing their values at 1, permits more flexibility in accounting for sibling dependence in detection (Hisakado et al. 2006; Martin et al. 2011). Preliminary model checking results supported this choice; the Bayesian p value when \(q=r=1\) was 0.04, and the Bayesian p value after we estimated q and r was 0.45.

3 Results

3.1 Simulation Study

The MCMC chains from each of our 100 model fits showed no lack of model convergence, supported by Gelman–Rubin diagnostics \(\le \) 1.01. More than 93 of 100 95% credible intervals of marginal posterior distributions of each model parameter covered the true values used to simulate data (Table 1), indicating that our model recovered true parameter values well (Little 2006). Posterior distributions had increased precision relative to prior distributions for all non-trivial simulated datasets (i.e., simulated datasets with at least one observation of an animal).

3.2 Mule Deer Application

The data collected from our field study are reported in Table 2. One female died prior to the birthing period and was removed from analysis and is not shown in Table 2. Of the 22 marked female adult deer, 10 were observed with two fawns at least once, 6 were observed with a maximum of 1 fawn, and 6 observed with zero fawns. Parameter estimates and 95% credible intervals of model 5 are reported in Table 3. The 95% credible interval for \(\beta _1\) (−2.21 to 3.68 ) overlapped zero, and survival estimates between the control group (no predator removal) and treatment group (predator removal) were similar (Table 3). Thus, either there was no treatment effect due to predator control or our data were insufficiently precise to detect a treatment effect of predator removal, if an effect was present.

Table 2 Data from aerial mule deer surveys

To follow-up on our inability to detect an effect of predator removal, we conducted a power analysis to estimate the number of adult animals that would be required to mark and relocate so that we could detect a 10% difference in survival between a treatment and control. We used the estimated median values reported in Table 3 for our parameter values to simulate the data, with the exception of \(\beta _1\), which we set to − 1.4. The value \(\beta _1=-1.4\) resulted in an approximately 10% lower survival rate for one group, compared to the other group. We found that if we relocated animals \(T=5\) times, we would need approximately 30 marked animals in each group (\(n=60\) total) to detect a difference in the treatment (Appendix 1 of ESM). We also examined how many radio-marked parents we would need to detect a 10% difference in survival using only \(T=3\) relocations, which was how many we used for our field study. We found the number of required animals was about 7 times more (\(n=425\)) to detect a difference between groups when \(T=3\) compared to when \(T=5\), suggesting a relatively higher information gain when increasing T versus n, for our system (Appendix 1 of ESM).

Table 3 Parameter estimates and 95% credible intervals of model 5 for the mule deer data shown in Table 2

4 Discussion

We presented a novel demographic model for estimating survival and reproduction of unmarked juvenile animals using aerial surveys that target marked adults with imperfect detection. Applying our methods to both simulated data and real data collected during a field study demonstrated that this framework is a promising tool for leveraging the dependency between juveniles and adults to estimate important juvenile demographic rates, while accounting for imperfect detection of both adults and juveniles. We provide an R function for model implementation and annotated code in Appendix 1 of ESM. Model estimates based on simulated data recovered true parameters well, with increased precision in posterior distributions compared to prior distributions, for all non-trivial simulated data sets.

There is a growing literature of methods for estimating wildlife population processes using data from unmarked individuals (e.g., MacKenzie et al. 2002; Royle 2004; Zipkin et al. 2014). However, information is lost when individuals are not marked, and this information loss can make reliable statistical modeling problematic without auxiliary data to help parameter estimation (Barker et al. 2018). The methods we present hybridize data on marked and unmarked individuals, using the auxiliary data from marked animals to help inform parameter estimation for unmarked individuals. For example, the N-mixture model of Royle (2004) assumes a closed population between replicate counts, and parameter estimates are sensitive to this assumption (Conn et al. 2018). The model we present does not assume juvenile abundance \(N_{i,t}\) is closed between subsequent surveys, only that it is non-increasing (e.g., no adoptions; equation 3). However, the relationship between parents and their offspring help to validate the closure assumption for many species (e.g., elephants, cranes, predators including wolves and fox, ungulates including deer and caribou, primates, marine mammals, bears, etc). If permanent juvenile emigration occurs during the study, survival estimates will reflect both mortality, and permanent emigration; often referred to as apparent survival (Cormack 1964; Jolly 1965; Seber 1965). Further, if juvenile adoptions or juvenile births occur after observations occur, detection probability will be under-estimated, and survival will be over-estimated. In our mule deer application, we designed our study to begin following the peak breading period to minimize the probability of new births after our initial survey. Further, many deer species are reported to be intolerant of nursing by non-offspring, and when non-offspring nursing occurs and is recognized by the doe, the non-offspring are driven away (Hirth 1985). However, although rare, exceptions do occur (see Nixon et al. 1990). Thus, investigators should design surveys to minimize potential violations of these assumptions. Alternative parameterizations to equation 3 to permit gains in juvenile abundance through time (through adoption, for example), analogous to the extensions of the N-mixture model by Dail and Madsen (2011) are also possible and have been considered in Schmidt et al. (2015).

Our field application focused on marked female mule deer and their unmarked juveniles. Estimates of mule deer juvenile survival have relied on capturing and radio-marking fawns (Hurley et al. 2014; Lomas and Bender 2007; Pojar and Bowden 2004; Unsworth et al. 1999). While radio-marking juveniles can provide valuable information on cause-specific mortality, it is often expensive and logistically challenging to collect. Alternatively, if adults have been marked, and they are being monitored synchronously to caring for juveniles, we can obtain survival rates at little or no additional cost. Further, by using aerial images to estimate juvenile survival, we can apply these methods to any spatial domain (including remote or dangerous areas) where adults can be captured and subsequently monitored with aircraft. We used aerial IR images to examine survival in remote areas in central Nevada, where ground-based field work can be challenging. Median juvenile survival rates were 0.98 (95% credible interval 0.84–1.0). Monteith et al. (2014) estimated juvenile survival rates during the first four weeks of life in the Sierra Nevada using 347 marked adult females, and 119 marked juvenile offspring. For comparison, they found week-two survival was 0.88 (\(\text {SE}=0.039\)), week-three survival was 0.96 (\(\text {SE}=0.026\)), and week 4 survival was 0.90 (\(\text {SE}=0.041\)). These estimates were similar, albeit slightly lower than our estimates. Our estimates incorporated survival information through August, a period when juvenile mule deer survival tends to be higher than in the first four weeks of life (Bishop et al. 2009; Heffelfinger et al. 2018). The period for our survival estimates was determined by our survey design, which can be tailored to any period of interest.

We examined differences in survival between an area that was undergoing active predator removal via trapping and aerial gunning and an area with no predator removal. We found no evidence that predator removal influenced juvenile survival. However, our sample was selected mainly to demonstrate proof-of-concept and was an artifact of another concurrent study. While we were able to successfully collect data to fit our model, we had a small sample size of 15 adults in the control site and 7 adults in the treatment site. It is likely that the power to detect an effect was low with our data. We provide R code for investigators interested in conducting power analyses to estimate the sample size required to detect differences in treatments, if differences exist (Appendix 1 of ESM). For our system, we estimated a sample of 30 marked animals in each group, and 5 replicate visits would be sufficient to detect a 10% difference in survival.

Our framework is the first to our knowledge that explicitly leverages the dependency between juveniles and adults to estimate juvenile survival, while accounting for parental detection probability. Despite using aerial IR imagery (Fig. 1) and using cues from radio-telemetry locations to locate marked adults, posterior median detection rate of adults was only 0.74. Thus, assuming perfect detection of adults may be misleading for many systems, especially when using aerial surveys in thick or dense terrain. This highlights the importance of using methods that explicitly account for adult and/or juvenile detection probability for designing efficient aerial survey methodology.

The modeling framework we developed is not exclusive to aerial survey data. For example, if adults can be uniquely identified, then the framework we present can be generally applied to data collected from camera traps, nest searches, brood surveys, or any other observation method relying on counting juveniles of a marked parent.