Use of stochastic patch occupancy models in the California red-legged frog for Bayesian inference regarding past events and future persistence.

Assessing causes of population decline is critically important to management of threatened species. Stochastic patch occupancy models (SPOMs) are popular tools for examining spatial and temporal dynamics of populations when presence-absence data in multiple habitat patches are available. We developed a Bayesian Markov chain method that extends existing SPOMs by focusing on past environmental changes that may have altered occupancy patterns prior to the beginning of data collection. Using occupancy data from 3 creeks, we applied the method to assess 2 hypothesized causes of population decline-in situ die-off and residual impact of past source population loss-in the California red-legged frog. Despite having no data for the 20-30 years between the hypothetical event leading to population decline and the first data collected, we were able to discriminate among hypotheses, finding evidence that in situ die-off increased in 2 of the creeks. Although the creeks had comparable numbers of occupied segments, owing to different extinction-colonization dynamics, our model predicted an 8-fold difference in persistence probabilities of their populations to 2030. Adding a source population led to a greater predicted persistence probability than did decreasing the in situ die-off, emphasizing that reversing the deleterious impacts of a disturbance may not be the most efficient management strategy. We expect our method will be useful for studying dynamics and evaluating management strategies of many species.

Introduction

Occupancy models, which consider observations of presence or absence of a species across habitat patches, have been used to test ecological hypotheses on metapopulations (Hanski 1994), species invasions (Yackulic et al. 2012), disease dynamics (Adams et al. 2010), species distributions (Gormley et al. 2011), population trends (With & King 1999), abiotic relationships (Cole & North 2014), and community‐level interactions (Welsh et al. 2006). Stochastic patch occupancy models (SPOMs) (Gyllenberg & Silvestrov 1994; Hanski 1994) are a family of occupancy models that describe transitions between occupancy states in terms of extinction and colonization. This class of models has advanced to accommodate imperfect detection, demographic dynamics, sparse data sets, and spatially explicit data. Some early SPOMs considered incidence‐function methods, in which the stationary probability of occupancy was used to infer extinction and colonization rates (Hanski 1994; ter Braak et al. 1998), based on the assumption they are time independent. To address this oversimplification, several researchers developed Markov chain models that permit time‐varying extinction and colonization probabilities dependent on values from other patches (O'Hara et al. 2002; ter Braak & Etienne 2003; Moilanen 2004) and that accommodate imperfect detection (MacKenzie et al. 2003; Johnson et al. 2013). A recent method with demographic dynamics is spatially explicit and allows imperfect detection (Sutherland et al. 2014).

Our study was motivated by the population decline of the U.S. federally threatened California red‐legged frog (Rana draytonii) and a desire to recover an R. draytonii metapopulation on Stanford University lands. The species has declined by 70% in its range (Hayes & Jennings 1988; Fisher & Shaffer 1996). On Stanford lands, R. draytonii was reported in 23 stream segments in 1997 and in 12 in 2012. Factors at multiple spatial scales may influence the populations: habitat loss (Davidson et al. 2002), predation by exotic species (Lawler et al. 1999), disease (Fisher et al. 2012), and climate change (Davidson et al. 2002).

We found existing SPOM methods imperfectly suited to our study system. Our data set had missing data in the sampling period 1997–2016, and did not include any data from before 1997, when events relevant to the decline likely took place. Many current spatially implicit methods use a Bayesian framework and allow years of missing data (Royle & Dorazio 2008; Fiske et al. 2011; Kéry & Schaub 2011; Bailey et al. 2014), but do not estimate parameters for those years. Although Markov methods can be used to infer parameters during missing years, they have not previously been designed to infer changes in parameter values during missing‐data periods (O'Hara et al. 2002; ter Braak & Etienne 2003; Risk et al. 2011).

We developed a Bayesian Markov chain method for inference under a novel SPOM that permits substantial missing data, handles imperfect detection, and explicitly models periods predating sampling, allowing temporal parameter changes owing to disturbance prior to sampling. The method enables inference of extinction and colonization rates and detection probabilities, hypothesis testing, and prediction of extinction risk. It entails a likelihood of model parameters given occupancy data and incorporates prior parameter information to obtain posterior probabilities and credible intervals. The framework uses inferred parameter distributions for probabilistic imputation of missing occupancy data and to estimate future probabilities of extinction. For the period before sampling, we compared 2 causes of decline in R. draytonii: increased in situ die‐off (e.g., due to disease or introduced predators) and source population loss (e.g., through habitat loss). We evaluated possible management actions and their impact on extinction risk. This example illustrates the potential of the new method for a variety of problems.

Methods

Data

From 1997 to 2016, up to 2 visual encounter surveys per year were completed in prespecified segments of 3 creeks known to contain California red‐legged frog (Fig. 1). The creeks are often partially dry during the summer dry season and are continuously wet during the winter wet season. The proportions of a creek that have pools, riffles, and runs vary over time and shift with flood events; the locations of these features also vary.

Figure 1

Locations of 3 creeks (Matadero, Deer, and San Francisquito) sampled for the California red‐legged frog on Stanford University lands.

Surveys focused on presence and absence of adult frogs, the life stage most closely linked to long‐term population persistence in r‐selected species such as R. draytonii (Biek et al. 2002; Vonesh & De la Cruz 2002). Summer surveys, when water levels are low, maximize detection probability of adult frogs and surveyor safety.

We designed our model for the annual life cycle of R. draytonii. Breeding and dispersal happen during the wet season. Breeding occurs in permanent or seasonal ponds or stream pools, then tadpoles metamorphose in subsequent months. Dispersal of metamorphosed frogs occurs both upstream and downstream. Mean dispersal distance is 50–500 m (Lannoo 2005), and dispersal can occur across interstitial habitats (Bulger et al. 2003). We assumed the population was closed during summer (Supporting Information); therefore, surveys that did not detect R. draytonii were regarded as detection failures if occupancy was detected in the same segment in other surveys of the same year.

Stream segments were numbered from upstream to downstream and had a mean length of 202 m (SD 41). There were 11 segments in Matadero Creek, 10 in Deer Creek, and 26 in San Francisquito Creek (Supporting Information). Missing data rates were 17% in Matadero Creek, 22% in Deer Creek, and 13% in San Francisquito Creek (Supporting Information).

General SPOM

In our model (Fig. 2a), we divided a linear habitat into N patches. The distance between patches i and j is denoted by d. In a discrete time period t, each patch is in 1 of 2 states, either occupied or unoccupied. Because our system has an annual cycle, we considered years rather than a generic time unit. We defined variables as follows: p is per‐patch detection probability of the species during survey j; J is number of surveys in patch i and year t; Y , , is detected occupancy (1, presence; 0, absence) in the jth survey of patch i in year t, and Y = (Y 1,1, , …, Y ,1, 1,2, , …, ) is the vector of observations of all patches in all surveys of year t.

Figure 2

(a) General stochastic patch occupancy model (SPOM) (t, year; z’ and z +1, hidden random variables describing occupancy states at the end of extinction and colonization phases, respectively; Y, observed occupancy in survey j of patch i in year t; p, detection probability; quantities on arrows, probabilities). (b) Hypothesis 1: In year t, in situ die‐off (e.g., due to infection) reduces the population size of each patch from K to K. (c) Hypothesis 2: In year t, a source population of size K at distance d from the creek is lost (e.g., due to anthropogenic disturbance). In (b) and (c), year considered for the initial condition is 1882. In (c), curved arrows indicate direction of flow.

In our model of sequential extinction‐colonization dynamics, we assumed occupancy measurement precedes extinction. First, patches are surveyed, providing Y . Following MacKenzie et al. (2003), survey results were assumed to depend on the detection probability p and the occupancy state z = (z 1, ,…,z , ). Conditional on presence of the species (z , = 1), the probability of observing presence in survey j of patch i in year t is p, and the probability of observing absence is 1 − p. The probability of Y , , = 1 given the species is absent (z , = 0) is 0, and the probability of Y , , = 0 given species absence is 1. We assumed independence of surveys, so Y has probability equal to the product of the probabilities of the Y , , across surveys j and patches i. We henceforth assumed p was constant across all j and equal to p for all surveys. The framework can be extended easily to allow the detection probability to vary spatially and temporally.

In the cycle corresponding to each discrete time step, after occupancy measurement, extinction occurs. The extinction phase represents the dry season in which some patches become empty. We denote by E the extinction probability of patch i (i.e., the probability the patch converts from occupied to unoccupied in the extinction phase of a discrete time unit) (Fig. 2a). Following Hanski and Ovaskainen (2000), we assumed extinction rates are inversely proportional to population size and colonized patches instantly reach carrying capacity (K) for patch i. Thus, E = e/K, where e is a global extinction parameter that is constant across patches; e but not E may exceed 1.

After the extinction phase of a time step, the colonization phase occurs, corresponding to the wet season in which patches are recolonized. We denote by C the colonization probability of patch i in year t (Fig. 2a). Following Hanski and Ovaskainen (2000), we assumed colonization rates are proportional to the total number of migrants entering from occupied patches. We also assumed an exponential dispersal kernel. Thus, the colonization rate of patch i at time t is C = c ∑ ≠ exp(‐αd)K, where c is the global colonization parameter (a constant across patches independent of i and t), α is the inverse of the mean dispersal distance of the species, and z’ is the occupancy of patch j after the extinction phase of year t. Although K values are unaffected by the extinction phase, C is indirectly affected because extinction leads to z’ = 0 for some segments, j, and thus decreases C. The c but not C may exceed 1.

Model likelihood computation is in Supporting Information. In the likelihood computation, the E and C are used to compute ɸ , the 2 × 2 transition matrix from all possible states in year t to states in year t+1, where entry ɸ represents the transition probability from state k to state l. The exact likelihood function of model parameters () given a series of detections measured from data Y 1997,…,Y 2016 is

In this equation, ɸ denotes the 1 × 2 probability vector of all 2 possible states in year t; q denotes the 2 × 1 column vector whose elements correspond to the probabilities of observing Y given each state; and D(q ) denotes the 2 × 2 diagonal matrix for which diagonal entries are the elements of q (Supporting Information). Thus, Eq. 1 matches Eq. 5 of MacKenzie et al. (2003), but our computation of ɸ is related to that of O'Hara et al. (2002) and ter Braak & Etienne (2003) instead of that from MacKenzie et al. (2003).

To reduce computation time when the number of segments N is large, we built an approximate likelihood function, (Supporting Information). The approximate likelihood function considers only the most likely occupancy states.

Hypothetical Causes of Population Decline

We built models corresponding to 2 hypothetical causes of population decline for R. draytonii: in situ die‐off and loss of a source population (Fig. 2b–c). For both models, we assumed an event changed the model parameters in the past, before the first survey that produced occupancy data. We assumed parameters were constant for a long time before the event of interest; thus, the initial occupancy state was assumed to have limited impact on occupancy at the time of the event. We also assumed parameters were constant after the event.

Under hypothesis 1 (in situ die‐off), population declines result from sudden mortality increases in all patches at time t, prior to the first sampling time. In situ die‐off could result from disease, introduced predators, or reductions in habitat quality or availability. We assumed population sizes were equal in all patches. The per‐patch population size before in situ die‐offs is labeled K and the subsequent size is denoted K. Because extinction and colonization are treated as functions of population size, we assumed in situ die offs increase population extinction rates and decrease colonization rates. See Supporting Information for derivation of the likelihood under hypothesis 1 from the general likelihood (Eq. 1).

Under hypothesis 2 (source population loss), habitat destruction or local extinction leads to loss of a neighboring source population (e.g., a pond not subject to seasonal disappearance) at time t, before the first sampling time. Stanford lands are bordered on 3 sides by urban developments that have increased in density and spatial scope over the last 50 years. We hypothesized that a source population of size K was formerly near the current habitat, at distance d, and that it became extinct at time t. We assumed this source population was simply present and not subject to extinction‐recolonization dynamics before the loss (e.g., a pond). We assumed equal population size for all patches other than the source population. The difference in hypothesis 2 compared with hypothesis 1 is that the past event is localized rather than occurring as a global event affecting all patches similarly. Under hypothesis 2, global parameters e and c do not change following source population loss. Within‐patch population sizes or extinction rates also do not change; these depend only on e and the population size. However, change occurs for the patch colonization rates—which depend on dispersal from the source population. The source‐population‐loss hypothesis encodes reduced colonization by preventing population rescue from the source population.

We computed the likelihood of the 2 hypotheses (see Supporting Information). With Θ denoting the parameters under hypothesis h, under hypothesis 1, Θ = (K), and, under hypothesis 2, Θ = (K). With t 0 denoting the initial year, the probability vector of all possible states in t 0 is . Because transition probabilities after the past event of interest do not depend on Θ , we denote the event date t or t by t. Then

Special cases of Eq. 2 corresponding to hypotheses 1 and 2 are in Supporting Information. The null hypothesis is a special case of hypothesis 1, with K = 1.

Bayesian Parameter Estimation

Estimation proceeded in 2 steps. First, using the whole occupancy data set, we inferred the detection probability p and the unknown parameters shared in both hypotheses: the mean dispersal distance 1/α and the global extinction and colonization parameters, e and c. The Bayesian estimation yielded posterior distributions of the parameters given the data set by multiplying the model likelihood (Eq. 1) by the prior distributions of model parameters (Supporting Information). We then used the mode of the joint posterior distribution to obtain maximum a posteriori estimates of , , and in each creek. We identified 95% credible intervals (CI) as the 2.5% and 97.5% quantiles of marginal posterior distributions. We used the posterior to perform missing data imputation (details in Supporting Information).

Second, we inferred the unknown parameters distinct to the 2 hypotheses: population size before infection K and infection timing t (hypothesis 1) and size of the lost habitat K, its distance to the creek d, and loss timing t (hypothesis 2). We multiplied the prior distribution of the parameters by the likelihood under the hypothesis (Eq. 2) to obtain the posterior distribution of Θ (Supporting Information). Maximum a posteriori parameter estimates and CI were computed as in the analysis used for both hypotheses. Our method provided a posterior distribution for the initial occupancy z 1882; however, this value is not of interest, and we integrated the joint posterior for the other parameters over all possible values of z 1882.

The numerical implementation of the inference method and an approximate method for the case of a large number of patches, which we used for San Francisquito Creek, are in Supporting Information. We implemented the method in software MIDASPOM (Supporting Information).

Parametrization

Model parameters and prior distributions are in Table 1 (details in Supporting Information). The mean dispersal distance 1/α was unknown and assumed equal for all creeks. The detection probability p and the global extinction and colonization parameters e and c were also unknown and differed by creek. The size of the populations under hypothesis 1 (K), source population size under hypothesis 2 (K), and timing of the hypothetical increase of in situ die‐off (t) and of the hypothetical source population loss (t) were also unknown, and we assumed events occurred between 1902 and 1982 (Padgett‐Flohr & Hopkins 2009; Fofonoff et al. 2017). Occupancies before t and t were unknown, and we assumed the species was present in 1882 (Hayes & Jennings 1988; Fellers 2005), so likelihood computation incorporated at least 20 years before t and t. The position of the hypothetical source population under hypothesis 2 was unknown. We assumed its distance d to the first creek segment was between 200 m and 4 km. Elements of the vector ɸ 1997 represent possible states z 1997, and their values provide prior probabilities of each state in 1997; we considered a Bernoulli prior.

Table 1

Summary of parameters of the stochastic patch occupancy model

Hypothesis and variable	Interpretation	Rangea	Value usedb	Reference
Parameters shared across hypothesesb
P	probability of detection	[0,1]	uniform prior	none
E	global extinction parameter	[0,1]	uniform prior	none
C	global colonization parameter	[0,1.5]	uniform prior	none
D	distance between consecutive segments (m)	[164,256] in MCc; [147,249] in DC; [152,446] in SFC	200 m	Supporting Information
1/α	mean annual dispersal distance	[50,500]	uniform prior	Bulger et al. 2003; Lannoo 2005; Fellers & Kleeman 2007
z 1997, i	occupancy state in patch i in 1997	[0,1]	Bernoulli prior	none
Hypothesis 1 In situ die‐off
tD	timing of infection or increased predation	[1902,1982]	uniform prior	Padgett‐Flohr & Hopkins 2009; Fofonoff et al. 2017
KD	per‐patch population size prior to infection	[0.1,100]	log‐uniform prior	none
Hypothesis 2 Source population loss
tL	timing of source population loss	[1902,1982]	uniform prior	Google Historical Imageryd
KL	population size of extinct population	[0.1,100]	log‐uniform prior	none
dL	distance to source population (m)	[200,4000]	uniform prior	none

Ranges correspond to the maximal and minimal values of the prior distributions detailed in Supporting Information.

Values of parameters with priors are inferred using our Bayesian method, and distributions are detailed in Supporting Information.

Abbreviations: MC, Matadero Creek; DC, Deer Creek; SFC, San Francisquito Creek.

Used to identify years when rapid expansion of human habitation in San Mateo and Santa Clara Counties around MC, DC, and SFC occurred.

We assumed all other parameters were known. Patches corresponded to creek segments, and distances d between segment pairs were assumed known from geographic data. Because the habitat is linear and distances between midpoints of consecutive segments were similar (Supporting Information), we assumed consecutive segments had fixed distance d. The relative carrying capacities K were treated as known from geographic data, and we assumed proportionality between patch area and population size (Hanski & Ovaskainen 2000). As a first approximation, we assumed K 1 = K 2 = … = K = 1 (Supporting Information).

Hypothesis Testing

To quantify fit of hypotheses to data, we performed Bayesian model selection. For hypotheses i and j, we computed Bayes factors by numerically integrating under each hypothesis the product of the model likelihood and the prior probabilities of parameters, considering all possible parameter values: where i and j equal 0 under the null hypothesis, 1 under hypothesis 1, and 2 under hypothesis 2, the likelihood was computed from Eq. 2, and the prior probabilities were computed from equations in Supporting Information. We interpreted B as follows (Jeffreys 1998): |log10(B)|<0.5, little support for either hypothesis; 0.5<|log10(B)|<1, substantial evidence for one hypothesis; and 1<|log10(B)|, strong evidence.

Extinction Probability Under Management Scenarios

We predicted future extinction probabilities in each creek under each of 4 management scenarios: no management, mortality reduction, habitat creation, and a combined scenario with both mortality reduction and habitat creation. In the mortality reduction scenario, local causes of mortality were controlled within every stream segment (e.g., non‐native predator removal) to increase carrying capacity K. This scenario is equivalent to reversing the deleterious effects of the population decline that occur under hypothesis 1. In the habitat creation scenario, a population of size K was restored at distance d from the creek. This scenario is equivalent to reversing the deleterious effects of the decline that occur with source population loss in hypothesis 2. The combined strategy both reduced local causes of mortality and created new habitat.

To implement these scenarios, we performed Monte Carlo simulations from the current state p 2016 with patch extinction (E i) and colonization (C i) probabilities computed from e and c values sampled from their joint posterior distribution (Supporting Information). For each creek, we recorded the proportion of simulations with all segments unoccupied (total extinction of the population within a creek) from 2016 to 2065. This quantity is the probability that extinction has occurred by a given year when applying the transition probabilities from the observed state in 2016 and assuming that occupancy data are missing in all future years (Supporting Information).

To account for the difficulty of implementing management plans, for each management scenario, we chose parameter values that corresponded to reversing some but not all effects of disturbance events (hypotheses 1 and 2; parameter estimates in Supporting Information). For the mortality‐reduction scenario above, this meant considering an increased local population size smaller than the estimated size before the event ( for hypothesis 1, from Fig. 4). We thus tested the impact of a small increase (K = 1 to 1.05) and a larger increase in population size (K = 1 to 1.25). For the habitat‐creation scenario, this meant considering an increased source population size smaller than the estimated size before the event ( for hypothesis 2, from Supporting Information). We thus tested the impact of the addition of a small population (K = 1) either near (200 m) or farther from (600 m) the creek.

Figure 4

Bayesian estimation of model parameters under the hypothesis of in situ die‐off of California red‐legged frogs (H1) in (a) Deer Creek and (b) San Francisquito Creek (K, population size before the event that increases in situ die‐off; gray, prior distribution of K [Supporting Information]; red, posterior conditional on the observed data [Supporting Information]; red dashed lines, mode of the posterior distribution used as a point estimate under hypothesis H1). Model parameters are defined in Table 1. Only parameter estimates under the hypotheses with the greatest support (based on the Bayes factor, see Supporting Information) are shown. Parameter estimates under the (less strongly supported) source population loss hypothesis (H2) and under hypothesis H1 for Matadero Creek are in Supporting Information.

Results

Estimation of Shared Model Parameters

We estimated the value of the mean dispersal distance at m (95% CI, 125–425) (Supporting Information). In addition, we estimated similar probabilities of detection in Matadero and Deer Creeks (Fig. 3a,c): in Matadero (95% CI, 0.69–0.82) and 0.75 in Deer (95% CI, 0.64–0.81). Values were lower in San Francisquito (Fig. 3e), (95% CI, 0.57–0.77).

Figure 3

In Matadero, Deer, and San Francisquito Creeks, Bayesian parameter estimation of (a, c, e) detection probability parameter and (b, d, f) joint distribution of extinction and colonization parameters. In (a, c, e), gray is prior distribution (Supporting Information), red is posterior distribution given observed data (from Eq. 1 for Matadero and Deer Creeks and from Supporting Information for San Francisquito Creek), and red dashed line is maximum a posteriori estimate. In (b, d, f), red is joint posterior distribution (posterior density scale 0–22) and red dashed lines are maximum a posteriori estimates (, respectively, extinction and colonization). An informative prior (Bernoulli, Supporting Information) was used for the missing data in the initial state of Deer Creek. Supporting Information contains the posterior distribution under an uninformative prior.

The extinction and colonization parameter estimates differed among creeks (Fig. 3b,d,f). In Matadero Creek, point estimates of both parameters were small, . The marginal posterior distributions were narrow, indicating the data set was informative about extinction and colonization dynamics (95% CI, 0.06–0.24 for e and 0.22–0.96 for c).

By contrast, point estimates in Deer Creek were large: and. The marginal posterior distributions of e (95% CI, 0.21–0.52) and c were also wider (95% CI, 0.75–1.87). The choice of either an informative or an uninformative prior for the missing data in the initial year (1998) had little effect on point estimates and credible intervals (Supporting Information).

In San Francisquito Creek, the extinction point estimate was the largest among all creeks, , and the colonization estimate was between those of Matadero and Deer Creeks, . The marginal posterior distributions of e and c were both moderately large (95% CI, 0.33–0.62 for e, 0.43–1.30 for c). Supporting Information provides an assessment of the accuracy of the estimation in Fig. 3.

The Bayes factor strongly supported the in situ die‐off hypothesis H1 in Deer Creek (log10 B 0,1 = −0.944, log10 B 0,2 = −0.387, log10 B 1,2 = 0.558) and San Francisquito Creek (log10 B0,1 = −44.616, log10 B 0,2 = 0.087, log10 B 1,2 = 44.680) (Supporting Information). In Matadero Creek, it did not reject the null hypothesis H1 (log10 B 0,1 = 0.058, log10 B 0,2 = −0.011, log10 B 1,2 = −0.069). Assessment of the accuracy of the hypothesis testing is in Supporting Information.

Parameter Estimation

For Deer Creek, the posterior probability of K plateaued above a threshold of 1.45. For San Francisquito Creek, it was highest for (Fig. 4). Current declines in Deer and San Francisquito Creeks, respectively, could therefore be explained by decreases of at least 31% and 98% of the population size in each segment due to increased in situ die‐off. The posterior distribution of the event timing was relatively flat for Deer Creek, indicating the data contained little information about this quantity (Supporting Information), and was largest in 1982 for San Francisquito Creek (Supporting Information).

For Matadero Creek, because H0 was not rejected, parameter estimation under H1 was not relevant (Supporting Information). Parameter estimation under H2 in the creeks is in Supporting Information.

The Deer Creek population was more likely to go extinct than the Matadero Creek population (8‐fold higher extinction probability by 2030, black line in Fig. 5a–b). The San Francisquito Creek population was probably already extirpated by 2008 (Supporting Information).

Figure 5

Cumulative probability of complete within‐creek extinction of California red‐legged frogs in 50 years computed using the maximum a posteriori estimates of extinction and colonization parameters: (a) Matadero, (b) Deer, and (c) San Francisquito Creeks. Maximum a posteriori extinction and colonization parameters are in Fig. 3.

All management interventions decreased the extinction probability in Deer and Matadero Creeks. In San Francisquito Creek, only source‐population creation restored the population and lowered the future extinction probability below 1. In the 3 creeks, source‐population creation 200 m from the creek generated the lowest extinction probability.

Discussion

We developed a method to use temporal patch occupancy data to infer extinction and colonization rates and to test hypotheses under a Bayesian framework and applied the method to 3 populations of R. draytonii. With the method, we were able to test for changes in parameter values that occurred before the first occupancy survey. It also enabled tests of hypotheses of past disturbance and estimates of future population trajectories.

A key innovation is that our method detects disturbances that predate data collection. Many spatially implicit SPOM methods (Royle & Dorazio 2008; Bailey et al. 2014) use a Bayesian framework, allow arbitrary missing data, and detect temporal parameter changes; however, they neither detect changes before the first survey nor test explicit hypotheses concerning such changes. Although the spatially explicit SPOMSIM (Moilanen 2004) and O'Hara et al.’s (2002) method permit years of missing data, they do not allow temporal parameter changes or missing data years before the first survey. The spatially explicit Bayesian method from ter Braak and Etienne (2003) allows years of missing data prior to data collection by setting a prior for the initial occupancy state, but it does not implement parameter changes prior to data collection. We extended their framework by using the whole data set to estimate current parameter values and the first year of data to estimate past values.

Past Influences on California Red‐Legged Frogs

Our method identified local variation in R. draytonii populations. Despite the proximity of 3 creeks, we found differences in extinction and colonization (Fig. 3). Local variation in parameter estimates might reflect differences in community structure, habitat quality, or disease dynamics. For example, non‐native predators are most abundant in San Francisquito Creek (A.E.L. and E.M.C. data).

We estimated the mean dispersal distance of the species at 175 m (Supporting Information). This estimate accords with movement studies in Point Reyes, California, which report a median of 185 m (Fellers & Kleeman 2007). However, this value may not be appropriate here because Point Reyes frogs breed in freshwater ponds rather than creeks.

We compared 2 past population‐decline scenarios: in situ die‐off and source‐population loss. An in situ die‐off starting between the 1960s and the early 1980s was the most likely explanation for decline in Deer and San Francisquito Creeks (Supporting Information). Several factors could contribute to in situ die‐off. Unfortunately, our model does not distinguish the effects of disease and introduced predators, phenomena that could have had similar impacts on in situ die‐off. Infection by Bd (Batrachochytrium dendrobatidis), a fungal pathogen, is widespread in amphibians in California (Padgett‐Flohr & Hopkins 2009) and worldwide (Fisher et al. 2012). Infections with Bd were observed in 1961 on the Stanford campus (Padgett‐Flohr & Hopkins 2009), but might have begun in California as early as the late 1800s (Huss et al. 2013). Bullfrogs (Lithobates catesbeianus), a non‐native predator of R. draytonii, were introduced into California around 1900 (Stebbins & McGinnis 2012), and signal crayfish (Pacifastacus leniusculus) were introduced in San Francisco watersheds as early as 1898 (Fofonoff et al. 2017). Non‐native predators present in Stanford creeks, including mosquitofish (Gambusia affinis) (Lawler et al. 1999), bullfrogs (Moyle 1973), and signal crayfish (Allan and Tennant 2000), have contributed to similar R. draytonii declines elsewhere in California. The trajectory of population decline suggests extinction debt, where local extinction occurs after substantial delay following habitat degradation (Kuussaari et al. 2009).

Predicted Dynamics and Conservation

Our method predicted future persistence probabilities for the R. draytonii populations. Interestingly, despite similar ratios of extinction and colonization rates in Matadero and Deer Creeks, stochastic dynamics produced a greater 50‐year persistence probability in Matadero Creek (Fig. 5). In Deer Creek, large estimated extinction and colonization rates magnified variability in occupancy. Thus, successive years of substantial occupancy reduction could lead to chance extinction. In 1 year, only 2 segments were occupied in Deer Creek. With extinction parameter e = 0.39, because extinction precedes colonization, the probability of complete extinction in the next year depended only on the extinction rate (0.392 = 0.1521). In Matadero Creek, however, smaller extinction and colonization parameters and consequent lower occupancy variation made successive years of substantial occupancy reduction unlikely. For example, in 3 recent instances, 3 segments were occupied (e = 0.12), and the extinction probability was a relatively safe 0.123 = 0.001728. Sensitivity to years with substantial occupancy reduction raises concerns for persistence under unpredictable climate change: a warming trend in California, where precipitation deficits and increased incidence of extreme drought (Diffenbaugh et al. 2015) have occurred, potentially threatens aquatic species (Meyer et al. 1999).

Our results suggest a preferred approach to maintenance of R. draytonii populations, predicting that establishment of a source population <200 m from Matadero and Deer Creeks would reduce extinction risk (Fig. 5). The greatest threat is stochasticity of extinction and colonization; adding a source would facilitate survival by population rescue. The Stanford Conservation Program has pursued this management approach in a similar project, creating ponds that support California tiger salamander (Ambystoma californiense) reproduction. Recovery results have also been seen in other similar projects (Petranka et al. 2007; Rannap et al. 2009). Perhaps surprisingly, among possible conservation actions, the preferred action (adding a source) did not most directly reverse the most likely past disturbance (increased in situ die‐off). Our results suggest prioritizing management actions for Deer Creek is likely to provide the most risk reduction for a given level of effort.

Extensions

Incorporating life‐stage or demographic data by modeling local population sizes as temporally varying rather than constant could enable modeling of disturbances affecting specific life stages. It may also enable modeling of phenological shifts in life‐history dynamics, such as differing dispersal by life stage. However, when data on population size are unavailable, as in our frog example, models integrating local population dynamics are not applicable (Sutherland et al. 2014). Another extension would allow patches to have different carrying capacities K, if estimates of these quantities are known. This change would simply rescale expressions for extinction and colonization rates E and C.

Our model permits nonlinear habitats and unequal patch sizes. Computations require only pairwise distances between populations and could therefore be applied to many geometries. Scenarios beyond increased in situ die‐off and source population loss can be examined. For example, reductions affecting only certain patches would reduce some population sizes K, cyclic droughts would periodically push certain population size values K i to 0, and so on. A covariate x that influences extinction or colonization can be accommodated by making the associated parameters in patch i at time t, e , and c , , functions of x , , for example with a logit link e , (x , ) = [exp(a 0+a 1 x , )]/[1+ exp(a 0+a 1 x , )]. Posterior distributions of parameters a 0 and a 1 can then be inferred.

In summary, our method enables assessments of influences on population persistence of complex scenarios, including local (e.g., adding new patches) or global actions (e.g., reducing predation or disease). It can both assess precise scenarios (e.g., specifying new habitat locations) and perform exploratory analyses (e.g., predictions across a range of possible locations). In our system, it showed that the best management action for reversing population decline was not the action that most directly reversed the original cause of decline. This observation is particularly important for management of species for which threats are hard to address and suggests that alternative strategies can improve persistence. We encourage use of our framework for advancement of ecological theory and management.

Supplementary methods (Appendices S1 to S5), parameter estimation (Appendix S6), hypothesis tests (Appendix S7), model selection based on AIC (Appendix S8), parameter estimation under exact and approximate algorithms (Appendix S9), parameter estimation of dispersal distance (Appendix S10), parameter estimation in Deer Creek with an uninformative prior for missing occupancy (Appendix S11), accuracy of parameter estimation from Fig. 3 (Appendix S12), accuracy of model comparisons (Appendix S13), estimation of timing of population reduction event (Appendix S14), estimation of model parameters under hypothesis 1 in Matadero Creek (Appendix S15) and hypothesis 2 (Appendix S16), estimated occupancy trajectories (Appendix S17), parameter estimation based on many‐patches approximate algorithm (Appendix S18), creek segment lengths (Appendix S19), segment occupancies (Appendix S20), and software MIDASPOM (Appendix S21) are available online.

Click here for additional data file.

Click here for additional data file.