Evaluating alternative distance criteria

Simulation framework. We developed a simulation framework to assess the performance of alternative distance criteria to correctly assign sightings to their respective true identities (IDs; Fig. 2). Schwartz et al. (2008) developed a computer algorithm to automate the application of the Knight et al. (1995) rule set consistent with its implementation. They then used location data of radiomarked females with cubs to simulate performance of the rule set under various hypothetical, but realistic, levels of “true” abundance of females with cubs. To accomplish the latter, they compiled a geospatial data layer with telemetry locations from multiple bear-years as if they had all been observed in a single year, and then randomly sampled from this superpopulation of observable bears. Live-trapping bears for radiomonitoring purposes is not feasible in some portions of the ecosystem; therefore, sets of known telemetry locations of females with cubs were placed on the data layer to populate areas in which few radiomarked females had been located but that were known to be occupied by adult female bears (Schwartz et al. 2008). The result was a representative distribution of bear locations for simulations to evaluate the Knight et al. (1995) rule set, with the goal of producing realistic inter-sighting distances and associated dates and times as crucial components of the rule set. They then took repeated samples (n = 500 simulations) of 10, 20, 40, 80, and 100 true females with cubs from this superpopulation to represent variability in samples obtained by chance through the sampling protocol.

Fig. 2.

Flow chart of analyses designed to evaluate and enhance techniques for monitoring the Yellowstone grizzly bear (Ursus arctos) population in the United States. See Results section for referenced figures and tables. See   Supplemental Materials for Fig. B1 (ursu-33-01-16_s01.pdf).

For this study, we built on the general approach of Schwartz et al. (2008). We used aerial telemetry locations (1 May–31 Aug; individual bears located every 10–14 days) and ground sightings (prior to 31 Aug) of radiocollared females with cubs collected annually during 1997–2019 ( Supplemental Materials, Fig. B.1 (ursu-33-01-16_s01.pdf)). This data set was more recent compared with Schwartz et al. (2008) and enabled us to evaluate potential changes over time. Following the approach of Schwartz et al. (2008), we created 1,000 simulated data sets with true population sizes of females with cubs at each of 5 plausible levels (N_true = 50, 60, 70, 80, and 90). We chose 50 as our minimum N_true to reflect the demographic recovery criterion of 48 unique females with cubs mentioned previously. We varied the total number of simulated sightings for each replicate as a ratio of N_true, based on empirical ratios of total sightings (n) and N_Chao2 estimates for the period 1997–2019 (total sightings:unique females with cubs [n/N_Chao2]; mean = 2.3, range = 1.5–3.2). For simulated sightings, we retained the empirical day, month, time, and coordinate values from the telemetry data.

For each replicate, we allowed only 1 sample-year (1 sample-year = 1 year of location data for 1 female with cubs) to be chosen for any female with multiple years of data to prevent unrealistic spatial overlap (Schwartz et al. 2008). Similarly, our location sample spanned >20 years, so spatial overlap among different individuals could occur that is unrealistic. For example, if a female died and her home range was later occupied by a different female, randomly selecting both individuals may create an unrealistic dyad for evaluation of distance criteria. Therefore, when selecting individuals, we required the activity center of a newly selected candidate female to be ≥1 km from any activity center of a female previously selected while still allowing 2 simulated females to have a high degree of spatial overlap.

Litter size is also used in the clustering algorithm to distinguish sightings of unique females. Thus, we randomly assigned litter size to the earliest sighting of each female using discrete inverse transformation sampling (Devroye 1986) of empirical litter size data for the period 1997–2019 (Haroldson et al. 2020). We then simulated changes in litter size caused by cub mortality by applying estimated daily cub survival rates (IGBST 2012) to the number of days between simulated sightings of the same female. We censored simulated sightings if complete litter loss occurred because actual counts do not include females without cubs.

In field conditions, when observed females with cubs are radiomarked and are individually identified with telemetry, this information is included in the clustering algorithm and increases algorithm accuracy because these individuals cannot be misidentified (Schwartz et al. 2008). To simulate collared females with cubs, we assigned an identifier to a proportion of the females as being radiomarked in each replicate, based on a random sample from the distribution of empirical radiomonitored females with cubs on an annual basis (1997–2019; range = 3–13/yr).

Our simulation framework was designed to directly compare a true number of sighted bears with the estimate of m, but was not designed to test the efficacy of the Chao2 adjustment (i.e., Thus, unsighted females were not simulated; therefore, inferences about N_Chao2 are based on the premise of correctly assigning f₁ and f₂ sighting frequencies (i.e., the simulated f₁ and f₂ counts correctly capturing the Chao2 adjustment; Keating et al. 2002, Cherry et al. 2007, Schwartz et al. 2008).

Analysis. To assess the accuracy of different distance criteria, we used the computer program of Schwartz et al. (2008) to cluster sightings into individuals, by varying the distance threshold from 12 km to 30 km in 2-km intervals (i.e., 10 distance criteria) and holding all other parameters the same, including setting the spatial extent to the Demographic Monitoring Area. We chose 12 km as a lower bound for the range of distance criteria, based on the original 15-km annual home-range diameter documented by Knight et al. (1995) and more recent findings from Bjornlie et al. (2014), indicating that female home ranges in areas with higher bear densities have decreased in size. This approach resulted in 50,000 output data sets, with 1,000 simulated data sets of females with cubs for each of the 10 distance criteria and the 5 levels (N_true = 50, 60, 70, 80, and 90) of females with cubs (5 population levels × 10 distance criteria × 1,000 replicates each = 50,000).

We evaluated classification performance and bias associated with the 10 distance criteria at the unique ID level and the sighting level, allowing us to examine how assignments influenced the 2 components of the Chao2 equation. Accurate counts of m influence the primary and largest component of the equation. Accurate assignment of f₁ and f₂ sighting frequencies to each female ID influences the Chao2 adjustment Chao 1989, Cherry et al. 2007). Rather than evaluate strict accuracy of assignments, which would not indicate whether false negatives or false positives were more common, we evaluated classification performance by considering precision and recall (Lever et al. 2016). At the unique ID-level, we compared the presence or absence of unique IDs of females with cubs in the simulation (true IDs) with modeled output (predicted IDs). This involved 3 distinct outcomes: 1) true positive (true ID correctly predicted); 2) false positive (ID erroneously predicted to be present when sightings of a single true ID are split into multiple IDs); and 3) false negative (failure to predict true ID because multiple IDs are combined into a single ID). The true number of unique IDs is the sum of those correctly classified (true positives), plus those that were missed (false negatives). Only when the number of false positives equals false negatives are the correct numbers of unique IDs predicted. We used the F_β score, which is an aggregative perfor mance metric of precision and recall The F_β score is bounded by 0 and 1, with higher values indicating better classification performance (Lever et al. 2016). The parameter β controls the balance of precision and recall, which we set to 1 for equal balance. At the sighting level, we assessed assignment of sightings to their respective true IDs by using a multiclass confusion matrix (Grandini et al. 2020). For multiclass classification where each observation can only be assigned to a single class label (i.e., female ID), average precision = recall = F_β score. Therefore, we used mean F_β scores to summarize classification performance at the sighting level.

To assess the likelihood of overestimation at the unique ID level, we calculated the proportion of simulations with positive bias >5 and >10% of N_true, the Chao2-adjustment value, and the resulting Chao2 estimates. Simulated data sets showed a strong positive correlation between the f₁:f₂ ratio and the Chao2 adjustment component in the Chao2 equation at all N_true levels (Spearman's rank r̄_s = 0.93, σ = 0.01). Therefore, we estimated the true and predicted adjustment component of the Chao2 equation to quantify bias at the sighting level because it is the direct application of f₁ and f₂ counts, and interpreta- tion in terms of Chao2 units is intuitive. We focused that analysis on the 3 top-performing distance criteria and, for comparison, the original 30-km criterion. We also focused on N_true levels of 60 and 70 unique females with cubs because empirical estimates of m based on the 30-km criterion and total sightings for 2001–2019, when linked to simulation results, suggest this range was most relevant to contemporary conditions ( Supplemental Materials, Appendix A (ursu-33-01-16_s01.pdf)).

Evaluating performance of generalized additive models

Simulation framework. We used another simulation framework, separate from the previous analysis, to create realistic variation in trends of annual N_Chao2 estimates (Fig. 2). Given that the 2016 Conservation Strategy specified a management objective that reflected the mean population size during the period 2002–2014 (Yellowstone Ecosystem Subcommittee 2016), we developed data sets to simulate a stable population experiencing a decline followed by a return to stability. We varied magnitude (10, 15, and 20%) and duration (5, 10, and 15 yrs) of the decline periods to reflect grizzly bear biology and the management approach. The combined duration × magnitude effect sizes correspond to constant population growth rates (λ) during decline years, ranging from 0.956 (20% decline over 5 yrs) to 0.993 (10% decline over 15 yrs;  Supplemental Materials, Table B.1 (ursu-33-01-16_s01.pdf)). We also included a null scenario of zero growth (λ = 1.0). Although we only simulated population declines to keep our analyses focused, we note that the ability to detect population increases is just as relevant to monitoring of the Yellowstone grizzly bear population. The flexibility of GAMs means that model performance would be the same and results would be equally applicable for equivalent scenarios of population increase.

To simulate annual N_Chao2 values, we added “residual-noise” to the deterministic trends shown in  Table B.1 (Supplemental Materials (ursu-33-01-16_s01.pdf)) with the goal of approximating process and sampling variance present in observed N_Chao2 estimates. We used the empirical residuals from a regression of N_Chao2 ∼ year from 2000 to 2018 data to parameterize residual noise and assumed a relatively stable true population during this period with noise equally distributed in positive and negative directions. We extracted the residuals from the regression and used their statistical properties (e.g., autocorrelation, standard deviation) to simulate an auto-regressive time series using the arima.sim() function in Program R (R Core Team 2019). We scaled simulated residuals by dividing by the deterministic value (mean value of empirical N_Chao2 estimates for 2000–2018; N_Chao2 = 55.8). This allowed us to use the same simulated time series of noise for each scenario, regardless of the deterministic values. Scaled residuals were subsequently “unscaled” by multiplying by the deterministic time series value and adding it to the deterministic value to create stochastic time series.

We simulated 1,000 replicate time series of 75 years for each of the 10 scenarios (i.e., 1 null and 9 treatment scenarios). This period length allowed for stabilization, or “burn-in,” time before the start of the decline period, and a postdecline stable period. We started declines in year 31 of the simulation. Stochastic noise of the simulated N_Chao2 values resulted in variation across simulations leading up to the decline period (i.e., sometimes higher than deterministic values, sometimes lower), which added a realistic variance component to the simulations. To account for the observed increase of the empirical N_Chao2 estimates for the Yellowstone grizzly bear population through the early 2000s (IGBST 2012), we set the first 10 years of each simulation to be an increasing linear trend, thus requiring the GAM models to make an initial “turn” from increasing to stable deterministic trends. Only contemporary trends are the target of our GAM application, so we chose a generic increase to challenge the model but were less concerned with exactly matching the empirical data for these first 10 years of the simulations. Thus, simulations started with a “burn-in” of 10 years of increase plus 20 years of stable population size, followed by a population decline lasting 5, 10, or 15 years and then a postdecline stable period of up to 40, 35, and 30 years, respectively.

Model parameterization. For each monitoring year of a simulation scenario, we fit a GAM to annual N_Chao2 estimates (N_Chao2t ):

where f is a smooth function of the covariate year from t = 1 to the current monitoring year and ε_t is a vector of error terms. We fit a similar model to 3-year moving averages of annual N_Chao2 estimates (). We fit models with the mgcv package (Wood 2004) in Program R and evaluated model performance using raw and 3-year simple moving average (x̄₃) of simulated N_Chao2 values. We chose to include the moving average based on exploratory work and previous research showing that any reduction in sampling variance would increase power to detect trends (Harris et al. 2007). Although use of simple moving averages lags the data by half the size of the sample window (e.g., 1.5 yrs in our application) and delays the onset of changes in the input signal, the reduc- tion in the signal-to-noise ratio generally outweighs this lag effect in model performance. Moving averages increase autocorrelation in the time series that could lead to overfitting if not accounted for; therefore, we modified the default GAM parameterization to protect against overfitting by upscaling the spline penalization. To provide a reasonable time-series for fitting models, we began model-fitting in simulation year 25 with 5 years of pre-impact before deterministic trends started. For each simulated monitoring year, we fit a GAM with N_Chao2 or its 3-year moving average as the response variable and year as the predictor variable. Use of fitted models followed existing monitoring protocols of using only the monitoring year, or last year of a fitted model, for interpretation and not back-correcting estimates of previous years as time advances and more data become available. To account for the presence of autocorrelation in the data and protect against overfitting, we increased the effective degrees of freedom penalty by 30% (γ = 1.3). This produced a smoother fit (Kim and Gu 2004; Wood 2006, 2017) and balanced overfitting protection while still allowing the smoother to respond nonlinearly to changes in trend. Failure to account for such dependencies in the N_Chao2 values could lead to overly complex model fitting and a greater probability of false-positive results (Simpson 2018). Following Wood (2011) and Simpson (2019), we used restricted maximum likelihood for parameter estimation. We set the smoother function to use univariate penalized cubic regression splines (Wood 2017).

Model outputs and inference. For each monitoring year (n = 50) and simulation replicate (n = 1,000), we fit GAMs and used the methods of Simpson (2018) to generate and store 1,000 posterior distributions of the smoothed N_Chao2 estimates and the first derivatives f′(year_t), or slopes, as a measure of the instantaneous rate of change ( Supplemental Materials, Appendix C (ursu-33-01-16_s01.pdf)). We used the variation of estimated slopes from the control (no decline) scenario simulations to optimize the α-level based on rates of false-positive events, defined by first derivatives being significantly different from zero.

We calculated model bias as the mean absolute error of smoothed estimates (i.e., predicted GAM values) for each monitoring year at the replicate level. We selected mean absolute error over the more conventional root mean squared error because it retains the directionality of bias (positive vs. negative). We calculated this metric as follows:

where fitted and deterministic are the median of the smoothed N_Chao2 posterior distribution and the simulated deterministic N_Chao2 values, respectively, associated with year t. We output 3 metrics associated with the GAM first derivatives (f′(year)) that might be used to interpret trend. The first was the median of the posterior distribution as the point estimate for trend. The second was the proportion of the posterior distribution that was <0, representing the probability of decline (pd). The third was “trend state,” for which we assigned a state of decline for years when (1 – α) confidence intervals of first derivatives were different from zero and a state of no decline for years when the confidence interval contained 0.

Quantifying results across replications is challenging because measures of central tendency are relevant to the annual and replicate level. At the annual level, results reflect the average dynamics for a given year but do not account for time dependency present within a time series of a single replicate. Therefore, we also considered replicate-level results and explicitly accounted for the time dependency of each entire time series (i.e., any trend among annual data within each replicate, representing how monitoring data are used in practice) by using an indicator variable for “trend state.” We defined detection of a decline event as ≥2 consecutive years in the decline state. For each scenario, we estimated the proportion of simulations where a decline event was detected, the lag (yrs) between the start of the decline and its detection, and the length of the detection event (yrs).

Evaluating alternative distance criteria

Summary of simulation data. The simulated sighting data to evaluate alternative distance criteria contained 1,139 locations of females with cubs, representing 117 unique bears ( Supplemental Materials, Fig. B.1 (ursu-33-01-16_s01.pdf)). For some individual females, we had multiple years of location data while accompanied by cubs, resulting in 154 sampling years. The number of simulated sightings per unique female varied from 1 to 36 (x̄ = 7.4, σ = 4.1). Median distance between simulated sightings within the same individual averaged 9.1 km (σ = 5.9 km), ranged from 0.2 to 37 km, and lacked evidence of directional trend during 1997–2019 (β = –0.05, P = 0.49, adjusted R² = –0.004). The median diameter of the smallest circle encompassing all simulated sightings for unique individuals was 16.4 km, with 86% of individuals' minimum diameter <30 km; no trend was evident over the period 1997–2019 (β = –0.24, P = 0.16, adjusted R² = 0.007).

The average distance from each individual's centroid location to that of their nearest neighbor's centroid decreased with increasing number of unique females (e.g., x̄₅₀ = 13 km, σ₅₀ = 8.6; x̄₉₀ = 9.2 km, σ₉₀ = 6.4; subscript represents N_true level).

The number of nearest neighbors with centroids within 30 km also increased with N_true levels, showing an 81% increase from 50 to 90 unique females (x̄₅₀ = 3.8, σ₅₀ = 2.3; x̄₉₀ = 6.9, σ₉₀ = 3.6). The spatial extent of the sampling frame (i.e., Demographic Monitoring Area) was fixed across all simulation replicates, so these patterns reflect increasing density of simulated females as N_true increases. Additional summary statistics are in  Supplemental Materials, Appendix D (ursu-33-01-16_s01.pdf).

Evaluation of alternative distance criteria. At the unique-ID level, mean F_β scores were highest at distance criteria between 14 and 18 km (Fig. 3 left panel), whereas mean bias was minimized between 12 and 16 km (Fig. 3 right panel). The top-performing distance criteria for both measures decreased with increasing N_true level, but those maximizing classification performance (Fig. 3 left panel) did not always match those minimizing bias (Fig. 3 right panel). However, the top-performing criteria were always within 1 distance increment (e.g., 14 vs. 16 km); differences in mean F_β scores were within 1% of each other and differences in mean bias were ≤5% of the N_true.

Fig. 3.

Classification performance at the unique ID-level shown by (left panel) F_β score and (right panel) predicted bias in the number of unique females (m bias) based on simulations applying varying distance criteria to the Knight et al. (1995) rule set to identify unique female grizzly bears (Ursus arctos) with cubs from sightings. For each N_true level, distance criteria range from 12 to 30 km in 2-km steps (indicated by color gradient and arranged from left to right). Each boxplot summarizes n = 1,000 simulated data sets.

At the sighting level, mean F_β scores were lower than the unique-ID level, although patterns related to distance criteria were similar (Fig. 4 left panel). Top-performing distance criteria ranged from 12 to 16 km for classification and bias, with smaller distance criteria performing better with increasing N_true. Differences between mean F_β of the top-performing distance criteria and its closest competitors were small (Fig. 4). Distance criteria showed distinct and consistent pattern across all N_true levels, with effects over the range of distance criteria (12 vs. 30 km) outweighing effects within distance criteria over the range of N_true (Fig. 4 right panel). Again, this assessment of bias in the Chao2 adjustment is reflective of the bias relative to the simulated frequencies (i.e., perfect clustering of all locations assigned to the correct unique ID), not bias associated with females not observed.

Fig. 4.

Classification performance at the sighting level based on simulations applying varying distance criteria to the Knight et al. (1995) rule set to identify unique female grizzly (Ursus arctos) bears with cubs from sighting. (left panel) F_β score and (right panel) predicted bias (predicted–simulated; expressed as no. of females with cubs) in the Chao2 adjustment of the Chao2 equation. For each N_true level, distance criteria range from 12 to 30 km in 2-km steps (indicated by color gradient and arranged from left to right). Each boxplot summarizes n = 1,000 simulated data sets based on micro-averaged F_β scores.

Correlations between the bias in m and the bias in Chao2 adjustment were moderate (r_s = 0.44–0.70) but strengthened with increasing numbers of simulated females with cubs (Fig. 5). Scatterplot patterns indicated the relationship between bias in m and the Chao2 adjustment were similar across different distance criteria within levels of N_true. However, distance criteria in the 12–16-km range best minimized bias of both m and the Chao2 adjustment (Fig. 5). The range of bias in the Chao2 adjustment for unbiased ranges of m highlights the additional challenge of simultaneously estimating the f₁ and f₂ sighting frequencies compared with only m. For example, for distances of 12 km to 16 km, even when m was predicted with reasonable accuracy (e.g., within ±2 females with cubs of the true value), although mean bias of f₁ and f₂ sighting frequencies was low, individual replicates varied from –10 to 13 for f₁ and –17 to 14 for f₂ (Fig. B.2). These biases were negatively correlated, and overestimation of f₁ generally corresponded to underestimates of f₂ and a positive bias in the adjustment component of N_Chao2, whereas underestimation of f₁ corresponded to overestimation of f₂ and a negative bias in the Chao2 adjustment (Fig. B.2). Despite the challenges of simultaneously reducing bias of m, f₁, and f₂ estimates at smaller distance criteria, improvements relative to the 30-km rule set were substantial. For example, even at the lowest N_true level of 50, where the 30-km distance criterion performed best, differences between 30 and 16 km were large: mean bias of m, f₁, and f₂ using the 16-km distance criterion were –0.6, –0.1, and –0.4, respectively, but for the 30-km criterion were –11.6, –8.7, and –5.3, respectively.

Fig. 5.

Relationships between bias (expressed as no. of unique females with cubs) in the parameter m and bias in the Chao2 adjustment (i.e., N_Chao2 – m) of the estimator based on simulations (n = 1,000 replicates for each combination of distance criterion and N_true level), applying varying distance criteria to the Knight et al. (1995) rule set to identify unique female grizzly bears (Ursus arctos) with cubs. Results are shown for distance criteria of 12, 14, and 16 km (columns) within each of 3 simulated levels of true females with cubs (N_true = 50, 70, and 90; rows). Blue contour lines represent 50th, 75th, and 90th isopleths, respectively.

The proportion of simulations with positive bias >5 or >10% of N_true indicated trade-offs among distance criteria for minimizing mean bias while also reducing risk of overestimation (Table 1). This pattern was similar for m, the Chao2 bias adjustment, and overall Chao2 estimates, but most distinct for the latter.

Table 1.

Measures of bias for top-ranking (12–16-km) and reference (30-km) distance criteria for simulations where N_true = 60 and 70 female grizzly bears (Ursus arctos) with cubs and total sightings within the empirical range of n ≤ 165. Results are based on 1,000 simulation replicates for each combination of distance criterion and N_true level using the Knight et al. (1995) rule set to identify unique females with cubs from sightings. (A) m bias (no. of unique females with cubs), and proportion of simulations >+5 and >+10% of N_true. (B) Chao2-adjustment bias, and proportion of simulations >+5 and >+10% of mean known Chao2 adjustment (simulation estimate). (C) Chao2 bias, and proportion of simulations >+5 and >+10% of mean known Chao2 (simulation estimate). The 5% adjustment of the known Chao2 was approximately 3.6 (N_true = 60) and 4.4 (N_true = 70).

Evaluating performance of generalized additive models

Model performance. Monitoring year estimates for the model under the null model scenario (i.e., no simulated decline) were unbiased, with >85% of monitoring years (n = 50,000) within 2 N_Chao2 units from the deterministic, or true, N_Chao2 value (mean absolute error = 0.029; σ = 1.33). The N_Chao2t model showed a slight positive bias under the null model scenario (mean absolute error = 0.484, σ = 1.33). These differences occurred mostly during the pre-impact phase and were associated with larger lag effects because of higher variance in the annual N_Chao2t values subsequent to the initial increase prior to stabilization. For the 9 treatment scenarios, the fitted bias varied as a function of the interaction of effect size and duration. General patterns in bias reflected the lag time required for models to distinguish declines from annual variation in N_Chao2. Smoothed estimates during the decline were positively biased and bias increased with the size and steepness of simulated declines. Similarly, during the postdecline stabilization period, there was a transition to a period of negative bias when responding to the change from decline to stabilization. Compared with the N_Chao2t models, variance reduction associated with the 3-year moving averages of the models resulted in less bias during impact phases and faster returns during the subsequent stable period with bias levels equivalent to the null model scenario (Fig. 6).

Fig. 6.

Mean absolute error (MAE) for N_Chao2t (annual; blue) and (3-yr moving average; red) fitted estimates of female grizzly bears (Ursus arctos) with cubs as a function of simulation year. Columns show magnitude of the impact (large = 20%, medium = 15%, small = 10% decline) and rows show duration of impact period (long = 15 yrs, medium = 10 yrs, short = 5 yrs). Dashed vertical lines indicate the start and end of the impact period and dashed horizontal line indicates reference bias of 0.

Trend detection. We present trend detection results only for the model because of its smaller bias compared with the model. For all scenarios, the annual probability of decline, pd, increased within the first 2–3 years of the decline period, indicating that existence of a declining trend was rapidly detected. Temporal dynamics based on median pd values closely tracked the different scenarios, with shorter durations and larger magnitudes resulting in faster shifts and larger probabilities. Peak values for annual medians ranged from 0.845 (decline = 10% over 15 yrs; λ = 0.993) to 0.999 (decline = 20% over 5 yrs; λ = 0.956). On average, median pd values reached maximum levels within 3 years of the end of the decline period for all durations of decline, except for the short duration of 5 years, which reached maximum values 1 or 2 years after completion of the decline period (Fig. 7). This pattern reflects that changes postdecline were more pronounced with shorter impact periods because it is difficult for models to fully capture these rapid dynamics.

Fig. 7.

Trend dynamics for first derivative (slope) posterior distributions of number of female grizzly bears (Ursus arctos) with cubs, for 9 treatment scenarios. Columns show magnitude of the impact (large = 20%, medium = 15%, small = 10% decline) and rows show duration of impact period (long = 15 yrs, medium = 10 yrs, short = 5 yrs). Black dashed lines indicate the start and end of the simulated decline periods. Density strips associated with each year reflect the distribution of posterior medians across replicates (n = 1,000). The width of each density strip reflects the average pd value, scaled such that pd = 1.0 is reflected by adjacent years having no space between their density strips. Color gradient indicates the proportion of simulations (n = 1,000/scenario) in a decline “state” where confidence intervals for ≥2 consecutive years do not contain zero.

At the annual level, all impact scenarios except for the most gradual decline (10% over 15 yrs; λ = 0.993) showed support for significant negative slopes when averaged across replicates (Fig. 7). Although the scenario of 10% decline over 15 years lacked power to detect slope significance, it is unlikely that trends would go undetected. For example, during the simulated decline phase, median posterior slope estimates were less than the previous year's estimates during 56% of simulation years, and the number of consecutive years under this pattern (current year's slope < previous year's slope) averaged 4.5 years. When coupled with inference on the probability of decline, pd, which averaged 0.72 during the impact period, the temporal trends provide substantial additional inference of changing conditions despite the lack of statistical significance.

Summarizing results at the replicate level using the state variable of decline versus no decline for each year of a time series, detection of simulated declines was high, with declines detected in >99.6% of replicates under the medium (15%) and large (20%) decline scenarios. For small magnitude scenarios (10%), declines were detected in 84.7% (15-yr duration) to 94.7% (5-yr duration) of replicates. The mean number of years from decline onset to year of first detection ranged from 3.7 (20% decline over 5 yrs) to 11.1 (10% decline over 15 yrs), and mean detected duration of events (i.e., consecutive years in detect state; range = 3.9–8.8 yrs) was correlated with interaction of decline duration × magnitude (Table 2). Patterns for detecting the return to stabilization postdecline were similar to decline detection and symmetrical around the peak support for decline for the 15- and 10-year decline scenarios (Fig. 7). Five-year scenarios showed slight asymmetry around the peak, with a longer and more linear return toward negligible levels. For all scenarios, rebounding trends were evident well before state transition from decline to no decline occurred, based on the relative change in pd and sustained increases in the median posterior distribution (Fig. 7).

Table 2.

Change detection metrics for 9 scenarios of decline for simulated time series of N_Chao2 estimates of female grizzly bears (Ursus arctos) with cubs, based on significance of first derivative of generalized additive models and 3-year simple moving averages (). We simulated 1,000 replicate time series for each scenario, each with a length of 75 years and with population decline starting in year 31 of the simulation. Detection of the decline event was defined as ≥2 consecutive years with first derivatives statistically different from 0. Mean and standard deviation (SD) for lag to detect reflect the number of years post–simulation decline before a detection and its variation. Mean length of detection events represent the mean number of consecutive years in each event and the mean slope estimate gives an indication of effect size.