Fuel models
Sarah Kincy
2025-09-12
Library and data:
library(ggplot2)
library(readr)
library(lme4)## Loading required package: Matrixlibrary(performance)
library(DHARMa)## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')library(see)
chpt1_data <- read_csv("maindatachpt1.csv")## Rows: 51 Columns: 20## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): position, season, patch.location, herbicide, burn.date
## dbl (14): block, pair, max.temp, ros, live.biomass, litter.biomass, dead.bi...
## time (1): time stamp
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Global ROS model
ros.global_model <- glmer(
ros ~ patch.location + season +
rel.humidity + mf.wind.sp +
(1 | block),
data = chpt1_data,
family = Gamma(link = "log")
)## boundary (singular) fit: see help('isSingular')summary(ros.global_model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: ros ~ patch.location + season + rel.humidity + mf.wind.sp + (1 |
## block)
## Data: chpt1_data
##
## AIC BIC logLik -2*log(L) df.resid
## -232.9 -217.4 124.5 -248.9 43
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.7521 -0.5508 -0.3753 0.1054 3.7305
##
## Random effects:
## Groups Name Variance Std.Dev.
## block (Intercept) 0.000 0.000
## Residual 1.747 1.322
## Number of obs: 51, groups: block, 5
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) -5.414e+00 9.606e-01 -5.636 1.74e-08 ***
## patch.locationinside 8.465e-02 4.445e-01 0.190 0.8490
## patch.locationoutside -6.128e-01 4.038e-01 -1.517 0.1291
## seasongrowing -9.782e-01 4.937e-01 -1.982 0.0475 *
## rel.humidity 8.962e-05 3.378e-02 0.003 0.9979
## mf.wind.sp 8.481e-01 3.959e-01 2.142 0.0322 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ptch.lctnn ptch.lctnt ssngrw rl.hmd
## ptch.lctnns -0.598
## ptch.lctnts -0.358 0.623
## seasongrwng 0.554 -0.430 -0.203
## rel.humidty -0.656 0.434 0.171 -0.760
## mf.wind.sp 0.170 -0.279 -0.136 0.543 -0.832
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')car::Anova(ros.global_model) ## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ros
## Chisq Df Pr(>Chisq)
## patch.location 4.4126 2 0.11011
## season 3.9265 1 0.04753 *
## rel.humidity 0.0000 1 0.99788
## mf.wind.sp 4.5882 1 0.03219 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sim_res <- simulateResiduals(fittedModel = ros.global_model, n = 1000)
plotQQunif(sim_res)
check_model(ros.global_model)
ROS global model but with no location:
ros.no.loc <- glmer(
ros ~ season +
rel.humidity + mf.wind.sp +
(1 | block),
data = chpt1_data,
family = Gamma(link = "log")
)## boundary (singular) fit: see help('isSingular')summary(ros.no.loc)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: ros ~ season + rel.humidity + mf.wind.sp + (1 | block)
## Data: chpt1_data
##
## AIC BIC logLik -2*log(L) df.resid
## -232.9 -221.3 122.4 -244.9 45
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.6694 -0.5026 -0.2894 -0.0066 3.9549
##
## Random effects:
## Groups Name Variance Std.Dev.
## block (Intercept) 0.000 0.000
## Residual 2.214 1.488
## Number of obs: 51, groups: block, 5
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) -5.24885 0.78160 -6.716 1.87e-11 ***
## seasongrowing -0.85842 0.44537 -1.927 0.0539 .
## rel.humidity -0.01096 0.03044 -0.360 0.7187
## mf.wind.sp 0.90913 0.38778 2.344 0.0191 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ssngrw rl.hmd
## seasongrwng 0.404
## rel.humidty -0.543 -0.680
## mf.wind.sp -0.008 0.455 -0.823
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')car::Anova(ros.no.loc) ## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ros
## Chisq Df Pr(>Chisq)
## season 3.7150 1 0.05392 .
## rel.humidity 0.1297 1 0.71872
## mf.wind.sp 5.4964 1 0.01906 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sim_res <- simulateResiduals(fittedModel = ros.no.loc, n = 1000)
plotQQunif(sim_res)
check_model(ros.no.loc)
Chi-squared likelihood test:
anova(ros.global_model, ros.no.loc, test = "Chisq")## Data: chpt1_data
## Models:
## ros.no.loc: ros ~ season + rel.humidity + mf.wind.sp + (1 | block)
## ros.global_model: ros ~ patch.location + season + rel.humidity + mf.wind.sp + (1 | block)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## ros.no.loc 6 -232.88 -221.29 122.44 -244.88
## ros.global_model 8 -232.90 -217.45 124.45 -248.90 4.0185 2 0.1341max temp global model
global.max.temp.model <- lmer(
max.temp ~ patch.location + season +
rel.humidity + mf.wind.sp +
(1 | pair) + (1 | block),
data = chpt1_data
)## boundary (singular) fit: see help('isSingular')summary(global.max.temp.model)## Linear mixed model fit by REML ['lmerMod']
## Formula: max.temp ~ patch.location + season + rel.humidity + mf.wind.sp +
## (1 | pair) + (1 | block)
## Data: chpt1_data
##
## REML criterion at convergence: 601.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1307 -0.6031 0.0554 0.5485 2.1127
##
## Random effects:
## Groups Name Variance Std.Dev.
## pair (Intercept) 2370 48.69
## block (Intercept) 0 0.00
## Residual 21361 146.15
## Number of obs: 51, groups: pair, 20; block, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 603.2589 125.5534 4.805
## patch.locationinside 130.8005 54.5356 2.398
## patch.locationoutside 17.2756 54.7028 0.316
## seasongrowing -80.0395 71.1418 -1.125
## rel.humidity 0.7413 4.7113 0.157
## mf.wind.sp -38.1243 58.4851 -0.652
##
## Correlation of Fixed Effects:
## (Intr) ptch.lctnn ptch.lctnt ssngrw rl.hmd
## ptch.lctnns -0.277
## ptch.lctnts -0.251 0.629
## seasongrwng 0.450 -0.101 -0.071
## rel.humidty -0.569 0.115 0.053 -0.720
## mf.wind.sp 0.052 -0.135 -0.074 0.492 -0.822
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')car::Anova(global.max.temp.model) ## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: max.temp
## Chisq Df Pr(>Chisq)
## patch.location 8.1128 2 0.01731 *
## season 1.2658 1 0.26056
## rel.humidity 0.0248 1 0.87498
## mf.wind.sp 0.4249 1 0.51449
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sim_res <- simulateResiduals(fittedModel = global.max.temp.model, n = 1000)
plotQQunif(sim_res)
check_model(global.max.temp.model)
Reduced max temp global model:
library(lme4)
reduced.max.temp.model <- lmer(
max.temp ~ season +
rel.humidity + mf.wind.sp +
(1 | pair) + (1 | block),
data = chpt1_data
)## boundary (singular) fit: see help('isSingular')summary(reduced.max.temp.model)## Linear mixed model fit by REML ['lmerMod']
## Formula: max.temp ~ season + rel.humidity + mf.wind.sp + (1 | pair) +
## (1 | block)
## Data: chpt1_data
##
## REML criterion at convergence: 628.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.33077 -0.63381 0.00008 0.50445 2.27857
##
## Random effects:
## Groups Name Variance Std.Dev.
## pair (Intercept) 0 0.0
## block (Intercept) 0 0.0
## Residual 26138 161.7
## Number of obs: 51, groups: pair, 20; block, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 661.4726 116.9938 5.654
## seasongrowing -74.1859 71.3413 -1.040
## rel.humidity -0.1488 4.9174 -0.030
## mf.wind.sp -24.0014 61.0997 -0.393
##
## Correlation of Fixed Effects:
## (Intr) ssngrw rl.hmd
## seasongrwng 0.454
## rel.humidty -0.567 -0.740
## mf.wind.sp 0.063 0.536 -0.847
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')car::Anova(reduced.max.temp.model) ## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: max.temp
## Chisq Df Pr(>Chisq)
## season 1.0813 1 0.2984
## rel.humidity 0.0009 1 0.9759
## mf.wind.sp 0.1543 1 0.6944sim_res <- simulateResiduals(fittedModel = reduced.max.temp.model, n = 1000)
plotQQunif(sim_res)
check_model(reduced.max.temp.model)
Chi-squared likelihood test:
anova(reduced.max.temp.model, global.max.temp.model)## refitting model(s) with ML (instead of REML)## Data: chpt1_data
## Models:
## reduced.max.temp.model: max.temp ~ season + rel.humidity + mf.wind.sp + (1 | pair) + (1 | block)
## global.max.temp.model: max.temp ~ patch.location + season + rel.humidity + mf.wind.sp + (1 | pair) + (1 | block)
## npar AIC BIC logLik -2*log(L) Chisq Df
## reduced.max.temp.model 7 673.29 686.82 -329.65 659.29
## global.max.temp.model 9 669.63 687.02 -325.82 651.63 7.6601 2
## Pr(>Chisq)
## reduced.max.temp.model
## global.max.temp.model 0.02171 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Fuel Models
#Histograms of Fuel Variables# Live.biomass:
ggplot(chpt1_data, aes(x = live.biomass)) +
geom_histogram(binwidth = 0.1, fill = "lightgreen", color = "black") +
labs(title = "Frequency of Live Biomass",
x = "Live Biomass",
y = "Frequency")## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).
Histogram of litter.biomass:
ggplot(chpt1_data, aes(x = litter.biomass)) +
geom_histogram(binwidth = 1, fill = "lightblue", color = "black") +
labs(title = "Frequency of Litter Biomass",
x = "Litter Biomass",
y = "Frequency")## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).
Histogram of Dead.biomass:
ggplot(chpt1_data, aes(x = dead.biomass)) +
geom_histogram(binwidth = 1, fill = "purple", color = "black") +
labs(title = "Frequency of Dead Biomass",
x = "Dead Biomass",
y = "Frequency")## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).
Histogram of Live fuel moisture:
ggplot(chpt1_data, aes(x = live.fmc)) +
geom_histogram(binwidth = 2, fill = "darkgreen", color = "black") +
labs(title = "Frequency of Live Fuel Moisture",
x = "Live Fuel Moisture (%)",
y = "Frequency")## Warning: Removed 18 rows containing non-finite outside the scale range
## (`stat_bin()`).
Histogram of Litter fuel moisture:
ggplot(chpt1_data, aes(x = litter.fmc)) +
geom_histogram(binwidth = 1, fill = "cyan", color = "black") +
labs(title = "Frequency of Litter Fuel Moisture",
x = "Litter Fuel Moisture (%)",
y = "Frequency")## Warning: Removed 16 rows containing non-finite outside the scale range
## (`stat_bin()`).
Histogram of Dead fuel moisture:
ggplot(chpt1_data, aes(x = dead.fmc)) +
geom_histogram(binwidth = 1, fill = "violet", color = "black") +
labs(title = "Frequency of Dead Fuel Moisture",
x = "Dead Fuel Moisture (%)",
y = "Frequency")## Warning: Removed 21 rows containing non-finite outside the scale range
## (`stat_bin()`).