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)
library(glmmTMB)
library(tweedie)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(sandwich) ##
## Attaching package: 'sandwich'## The following objects are masked from 'package:glmmTMB':
##
## meatHC, sandwichlibrary(lmtest) ## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(merDeriv) ## Loading required package: nonnest2## This is nonnest2 0.5-8.
## nonnest2 has not been tested with all combinations of supported model classes.## Loading required package: lavaan## This is lavaan 0.6-20
## lavaan is FREE software! Please report any bugs.library(clubSandwich)## Registered S3 methods overwritten by 'clubSandwich':
## method from
## bread.lmerMod merDeriv
## bread.mlm sandwichlibrary(emmeans)## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'library(GLMMadaptive)##
## Attaching package: 'GLMMadaptive'## The following objects are masked from 'package:lavaan':
##
## anova, coef, fitted, nobs, predict, residuals## The following object is masked from 'package:lme4':
##
## negative.binomialchpt1_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.chpt1_data <- chpt1_data %>%
mutate(
dead.biomass = ifelse(dead.biomass == 0, 0.0001, dead.biomass),
live.biomass = ifelse(live.biomass == 0, 0.0001, live.biomass)
)Global ROS model
ros.global_model <- glmer(
ros ~ patch.location + season +
rel.humidity + mf.wind.sp + herbicide +
(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 + herbicide +
## (1 | block)
## Data: chpt1_data
##
## AIC BIC logLik -2*log(L) df.resid
## -231.0 -213.6 124.5 -249.0 42
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.7500 -0.5454 -0.3882 0.0913 3.9357
##
## Random effects:
## Groups Name Variance Std.Dev.
## block (Intercept) 1.325e-08 0.0001151
## Residual 1.756e+00 1.3251306
## Number of obs: 51, groups: block, 5
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) -5.210421 1.145500 -4.549 5.4e-06 ***
## patch.locationinside 0.087978 0.444141 0.198 0.8430
## patch.locationoutside -0.656241 0.424994 -1.544 0.1226
## seasongrowing -0.976346 0.493044 -1.980 0.0477 *
## rel.humidity -0.002232 0.034566 -0.065 0.9485
## mf.wind.sp 0.834779 0.397051 2.102 0.0355 *
## herbicideyes -0.157157 0.469621 -0.335 0.7379
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ptch.lctnn ptch.lctnt ssngrw rl.hmd mf.wn.
## ptch.lctnns -0.497
## ptch.lctnts -0.449 0.581
## seasongrwng 0.473 -0.427 -0.193
## rel.humidty -0.653 0.419 0.216 -0.747
## mf.wind.sp 0.091 -0.274 -0.093 0.539 -0.790
## herbicideys -0.542 -0.017 0.317 -0.015 0.208 0.096
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')sim_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 + herbicide +
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 + herbicide + rel.humidity + mf.wind.sp + (1 | block)
## Data: chpt1_data
##
## AIC BIC logLik -2*log(L) df.resid
## -231.0 -217.5 122.5 -245.0 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.6739 -0.5115 -0.2860 -0.0306 4.0561
##
## Random effects:
## Groups Name Variance Std.Dev.
## block (Intercept) 9.943e-10 3.153e-05
## Residual 2.184e+00 1.478e+00
## Number of obs: 51, groups: block, 5
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) -5.492801 1.005351 -5.464 4.67e-08 ***
## seasongrowing -0.867132 0.445701 -1.946 0.0517 .
## herbicideyes 0.170150 0.450562 0.378 0.7057
## rel.humidity -0.007519 0.031638 -0.238 0.8121
## mf.wind.sp 0.921376 0.388877 2.369 0.0178 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ssngrw hrbcdy rl.hmd
## seasongrwng 0.345
## herbicideys -0.627 -0.049
## rel.humidty -0.580 -0.666 0.277
## mf.wind.sp -0.064 0.447 0.089 -0.761
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')sim_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 + herbicide + rel.humidity + mf.wind.sp + (1 | block)
## ros.global_model: ros ~ patch.location + season + rel.humidity + mf.wind.sp + herbicide + (1 | block)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## ros.no.loc 7 -231.03 -217.50 122.51 -245.03
## ros.global_model 9 -231.01 -213.62 124.51 -249.01 3.9848 2 0.1364max temp global model
global.max.temp.model <- lmer(
max.temp ~ patch.location + season + herbicide +
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 + herbicide + rel.humidity +
## mf.wind.sp + (1 | pair) + (1 | block)
## Data: chpt1_data
##
## REML criterion at convergence: 591
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.03885 -0.57043 0.04223 0.54815 2.01876
##
## Random effects:
## Groups Name Variance Std.Dev.
## pair (Intercept) 3465 58.87
## block (Intercept) 0 0.00
## Residual 20829 144.32
## Number of obs: 51, groups: pair, 20; block, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 633.9914 137.0048 4.628
## patch.locationinside 123.3802 55.4542 2.225
## patch.locationoutside -2.9217 62.0186 -0.047
## seasongrowing -70.6592 73.5358 -0.961
## herbicideyes -42.7880 61.8966 -0.691
## rel.humidity -0.1268 4.8583 -0.026
## mf.wind.sp -28.7496 59.9964 -0.479
##
## Correlation of Fixed Effects:
## (Intr) ptch.lctnn ptch.lctnt ssngrw hrbcdy rl.hmd
## ptch.lctnns -0.316
## ptch.lctnts -0.353 0.646
## seasongrwng 0.457 -0.123 -0.119
## herbicideys -0.314 0.221 0.486 -0.120
## rel.humidty -0.597 0.158 0.150 -0.714 0.215
## mf.wind.sp 0.090 -0.174 -0.158 0.480 -0.196 -0.817
## 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.7291 2 0.01272 *
## season 0.9233 1 0.33661
## herbicide 0.4779 1 0.48939
## rel.humidity 0.0007 1 0.97917
## mf.wind.sp 0.2296 1 0.63180
## ---
## 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 + herbicide +
(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 + herbicide + (1 |
## pair) + (1 | block)
## Data: chpt1_data
##
## REML criterion at convergence: 618.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.31977 -0.64950 0.02826 0.52951 2.34304
##
## Random effects:
## Groups Name Variance Std.Dev.
## pair (Intercept) 0 0.0
## block (Intercept) 0 0.0
## Residual 26651 163.3
## Number of obs: 51, groups: pair, 20; block, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 669.3434 120.8710 5.538
## seasongrowing -72.4103 72.2688 -1.002
## rel.humidity -0.4044 5.0344 -0.080
## mf.wind.sp -21.3353 62.3013 -0.342
## herbicideyes -16.8612 54.7618 -0.308
##
## Correlation of Fixed Effects:
## (Intr) ssngrw rl.hmd mf.wn.
## seasongrwng 0.459
## rel.humidty -0.581 -0.741
## mf.wind.sp 0.091 0.540 -0.851
## herbicideys -0.211 -0.080 0.165 -0.139
## 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.0039 1 0.3164
## rel.humidity 0.0065 1 0.9360
## mf.wind.sp 0.1173 1 0.7320
## herbicide 0.0948 1 0.7582sim_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 + herbicide + (1 | pair) + (1 | block)
## global.max.temp.model: max.temp ~ patch.location + season + herbicide + rel.humidity + mf.wind.sp + (1 | pair) + (1 | block)
## npar AIC BIC logLik -2*log(L) Chisq Df
## reduced.max.temp.model 8 675.19 690.64 -329.59 659.19
## global.max.temp.model 10 671.22 690.54 -325.61 651.22 7.9702 2
## Pr(>Chisq)
## reduced.max.temp.model
## global.max.temp.model 0.01859 *
## ---
## 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()`).
Live biomass global model:
chpt1_livenoNA <- chpt1_data %>%
filter(!is.na(live.biomass))
live.biomass.global.model <- glm(
live.biomass ~ season + patch.location,
data = chpt1_livenoNA,
family = Gamma(link = "log")
)
summary(live.biomass.global.model)##
## Call:
## glm(formula = live.biomass ~ season + patch.location, family = Gamma(link = "log"),
## data = chpt1_livenoNA)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7182 0.3713 1.934 0.061 .
## seasongrowing -0.1843 0.4127 -0.447 0.658
## patch.locationinside -0.5853 0.4704 -1.244 0.221
## patch.locationoutside -0.2740 0.4704 -0.582 0.564
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 1.427309)
##
## Null deviance: 110.47 on 39 degrees of freedom
## Residual deviance: 107.99 on 36 degrees of freedom
## AIC: 100.11
##
## Number of Fisher Scoring iterations: 14car::Anova(live.biomass.global.model) ## Analysis of Deviance Table (Type II tests)
##
## Response: live.biomass
## LR Chisq Df Pr(>Chisq)
## season 0.18533 1 0.6668
## patch.location 1.55176 2 0.4603sim_res <- simulateResiduals(fittedModel = live.biomass.global.model, n = 1000)
plot(sim_res)
plotQQunif(sim_res)
check_model(live.biomass.global.model)
live.biomass.vcov <- sandwich(live.biomass.global.model)
coeftest(live.biomass.global.model, vcov. = live.biomass.vcov)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.71822 0.41384 1.7355 0.08265 .
## seasongrowing -0.18429 0.37530 -0.4910 0.62339
## patch.locationinside -0.58525 0.47852 -1.2231 0.22131
## patch.locationoutside -0.27396 0.45230 -0.6057 0.54471
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Histogram 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()`).
Litter biomass model:
chpt1_litternoNA <- chpt1_data %>%
filter(!is.na(litter.biomass))
litter.biomass.global.model <- glmer(
litter.biomass ~ season + patch.location +
(1 | block) + (1 | pair),
data = chpt1_litternoNA,
family = Gamma(link = "log"),
control = glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 1e5))
)
summary(litter.biomass.global.model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: litter.biomass ~ season + patch.location + (1 | block) + (1 |
## pair)
## Data: chpt1_litternoNA
## Control: glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 1e+05))
##
## AIC BIC logLik -2*log(L) df.resid
## 336.6 348.4 -161.3 322.6 33
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3244 -0.5939 -0.1774 0.3039 2.3226
##
## Random effects:
## Groups Name Variance Std.Dev.
## pair (Intercept) 0.41969 0.6478
## block (Intercept) 0.03398 0.1843
## Residual 0.48656 0.6975
## Number of obs: 40, groups: pair, 14; block, 4
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 2.84319 0.43493 6.537 6.27e-11 ***
## seasongrowing -0.03515 0.66695 -0.053 0.958
## patch.locationinside -0.11391 0.26553 -0.429 0.668
## patch.locationoutside 0.30009 0.25954 1.156 0.248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ssngrw ptch.lctnn
## seasongrwng -0.578
## ptch.lctnns -0.319 -0.010
## ptch.lctnts -0.304 -0.032 0.522reduced_patch <- update(litter.biomass.global.model, . ~ . - patch.location)
anova(reduced_patch, litter.biomass.global.model, test = "Chisq")## Data: chpt1_litternoNA
## Models:
## reduced_patch: litter.biomass ~ season + (1 | block) + (1 | pair)
## litter.biomass.global.model: litter.biomass ~ season + patch.location + (1 | block) + (1 | pair)
## npar AIC BIC logLik -2*log(L) Chisq Df
## reduced_patch 5 335.33 343.77 -162.66 325.33
## litter.biomass.global.model 7 336.57 348.39 -161.28 322.57 2.763 2
## Pr(>Chisq)
## reduced_patch
## litter.biomass.global.model 0.2512reduced_season <- update(litter.biomass.global.model, . ~ . - season)
anova(reduced_season, litter.biomass.global.model, test = "Chisq")## Data: chpt1_litternoNA
## Models:
## reduced_season: litter.biomass ~ patch.location + (1 | block) + (1 | pair)
## litter.biomass.global.model: litter.biomass ~ season + patch.location + (1 | block) + (1 | pair)
## npar AIC BIC logLik -2*log(L) Chisq Df
## reduced_season 6 334.57 344.70 -161.28 322.57
## litter.biomass.global.model 7 336.57 348.39 -161.28 322.57 0.0033 1
## Pr(>Chisq)
## reduced_season
## litter.biomass.global.model 0.9544sim_res <- simulateResiduals(fittedModel = litter.biomass.global.model, n = 1000)
plot(sim_res)
plotQQunif(sim_res)
check_model(litter.biomass.global.model)
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()`).
Dead biomass model:
chpt1_deadnoNA <- chpt1_data %>%
filter(!is.na(dead.biomass))
dead.biomass.global.model <- glm(
dead.biomass ~ season + patch.location,
data = chpt1_deadnoNA,
family = Gamma(link = "log"),
)
summary(dead.biomass.global.model)##
## Call:
## glm(formula = dead.biomass ~ season + patch.location, family = Gamma(link = "log"),
## data = chpt1_deadnoNA)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.5657 0.3173 8.086 1.3e-09 ***
## seasongrowing -0.9459 0.3526 -2.682 0.011 *
## patch.locationinside 0.6462 0.4020 1.608 0.117
## patch.locationoutside -0.1827 0.4020 -0.454 0.652
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 1.042267)
##
## Null deviance: 90.499 on 39 degrees of freedom
## Residual deviance: 76.615 on 36 degrees of freedom
## AIC: 280.25
##
## Number of Fisher Scoring iterations: 17car::Anova(dead.biomass.global.model) ## Analysis of Deviance Table (Type II tests)
##
## Response: dead.biomass
## LR Chisq Df Pr(>Chisq)
## season 5.8942 1 0.01519 *
## patch.location 5.0514 2 0.08000 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sim_res <- simulateResiduals(fittedModel = dead.biomass.global.model, n = 1000)
plot(sim_res)
plotQQunif(sim_res)
check_model(dead.biomass.global.model)
dead.biomass.vcov <- sandwich(dead.biomass.global.model)
coeftest(dead.biomass.global.model, vcov. = dead.biomass.vcov)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.56574 0.26927 9.5285 < 2e-16 ***
## seasongrowing -0.94593 0.41372 -2.2864 0.02223 *
## patch.locationinside 0.64620 0.28817 2.2424 0.02494 *
## patch.locationoutside -0.18268 0.42465 -0.4302 0.66705
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Dead biomass emmeans:
emm_patch <- emmeans(dead.biomass.global.model, ~ patch.location)
pairs_patch <- pairs(emm_patch, adjust = "tukey")
pairs_patch## contrast estimate SE df t.ratio p.value
## edge - inside -0.646 0.402 36 -1.608 0.2555
## edge - outside 0.183 0.402 36 0.454 0.8927
## inside - outside 0.829 0.386 36 2.148 0.0944
##
## Results are averaged over the levels of: season
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimatespairs_df <- as.data.frame(pairs_patch)
head(pairs_df)## contrast estimate SE df t.ratio p.value
## edge - inside -0.6461973 0.4019767 36 -1.608 0.2555
## edge - outside 0.1826815 0.4019767 36 0.454 0.8927
## inside - outside 0.8288788 0.3858696 36 2.148 0.0944
##
## Results are averaged over the levels of: season
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimatesHistogram 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()`).
live fuel moisture models:
chpt1_fmc_noNA <- chpt1_data %>% filter(!is.na(live.fmc))
live.fmc.model <- lmer(
live.fmc ~ season + patch.location + (1 | block) + (1 | pair),
data = chpt1_fmc_noNA
)
summary(live.fmc.model)## Linear mixed model fit by REML ['lmerMod']
## Formula: live.fmc ~ season + patch.location + (1 | block) + (1 | pair)
## Data: chpt1_fmc_noNA
##
## REML criterion at convergence: 317.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.87271 -0.49532 0.05793 0.59459 2.16896
##
## Random effects:
## Groups Name Variance Std.Dev.
## pair (Intercept) 1281.4 35.80
## block (Intercept) 172.5 13.13
## Residual 1558.0 39.47
## Number of obs: 33, groups: pair, 13; block, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 114.034 20.179 5.651
## seasongrowing 26.270 30.808 0.853
## patch.locationinside -3.626 18.131 -0.200
## patch.locationoutside -12.274 18.086 -0.679
##
## Correlation of Fixed Effects:
## (Intr) ssngrw ptch.lctnn
## seasongrwng -0.421
## ptch.lctnns -0.496 -0.058
## ptch.lctnts -0.517 -0.015 0.592car::Anova(live.fmc.model) ## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: live.fmc
## Chisq Df Pr(>Chisq)
## season 0.7271 1 0.3938
## patch.location 0.5234 2 0.7697sim_res <- simulateResiduals(fittedModel = live.fmc.model, n = 1000)
plot(sim_res)
plotQQunif(sim_res)
check_model(live.fmc.model)
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()`).
Litter fuel moisture model:
chpt1_litter_noNA <- chpt1_data %>% filter(!is.na(litter.fmc))
#litter.fmc.model <- glmer(
# litter.fmc ~ season + patch.location +
# (1 | pair),
# data = chpt1_litternoNA,
# family = Gamma(link = "log"),
# control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e5))
#)
library(geepack)
gee_model.litter.fmc <- geeglm(
litter.fmc ~ season + patch.location,
id = pair,
data = chpt1_litternoNA,
family = Gamma(link = "log"),
corstr = "exchangeable"
)
summary(gee_model.litter.fmc)##
## Call:
## geeglm(formula = litter.fmc ~ season + patch.location, family = Gamma(link = "log"),
## data = chpt1_litternoNA, id = pair, corstr = "exchangeable")
##
## Coefficients:
## Estimate Std.err Wald Pr(>|W|)
## (Intercept) 2.0452 0.1813 127.247 < 2e-16 ***
## seasongrowing 0.9198 0.2955 9.686 0.00186 **
## patch.locationinside -0.3341 0.1683 3.939 0.04719 *
## patch.locationoutside -0.1137 0.1814 0.393 0.53064
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation structure = exchangeable
## Estimated Scale Parameters:
##
## Estimate Std.err
## (Intercept) 0.2677 0.06141
## Link = identity
##
## Estimated Correlation Parameters:
## Estimate Std.err
## alpha 0.4345 0.2057
## Number of clusters: 13 Maximum cluster size: 3resid_pearson <- residuals(gee_model.litter.fmc, type = "pearson")
fitted_vals <- fitted(gee_model.litter.fmc)
plot(fitted_vals, resid_pearson,
xlab = "Fitted values",
ylab = "Pearson residuals",
main = "GEE: Residuals vs Fitted")
abline(h = 0, col = "red", lty = 2)
qqnorm(resid_pearson)
qqline(resid_pearson, col = "red")
#summary(litter.fmc.model)
#full <- litter.fmc.model
#reduced_patch <- update(full, . ~ . - patch.location)
#anova(reduced_patch, full, test = "Chisq")
#reduced_season <- update(full, . ~ . - season)
#anova(reduced_season, full, test = "Chisq")
#sim_res <- simulateResiduals(fittedModel = gee_model.litter.fmc, n = 1000)
#plot(sim_res)
#plotQQunif(sim_res)
check_model(gee_model.litter.fmc)## Converting missing values (`NA`) into regular values currently not
## possible for variables of class `NULL`.## Cannot simulate residuals for models of class `geeglm`. Please try
## `check_model(..., residual_type = "normal")` instead.
Multiple comparisons test for patch.location term in litter fuel
moisture model:
emm_patch <- emmeans(gee_model.litter.fmc, ~ patch.location)
pairs_patch <- pairs(emm_patch, adjust = "tukey")
pairs_patch## contrast estimate SE df t.ratio p.value
## edge - inside 0.334 0.168 31 1.985 0.1329
## edge - outside 0.114 0.181 31 0.627 0.8065
## inside - outside -0.220 0.136 31 -1.617 0.2538
##
## Results are averaged over the levels of: season
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimatessummary(pairs_patch)## contrast estimate SE df t.ratio p.value
## edge - inside 0.334 0.168 31 1.985 0.1329
## edge - outside 0.114 0.181 31 0.627 0.8065
## inside - outside -0.220 0.136 31 -1.617 0.2538
##
## Results are averaged over the levels of: season
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimatesHistogram 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()`).
Dead fuel moisture model:
chpt1_dead_noNA <- chpt1_data %>% filter(!is.na(dead.fmc))
gee_dead_fmc <- geeglm(
dead.fmc ~ season + patch.location,
id = pair,
data = chpt1_litternoNA,
family = Gamma(link = "log"),
corstr = "exchangeable"
)
summary(gee_dead_fmc)##
## Call:
## geeglm(formula = dead.fmc ~ season + patch.location, family = Gamma(link = "log"),
## data = chpt1_litternoNA, id = pair, corstr = "exchangeable")
##
## Coefficients:
## Estimate Std.err Wald Pr(>|W|)
## (Intercept) 1.845 0.283 42.63 6.6e-11 ***
## seasongrowing 1.104 0.349 10.02 0.0016 **
## patch.locationinside -0.145 0.329 0.20 0.6588
## patch.locationoutside 0.538 0.381 1.99 0.1579
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation structure = exchangeable
## Estimated Scale Parameters:
##
## Estimate Std.err
## (Intercept) 0.709 0.279
## Link = identity
##
## Estimated Correlation Parameters:
## Estimate Std.err
## alpha 0.0397 0.156
## Number of clusters: 14 Maximum cluster size: 3emm_patch <- emmeans(gee_dead_fmc, ~ patch.location)
pairs_patch <- pairs(emm_patch, adjust = "tukey")
pairs_patch## contrast estimate SE df t.ratio p.value
## edge - inside 0.145 0.329 26 0.442 0.8990
## edge - outside -0.538 0.381 26 -1.412 0.3490
## inside - outside -0.683 0.402 26 -1.698 0.2250
##
## Results are averaged over the levels of: season
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimatessummary(pairs_patch)## contrast estimate SE df t.ratio p.value
## edge - inside 0.145 0.329 26 0.442 0.8990
## edge - outside -0.538 0.381 26 -1.412 0.3490
## inside - outside -0.683 0.402 26 -1.698 0.2250
##
## Results are averaged over the levels of: season
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimatesresid_pearson <- residuals(gee_dead_fmc, type = "pearson")
fitted_vals <- fitted(gee_dead_fmc)
plot(fitted_vals, resid_pearson,
xlab = "Fitted values",
ylab = "Pearson residuals",
main = "GEE: Residuals vs Fitted")
abline(h = 0, col = "red", lty = 2)
qqnorm(resid_pearson)
qqline(resid_pearson, col = "red")
Multiple comparisons test for patch.location term in dead fuel moisture
model:
#library(emmeans)
#emm_patch <- emmeans(dead.fmc.model, ~ patch.location | season)
#pairs_patch <- pairs(emm_patch, adjust = "tukey")
#pairs_patch
#pairs_df <- as.data.frame(pairs_patch)
#head(pairs_df)