Spatial Data Compilation

We compiled multiple spatially and temporally explicit remotely sensed data layers to evaluate the effects of environmental characteristics on sage-grouse habitat selection and demographic performance across reproductive life stages and seasons (table 1). We used data from the year each radio-marked bird was monitored for temporally varying covariates. Because animals commonly select habitat in a scale-dependent manner (Boyce, 2006), we evaluated the effects of land cover variables at multiple spatial scales. We first calculated the proportion of each land cover covariate or density of linear features within circular moving windows with radii representing the averages of minimum (167.9 m), mean (439.5 m), and maximum (1,451.7 m) daily distances typically traveled by sage-grouse (Coates and others, 2016a). Movement patterns may differ during reproductive life stages, when a female is associated with either a single nest site or a brood with limited mobility. Therefore, we included additional spatial scales for nesting and brood-rearing analyses, including a finer scale (75 m), representing the core (that is, most consistent) area used by nesting females during incubation recesses, and an intermediate scale (260 m), representing the average distance moved during incubation recesses (Dudko and others, 2019). An additional scale captured the median movement distance of a female with a brood during either the early (260 m) or late (370 m) brood-rearing period (Coates and others, 2024a). We transformed distance-based predictors using an exponential decay function where e=exp(−d/α) and α represented the mean value at all locations for a given analysis (Coates and others, 2016a), which allowed for the effect to decay with increasing distance. We scaled and centered all variables before analysis.

Analysis

Our analysis assessed both habitat selection and survival for three seasons (spring, summer, and winter) and two reproductive life stages, nesting and brood-rearing, with the latter stage further divided into early (less than or equal to 21 days post-hatch) and late (greater than 21 days post-hatch) brood-rearing periods (Blomberg and others, 2014). We defined spring as 16 March to 30 June, summer as 1 July to 15 October, and winter as 16 October to 15 March. For habitat selection analyses, we used resource selection functions (RSFs) to compare used and available points (Manly and others, 2002), and our selection model (eq. 1) took the following form:

where Telemetry locations collected from 2003 to 2019 were considered used points (Y=1). For reproductive life stage models, we sampled available points within the 90th percentile of the distance traveled from leks or an individual’s nest for nest and brood selection analyses, respectively. The 90th percentile was used as a cutoff to capture sage-grouse movement patterns while also excluding extreme values that had the potential to bias the availability sample. For seasonal models, we sampled available locations within the 99-percent kernel utilization distribution calculated using all locations for a given season. We sampled available points at a 5:1 available:used ratio to capture heterogeneity in the ecosystem for contrasting use relative to availability. Distance to lek was included as a confounder covariate in the nest and seasonal models and distance to nest was included in the brood analyses, which accounted for sage-grouse movement behavior often clustered near leks that could otherwise produce spurious associations with other environmental covariates. We used the logistic regression function defined above (eq. 1) to estimate coefficients, discarding the intercept and applying the exponential link function to generate the RSF (Johnson and others, 2006; McDonald, 2013; Northrup and others, 2013).

We used logistic-exposure models in a Bayesian framework to model survival for each reproductive life stage (Shaffer, 2004; Sinnott and others, 2022). We calculated daily survival (eq. 2) using the following form:

where We assumed the survival probability of nest or brood h over interval i followed a Bernoulli distribution, $y_{h,i}~ Bernoulli\left({\theta }_{h,i}\right)$, where y_h,i=1, if the nest or brood survived the interval and y_h,i=0, if the nest or brood failed. The probability of the nest or brood surviving the interval, ${\theta }_{h,i}$, equaled $D{S}_{h,i}^t$ or the daily survival probability (DS) for nest or brood h raised to the length (time=t) of interval i. We structured encounter histories to represent the interval between consecutive relocations of a nest or brood and included the day of entry and exit, the length of the interval, the fate of the nest or brood, and the habitat features measured at the beginning of each interval (Coates and others, 2024a).

A single “best” scale for all habitat variables rarely exists because relationships between wildlife and habitat are typically scale-dependent (Stuber and Fontaine, 2019). To address this, we used Bayesian latent indicator variable scale selection (BLISS; Stuber and others, 2017) to evaluate the most influential scale among each group of variables (see the “Category” column in table 1). BLISS estimates latent scale indicator variables with reversible-jump Markov chain Monte Carlo (MCMC) sampling to evaluate the scale with highest statistical support within a group, without issues of collinearity (Stuber and others, 2017). We grouped variables into 12 categories to represent similar features and correlated variables (|r|≥0.5), such that no variables between groups were highly correlated while co-occurring in the model (table 1). Full results of variable selection analyses are reported in appendix 1.

We fit a final model to include only the most influential variables from each group to allow for straightforward inferences. To evaluate statistical support, we calculated the probability of direction (P(|β|>0; Makowski and others, 2019) and considered probabilities ≥0.85 and ≥0.95 to represent moderate and strong evidence of effects, respectively. We also calculated relative selection strengths (RSS), which provide an estimate of the relative intensity of two locations differing by one standard deviation and assuming that both locations were equally available (Avgar and others, 2017). Values greater than one represented a positive effect on habitat selection, whereas values less than one represented a negative effect.

All models were fit using MCMC simulations with Just Another Gibbs Sampler (JAGS version 4.3.0, mcmc-jags.sourceforge.net, accessed March 2021) in the “R2Jags’ package (Su and Yajima, 2015). We included vague normal priors for random effects and their measures of error (Kéry, 2010), and to prevent overfitting, we specified Lasso (that is, Laplace) prior distributions for each habitat covariate specifying an uninformative hyperprior for the tuning parameter lambda (Park and Casella, 2008; Hooten and Hobbs, 2015). We ran each model for 30,000 iterations with a thinning factor of 5, discarded the first 20,000 samples, and made inferences based on the remaining 6,000 samples from 3 independent MCMC chains. To evaluate convergence, we assessed MCMC chain mixing visually and based on Gelman-Rubin convergence statistics (<1.1; Gelman and Hill, 2006).

Model Validation

We used independent testing data from 50 individual birds or broods randomly withheld to validate each model’s predictive ability. To validate selection models, we used cross-validation methods for RSFs (Johnson and others, 2006) and used-habitat calibration (UHC) plots (Fieberg and others, 2018). The UHC plots compare observed distributions for each habitat variable to the expected distribution based on the model (Fieberg and others, 2018) and use the fitted RSF to generate a predictive distribution of used habitat that is compared to the independent testing data for each habitat variable. To validate survival models, we generated survival predictions based on the posterior distributions from the final model. We then created predicted survival curves that were compared to the observed survival curves from the independent testing data and calculated a post-hoc Bayesian predictive P-value (Gelman and others, 2013), with values near 0 or 1 indicating poor model fit (Schmidt and others, 2010).

Selection and Survival Mapping

We mapped predicted selection and survival for each reproductive life stage and season across the Bi-State region using estimates from the final models and the spatial layers described above (Coates and others, 2024a). We used the median value of the posterior distribution and applied those in the model equations outlined above, where the matrix X represented the raster values for each covariate across each 900-square meter (m²) pixel. We used raster values from 2021 for time-varying covariates because that represented the most recent nadir in population abundance based on an integrated population model (see the “Objective 2” section). Distance to lek or nest was not included in the habitat selection maps because we were interested in the spatial distribution of areas suitable for sage-grouse based on underlying habitat characteristics, whereas distance to lek would provide information on the distribution of occupied habitat that was addressed through integration with an abundance and space use index (ASUI; see below) following Coates and others (2016a).

To create selection maps, we transformed estimates using a habitat selection index (HSI), where $HSI= \frac{w\left(x\right)}{1+w\left(x\right)}$ and indicates habitat use proportional to availability on a scale of 0‒1 (Coates and others, 2016a). We then implemented the isopleth method at used locations with cutoff values at the 50th, 25th, and 5th percentiles to categorize the continuous selection surface into four categories (non-habitat, low, moderate, and high selection; Doherty and others, 2016; O’Neil and others, 2020). For survival maps, we exponentiated values of daily survival based on the number of days in either the nesting or brood-rearing periods (nesting=38 days, early brood-rearing=21 days, late brood-rearing=28 days) to calculate cumulative survival. We then categorized the continuous survival surfaces into four categories based on the distribution of values at failed and successful points following O’Neil and others (2020).

We then calculated reproductive selection and survival indices using the continuous selection and survival surfaces for each reproductive life stage to represent the entire reproductive period. We first relativized the habitat selection indices for each life stage by dividing by the maximum value in that life stage (Coates and others, 2019), thus weighting all reproductive life stages equally, and then multiplied the three relativized surfaces together. We categorized the reproductive selection surface into four categories as described for the individual selection layers. To calculate the survival index, we did not relativize values because the outputs represented true survival probabilities. We therefore multiplied the exponentiated survival surfaces during each life stage to create the reproductive survival index, which we categorized using the same methods for individual survival maps. We did not calculate the survival index for pixels that were considered non-habitat based on the reproductive selection index.

For each reproductive life stage, we created a ranked index that represented the overlap between the categorized selection and survival maps for each life stage and the reproductive indices. Pixels that had both high selection and survival constituted the highest habitat rank and were assumed to represent productive habitat for that life stage. In contrast, pixels that had high selection but low predicted survival were assigned the lowest rank because they could represent maladaptive habitat selection, potentially contributing to ecological traps.

We also created maps of example habitat management categories as an example for managers, which represented the overlap between the reproductive selection index, reproductive source habitat, and an ASUI that was calculated from breeding sage-grouse (lek count) surveys (Coates and others, 2016a) and allowed us to delineate high-quality habitat occupied by sage-grouse. We included the ASUI to represent the distribution of habitat that was predicted to be occupied by sage-grouse after accounting for the configuration of leks, the distance to leks, and the predicted abundance at each lek over time. We combined the ASUI with the reproductive selection map and identified source habitat (high selection and high survival in any reproductive life stage) to differentiate between likely occupied habitat and potential habitat that had zero or low occupancy. We delineated four example habitat management categories: (1) priority+, assumed to represent important source habitat for each reproductive life stage supporting both high selection and high survival within high use areas based on the ASUI; (2) priority areas that had either predicted high, moderate, or low selection (but excluded non-habitat) within high use areas based on the ASUI or were source areas (that is, high selection and high survival) outside of high use areas; (3) general areas that represented either the intersection between high selection and low to no use from the ASUI or non-habitat within high use areas based on the ASUI; and (4) other areas where there was moderate selection combined with low to no use from the ASUI.

Finally, we created separate habitat layers for each population nadir (1995, 2001, 2008, and 2021) estimated from an integrated population model (see the “Objective 2” section) for each life stage and season using time-varying remotely sensed vegetation cover layers to capture changes in habitat over time (Coates and others, 2024b). Time-varying covariates included sagebrush cover, shrub cover, herbaceous cover, bare ground, annual grass cover, and cumulative burned area (table 1). We also calculated composite selection, survival, and habitat suitability indices to capture overall changes in time. To create the annual composite selection index, we first calculated seasonal composite selection layers, where the spring composite selection index was a combination of seasonal spring, nest, and early brood-rearing selection layers, the summer composite selection index included seasonal summer and late brood-rearing selection layers, and the winter composite selection index only included the seasonal winter selection layer. To calculate the seasonal composite selection indices, we first relativized each layer by dividing by its maximum value and then added the individual selection components together. We then relativized the seasonal composite selection indices, added them together, and scaled the resulting layers between 0 and 1 to calculate the annual composite selection index. To create the annual composite survival index, we followed the same procedure but only created seasonal composite survival indices for spring, which included nest and early brood-rearing survival layers, and summer, which only included the late brood-rearing survival layer. To create the annual composite habitat suitability index, we first added the seasonal composite selection and survival indices together and relativized them to create seasonal composite suitability indices for spring, summer, and winter. We then added the three seasonal composite suitability layers together and scaled the resulting layer between 0 and 1 to calculate the annual composite habitat suitability index. For each selection layer (spring, summer, winter, nest, early brood, and late brood), survival layer (nest, early brood, and late brood), and the composite indices (selection, survival, and suitability), we calculated two metrics: (1) the percentage change in area excluding the lowest category (that is, excluding non-habitat for selection layers and very low survival for survival layers) based on the categorized surface from three previous population nadirs (1995, 2001, and 2008) to the most recent nadir (2021) based on an integrated population model (see the “Objective 2” section) and (2) the percentage change in the median value of the continuous surface between the three previous population nadirs and the most recent nadir. All values are reported in appendix 2. We also calculated the percentage change in distribution between previous population nadirs and the most recent nadir, where distribution was calculated as the intersection between the annual composite selection index after excluding non-habitat and the abundance and space use index, with values reported in appendix 2.

Nest Selection and Survival

During the nesting life stage, the relative probability of selection increased with taller sagebrush and at higher elevations and decreased with greater herbaceous cover, pinyon-juniper cover, bare ground, burned area, streams, springs, saline lakes, and in areas with south-southwesterly aspects (fig. 2; table 2). Model validation indicated that the nest selection model had reasonable out-of-sample predictive capabilities (Spearman’s rank ρ=0.98, R²=0.85, β_predict=0.64), and UHC plots demonstrated that model-predicted habitat was consistent with observed nest locations.

Survival of nests increased with more burned area, shrubs, more topographic roughness, and at higher elevation (fig. 3; table 2). Although the positive effects of burned area and shrubs may appear contradictory, they were operating at different spatial scales and a positive effect of burned area at 1,451 m does not necessary imply less shrub cover within 75 m of a specific location. In addition, selection and survival were misaligned for burned areas with apparent avoidance but higher survival in burned areas. However, only 2 percent of nests were within a fire perimeter, so the moderate positive relationship with burned areas was interpreted with caution. Model validation indicated that predicted survival curves were consistent with testing data (Bayesian P-value=0.55).

Early Brood Selection and Survival

During the early brood-rearing period, the relative probability of selection increased with taller sagebrush, more bare ground, and higher elevations, but decreased with more wet meadows, pinyon-juniper cover, more cumulative burned area, more perennial streams, greater topographic roughness, and higher heat load index (fig. 2; table 3). Model validation indicated that the early brood selection model had strong out-of-sample predictive capabilities (Spearman’s rank ρ=0.99, R²=0.96, β_predict=1.00), and UHC plots showed that model-predicted habitat was consistent with observed early brood locations.

During the early brood-rearing phase, survival was predicted to increase with more burned area and closer to conifer cover and decrease with more shrubs and more perennial streams (fig. 3; table 3). Model validation indicated that predicted survival curves were consistent with testing data (Bayesian P-value=0.55). Selection and survival only differed for shrubs, burned area, and conifer cover, with moderate or strong effects. Early broods tended to select areas with taller shrubs, but experienced higher survival in areas with lower shrub density. Early broods also tended to select areas with less burned area closer to conifer cover, which was associated with lower survival, although there was high uncertainty around the effect of conifer cover on early brood survival.

Late Brood Selection and Survival

During the late brood-rearing period, the relative probability of selection was higher with taller sagebrush, more herbaceous cover, and at high elevations; however, it was lower with more pinyon-juniper cover, bare ground, cumulative burned area, perennial streams, saline lakes, steeper slopes, and south-southwesterly aspects (fig. 2; table 4). Model validation indicated that our late brood selection model had reasonable out-of-sample predictive capabilities (Spearman’s rank ρ=0.95, R²=0.62, β_predict=0.60), and UHC plots demonstrated that model-predicted habitat was consistent with observed late brood locations.

Survival during the late brood-rearing period increased with more intermittent streams, more springs, steeper slopes, and higher values of the compound topographic index (fig. 3; table 4). Model validation indicated that predicted survival curves were consistent with testing data (Bayesian P-value=0.55).

Selection and survival differed for two variables with moderate or strong effects during the late brood-rearing period: streams and slope. Selection was lower but survival higher in places with more intermittent or perennial streams. Selection was also higher but survival lower with steeper slopes.

Spring Selection

The relative probability of selection during the spring increased with more sagebrush cover, annual grass cover, and at higher elevations (fig. 4; table 5). The relative probability of selection decreased with greater pinyon-juniper cover, more wet meadows, more cumulative burned area, more springs and streams, greater topographic roughness, more south-southwesterly aspects, and more saline lakes (fig. 4; table 5). Our spring selection model had reasonable out-of-sample predictive capabilities (Spearman’s rank ρ=0.99, R²=0.85, β_predict=0.79), and UHC plots demonstrated that model-predicted habitat was consistent with observed spring locations.

Summer Selection

The relative probability of selection during the summer increased with more sagebrush cover, more wet meadows, at higher elevations, and closer to perennial streams (fig. 4; table 5). During the summer, the relative probability of selection decreased with more annual grass cover, more pinyon-juniper cover, more cumulative burned area, more bare ground, higher heat load index, greater topographic roughness, and more saline lakes (fig. 4; table 5). Our summer selection model had reasonable out-of-sample predictive capabilities (Spearman’s rank ρ=0.98, R²=0.85, β_predict=0.83), and UHC plots demonstrated that model-predicted habitat was consistent with observed summer locations.

Winter Selection

The relative probability of selection during the winter increased with greater sagebrush height, more bare ground, more intermittent streams, higher elevations, and greater proximity to wet meadows (fig. 4; table 5). During the winter, the relative probability of selection decreased with more annual grass cover, more pinyon-juniper cover, more cumulative burned area, greater topographic roughness, more south-southwesterly aspects, more saline lakes, and greater proximity to springs (fig. 4; table 5). Our winter selection model had reasonable out-of-sample predictive capabilities (Spearman’s rank ρ=0.98, R²=0.84, β_predict=0.86), and UHC plots demonstrated that model-predicted habitat was consistent with observed winter locations.

Mapping

We mapped habitat selection (fig. 5) and survival (fig. 6) for reproductive life stages and three seasons (fig. 7), and we created combined selection-survival maps for the reproductive life stages (fig. 8). For the reproductive period, we calculated composite selection and survival indices across the Bi-State region by combining the selection and survival maps for the three reproductive life stages. Based on the composite reproductive selection index, 5.3 percent of the Bi-State region was classified as high selection (844 km²), 15.2 percent as moderate (2,419 km²), 41.2 percent as low selection (6,536 km²), and 38.3 percent as non-habitat (6,325 km²; fig. 5). Of the area classified as habitat, a total of 17.2 percent of the Bi-State region was classified as high survival (1,684 km²), 19.9 percent as moderate (1,947 km²), 14.8 percent as low (1,445 km²), and 48.1 percent as very low (4,621 km²; fig. 6). Based on the ranked index that combined selection and survival maps during reproductive life stages, productive habitat (high or moderate selection paired with high or moderate survival) represented 14.1 percent (1,375 km²), whereas areas with potentially maladaptive selection (high or moderate selection paired with low or extremely low survival) represented 19.3 percent of the Bi-State DPS (1,888 km²; fig. 8). We also delineated example habitat management categories as an example for managers using the reproductive selection map and the ASUI, with priority+ representing 3.3 percent (559 km²), priority areas representing 16.8 percent (2,829 km²), general areas representing 7.4 percent (1,243 km²), and other areas representing 6.1 percent (1,033 km²) of the Bi-State region (fig. 9).

Data Compilation

We cleaned and filtered lek count data following O’Donnell and others (2021) and Coates and others (2021), restricting analyses to leks for which there were sufficient data to reliably estimate trends. Sage-grouse populations demonstrate cyclic patterns in abundance, and estimates of population change can be sensitive to where the start and end points occur relative to population cycles (Garton and others, 2011; Edmunds and others, 2018; Coates and others, 2021). Therefore, we used lek counts spanning 1960–2023 to identify nadirs in population abundance. Analyses were necessarily restricted to data from known leks and unknown leks may occur throughout the Bi-State DPS, particularly in the White Mountains which are difficult to access reliably in the spring. However, we accounted for unknown leks by applying a correction factor for lek detection when calculating total abundance in the analyses described below.

Following Coates and others (2021) and O’Donnell and others (2022), we assigned all modeled leks and radio-marked sage-grouse to a subpopulation based on their location and capture location, respectively. Subpopulations within the Bi-State DPS were developed using spatial clustering algorithms delineating the geographic structure of sage-grouse distributions based on the spatial arrangement of leks and barriers to connectivity (O’Donnell and others, 2022). Subpopulations were assumed to represent closed population units, wherein immigration and emigration were negligible and changes in abundance were more likely to be driven by changes in demographic rates than migration. The Bi-State DPS included 11 subpopulations, and we had sufficient data collected from leks and individually marked birds within 9 of those subpopulations (fig. 1).

Methods: Population Modeling and Estimated Abundance

To model trends in population abundance, we used two models that spanned different time frames to align with and make best use of available data sources. Because reliable population counts were available starting in 1960, we used a population state-space model (SSM; Royle and Dorazio, 2009; Kéry and Schaub, 2012) to estimate the population rate of change ($\widehat{\lambda }$) over six distinct nadir-to-nadir cycles: long-term (1969–2019), medium/long-term (1978–2019), medium-term (1983–2019), short/medium-term (1995–2019), short-term (2002–19), and recent-term (2008–19) following Coates and others (2023). Consistent demographic data from marked birds was available starting in 2003 and we therefore also used an integrated population model, which combines multiple sources of data and can provide increased precision over models that only consider population counts (Schaub and Abadi, 2011; Kéry and Schaub, 2012; Coates and others, 2019), to estimate $\widehat{\lambda }$ over the two most recent nadir-to-nadir cycles. Nadir-to-nadir cycles differed slightly for the nadirs estimated using the integrated population model (IPM): short-term (2001–21) and recent-term (2008–21).

We used SSMs (Royle and Dorazio, 2009; Kéry and Schaub, 2012) within a Bayesian hierarchical framework to estimate $\widehat{\lambda }$ at each lek from 1969 to 2019. State-space models separate variance among the observation process and the state process, thus accounting for observation error (such as, imperfect detection). In addition, SSMs are well-suited for modeling time series data (Kéry and Schaub, 2012; Aeberhard and others, 2018) and have previously been used to model range-wide population trends of sage-grouse (Coates and others, 2021). We modeled the number of males per lek (l) and time (t) as a function of the number of males the prior year and the intrinsic rate of change (r) using the following equations (eqs. 3–10):

where The observation process (eq. 11) in our model represented the relationship between the observed data y_l,t, and the latent state estimates (N_l,t), where observed data were assumed to arise from the following count distribution: where

Estimated abundance (${\widehat{N}}_{i,t}$) for subpopulation i in time t was calculated as the sum of $\widehat{N}$ at all leks within subpopulation i in time t after correcting for lek detection (probability of 0.95), average detection bias in sightability (probability of 0.86; Coates and others, 2019) and average fluctuations in attendance by males (probability of 0.85; Wann and others, 2019). In addition, we included a correction on lek counts to estimate the number of females in the population using sex ratios ($\text{Normal}\left[2.04,11.11\right]$; Braun and others, 2015), which ultimately provided an estimate of the total abundance per subpopulation (Coates and others, 2024b). We also summed the annual abundance estimates across subpopulations to calculate an annual estimate of $\widehat{N}$ across the Bi-State DPS and estimate regional $\widehat{\lambda }$. We summarized all posterior distributions for parameters using the median and 95-percent CRI. We used the “rjags” package (Plummer, 2016) to fit the SSM. We ran the model on three chains of 100,000 iterations each, discarded the first 50,000 iterations from each chain as burn-in, and used a thinning rate set to 10, retaining 15,000 posterior samples for all parameters monitored. We evaluated chain convergence based on a potential scale reduction factor value less than 1.2 (Gelman and Rubin, 1992; Brooks and Gelman, 1998).

For nadir-to-nadir cycles for which we had consistent demographic data, we used an IPM within a Bayesian hierarchical framework, which combined multiple sources of data, including both lek counts and annual demographic data from individually marked birds, to estimate $\widehat{\lambda }$ (Schaub and Abadi, 2011; Kéry and Schaub, 2012; Coates and others, 2019). Integrated population models have previously been used in the Bi-State DPS to evaluate population trends and the effects of climate on population processes (Coates and others, 2018, 2019). Within the IPM, the state process was structured as a stochastic demographic matrix model, which included two age classes (yearling and adult) and subcomponent models for individual vital rates (for example, annual survival and recruitment). Lek counts informed the observation process of the IPM and we generated estimates of the finite rate of population growth using a joint likelihood reflecting annual survival and recruitment. Our modeling framework also allowed for the inclusion of random effects for subpopulation and year, which further refined and facilitated variability in the demographic process models over space and time.

We separated the demographic process into annual survival and fecundity, which in turn were divided into multiple subcomponents (Coates and others, 2019). We assumed that migration among subpopulations was negligible (O’Donnell and others, 2022), but sage-grouse frequently move among leks within a subpopulation. Therefore, we assumed that leks within the same subpopulation shared similar demographic rates.

We estimated adult and yearling survival using frailty models (Halstead and others, 2012), wherein we assumed a constant hazard rate and used discrete monthly intervals created from live-dead encounter histories. Individuals with unknown final fates were right-censored, and sage-grouse were graduated to the next age class if they were alive through the month of March of the subsequent year. Monthly unit hazards (UH, eq. 12) were specified using the following model structure:

where

Following the modeling process, we derived the annual survival parameter (s; eq. 18) specific to age class, subpopulation, and year using the following equations:

where

Inferences of survival models were made from VHF-marked sage-grouse due to reduced survival probabilities of GPS-marked sage-grouse (Severson and others, 2019). Juvenile sage-grouse (post-fledging, greater than 35 days old and less than 1 year old) were radio-marked in the fall, and we modeled juvenile survival separately following the same approach for adults and yearlings, less the age effect.

Telemetry data were insufficient to estimate parameters for nesting propensity, so we used estimates from Taylor and others (2012) of 0.96 (95-percent confidence interval [CI]: 0.94,0.97) and 0.89 (95-percent CI: 0.87,0.91) as informative priors for adults and yearlings, respectively. To estimate clutch size of first (cl1) and second (cl2) nests, we used a log-linear model. We assumed clutch size arose from a Poisson distribution (eq. 20) taking the form:

where

Thus, the log of the expected clutch size (μ_cl) was a linear combination of nest attempt (β_attempt), age (β_age), subpopulation (α_i), year (γ_j), and subpopulation-by-year (${\zeta }_{ij}$) effects, with intercept ${\beta }_0$ (eqs. 20–27). Subscript c represented clutch number (that is, first or second), and all other subscripts followed the adult survival equations above (eqs. 12–17).

Nest survival for first (ns1) and second (ns2) nests was modeled using the frailty models described above for adult survival (eqs. 12–19) at discrete daily intervals (t=1,2, …, T) to estimate cumulative survival (T=38) across both the laying and incubation phases. We included a fixed effect for nest attempt and an additional random effect for individual hen to account for individuals with multiple nests within or across seasons.

Egg hatchability (hatch) was modeled with data from successful nests (one or more egg hatched) using a binomial distribution (eq. 28) that took the following form:

where

Through the logit-link, p_hatch,ahij was modeled as a linear function of age (${\beta }_{age}$), subpopulation ($a_i$), year (${\gamma }_j$), and subpopulation-by-year (${\zeta }_{ij}$) effects, with intercept ${\beta }_0$ (eqs. 29–34).

The probability of a second nest attempt (np2) was modeled using a Binomial distribution (eq. 35) and took the form:

where In this model, y_np2,aij represents the number of renests and p_np2 (probability of renesting) was modeled through the logit link as a linear function of age (${\beta }_{age}$), subpopulation ($a_i$), year (${\gamma }_j$), and subpopulation-by-year (${\zeta }_{ij}$) effects, with intercept ${\beta }_0$ (eqs. 36–41).

We derived chick survival (cs) probabilities from two brood counts with time intervals that varied across the study period. The maximum number of days that broods were observed was either 50 days (2003–05, 2012–19), 35 days (2007–09, 2021–22), or 28 days (2010–11). To account for these differences, we applied a correction factor to the numerator in the inverse-logit function based on the number of days d in the observation period for year j ($\frac{d_j}{35}$). We then modeled chick survival from the brood count data (eq. 42) using the following equation:

where

Parameter terms and their definitions (eqs. 43–48) follow the same structure as for second nest propensity. The initial brood size at hatch served as the number of trials and the number of chicks that survived to the end of the observation period represented the number of successes.

We defined recruitment as the number of females added to the population on an annual time step from reproduction. For each age class a, year j, and subpopulation i, we multiplied the probabilities of all the demographic subcomponent models and multiplied that result by 0.5, assuming a balanced male:female ratio to estimate the number of females recruited into the population the following breeding season. Thus, recruitment (R_aij; eq. 49) was specified using the following equation:

where

To restore a critically small subpopulation (A-006), both male and female sage-grouse were translocated from 2017 to 2022, with all translocated grouse explicitly accounted for in the IPM. We estimated separate retention rates for translocated males and females by calculating the proportion of individuals that were present at the release site greater than 30 days after translocation. Within the IPM, males and females were removed from the capture subpopulation and added to the release subpopulation, with an associated probability of retention, in the state process functions (Kéry and Schaub, 2012). Translocated males were only included in subsequent lek counts if they remained and joined the subpopulation.

Within the IPM, an underlying SSM (Kéry and Schaub, 2012) estimated observation error using the observed lek count time series data (see the “Field Methods” section above) coupled with latent (estimated) total population size and annual population rates of change ($\widehat{\lambda }$). We applied the same corrections for sightability, lek attendance, and sex ratios described for the SSM above, which ultimately provided an estimate of the total abundance per subpopulation. Lek counts (y) informed latent abundance estimates (eqs. 50–51) using the following equation:

where

We then used joint likelihoods, where lek counts informed the observation process and demographic data informed the state process (Kéry and Schaub, 2012; Coates and others, 201832). Including random effects allowed us to estimate abundance ($\widehat{N}$) through time and across subpopulations, which were used to derive the estimated population rate of change, ${\widehat{\lambda }}_{ij}$, (eq. 52) using the following form:

where

Similar to the SSM described above, we summed the annual abundance estimates across the subpopulations to calculate an annual estimate of $\widehat{N}$ across the Bi-State DPS and estimate regional $\widehat{\lambda }$. We summarized all posterior distributions for parameters using the median and 95-percent CRI. We used the “rjags” package (Plummer, 2016) to fit the IPM. We ran the model on three chains of 100,000 iterations each, discarded the first 50,000 iterations from each chain as burn-in, and used a thinning rate set to 10, retaining 15,000 posterior samples for all parameters monitored. We evaluated chain convergence based on a potential scale reduction factor value less than 1.2 (Gelman and Rubin, 1992; Brooks and Gelman, 1998).

Distributional Areas and Modeling Changes in Sage-Grouse Distribution

Using standardized kernel point density models (Doherty and others, 2016), we created a probability density function across the Bi-State DPS using lek locations. Kernel estimators are appropriate for this type of analysis because they are non-parametric analyses that allow probabilistic inference of animal utilization based on density of geographic point locations. Although location data from individual animals are typically used to approximate probability density functions, point locations of breeding sites, such as leks, can provide an index of population-level density and distribution when abundance data are associated with each lek (Coates and others, 2016a; Doherty and others, 2016). Because we were interested in broad-scale population distributions, we estimated two probability density functions using separate bandwidths (Gitzen and others, 2006), specifically 6.4 and 10.6 kilometers (km), which represented estimates of breeding season and year-round population-level distribution patterns of radio-marked sage-grouse relative to lek sites within the Bi-State DPS (Coates and others, 2013). We compiled median $\widehat{N}$ for each lek in each year using estimates from the SSM described above (Coates and others, 2024b). We weighted each lek location by its corresponding median $\widehat{N}$ for each year, allowing for variation in spatial distribution based on variation in abundance among leks over time. The two probability density surfaces were then standardized and averaged for each pixel to account for both breeding areas and other seasonally used areas.

To estimate boundaries of core and overall sage-grouse distribution, we calculated 50-percent and 99-percent isopleths of the probability density function, respectively. We then removed non-habitat areas within these isopleths by intersecting them with areas of sage-grouse habitat defined in the “Objective 1” section to minimize inclusion of landcover types rarely used by sage-grouse (for example, salt flats, open water). We used the annual composite selection layer to represent habitat and considered non-habitat to be anything below the 5th percentile of values at all used locations (see the “Objective 1” section for full details). We calculated separate habitat layers for each population nadir estimated using the SSM to capture changes in habitat across time back to 1995, which was the earliest nadir for which time-varying remotely sensed layers were available. The resulting distributional layers created by intersecting the isopleths with the habitat layer represent occupied sage-grouse habitat during each time period and were scaled from 0 to 1. We then calculated two metrics of spatial distribution, the total area and proportional volume of the distributional raster, within each subpopulation, as well as the Bi-State DPS as a whole, in each year. We calculated area as the number of pixels in the resulting raster layers multiplied by the size of the pixel, and we calculated proportional volume for each year by dividing the raster representing area by the sum of all the pixels in that same raster, which incorporates the height of the distributional layer in addition to the spatial area. To validate the distributional areas, we calculated the proportion of modeled leks (94.8 percent) that fell within the 95-percent isopleth, indicating that sage-grouse space use patterns in this region continued to meet the assumptions used to generate distributional areas.

To evaluate changes in sage-grouse distribution across the study period (1995–2019) for each subpopulation and the Bi-State DPS as a whole, we used linear mixed models in a Bayesian framework. We analyzed subpopulation-specific trends in total occupied area on a common scale by standardizing the annual area for each subpopulation around its mean and standard deviation for the series. This allowed trend estimates to be comparable across subpopulations. Proportional volume did not require standardization because it was bounded between zero and one. By including both random intercepts, ${\alpha }_i$, and slopes, ${\beta }_i$, for each subpopulation i, we allowed trends over time to vary across subpopulations. We estimated the trend using the marginal distribution of the time effect for each subpopulation, where a negative or positive sign would indicate a declining or increasing trend, respectively. We assumed random intercepts and slopes arose from Normal distributions with means of ${\mu }_{\alpha }$ and ${\mu }_{\beta }$ and variances of ${\sigma }_{\alpha }^2$ and ${\sigma }_{\beta }^2$ corresponding to subpopulation intercept and slope, respectively. For inferences at the Bi-State DPS-wide scale, we used the hyperparameters of the random slopes and we used differences in distribution between the initial and final years to estimate net gain and loss for each subpopulation and the Bi-State DPS as a whole. We report estimates of net gain/loss (km²), rate of gain/loss (${\beta }_i$), and probability (P) of $\left|{\beta }_i\right|\text{>0}$ for each subpopulation and the Bi-State DPS in each time period based on predicted posterior distributions. To evaluate if distributional trends deviated from abundance trends, we compared trends in both area and proportional volume for each subpopulation to trends in abundance based on estimates from the SSM (see the “Methods: Population Modeling and Estimated Abundance” section).

Results: Population Modeling and Estimated Abundance

We compiled 1,720 counts from 58 leks across the Bi-State DPS to inform the observation processes of the SSM and IPM. To inform the state process of the IPM, we compiled individual-based life history data from 698 marked sage-grouse. Sample sizes varied for the estimation of specific vital rates (survival: n=698, nest survival: n=715, clutch size: n=494, hatchability: n=361, second nest propensity: n=283, chick survival: n=349).

Using an SSM, we evaluated changes in population abundance over six different temporal scales encompassing six complete population cycles: long-term (1969–2019), medium/long-term (1978–2019), medium-term (1983–2019), short/medium-term (1995–2019), short-term (2002–19), and recent-term (2008–19). Based on estimates from the state-space model the Bi-State DPS demonstrated evidence of declines of ~1.2–2.5 percent annually over the long-term ($\widehat{\lambda }$=0.979, 95-percent CRI: 0.968, 0.987), medium/long-term ($\widehat{\lambda }$=0.977, 95-percent CRI: 0.963, 0.987), medium-term ($\widehat{\lambda }$=0.985, 95-percent CRI: 0.969, 0.995), short/medium-term ($\widehat{\lambda }$=0.987, 95-percent CRI: 0.970, 0.999), short-term ($\widehat{\lambda }$=0.975, 95-percent CRI: 0.963, 0.985), and recent-term ($\widehat{\lambda }$=0.988, 95-percent CRI: 0.973, 1.001), which translated to estimated declines of 65.9 percent since the nadir in 1969, 62.0 percent since the nadir in 1978, 42.7 percent since the nadir in 1983, 26.7 percent since the nadir in 1995, 34.5 percent since the nadir in 2002, and 12.9 percent since the nadir in 2008 (table 6). Estimated trends over the short-term ($\widehat{\lambda }$=0.980, 95-percent CRI: 0.971, 0.990) and recent-term ($\widehat{\lambda }$=0.985, 95-percent CRI: 0.973, 0.997) from the integrated population model were similar to estimates from the SSM above and previous studies (fig. 10; Coates and others, 2021), although there were slight differences in the estimated nadir years and credible intervals were smaller due to the inclusion of demographic data. After accounting for variation in male lek attendance, sightability, and sex ratios, estimated population abundance across the Bi-State DPS from the SSM was 7,480 (95-percent CRI: 4707,13168) in 1969, 6,707 (95-percent CRI: 4214,12502) in 1978, 4,456 (95-percent CRI: 2841,8258) in 1983, 3,482 (95-percent CRI: 2414,5427) in 1995, 3,898 (95-percent CRI: 2955,5232) in 2002, 2,929 (95-percent CRI: 2261,3742) in 2008, and 2,551 (95-percent CRI: 1984,3194) in 2019.

Most of the Bi-State DPS was concentrated in three subpopulations (A-003, A-004, and A-007), considered core subpopulations, and the remaining subpopulations were considered peripheral subpopulations (fig. 1B). Trends in abundance based on the SSM varied by subpopulation (table 6). However, only one core (A-003) and one peripheral subpopulation (A-002) had positive population growth over the medium- and short/medium-term, and one core (A-004) and one peripheral population (A-005) had positive population growth over the recent-term (table 6). All other subpopulations demonstrated evidence of negative population growth across all time periods, although credible intervals did overlap one for some subpopulations in some time periods (table 6). However, some subpopulations, particularly the subpopulation in the White Mountains (A-001), likely contain unknown leks and are difficult to access during the lekking season, which could affect estimates if trends at unknown leks or at leks lacking sufficient data due to access differ from trends at modeled leks. Declines were substantial in some subpopulations, with 5 subpopulations declining more than 50 percent and 3 peripheral subpopulations (A-001, A-006, and A-011) declining more than 75 percent since 1995.

Distributional Areas and Changes in Sage-Grouse Distribution

The distribution of occupied sage-grouse habitat declined in the Bi-State DPS (fig. 11) and shifted among subpopulations over 25 years and 3 population cycles. Whereas the area of two core subpopulations (A-003 and A-004) and one peripheral subpopulation (A-002) and the volume of one core subpopulation (A-003) increased, the area and volume of six and eight subpopulations, respectively, decreased, with the evidence for increasing or decreasing area and volume strong across most subpopulations (tables 7, 8; figs. 12, 13). Annual rates of loss of area and volume were most substantial in one core subpopulation (A-007; figs. 12, 13; tables 7, 8) and three peripheral subpopulations (A-005, A-006, and A-009). Gains in three subpopulations (fig. 12) helped offset losses in the other subpopulations over time, but population-level trends for the Bi-State DPS were negative for both area (β= −0.02; 95-percent CRI: −0.08, 0.04) and volume (β= −0.03; 95-percent CRI: −0.06, 0.002) with 77- and 97-percent probabilities of contraction in area and volume, respectively. Overall, even though there was not strong evidence for a decline in area due to variation in trends among subpopulations, we estimated the Bi-State DPS lost ~156 km² of occupied area between 1995 and 2019.

The core distributional area (that is, the 50-percent isopleth) was comprised of three subpopulations (A-003, A-004, and A-007) over the 25-year period. During that time, two subpopulations (A-004 and A-007) contracted whereas the third (A-003) expanded (table 9; fig. 14). There was no significant change in the overall core distribution (β= −0.02; 95-percent CRI: −0.91, 0.96, P(|β|>0)=0.57), although there was high uncertainty.

The direction of trends in distribution and abundance were consistent for five subpopulations but diverged for four subpopulations (fig. 12). Increases in area and volume were typically greater in populations that were also experiencing increases in abundances, whereas decreases in population abundance were greater or similar in populations experiencing decreases in area and volume. Overall, this indicates important variation among subpopulations, with increases in both population abundance and distribution in expanding subpopulations generally being insufficient to offset decreases in subpopulations declining in both abundance and distribution.

Conservation Efforts Database

The USFWS, USGS, and the Great Northern Landscape Conservation Cooperative (GNLCC) developed the CED (U.S. Geological Survey, 2023) to track conservation efforts for sage-grouse across broad regions managed by multiple entities (public and private). To that end, the CED includes conservation efforts carried out from 1960 to present (2024) across the range of sage-grouse and is continuously updated as additional conservation efforts are implemented. Advantageously, the CED facilitates various agencies or organizations to enter data records of conservation actions, as well as supporting documents or reports, to create a centralized, comprehensive database of conservation efforts undertaken within the sagebrush ecosystem.

We queried the CED and generated summaries of all conservation actions and treatments completed on public and private lands within the Bi-State region between 2012 and 2019. These years were chosen because they aligned with years of robust lek count data collection, which included periods before and after efforts were implemented. We recorded information about the type of conservation effort, the year the effort was implemented, the size of the treatment area, and method of treatment if applicable (for example, mechanical versus chemical removal of invasive weeds). We considered the following conservation efforts in the Bi-State region and recorded in the CED: sagebrush restoration, fence modification (fence removal, fence letdown, or fence marking), road closure, conifer removal, invasive weed reduction, and wet meadow restoration (fig. 1C; table 10). We excluded some conservation efforts (CEs) that are included in the CED (for example, energy, fire, and herd management), due to inadequate sample size or lack of details about the specific treatment performed.

Modeling Rates of Population Change in Abundance

We obtained Bi-State DPS lek count data spanning 1990–2021 following the same procedures described for the “Objective 2” section. We retained observations that were recorded during the season of peak lek attendance (March 1 to May 31; Emmons and Braun, 1984; Wann and others, 2019), carried out within 30 minutes before and 90 minutes after local sunrise (Jenni and Hartzler, 1978) and were collected using one of the following standardized survey methods: ground-based lek counts, aerial helicopter-based counts, aerial fixed-wing aircraft-based counts, aerial helicopter-mounted HD/IR camera counts, aerial fixed-wing aircraft-mounted HD/IR camera counts, unknown aerial based counts, ground-based lek count routes, and unknown ground-based counts. We assumed counts with missing date, time, or survey methods were done within daily and seasonal windows and followed an approved survey method. We only included the maximum count for each year when multiple, within-year counts were done. Finally, we did not include leks that were considered active (greater than or equal to two males present/visit) in fewer than 5 years during the inferential period (1990–2021), based on findings from Coates and others (2020) that sporadic time series data can generate biased estimates of the intrinsic rate of population change ($\widehat{r}$).

We used the software package R (R Core Team, 2022) with the JAGS library utilizing the “rjags” package (Plummer, 2016) to fit a SSM on the filtered lek count data. State-space models produce estimates of $\widehat{r}$ that represent trends in the sampled population when count data are collected with imperfect, but relatively constant detection rates across time (Kéry and Schaub, 2012; Monroe and others, 2019). Additionally, SSMs consist of a set of linked likelihoods, which allows for parameters to be fit according to their effect on process (that is, environmental) variance or observation error (Kéry and Schaub, 2012). The process component of our model (eq. 53) incorporated spatiotemporal relationships that linked population dynamics across fine (individual lek), intermediate (neighborhood; a cluster of leks), and coarse (Bi-State DPS) spatial scales (O’Donnell and others, 2019; Coates and others, 202130):

where

Here ${\widehat{N}}_{l,t}$ represents the latent state (that is, population size) of lek l in year t, and the intrinsic rate of change in abundance (${\widehat{r}}_{l,t}$), describes changes in log-population size from one year to the next (that is, from t to t+1). We used an exponential growth model, which uses a log-scale for abundance and rate of change parameters (Lande and others, 2003) because of the stochasticity of sage-grouse population dynamics (Moynahan and others, 2007; Fedy and Aldridge, 2011; Fedy and Doherty, 2011; Dinkins and others, 2016). We modeled ${\widehat{r}}_{l,t}$ as realizations of a normally distributed random process where population dynamics were propagated across progressively larger spatial extents:

where

The mean (${\widehat{r}}_{n,t}$) represents the annual rate of change within intermediate population levels (neighborhood n) within which each lek (l) is spatially nested. The standard deviation (${\widehat{\sigma }}_{r_l}$) quantifies inter-annual, environmental variability affecting lek-level rates of change. We chose a weakly informative prior to represent ${\widehat{\sigma }}_{r_l}$ (eq. 55). Intrinsic growth rates within intermediate population levels (${\widehat{r}}_{n,t}$) were modeled as normally distributed realizations with mean (${\widehat{r}}_{c,t}$) and standard deviation (${\widehat{\sigma }}_{r_n}$) parameters representing central tendencies of the Bi-State DPS (c) and inter-annual, environmental variability at the intermediate population level, respectively. Bi-State DPS rates of change were described using a similar specification, but with the mean (${\widehat{\mu }}_{r_C}$) parameter representing the long-term average intrinsic rate of population change of the Bi-State DPS. We assigned a vague, uniform prior for ${\widehat{\mu }}_{r_C}$ (eq. 60), which allows for percentage changes ranging from approximately −10 to 10 percent, which is beyond estimates reported for long-term sage-grouse population performance range-wide (Garton and others, 2011), and within the Bi-State DPS (Coates and others, 2018).The final parameter within our model (${\widehat{\sigma }}_r$) represented a measure of inter-annual, environmental variability affecting population rates of change within the distributional extent of the species’ range within the United States (Coates and others, 2021).

The observation process of our model mapped the latent state of abundance onto the observed data ($y_{l,t}$), which were maximum counts of male sage-grouse (y) attending leks (l) in each year (t; eq. 61).

where

Counts were modeled using a Poisson distribution, which assumes increasing variance with increasing individuals counted. We specified a vague prior with lower and upper bounds of 2 and 500, respectively, to describe initial populations for each lek (l) during year one (t=1) (eq. 62).

We ran the model on three chains of 100,000 iterations each, discarded the first 50,000 iterations from each chain as burn-in, and used a thinning rate set to 50, retaining 3,000 posterior samples for all parameters monitored. We evaluated chain convergence based on a potential scale reduction factor value less than 1.2 (Gelman and Rubin, 1992; Brooks and Gelman, 1998).

Evaluating Conservation Action Effectiveness

We used cross-scale comparisons of sage-grouse population performance before and after treatments were applied (BACI study design) to evaluate the effectiveness of CEs. Using the posterior distributions of abundance from the SSM as a performance metric, we partitioned that index into one of two treatment (control or impact) and period (before or after) categories. We split treatment periods at the year in which the CE occurred. Leks were assigned to the impact group if they were within 5 km (Coates and others, 2013) of at least one CE and were considered active (recording greater than two males during at least one count) during at least 1 year during the post-treatment period. By requiring that leks remained active for at least 1 post-treatment year, we ensured that CEs analyzed here had the opportunity to affect population performance. Lek count data were collected during the spring, so we assigned abundance estimates from the same year as the CE to the period before implementation because it was likely that CE implementation occurred after those counts. We lacked a conventional control and so used all leks within the DPS that were not within 5 km of the CE being assessed as the spatial reference. Finally, we only considered CEs that were implemented within a 5-km radius of greater than or equal to three leks. The number of leks used to evaluate the effectiveness of specific CEs are provided in table 10.

We incorporated a paired-series component (BACIPS) into the traditional BACI methodology that accounted for time-dependent changes in our response variable. Both traditional BACI and BACIPS analyses assume that the system being measured responds in an immediate and constant manner. This assumption is observed during steps that reduce the temporal variation of the time series (Thiault and others, 2017) into four categories (before-control, before-impact, after-control, after-impact). Although this may be appropriate in some cases, the effect of CEs on sage-grouse population dynamics is expected to be more complex, with lagged and potentially non-linear abundance responses. The PC BACIPS framework (Thiault and others, 2017) uses multi-model inference to evaluate alternative functional relationships for the relativized BACI effect. We used PC BACIPS methods to investigate whether CEs implemented within the Bi-State region had effects on sage-grouse populations that were either (1) immediate and constant (step change; classic BACIPS model), (2) linear (the difference between control and affected area changes at a constant rate through time), (3) asymptotic (the rate of change monotonically decreases to zero as time passes), or (4) sigmoidal (the rate of change initially increases over time but then decreases to zero). Conservation efforts were implemented at different times for different leks, but our interest was in the general effectiveness of CEs rather than on any one population’s (that is, lek) response. Therefore, we standardized lek-level time series such that CE effectiveness was measured in years since CE implementation. We first used a standard BACIPS analysis applied to posterior parameter distributions (Conner and others, 2016) for each lek but did not average over the post-treatment intervals. We calculated R-ratios of abundance during the pre-treatment period that represented long-term (9-year) averages from immediately before CE implementation (eq. 63):

where

However, we calculated R-ratios of abundance during the post-treatment period using an annual form (eq. 64):

where Any single lek was monitored following CE implementation for a maximum of 9 years (initial year plus 8 lag years). This length of time informed the duration of the pre-treatment period and marked the total number of years captured in the BACIPS ratio (eq. 65): where

We then standardized the temporal indices of the BACIPS ratios by subtracting the year of CE implementation, which allowed us to aggregate lek-level BACIPS ratios using intra-interval averages. The PC BACIPS analysis requires a complete time series including pre- and post-treatment periods for each treatment (impact and control), so we decomposed the time-standardized BACIPS ratios into the individual time series components. For the control group, we set the time series (9 years before and 9 years after) to a fixed value equaling one for simplicity and we assigned the same value to the before period for the impact time series. As a result, we could use the time-standardized BACIPS ratios as the after-impact time series because all other denominators had been cancelled out. To include the uncertainty in the SSM outputs (that is, posterior draws of abundance parameters), we iteratively applied the PC BACIPS method to a single set of BACIPS, recording the top-performing model and corresponding parameter estimates during each iteration. We then summed the number of times each model (that is, step, linear, sigmoid, or asymptotic) was scored as the top model and chose the model with the greatest number of overall selections for inference and line fitting. We based inference and line fit on the mean estimate of model parameters from the iterations corresponding to the overall, top-performing model. We also calculated traditional BACI ratios in addition to the PC BACIPS analysis to identify instances where the traditional BACI method underestimated current (2021) and future effects following CE implementation. Finally, to use as much information from the CED as possible, we implemented the same steps as described above but summarized each unique CE independently. This part of the analysis could be used by land managers to prioritize specific types of CEs that could be implemented in the future based on their relative effects on sage-grouse population performance.