Estimating arthropod survival probability from field counts: a case study with monarch butterflies

Survival probability is fundamental for understanding population dynamics. Methods for estimating survival probability from field data typically require marking individuals, but marking methods are not possible for arthropod species that molt their exoskeleton between life stages. We developed a novel Bayesian state‐space model to estimate arthropod larval survival probability from stage‐structured count data. We performed simulation studies to evaluate estimation bias due to detection probability, individual variation in stage duration, and study design (sampling frequency and sample size). Estimation of cumulative survival probability from oviposition to pupation was robust to potential sources of bias. Our simulations also provide guidance for designing field studies with minimal bias. We applied the model to the monarch butterfly (Danaus plexippus), a declining species in North America for which conservation programs are being implemented. We estimated cumulative survival from egg to pupation from monarch counts conducted at 18 field sites in three landcover types in Iowa, USA, and Ontario, Canada: road right‐of‐ways, natural habitats (gardens and restored meadows), and agricultural field borders. Mean predicted survival probability across all landcover types was 0.014 (95% CI: 0.004–0.024), four times lower than previously published estimates using an ad hoc estimator. Estimated survival probability ranged from 0.002 (95% CI: 7.0E−7 to 0.034) to 0.058 (95% CI: 0.013–0.113) at individual sites. Among landcover types, agricultural field borders in Ontario had the highest estimated survival probability (0.025 with 95% CI: 0.008–0.043) and natural areas had the lowest estimated survival probability (0.008 with 95% CI: 0.009–0.024). Monarch production was estimated as adults produced per milkweed stem by multiplying survival probabilities by eggs per milkweed at these sites. Monarch production ranged from 1.0 (standard deviation [SD] = 0.68) adult in Ontario natural areas in 2016 to 29.0 (SD = 10.42) adults in Ontario agricultural borders in 2015 per 6809 milkweed stems. Survival estimates are critical to monarch population modeling and habitat restoration efforts. Our model is a significant advance in estimating survival probability for monarch butterflies and can be readily adapted to other arthropod species with stage‐structured life histories.


INTRODUCTION
Survival probability is a fundamental demographic vital rate that is an important parameter in population models used to support research, management, and conservation. Survival probabilities are needed to parameterize stage-based matrix population models (Caswell 2001), to identify risks to population viability (Crouse et al. 1987), and to propose actions to reach management objectives (Crowder et al. 1994). Estimation of survival probabilities for most vertebrates depends on marking live individuals and later re-encountering live or dead individuals under a variety of sampling designs (Cormack 1964, Jolly 1965, Seber 1965, Lebreton et al. 1992, Williams et al. 2002. In contrast to vertebrates, survival probability estimation for arthropods presents unique field sampling and statistical challenges. Mark-recapture methods are ineffective for immature stages because arthropods molt their exoskeleton between life stages. In many insects, survival from egg to adult is often less than five percent, which may require that hundreds of eggs be followed to obtain an adequate sample size of surviving adults (Wood 1994, De Anda and. Furthermore, frequentist (non-Bayesian) statistical methods may have difficulty with numerical estimation of very small probabilities, especially when sample sizes are small (King et al. 2009). Thus, innovative Bayesian analytical approaches that can use field counts of unmarked individuals are urgently needed, particularly given the number of arthropods at risk of extinction (Dirzo et al. 2014).
Organisms with progressive life stages can occur in multiple cohorts as eggs are laid each day. Field counts of such populations typically involve periodic counts of individuals from each life stage at a field site. Such data have been called stage-frequency data (Hoeting et al. 2003), stagestructured data (Manly 1990), or cohort data (de Valpine and Knape 2015). Several approaches have been proposed to derive statistical inference from such data. Pinder et al. (1978) and Oberhauser et al. (2001) estimated the shape of the survivorship curve from egg to pupation, but this approach ignores the stage structure of the population. Manly (1990) reviewed methods through the 1980s, but the intractability of the problem precluded substantial development until computational capacity increased. Hoeting et al. (2003) developed method-of-moments estimators based on the Laplace transform and associated variance estimators to bypass computational challenges with maximum-likelihood estimation. More recently, computer-intensive maximum-likelihood estimation algorithms have been reported (de Valpine and Knape 2015). Increased computational capacity also brought with it the advent of Bayesian approaches (Knape and de Valpine 2016). Although the model of Knape and de Valpine (2016) is useful for abundant species that are easily sampled, it may not be suitable for lowdensity populations, such as many species at risk (Johst et al. 2006).
Uncertainty in monarch butterfly (Danaus plexippus) population models is illustrative of the need for enhanced methods to estimate larval survival probabilities. Monarch butterfly populations in eastern North America have declined by 80% over the past two decades (Brower et al. 2012, Semmens et al. 2016. Monarch butterflies are a multivoltine species with five to six generations per year (Batalden et al. 2007). Monarchs migrate from Mexico to the United States and Canada each spring to reproduce (Malcolm et al. 1993, Miller et al. 2012). Most of the eggs laid by the overwintering generation are deposited in Texas and Oklahoma, USA, as they travel north. Successive generations lay eggs across the eastern United States, east of Rocky Mountains and north to southern Canada, which is the range limit of milkweed species (Asclepias ssp.), their obligate host plant (Flockhart et al. 2013). This colonization pattern results in stochastic spatiotemporal variation in population density (Prysby and Oberhauser 2004) that is well documented in egg and larval counts performed by citizen-science programs (Ries and Oberhauser 2015).
Monarchs have an egg stage, five instars, a pupal stage, and an adult stage. At typical summer temperatures, eggs hatch after 3-4 d (Zalucki 1982) and larval development takes 9-14 d (Oberhauser 2004). Field surveys for immature monarchs consist of repeated egg and instar-specific counts in fixed habitat patches. Pupae cannot be reliably counted because they are well camouflaged and typically not located near host plants. Surveys are typically conducted once per week (Prysby and Oberhauser 2004, ❖ www.esajournals.org 2 April 2020 ❖ Volume 11(4) ❖ Article e03082 Pitman et al. 2018). Survival estimates based on such surveys using an ad hoc estimator suggest that survival probabilities may be decreasing over time (Nail et al. 2015). Continental-scale population models of monarch butterflies suggest larval survival probabilities are an important component of population growth rate and critical to conservation of this species (Flockhart et al. 2015, Oberhauser et al. 2017. While survival probabilities of egg and larval stages are critical for understanding monarch population dynamics, limited techniques are available to estimate survival probabilities of these life stages (Grant and Bradbury 2019). In some cases, researchers have attempted to follow individual monarch eggs and larvae (De Anda andOberhauser 2015, Myers et al. 2019). However, after the third instar, larvae typically abandon their natal host plant and cannot be reliably followed (De Anda and Oberhauser 2015). To date, most estimation of survival probabilities has relied on the ratio of counts of the last larval stage (fifth instar) to egg counts (Oberhauser et al. 2001, Nail et al. 2015. These estimates do not explicitly model temperature-dependent larval development, latent cohort structure present in field counts, nor detection bias. We developed a Bayesian state-space model to estimate survival of monarch eggs and larvae from field counts. This model takes advantage of deterministic stage duration in arthropods and uses ambient temperatures from field sites to predict stage duration for each cohort (Zalucki 1982). We tested our model for potential bias in parameter estimates using simulation studies. Hypothesized sources of bias include variable detection probability among stages and individual variation in stage duration (de Valpine et al. 2014). We then applied our model to field data from Iowa, USA, and Ontario, Canada. Our model provides the first rigorous estimates of survival probability for monarch butterflies and can be readily adapted to many other arthropods for which stage duration can be predicted with ancillary information.

Model notation
Our model applies to organisms that progress through successive life stages, occur in multiple cohorts, and for which marking individuals is not possible or logistically feasible. Stages are denoted i = 1, 2, . . ., I. The study continues for j = 1, 2, . . ., M days. Time intervals other than days may be used, if appropriate. Daily temperatures at the study site are used to predict stage duration and are denoted t = (t 1 , . . ., t M ). Other data may be used to predict stage duration, if appropriate. Field data consist of counts, y ij , of individuals in each stage i on each day j. Counts do not need to be conducted on each day of the study.

Assumptions
We describe seven assumptions for this model. Researchers should carefully consider these assumptions when designing a study. The first assumption is that differences in cohort populations between time periods are due to mortality, which is the basis for inference in this model. The model controls for latent cohort population structure within observed counts and random variation in counts so as to estimate the differences in cohort populations each day. Variation in detection probability between stages can violate this assumption because counts may be relatively lower for stages that are more difficult to detect. Emigration from the survey site or immigration into the survey site could also violate this assumption. If emigration and immigration are equal and survival probabilities are the same for individuals who immigrate and emigrate, this assumption is not violated. We employ simulation studies to address violations of this first assumption. The second assumption is that stage durations are the same for each individual within a cohort. This assumption may cause bias in parameter estimates if there is individual variation in stage duration (de Valpine et al. 2014). We perform simulation studies to address potential bias due to this assumption. The third assumption is that counts are Poisson-distributed, a common assumption for count data. The fourth, fifth, and sixth assumptions are that survival probabilities are the same for all individuals within a stage, eggs are laid instantaneously each day at the same time, and sampling occurs instantaneously immediately after eggs are laid. These are typical, though often implicit, assumptions for mark-recapture and other ecological studies. These assumptions should be addressed in field ❖ www.esajournals.org 3 April 2020 ❖ Volume 11(4) ❖ Article e03082 study designs to minimize sources of variability in survival and underscore the importance of conducting surveys at the same time each day. The seventh assumption addresses survival probability between stages. For stage-specific survival probabilities, the question arises: For the time period over which a transition to a new stage occurs, should the survival probability be assigned to the pre-transition stage or post-transition stage? For example, if a larva is found as a second instar on day 5 and third instar on day 6, should the survival probability between day 5 and day 6 be modeled as the second or third instar survival probability? In our model, the survival probability is assigned to the second instar. Put another way, we assume the individual is a second instar until immediately before sampling. It then instantaneously molts to a third instar and is available to be detected. Thus, the seventh assumption of the model is that for the time interval between two life stages, survival probability in the model is assigned to the first stage. A similar issue occurs with mark-recapture multistate models (Brownie et al. 1993), and we employ the same solution.

A Bayesian state-space model
We used a Bayesian state-space framework to model the relationship between field counts and estimated parameters (Royle and Dorazio 2008, K ery and Schaub 2012, Newman et al. 2014. In this modeling framework, the data-generating observation process and associated measurement errors are distinct from the unobservable state processes and associated process variation. Thus, a state-space model is comprised of an observation model and a state model, and the observation model is conditioned on the state model. These models have also been called Bayesian hierarchical models because of the hierarchical relationship between the observation and state models (Royle and Dorazio 2008). This framework allows for better control of observation or measurement error that is a common concern in ecological studies (Royle and Dorazio 2008). In our observation model, the counts are modeled as drawn from a Poisson distribution with expected count ' ij : It is possible that ' ij = 0, rendering the Poisson a degenerate distribution. In such cases, the likelihood function is defined as Expected field counts are conditional on state processes driven by three types of parameters: days spent in each life stage (stage duration); number of eggs laid in each cohort, b k , for k in 1, . . ., M; and daily survival probability for each stage, p i , for i in 1, . . ., I.
Stage duration is calculated using daily temperature data, t. By calculating stage duration from ancillary data and fixing it in the likelihood, greater statistical power is available to estimate b k and p i . The accumulated degree days for each cohort k up until day j is calculated asD jk ¼ P j h¼k t h . The stage s jk for cohort k on day j is determined by D jk , following the thresholds for stage transitions reported by Zalucki (1982): if D jk 2 77:3; 105:1 ½ Þ 4 if D jk 2 105:1; 129:6 ½ Þ 5 if D jk 2 129:6; 165:3 ½ Þ 6 if D jk 2 165:3; 231:9 ½ Þ 7 i fD jk ! 231:9 Stage duration, c ik , is the number of days in which cohort k is in stage i: The expected counts, ' ij , are the sum of cohort populations that are in stage i on day j. The contribution to this sum from each cohort is the number of eggs originally laid in that cohort, b k , multiplied by the cumulative survival probability up until current day j. Thus, the expected count ' ij for stage i on day j is where the indicator function I takes the value of 1 if s jk = i, and 0 otherwise; p ik is the survival probability for stages 1 through i -1; and p i is the survival probability for the current stage i raised to the power of how many days the cohort has been in the stage. The cumulative survival probability up until current day j for cohort k is the product of the ❖ www.esajournals.org 4 April 2020 ❖ Volume 11(4) ❖ Article e03082 daily survival probabilities over i stages to day j (p ik Á p jÀkÀb ik i in Eq. 4). The survival probability for stages 1 through i-1 for the current stage i is where h are stages, p h is the daily survival rate for stage h, and c hk is the stage duration for stage h of cohort k. The survival probability for the current stage is p jÀkÀb ik ; where b ik is the number of days cohort k has been alive before stage i: b ik ¼ P iÀ1 h¼1 c hk : b k is modeled as drawn from a Poisson distribution with mean k: b k~P oisson(k), and a Gamma prior distribution on k: k~Gamma(0.001, 0.001). Daily survival probabilities, p i , are given a uniform prior distribution (p i~U nif(0,1)).
Although the details are complex, the overall structure of the state model is relatively simple. The expected field counts (' ij ) are the total number of individuals counted from all cohorts that are present in stage i on day j. However, these cohort populations are latent unobservable parameters for which we need more information to estimate. Ancillary information on stage duration from field data is used to determine which cohorts are in stage i on day j. For monarch butterflies, the individuals in stage i on day j are determined by temperature (Zalucki 1982); hence, stage durations can be calculated for each cohort based on daily mean temperatures at field sites and the accumulated degree days over time. Predicting stage duration from temperature data allows better estimation of survival probabilities with smaller sample sizes, thereby overcoming a need for large sample sizes to estimate stage duration and survival probability concurrently (Knape and de Valpine 2016).
Cumulative survival probability per stage, p i,c , and the cumulative survival probability from oviposition to pupation, p c , can be derived from p i . We calculate stage survival probability and cumulative survival probability using mean stage durations across cohorts. For example, cumulative egg survival probability is calculated as, where p 1 is daily egg survival probability, and c 1 is mean stage p 1;c ¼ p c 1 1 duration for eggs in the study. Cumulative survival probability from oviposition to pupation, p c , is calculated as Our model estimates survival from one site, where a site is defined by the researcher to meet the assumptions of the model and their research objectives. The complexity and assumptions of our model render inclusion of multiple sites in one analysis difficult, if not impossible. Consequently, survival probabilities from different sites are best analyzed together in a two-stage analysis or meta-analysis using standard Bayesian regression methods.

Simulation studies
Several biological factors such as detection probability (Royle and Dorazio 2008) and individual variation in stage duration (de Valpine et al. 2014) may introduce bias into survival estimates from stage-frequency data. In a Bayesian context, we desire the median of the posterior distribution to have minimal difference from the true value of the parameter (Gelman et al. 2013). If such bias cannot be controlled through study design or statistical modeling, we want to understand the magnitude of the bias and the effect it may have on ecological inferences.
To better understand bias introduced by factors such as detection probability and individual variation in stage duration, we simulated monarch butterfly count data using published survival probabilities (Oberhauser et al. 2001, Flockhart et al. 2015. We then analyzed the simulated data and compared the estimated parameters with the true parameters used to simulate the data. We begin with a simple test case scenario (Scenario 1 in Table 1) and build realistic complexity onto the test case in subsequent scenarios. For Scenario 1, true values of daily eggs laid for ten days are chosen to replicate the first wave of egglaying monarchs passing through breeding habitat during spring migration: b 1 , . . ., b 10 = {50, 100, 200, 500, 750, 500, 200, 100, 50, 25}. The true daily survival probabilities were calculated from Oberhauser et al. (2001) following Flockhart et al. (2015; see Table 2).
The simple scenario provides a baseline for comparison to more complicated scenarios to assess model robustness and provide guidance for study design (Scenarios 2-6 in Table 1). To test the second assumption of our model-stage durations are the same for each individual within a cohort-we simulated cohort data where each individual monarch had a different stage duration (Scenario 2). Stage duration was taken from a truncated normal distribution based on the ❖ www.esajournals.org 5 April 2020 ❖ Volume 11(4) ❖ Article e03082 empirical results of Zalucki (1982). To test the effect of stage-varying detection probability, a possible violation of the first assumption of our model, we simulated data where detection probability varied among stages (Scenario 3). To provide practical advice for study design, we simulated data with varying sampling intervals and population sizes, with and without individual variation in stage duration (Scenarios 4 and 5). Finally, we simulated data similar to our example datasets to determine possible bias in survival estimates (Scenario 6). R code for simulations is available at https://github.com/tgrant7. Further details on methods and results for Scenarios 2-6 are provided in Appendix S1. For all simulation scenarios, we simulated count data with known parameter values (h ¼ b 1 ; . . .; b M ; p 1 ; . . .; p 6 Þ for 100 datasets and fit the model to the data to estimate the marginal posterior distributionsĥ ¼b 1 ; . . .;b M ;p 1 ; . . .;p 6 for each dataset. The mean median of the posterior distributions ( h) was compared with the true parameter value (h) for bias calculations. We calculated absolute bias as h À h and percent bias as h À h =h Â 100. We define minor bias as a bias of <10%, moderate bias as 10-50%, and severe bias as >50%.
Simulations and analysis were performed in R (R Core Team 2016). For analysis, we employed Markov chain Monte Carlo (MCMC) methods in a Bayesian framework to estimate parameters using JAGS (Plummer 2003) and the R package runjags (Denwood 2016). We used non-informative, flat priors: Uniform (0,1) on p parameters and Gamma (0.001,0.001) on k parameters. We sampled three independent Markov chains for 10,000 iterations after discarding the first 4000 burn-in iterations and 1000 adaptive iterations. Convergence was evaluated by visually inspecting chains for convergence and ensuring the potential scale reduction factor,R; was <1.05 (Gelman et al. 2013). To reduceR, additional iterations were added if necessary. Parameters without a sufficiently smallR were not used for inference.

Estimating monarch butterfly survival probabilities
We conducted monarch butterfly field counts in relatively large areas to meet the assumptions of the Poisson distribution. Effort was kept constant for each survey day to ensure that differences in observed counts are due to mortality and not a result of sampling effort. Because the population consisted of multiple cohorts developing at different rates, observed counts Notes: The number of days in the study, M, is the number of days that eggs could have been laid, and cohorts are the number of days eggs were laid. Sampling frequency is the days of the study when counts were collected. Sampling frequency was daily, once weekly, twice weekly, and two consecutive days weekly. In all scenarios, mean daily temperature was assumed to be a constant at 21°C. consisted of individuals from several cohorts. An example of a monarch butterfly multi-cohort population over time is provided in Fig. 1. The latent population in each stage on each day is the sum of the individual cohort populations in that stage on that day. We estimated monarch survival probabilities from 19 field sites in Iowa, USA, and Ontario, Canada, in four study areas: Iowa road right-ofways (ROWs), Ontario ROWs, Ontario agricultural field borders, and Ontario natural areas. Sites chosen for analysis had at least one 5th instar observation. Iowa ROWs were two sites surveyed in 2015 along rural roadsides bordered by corn and soybean fields in Story and Boone Counties, Iowa (R. Bitzer et al., unpublished data). The 17 sites in Ontario surveyed in 2015 and 2016 are described in Pitman et al. (2018) and include five ROW sites, six agricultural field border sites, and six natural area sites. Agricultural border sites consisted of the first three rows of a herbicide-treated corn or soybean. Right-ofway sites consisted of the vegetation along roads. Natural areas consisted of restored meadow sites or private gardens. Counts were made approximately once per week starting in mid-July and ending in late August, with the earliest survey date being 13 July and the latest being 31 August. Among the 19 sites, the number of egg and fifth instar observations per site ranged from 27 to Cohort 1 is laid on day 1, cohort 2 is laid on day 2, and so on. Darkest gray corresponds to egg stage (Stage 1), progressively lighter shades of gray correspond to instars, and the lightest shade of gray corresponds to the fifth instar (Stage 6). The total population of stage i observed on day j is the sum of cohort populations. For example, the total population in Stage 2 on day 6 is 40, comprised of 10 individuals from cohort 2 and 30 individuals from cohort 3. ❖ www.esajournals.org 7 April 2020 ❖ Volume 11(4) ❖ Article e03082 627, and 1 to 23, respectively. Daily mean temperature data were obtained from the Ames 8 WSW weather station from the Climate Data Online portal (NOAA 2018) and the Delhi CS weather station in the Historical Climate Data Web portal (Government of Canada 2018). Daily mean temperature was calculated as the mean of high and low temperatures for each day (Arguez et al. 2012). The mean daily temperature was 21.8°C (range 15.3-28.3) and 20.9 (range 13.0-26.8) for Iowa and Ontario sites, respectively. Analysis of field data followed the same methods as for analysis of simulated data, except we sampled seven independent Markov chains, after discarding the first 4000 burn-in iterations and 1000 adaptive iterations, until convergence was reached. Convergence was evaluated by visually inspecting chains for convergence and ensuring the potential scale reduction factor,R, was less than 1.05 (Gelman et al. 2013). Additionally, 16 leading days (the maximum time from egg to pupation) are added to the study period to allow the model to estimate eggs laid before surveys began.
In our meta-analysis, we used survival probabilities from 18 sites (one site would not converge) as the response variable in a Bayesian regression. The cumulative survival rates were treated as drawn from a normal distribution p ci $ N l; r 2 À Á ; where = a + d 1 * X 1i + d 2 * X 2i + d 3 * X 3i , and X 1 , X 2 , and X 3 are indicator variables corresponding to landcover types. This analysis was also conducted in R (R Core Team 2016) using JAGS (Plummer 2003) and the R package runjags (Denwood 2016).
Lastly, we calculated the monarch production as adults produced per milkweed stem by multiplying cumulative egg, and larval and pupal survival probabilities by the number of eggs counted per milkweed stem. We used egg and larval survival probabilities estimated with our model and the pupal survival probability of 0.76 reported by Nail et al. (2015). Our objective was to estimate the relative number of adult monarchs produced from different landcover types over one year in each study site. Thus, we summed the number of eggs counted in each landcover type over the summer as an index of the total number of eggs laid over the season.

Simulation studies
Bias was minor for all parameter estimates in the Scenario 1 simulated dataset (Table 3). Cumulative survival probability from oviposition to pupation, p c , also had low bias and high precision. The greatest percent bias (~2-3%) occurred in estimation of second instar daily (p 3 ), fourth instar daily (p 5 ), second instar stage (p 3,c ), fourth instar stage (p 5,c ), and fifth instar stage (p 6, c ) survival probabilities.
Percent bias was moderate to severe for nearly all parameters when individual variation in stage durations was added; however, absolute bias for cumulative survival (p c ) was only 0.023 (Scenario 2; Appendix S1: Table S1). Given the bias induced in other parameters by individual variation in stage duration, we report primarily on p c for the remaining simulations. Bias induced by detection probability had a very different effect on parameter estimates depending on whether eggs or fifth instars were most detectable (Scenario 3; Appendix S1: Table S3). When fifth instars were assumed to have lower detectability than eggs, bias was minor in all parameter estimates except p c , which had moderate negative bias. When eggs had lower detectability than fifth instars, there was a moderate positive bias in p c .
Scenario 4 simulations showed inherent model bias caused by study design without individual variation in stage duration. Under these conditions, bias was minimized with daily sampling or twice-weekly sampling (e.g., every Tuesday and Friday); bias was greatest for once-weekly sampling (Appendix S1: Fig. S1). When individual variation in stage duration was included to provide a more realistic dataset (Scenario 5), the pattern was reversed, with bias minimized with weekly sampling and greatest with daily sampling (Appendix S1: Fig. S3). Individual variation in stage duration decreased precision (Appendix S1: Fig. S2; Appendix S1: Fig. S4). Precision was increased in all cases by increasing the sample size of counts.
For scenario 6, a dataset was simulated to represent the average case where the true value of cumulative larval survival was p c = 0.014. The daily p values were based on the Iowa ROW2 site. The mean estimate based on the simulated ❖ www.esajournals.org 8 April 2020 ❖ Volume 11(4) ❖ Article e03082 data wasp c ¼ 0:017 (standard deviation [SD] = 0.007). Mean bias was 0.003, and percent bias was 21%.

Analysis of monarch butterfly field counts
Convergence and parameter estimation were achieved when fitting the model to data from 18 of the 19 sites; data from an agricultural field border site would not converge, apparently because of an identifiability problem with some of the daily survival parameters. Estimated cumulative survival probability from egg to pupation for the remaining 18 sites ranged from 0.002 (95% CI: 7.0EÀ7 to 0.034) to 0.058 (95% CI: 0.013-0.113). Survival estimates for individual sites are reported in Appendix S2, where they are also compared with survival estimates from the ratio estimator.
Monarch production varied by year, landcover type, and location. For easier comprehension, monarch production was scaled to adults produced per 6809 milkweed stems (the number of milkweed needed to produce one adult monarch in the lowest production landcover type). Adult production varied from 1 (SD = 0.68) in Ontario ROWs in 2016 to 29.0 (SD = 10.42) in Ontario agricultural field borders in 2015 (Fig. 3). In Notes: Mean and standard deviation (SD) of estimated parameter values ( h) with mean 95% credible interval (CI) width from analysis of the simulated datasets. Bias ( h À h) and percent bias h À h =h Â 100 of estimated parameters. Parameters: b k is number of eggs laid per day, p 1 À p 6 are daily survival probabilities for each stage (egg through fifth instar), p 1,c À p 6;c are cumulative survival probabilities per stage, and p c is cumulative survival probability from oviposition to pupation. For these simulated datasets, temperature was held constant at 21°C for each day of the study.
❖ www.esajournals.org 9 April 2020 ❖ Volume 11(4) ❖ Article e03082 Ontario, 2016 had higher production than 2015, and agricultural field borders had highest production, while ROWs had lowest production. Iowa ROWs had production similar to Ontario agricultural field borders (Fig. 3).

DISCUSSION
Arthropod populations are declining worldwide, and many species with threatened or endangered status are data-deficient (Dirzo et al. 2014), lending urgency to the need for improved methods to estimate population vital rates. We present a new Bayesian state-space model to estimate survival probabilities from field counts of stage-structured arthropod populations. The model takes advantage of known stage durations, which can be estimated with ancillary degree-day information, and does not require large sample sizes. These characteristics make the approach suitable for a large number of arthropods routinely monitored in the field. Using this model, we obtained the first estimates of monarch butterfly survival from egg to pupation based on a standard statistical theory with well-defined assumptions, known bias properties, and measures of precision. We found that monarch survival and production vary across space, time, and landcover type. While our work was developed to support monarch butterfly conservation research, the only substantial change needed to adapt the model for other arthropod species would be adjustments for the number of stages and their durations.
New statistical models should be evaluated for bias and precision under various conditions. Statistical bias is introduced with measurement error caused by ecological processes. We addressed many of the potential biases reported over decades of population ecology research to provide a picture of estimation performance under realistic conditions. We addressed two major potential biases: detection probability and individual variation in survival probability.
A detection probability different from one is well known to cause bias in survival estimation (Royle and Dorazio 2008) and induce bias in closed mark-recapture estimators. Quantitative ecologists have spent substantial effort in developing estimators to control for detection probability (Pledger 2000), and this is an ongoing area of research (Morgan and Ridout 2008). The detection probability in Poisson models is different than in typical mark-recapture models. ❖ www.esajournals.org 10 April 2020 ❖ Volume 11(4) ❖ Article e03082 Detection probability in mark-recapture models is defined as the probability that an individual is detected/captured in a particular sampling occasion, given it is present (Lebreton et al. 1992). At the population scale, this is realized as the proportion of the population that is detected. We have modeled counts using the Poisson distribution where the total population size is unknown. The relationship between the rate parameter of the Poisson distribution (k) and detection probability is different than the relationship between population size and detection probability in mark-recapture models. The rate parameter for our model (the expected number of individuals in stage i observed on day j) is a function of population size, survey effort, and per-individual detection probability. However, because the population size is never known and does not enter directly into survival estimation, the per-individual detection probability is not relevant. What is needed is relative detection probability: the number of individuals actually detected relative to the number that would be detected if detection probability was equal across all stages. In practice, the most detectable stage then has a relative detection probability of 1, and other stages have a relative detection probability between 0 and 1. Detection bias can be caused by inherent variation in detectability among stages and also by variability in observer detection skill. Empirical studies are needed to determine these detection probabilities. Detection probability could then be incorporated into our model as an assumed value, in manner similar to the temperature data. While it would be preferable to estimate detection probability from the data, that is not likely possible. N-mixture models have been used to estimate detection probability from counts when estimating abundance (Barker et al. 2018), but these models require that abundance stay constant over several surveys. Many arthropods, such as monarch butterflies, move through stages in a single day (Zalucki 1982); thus, the logistics of multiple surveys per stage are formidable. Our simulations demonstrate that absolute bias due to detection probability is likely small and the survival estimates are manyfold more accurate than previously available estimates, but incorporating detection probability will likely improve estimates. ❖ www.esajournals.org 11 April 2020 ❖ Volume 11(4) ❖ Article e03082 Individual variation in stage duration could potentially cause severe bias in survival estimates (de Valpine et al. 2014). Our model can be classified as a delay-differential model (de Valpine et al. 2014), and a central critique of delaydifferential models is that individual variation in stage duration is not represented, which can induce bias (de Valpine et al. 2014). While this bias cannot be entirely eliminated, we determined the magnitude of bias and the limitations it places on inference. For monarch butterflies, absolute bias is low for cumulative survival probability and does not seriously hamper inference, especially relative to other available methods. de Valpine et al. (2014) also noted that correlated stage durations, wherein individuals with short or long stage durations in one stage also have correspondingly short or long stage durations in other stages, can also have an effect on parameter estimates. Currently, there are no field data to determine whether or not there is correlation in stage duration for individual monarchs in the field. Such data are needed to evaluate any potential effect.
Ecological factors such as detection probability and individual variation in stage duration are nuisances that ecologists must grapple with in any survival studies. The fundamentals of the model we have presented are sound, as shown in Scenario 1 and Table 1. By providing realistic simulations to determine bias, we have provided an objective picture of model performance under realistic conditions. Because data collection processes are imperfect, it is important to understand how ecological processes affect data collection and influence model estimates.
We performed several simulations to evaluate model performance under a variety of realistic field sampling scenarios. Sampling frequency and sample size, in addition to detection probability and individual variation in stage duration, can induce bias to greater or lesser degrees. Given the bias that may be induced from various sources, our model is best suited for estimating cumulative monarch larval survival probability, p c . Bias in p c can be minimized by careful study design using the simulation tools we have developed. Under realistic conditions of individual variation in stage duration, once-weekly or twice-weekly sampling had the least bias in p c , with once-weekly sampling having no bias for some sample sizes (Appendix S1: Fig. S3). If researchers require estimates of daily and stagespecific survival probabilities for monarch larvae, simulations should be conducted to determine potential bias in these parameters under different potential study designs because of the potential estimation bias for these parameters.
We assumed that the studies of Zalucki (1982) accurately reflect the development rates of wild monarchs. Zalucki (1982) fed his monarch larvae on Gomphocarpus fruticosa (syn. Asclepias fruticosa), a species not present in North America. It is possible that different species of milkweed or differing quality of milkweed could affect development rates. Further empirical research into developmental rates on different milkweed species would be a valuable improvement to the model. The model can be easily updated by changing the number of degree days to transition between stages.
Our cumulative survival estimates were lower than previous estimates and varied among landcover type, location, and year. Overall, cumulative monarch larval survival (p c = 0.014) was four times lower than a previously published estimate that used the ratio estimator on 18 yr of data (p c = 0.058; Nail et al. 2015). Our results are, however, consistent with an empirical study that found cumulative survival of 0.017 through the third instar (De Anda and Oberhauser 2015). Our survival estimate for individual sites was always substantially lower than ad hoc estimates using the ratio estimator (Appendix S2: Table S1). These results suggest that previous estimates of monarch larval survival based on ratios of eggs to fifth instars are biased high. This may not be surprising since the ratio estimator does not account for temperature-driven stage durations or separate cohorts. The ratio estimator has not been tested for bias and does not provide measures of precision, such as CIs. The assumptions of the ratio estimator, in terms of cohorts and developmental rates, are also not known. Nevertheless, the analysis of Nail et al. (2015) covered an extensive geographical area and time period and the relative survival trends are likely real.
Previous research suggests that larval survival probabilities may differ among landcover or habitat types (Oberhauser et al. 2001, Nail et al. 2015, Haan and Landis 2019b, Myers et al. 2019).
❖ www.esajournals.org 12 April 2020 ❖ Volume 11(4) ❖ Article e03082 Our improved estimator found that larval survival was higher in agricultural borders, consistent with Oberhauser et al. (2001) and Myers et al. (2019), and lower in right-of-ways and natural areas. Variation in survival across landcover types, in concert with field studies that show variation in egg density among landcover types (Oberhauser et al. 2001, Pleasants and Oberhauser 2013, Pitman et al. 2018) and variation in milkweed density (Hartzler and Buhler 2000), suggests that productivity varies substantially across space and time. Predation rates may be an important factor causing survival to vary among landcover types. Several authors have found that arthropod predators are more abundant in grasslands (Werling et al. 2011, Pitman et al. 2018, Haan and Landis 2019a, though some authors have not found varying predation rates despite higher observed numbers of predators (Pitman et al. 2018). We found that monarch production was consistently highest in Ontario in agricultural field borders, lowest in right-of-ways, with natural areas in between. However, right-ofways in Iowa had monarch production levels similar to agricultural field borders in Ontario.
Our results add to the evidence that survival and monarch production may vary between habitat types and geographical areas, which may have implications for monarch conservation and milkweed restoration planning. Nail et al. (2015) used their average cumulative survival probability of 0.058 to estimate the number of milkweed stems needed to produce one adult monarch and thus how many milkweed stems may be needed to increase monarch production in the U.S. summer breeding range. By multiplying an estimated survival of 0.058 by the average seasonal egg density of 0.60 and taking the inverse, Nail et al. (2015) calculated that one adult was produced from 29 milkweed stems, on average. We estimated one adult produced per 340-6809 milkweed stems. The difference between our estimates and the Nail et al. (2015) estimate is a function of lower estimated survival and lower egg density in the sites we sampled. Eggs per milkweed stem can vary substantially from year to year. Our egg densities ranged from 0.015 to 0.43. It seems unlikely that egg density was at carrying capacity in any of our study years; thus, we cannot conclude that 340-6809 stems are needed to produce one adult.
We believe the calculation of stems needed per adult of Nail et al. (2015) should be used with care. Simply because 29 or 6809 milkweed were observed to produce one adult, it does not mean that 29 or 6809 milkweed are necessary to produce one adult; it is possible that fewer are necessary. The minimum number of milkweed that produced one adult over time is likely a better measure of the number of milkweed needed, rather than the average number of milkweed. Future calculations of the number of milkweed needed to increase milkweed production should incorporate our improved survival estimates. Maximizing monarch production likely involves not only milkweed numbers and density, but also the spatial arrangement of milkweed on the landscape (Grant et al. 2018).
Our model provides a framework to conduct applied analyses of spatiotemporal larval arthropod survival. At the patch scale, larval survival estimates in different habitat types will help evaluate the hypothesis that survival rates are lower in natural grasslands compared with agricultural fields (Pitman et al. 2018). This model can support development of habitat restoration strategies that increase monarch production at the landscape-scape scale (Grant et al. 2018) and allow integration with conservation of other species as well. At continental scales, analyses could test life history trade-offs between migration and recruitment across space and time (McKinnon et al. 2010, Chapman et al. 2012. Monarch eggs and larvae have been counted at hundreds of sites across North America by citizen scientists for the Monarch Larval Monitoring Program (Ries and Oberhauser 2015). Analysis of large datasets would allow consideration of the local habitat conditions such as milkweed patch size and milkweed density (Pitman et al. 2018) and habitat configuration at the landscape scale (Grant et al. 2018) that affect egg density. The model described here can likely be used with these aggregate data if they can be vetted such that the model assumptions are met. Applying these survival estimates has implications in quantifying regional productivity to predict monarch overwintering population size (Flockhart et al. 2017) and for national population models that seek to develop monarch conservation plans focused on habitat restoration at continental scales (Flockhart et al. 2015, Oberhauser ❖ www.esajournals.org 13 April 2020 ❖ Volume 11(4) ❖ Article e03082 et al. 2017, Thogmartin et al. 2017). Similar advances in knowledge are possible for other species with application of improved methods for estimating survival probabilities.