Burn severity measurements and analysis
We measured burn severity using the composite burn severity index (CBI; Key and Benson 2006) modified to accommodate our study area and objectives. We calculated a CBI value for each survey point (0–3 range) representing the mean of up to 11 components quantifying aspects of canopy structure, understory cover, and downed woody fuels (for details, see Appendix A). We derived components from field measurements of these attributes before fire (2002–2003), after prescribed fire (2004–2007), and after wildfire (2008–2010). Components represented either changes in these attributes from before to after disturbance, or aspects of burn severity apparent after disturbance (e.g., extent of char). We only measured points within burned units and assumed CBI = 0 for units that were not burned during our study. Others describe in detail how CBI values correspond with changes in various aspects of vegetation structure (Key and Benson 2006; Saab et al. 2006). In short, CBI = 0, 0 < CBI < 1.25, 1.25 < CBI < 2.25, and CBI > 2.25 are interpretable as unburned, low severity, moderate severity, and high severity, respectively. In general, low severity fire primarily affects understory vegetation with minimal canopy mortality (< 40%), whereas high severity fire results in much greater canopy mortality (> 70 %).
We quantified prescribed fire CBI (hereafter CBIPF) using environmental data collected before (2002–2003) versus after (2004–2007) prescribed fire, and wildfire CBI (hereafter CBIWF) with data from immediately before (2004–2007) versus after wildfire (2008–2010). Unfortunately, wildfire burned one treated study unit (Fitsum Creek) before we could measure post-treatment prescribed fire conditions (2004–2007). For this unit, we imputed CBIWF by 1) calculating CBITotal representing overall burn severity (i.e., changes from 2002–2003 to 2008–2010), 2) regressing CBIWF as a linear function of CBITotal at units where both were available (Buckhorn, Dutch Oven, Williams, and Deadman), and 3) using the resulting regression model (\({CBI}_{WF}={\beta }_{0}+{\beta }_{1}\times {CBI}_{Total}\), with estimates β0 [s.e.] = -0.04623 [0.07653] and β1 = 0.98415 [0.04271]) to impute missing data. As a covariate of occupancy, we imputed missing CBIWF values using a normally distributed prior with mean and SD representing model-predicted CBIWF.
Data for calculating CBIPF were relatively limited, so we did not use CBIPF as a covariate of occupancy. Rather, we modeled occupancy with a categorical treatment effect (TRTPF = 0 or 1 for survey points in untreated versus treated units, respectively). We then summarized CBIPF values where available (Dutch Oven, Parks Creek) to inform inference and compare with CBIWF. We compared CBIWF between treated and untreated points within treatment-control unit pairs and related CBIWF with CBIPF where available to evaluate the effect of prescribed fire treatments on subsequent wildfire severity.
Occupancy models
We analyzed avian relationships with prescribed fire using community occupancy models formulated within a Bayesian hierarchical modeling framework (Dorazio et al. 2006; Russell et al. 2009). Occupancy models leverage repeat-survey data to estimate species detectability (p) conditional upon occupancy (species occurrence within a specified time period), allowing unbiased estimation of occupancy probabilities (ψ) given sufficient data and adherence to model assumptions (MacKenzie et al. 2002; MacKenzie et al. 2006). We assumed that the occupancy states of species could change among years, but not between visits within a year. We estimated species-specific parameters as random variables governed by community-level parameters. The use of a common distribution among species improves the precision of species-specific parameter estimates, particularly for rare species (Dorazio et al. 2006, Russell et al. 2009). We excluded raptors, owls, and grouse from analysis because they were not readily detectable with our survey methods. Additionally, we only included species that bred in our study areas. For mobile animals such as birds, detectability (p) estimated with surveys repeated over a season quantifies both within-season movement and the observation process (i.e., availability and perceptibility; sensu Chandler and Royle 2013; Amundson et al. 2014). In principle, occupancy probabilities thereby estimated the probability of a surveyed point intersecting ≥ 1 home range for a given species (Efford and Dawson 2012; Latif et al. 2016a).
We compiled a 3-dimensional data matrix y, where element yijt was the sum of binary indicators for species detection (Sanderlin et al. 2014). Given a binary indicator xijkt = 1, we detected species i (i = 1,…,N) at survey point j (j = 1,…,J) during visit k (k = 1,…,K) in year t (t = 1,…,T; T = 4). Because we did not have covariates that differed for detection between visits, we analyzed the sum of all binary detections for species i over all visits at each survey point j in year t, where \({y}_{ijt}= \sum _{k=1}^{K}{x}_{ijkt}\) and yijt ϵ {0,1,…,K}. We modeled these data given probability of detection pi, and occupancy latent state zijt using a Bernoulli distribution with probability of success pi × zijt:
where the latent variable zijt for occupancy given probability of occupancy ψijt was modeled as:
We analyzed changes in species occupancy patterns using a model that leveraged our before-after, control-impact (BACI) sampling for examining disturbance effects (Popescu et al. 2012; Russell et al. 2015). For prescribed fire effects, we modeled occupancy (ψijt) as a function of treatment (TrtPF,j), period (Perjt = 0 or 1 for before or after survey point j was treated, respectively), and the interaction of the two (Trtj × Perjt). Thus,
\(\text{l}\text{o}\text{g}\text{i}\text{t}\left({\psi }_{ijt}\right)={\beta }_{0,i}+{\beta }_{{Per}_{PF}, i}\times {Per}_{PF,t}+{\beta }_{{Trt}_{PF}, i}\times {Trt}_{PF,j}+{\beta }_{{Per}_{PF}\times {Trt}_{PF}, i}\times {Per}_{PF,t}\times {Trt}_{PF,j}\) (Eq. 3),
where \({\beta }_{0,i}\) is the intercept and \({\beta }_{{Per}_{PF}, i}\), \({\beta }_{{Trt}_{PF}, i}\), and \({\beta }_{{Per}_{PF}\times {Trt}_{PF}, i}\) describe the additive and interactive effects of covariates \({Per}_{PF,t}\) and \({Trt}_{PF,j}\) on occupancy of species i at survey point j in year t. We restricted analysis of prescribed fire effects to data collected before wildfire (2002–2007). For wildfire effects, we analyzed data collected 2 years before and 3 years after wildfire (2006–2010) using two models. The first model analyzed overall wildfire effects:
\(\text{l}\text{o}\text{g}\text{i}\text{t}\left({\psi }_{ijt}\right)={\beta }_{0,i}+{\beta }_{{Per}_{WF}, i}\times {Per}_{WF,t}+{\beta }_{{CBI}_{WF}, i}\times {CBI}_{WF,j}+{\beta }_{{Per}_{WF}\times {CBI}_{WF}, i}\times {Per}_{WF,t}\times {CBI}_{WF,j}\) (Eq. 4).
The second model analyzed differences in wildfire effects between units treated versus untreated with prescribed fire:
\(\text{l}\text{o}\text{g}\text{i}\text{t}\left({\psi }_{ijt}\right)={\beta }_{0,i}+{\beta }_{{Per}_{WF}, i}\times {Per}_{WF,t}+{\beta }_{{Trt}_{PF}, i}\times {Trt}_{PF,t}+{\beta }_{{CBI}_{WF}, i}\times {CBI}_{WF,j}+{\beta }_{{Per}_{WF}\times {CBI}_{WF}, i}\times {Per}_{WF,t}\times {CBI}_{WF,j}+{\beta }_{{Trt}_{PF}\times {CBI}_{WF}, i}\times {Trt}_{PF,j}\times {CBI}_{WF,j}+{\beta }_{{Per}_{WF}\times {Trt}_{PF}\times {CBI}_{WF}, i}\times {Per}_{WF,t}\times {Trt}_{PF,j}\times {CBI}_{WF,j}\) (Eq. 5).
As in Eq. 3, \({\beta }_{0,i}\) is the intercept and all remaining β parameters describe additive and interactive effects of covariates on avian occupancy in Equations 4 and 5. All estimated parameters were species-specific normal random effects. For numerical purposes, we rescaled CBIWF values to mean = 0 and SD = 1 prior to analysis.
For all three models above (Equations 3–5), we drew inference of disturbance (prescribed fire or wildfire) effects from the extent to which occupancy shifted towards or away from burned (or unburned) survey points following disturbance. Interaction parameters in Equations 3, 4, and 5 quantified these shifts, whereas additive parameters controlled for potentially confounding environmental variation among survey points and time periods (Popescu et al. 2012). We considered statistically supported interaction parameters (90% Bayesian credible interval [BCI] excluded zero) strong evidence for disturbance effects.
We used one additional model to analyze annual changes in bird occupancy and time-dependent disturbance effects with all available data (2002–2010). This model included a random year effect and year-specific prescribed fire and wildfire effects:
\(\text{l}\text{o}\text{g}\text{i}\text{t}\left({\psi }_{ijt}\right)={\beta }_{0,it}+{\beta }_{{Trt}_{PF}, i{t}_{PF}}\times {Trt}_{PF,j}+{\beta }_{{CBI}_{WF}, i{t}_{WF}}\times {CBI}_{WF,j}\) (Eq. 6).
The intercept, \({\beta }_{0,it}\), varied with species and year according to nested normal random effects (year within species). Prescribed fire effects (\({\beta }_{{Trt}_{PF}, i{t}_{PF}}\)) were estimated separately for 4 distinct time periods, pre-treatment (tPF = 0) and 1–3 years post-treatment (tPF = 1–3, respectively). Similarly, wildfire effects (\({\beta }_{{CBI}_{WF}, i{t}_{WF}}\)) were estimated for 4 time periods, pre-fire (2006–2007; tWF = 0) and 1–3 years post-fire (2008–2010; tWF = 1–3, respectively). \({\beta }_{{Trt}_{PF}, i{t}_{PF}}\) was not estimated for 2008–2010 and \({\beta }_{{CBI}_{WF}, i{t}_{WF}}\) was not estimated for 2002–2005 for comparability with other models (see Equations 3–5). We used this model to look for time-dependencies in disturbance effects (i.e., where 95% BCIs for \({\beta }_{{Trt}_{PF}, i\left({t}_{PF}ϵ\left\{1, 2, 3\right\}\right)}-{\beta }_{{Trt}_{PF}, i\left({t}_{PF}=0\right)}\) or \({\beta }_{{CBI}_{WF}, i\left({t}_{WF}ϵ\left\{1, 2, 3\right\}\right)}-{\beta }_{{CBI}_{WF}, i\left({t}_{WF}=0\right)}\) excluded zero, Eq. 6). Additionally, we scanned yearly occupancy estimates for surveyed sites (\({\psi \text{'}}_{t}={\sum }_{j=1}^{J}{z}_{ijt}/J\)) to identify notable changes among pre-treatment (2002–2003), post-treatment (2004–2007), and post-wildfire (2008–2010) periods. All surveyed sites were burned by wildfire to some degree (min CBI = 0.39, see Results), so we expected some changes in overall occupancy for species with similar responses to low- versus high-severity wildfire. We considered inference from changes in overall occupancy weaker, however, because estimates of these changes did not control for potentially confounding factors as did shifts in occupancy with respect to CBI (see above).
In addition to species-specific relationships, we plotted emergent changes between species richness with treatment condition. We estimated species richness (Njt) at each survey point j and year t: \({N}_{jt}=\sum _{i=1}^{\text{m}\text{a}\text{x}\left(i\right)}{z}_{ijt}\). Similar to some (Russell et al. 2009, Latif et al. 2016) and unlike others (Dorazio et al. 2006, Kéry et al. 2009), we did not augment data to represent unobserved species, so community-level inferences were restricted to the subset of species observed at least once during our studies.
We modeled detectability as a species-specific normal random effect b0,i:
\(\text{l}\text{o}\text{g}\text{i}\text{t}\left({p}_{i}\right)={b}_{0,i}\) (Eq. 7),
where pi is the probability of detecting species i when surveying a given survey point in a given year when the species was present (i.e., ≥ 1 home range intersected the 100m point neighborhood). Unlike others (Russell et al. 2015), we did not consider treatment effects on detectability. Estimated effects on detectability from preliminary analyses were imprecise (credible intervals overlapped 0 for all species) and model convergence was difficult to achieve, suggesting the additional complexity strained limits of the data (Q. Latif unpublished data). We therefore only modeled heterogeneity in detectability among species and assumed detectability did not change with treatment condition. We modeled heterogeneity among species using a correlation term (ρ) between species intercepts of detection probability (b0,i) with occupancy probability (β0,i) (Dorazio and Royle 2005, Kéry et al. 2009).
We sampled posterior parameter distributions for all models using JAGS v4 (Plummer 2003) programmed from R (R Core Team 2013, Su and Yajima 2014). We used independent non-informative priors for all parameters (for priors, see Appendix B). For each model, we ran 6 parallel MCMC chains of length 100,000 it, burn-in 10,000 it, and thinning 100 it to sample posterior distributions. We verified that neffective ≥ 100 and \(\widehat{R}\) ≤ 1.1 for all parameters (Gelman and Hill 2007). We examined model goodness-of-fit (GOF) using posterior predictive testing (Gelman and Hill 2007). Specifically, we calculated a Bayesian p-value representing the proportion of simulated datasets drawn from model posterior predictive distributions with deviance higher than deviance for observed datasets from each location, where p < 0.05 or p > 0.95 constituted evidence for lack of fit.