Take-aways
- Burn severity in vegetation increases with fuel load.
- There is no statistically-signficant relationship between substrate severity and fuel load.
- Other explanatory variables might include:
- Fire weather (RH, wind speed)
- Fuel (overall moisture, moisture by class + live:dead ratio, structure)
Data distributions
The Gamma distribution fits the severity data well.
Distribution of severity data as proportion of substrate (Sub) and vegetation (Veg) consumed by fire. Yellow line fits Gamma distribution.
Distribution of severity data as proportion of substrate (Sub) and vegetation (Veg) by location.
Distribution of plot-level fuel load data.
Relationship between severity and fuel load
- There is a statistically-significant relationship identifying fuel load as a driver of vegetation burn severity, but not substrate.
Scatterplot of burn severity against fuel load. Note that the severity axis is reversed to reflect 1 = highest severity and 5 = lowest severity.
Generalized linear mixed-effect regression models based on the Gamma distribution. mod1 fits burn severity against fuel load; mod0 fits null, intercept-only models.
Results of GLMM testing relationship between substrate severity and fuel load.
| Df | AIC | BIC | logLik | deviance | Chisq | Chi Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| sub0 | 3 | 167.9 | 176.7 | -80.93 | 161.9 | NA | NA | NA |
| sub1 | 4 | 169.5 | 181.3 | -80.75 | 161.5 | 0.3621 | 1 | 0.5473 |
Results of GLMM testing relationship between vegetation severity and fuel load.
| Df | AIC | BIC | logLik | deviance | Chisq | Chi Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| veg0 | 3 | 204.4 | 213.2 | -99.18 | 198.4 | NA | NA | NA |
| veg1 | 4 | 201.8 | 213.6 | -96.92 | 193.8 | 4.526 | 1 | 0.03338 |
| term | coef | lwr | upr |
|---|---|---|---|
| (Intercept) | 2.966 | 2.412 | 3.52 |
| AvgBiomass | 0.1122 | 0.2139 | 0.01051 |
Data & script
Data previews
severity## # A tibble: 140 x 8
## repID Sub Veg location site AvgBiomass tsf BurnYear
## <chr> <dbl> <dbl> <chr> <chr> <dbl> <chr> <chr>
## 1 GWCA_1 2.28 2.25 GWCA 01 3.30 2 2008
## 2 GWCA_2 3.06 1.94 GWCA 01 1.94 2 2010
## 3 GWCA_4 2.53 1.62 GWCA 01 3.07 2.4 2014
## 4 GWCA_5 3.97 4.03 GWCA 02 1.79 2 2008
## 5 GWCA_6 3.19 2.12 GWCA 02 2.03 2 2010
## 6 GWCA_8 2.72 2.16 GWCA 02 2.72 2.4 2014
## 7 GWCA_9 3.03 2.91 GWCA 03 2.49 3 2007
## 8 GWCA_10 3.22 2.41 GWCA 03 1.02 2 2010
## 9 GWCA_12 2.28 1.31 GWCA 03 3.31 2.4 2014
## 10 GWCA_13 3.09 2.66 GWCA 04 3.07 2 2008
## # … with 130 more rowsbiomass## # A tibble: 289 x 6
## location site repID AvgBiomass tsf BurnYear
## <chr> <chr> <chr> <dbl> <chr> <chr>
## 1 GWCA 01 GWCA_1 3.30 2 2008
## 2 GWCA 01 GWCA_2 1.94 2 2010
## 3 GWCA 01 GWCA_3 3.06 2 2012
## 4 GWCA 01 GWCA_4 3.07 2.4 2014
## 5 GWCA 02 GWCA_5 1.79 2 2008
## 6 GWCA 02 GWCA_6 2.03 2 2010
## 7 GWCA 02 GWCA_7 1.61 2 2012
## 8 GWCA 02 GWCA_8 2.72 2.4 2014
## 9 GWCA 03 GWCA_9 2.49 3 2007
## 10 GWCA 03 GWCA_10 1.02 2 2010
## # … with 279 more rowsR script
# Packages
pacman::p_load(tidyverse, magrittr)
#
severity %<>% filter(complete.cases(.))
# Fit gamma distributions
# MASS::fitdistr(severity$Sub+0.001, "Gamma")
# MASS::fitdistr(severity$Veg+0.001, "Gamma")
severity %<>% pivot_longer(cols = c(Sub, Veg),
names_to = "response",
values_to = "severity")
sev_gam <-
severity %>%
group_by(response) %>%
do(x = .$severity,
gamma = case_when(
.$response == "Sub" ~ dgamma(.$severity,
shape = 32.1,
rate = 11.1),
.$response == "Veg" ~ dgamma(.$severity,
shape = 13.4,
rate = 6)
) ) %>%
unnest(c(x, gamma))
# GGplot severity distribution
# Overall
severity %>%
ggplot() + theme_bw(16) +
geom_histogram(aes(x=severity,
y=..density..),
binwidth = 0.1,
fill=cbPal5[1],
alpha = 0.9) +
geom_line(data=sev_gam,
aes(x=x, y=gamma),
lwd = 2, alpha = 0.75,
color = cbPal5[2]) +
facet_wrap(~ response)
# By locations
ggplot(severity) + theme_bw(16) +
geom_histogram(aes(x=severity, y=..density..),
binwidth = 0.1,
fill=cbPal5[1], alpha = 0.9) +
facet_grid(location ~ response)
# Biomass distribution
ggplot(biomass) + theme_bw(16) +
geom_histogram(aes(x=AvgBiomass, y=..density..),
binwidth = 0.1,
fill=cbPal5[1], alpha = 0.9)
# Combine datasets
severity %<>%
pivot_wider(names_from = response,
values_from = severity,
id_cols = repID) %>%
inner_join(biomass)
# Severity vs. fuel load scatterplot
severity %>%
pivot_longer(cols = c("Sub", "Veg"),
names_to = 'response',
values_to = 'severity') %>%
ggplot() + theme_bw() +
geom_smooth(aes(x=AvgBiomass, y = severity,
color = location),
se = F, method = "lm") +
geom_point(aes(x=AvgBiomass, y = severity,
fill = location),
pch = 21, col="black",
stroke = 1.5, alpha = 0.75) +
scale_y_reverse() +
scale_fill_manual(values = cbPal5) +
scale_color_manual(values = cbPal5) +
facet_wrap(~ response) +
labs(x = "Fuel load",
y = "Severity ")
# GLMMs
sub0 <- lme4::glmer(Sub ~ 1 + (1|location),
family = Gamma(link = 'log'), severity)
sub1 <- lme4::glmer(Sub ~ AvgBiomass + (1|location),
family = Gamma(link = 'log'), severity)
veg0 <- lme4::glmer(Veg ~ 1 + (1|location),
family = Gamma(link = 'identity'), severity)
veg1 <- lme4::glmer(Veg ~ AvgBiomass + (1|location),
family = Gamma(link = 'identity'), severity)
# Model comparisons
anova(sub0, sub1)
anova(veg0, veg1)
# Coefficients
lme4::fixef(veg1)
confint(veg1, parm="beta_", method="Wald")