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.