Pyrodiversity drivers
We assessed the hypothesized drivers of climate, topography and human influence on pyrodiversity using a 1) pyrodiversity model and a 2) burn activity model. These models represent direct and indirect (burn activity-mediated) effects on pyrodiversity, respectively (Fig. 1c). We model direct effects on pyrodiversity as:
[Eq. 1]
\begin{equation} \begin{matrix}\text{pyrodiversit}y_{i,j}\sim&Beta({\overline{P}}_{i,j},\theta)\\ logit({\overline{P}}_{i,j})=&\alpha_{0}+\alpha_{j}+\\ &\beta_{\text{AET}}*X_{1,i}+\beta_{\text{CWD}}*X_{2,i}+\beta_{AET*CWD}*X_{1,i}*X_{2,i}+\\ &\beta_{\text{elev}}*X_{3,i}+\beta_{\text{rough}}*X_{4,i}+\beta_{elev*rough}*X_{3,i}*X_{4,i}+\\ &\beta_{\text{pop.den}}*X_{5,i}+\beta_{\text{wild}}*X_{6,i}+\\ &\beta_{\text{prop.burn}}*X_{7,i}+\beta_{\text{prop.bur}n^{2}}*X_{8,i}\\ \alpha_{j}\sim&Normal(0,\sigma_{HUC4})\\ \end{matrix}\nonumber \\ \end{equation}
Where actual evapotranspiration (\(\beta_{\text{AET}}\)), accumulated climate water deficit (\(\beta_{\text{CWD}}\)) and their interaction (\(\beta_{AET*CWD}\)) are estimates of climatic effects. Elevation (\(\beta_{\text{elev}}\)), roughness (\(\beta_{\text{rough}}\)), and their interaction (\(\beta_{elev*rough})\) are topographic effects. Population density (\(\beta_{\text{pop.den}}\)) and proportion of watershed i ’s land area in wilderness (\(\beta_{\text{wild}}\)) are surrogates for human influence. We hypothesize much of the effects of these ultimate drivers are mediated by burn activity, here represented by proportion of flammable area burned between 1985 and 2018 (\(\beta_{\text{prop.burn}}\)) and its quadratic (\(\beta_{\text{prop.bur}n^{2}}\)). This metric is cumulative and can exceed 1 in the case of multiple burns in the same area. Our sample uniti are forested watersheds delineated by the 10-digit Hydrologic Unit Code (HUC10; median area = 755 km2). We account for spatial structuring of these units by including larger watershedj (2-digit HUC2s; median area = 437,000 km2), within which HUC10s are nested, as varying random intercepts. In total we assessed 1971 watersheds and 3306fires.
To quantify indirect effects on pyrodiversity, we modeled proportion burned area as a function of the same climate, topographic, and human influence variables as in equation 1, excluding\(\beta_{\text{prop.burn}}\) and \(\beta_{\text{prop.bur}n^{2}}\) (Eq. S1). This burn activity model is linked with the pyrodiversity model via a Bayesian multivariate and multi-level model using the brms and rstan packages in R (Bürkner 2017; Stan Development Team 2018; R Core Team 2019). The multivariate model allows us to predict direct and indirect effects of the ultimate drivers and quantify their combined effect while properly propagating uncertainty through the model chain. In this way marginal effects are estimated by first fitting the burn activity model to generate a posterior distribution of proportion burned area and subsequently incorporating this full distribution as predictors of\(\beta_{\text{prop.burn}}\) and \(\beta_{\text{prop.bur}n^{2}}\) in the pyrodiversity model. Model code, data, and additional methodological details can be found in the supplementary material.