chpt 2 models finalized
Sarah Kincy
2026-05-12
reading in data:
library(readr)
chpt2_r_data <- read_csv("~/Downloads/chapter 2 /chpt2.r.data.csv",
col_types = cols(block = col_character(),
patch = col_character(), position = col_character(),
location = col_factor(), # ← add this
season = col_factor(), # ← add this
max.temp = col_number(), ros = col_number(),
rcd10 = col_number(), rcd20 = col_number(),
mortality1 = col_number(), mortality2 = col_number(),
mortality3 = col_number(), mortality10 = col_number(),
mortality20 = col_number(), mortality100 = col_number(),
mortality200 = col_number(), green.needles = col_number(),
tree.char = col_number(), tree.scorch = col_number(),
rcd1.avg = col_number(), rcd2.avg = col_number(),
rcd3.avg = col_number(), rcd4.avg = col_number(),
difn = col_number()))## New names:
## • `` -> `...27`
## • `` -> `...28`
## • `` -> `...29`## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
## dat <- vroom(...)
## problems(dat)chpt2_r_data <- chpt2_r_data[-(37:52), ]
chpt2_r_data <- chpt2_r_data[, !(names(chpt2_r_data) %in% c("...27", "...28", "...29"))]Library:
library(glmmTMB)
library(DHARMa)## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')library(performance)
library(emmeans)## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'library(multcomp)## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyserlibrary(viridis)## Loading required package: viridisLitelibrary(ggplot2)
library(car)## Loading required package: carDataModel 1: growth of longleaf seedlings pre-fire
model1 <- glmmTMB(
rcd10 ~ location + difn + (1 | patch) + (1 | block),
data = chpt2_r_data,
family = Gamma(link = "log")
)
summary(model1)## Family: Gamma ( log )
## Formula: rcd10 ~ location + difn + (1 | patch) + (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## 89.3 100.4 -37.7 75.3 29
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 0.004362 0.06604
## block (Intercept) 0.165233 0.40649
## Number of obs: 36, groups: patch, 12; block, 4
##
## Dispersion estimate for Gamma family (sigma^2): 0.0834
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.344923 0.377524 0.914 0.3609
## locationedge 0.054709 0.118409 0.462 0.6441
## locationoutside 0.255773 0.120000 2.131 0.0331 *
## difn 0.006028 0.004723 1.276 0.2019
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(model1, type = "II")## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: rcd10
## Chisq Df Pr(>Chisq)
## location 5.0190 2 0.08131 .
## difn 1.6289 1 0.20186
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm_loc1 <- emmeans(model1, ~ location, type = "response")
pairs(emm_loc1, adjust = "holm")## contrast ratio SE df null z.ratio p.value
## inside / edge 0.947 0.1120 Inf 1 -0.462 0.6441
## inside / outside 0.774 0.0929 Inf 1 -2.131 0.0992
## edge / outside 0.818 0.0977 Inf 1 -1.682 0.1849
##
## P value adjustment: holm method for 3 tests
## Tests are performed on the log scaleemm_loc1## location response SE df asymp.LCL asymp.UCL
## inside 2.07 0.458 Inf 1.34 3.19
## edge 2.18 0.483 Inf 1.41 3.37
## outside 2.67 0.591 Inf 1.73 4.12
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scalecld(emm_loc1, adjust = "holm", Letters = letters)## location response SE df asymp.LCL asymp.UCL .group
## inside 2.07 0.458 Inf 1.21 3.51 a
## edge 2.18 0.483 Inf 1.28 3.70 a
## outside 2.67 0.591 Inf 1.57 4.54 a
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## Intervals are back-transformed from the log scale
## P value adjustment: holm method for 3 tests
## Tests are performed on the log scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.#sim_res1 <- simulateResiduals(model1)
#plot(sim_res1)
##testDispersion(sim_res1)
#testZeroInflation(sim_res1)
#testOutliers(sim_res1)
#check_model(model1)model 1 figure
# Estimated marginal means
emm_loc1 <- emmeans(model1, ~ location, type = "response")
# CLD (still created in case you need it elsewhere, but NOT plotted)
cld_loc1 <- cld(emm_loc1, adjust = "holm", Letters = letters)
cld_df1 <- as.data.frame(cld_loc1)
cld_df1$.group <- gsub(" ", "", cld_df1$.group)
# Colors
my_colors <- viridis(100, option = "D")[c(10, 40, 70, 90)]
# Plot
ggplot(cld_df1, aes(x = location, y = response, fill = location)) +
geom_col() +
geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.2) +
annotate("text",
x = -Inf,
y = Inf,
label = "χ² = 5.019 df = 2 p = 0.0813",
hjust = -0.1,
vjust = 1.1,
size = 5) +
scale_fill_manual(values = my_colors) +
labs(
x = "Patch Location",
y = "Change in RCD (mm)"
) +
coord_cartesian(clip = "off") +
theme_classic() +
theme(
legend.position = "none",
axis.title = element_text(size = 16),
axis.text = element_text(size = 14),
plot.margin = margin(10, 25, 10, 25)
)
Model 2: mortality of longleaf seedlings pre-fire (binary)
chpt2_r_data$mort_bin <- ifelse(chpt2_r_data$mort.1 == "yes", 1, 0)
model2 <- glmmTMB(
mort_bin ~ location + difn + (1 | patch) + (1 | block),
data = chpt2_r_data,
family = binomial(link = "logit")
)
summary(model2)## Family: binomial ( logit )
## Formula: mort_bin ~ location + difn + (1 | patch) + (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## 50.1 59.6 -19.1 38.1 30
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 0.2021 0.4496
## block (Intercept) 0.5298 0.7279
## Number of obs: 36, groups: patch, 12; block, 4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.27692 3.04316 -1.077 0.282
## locationedge 1.06164 1.09437 0.970 0.332
## locationoutside 0.68888 1.12102 0.615 0.539
## difn 0.02420 0.04006 0.604 0.546Anova(model2, type = "II")## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: mort_bin
## Chisq Df Pr(>Chisq)
## location 0.9424 2 0.6243
## difn 0.3647 1 0.5459emm_loc2 <- emmeans(model2, ~ location, type = "response")
pairs(emm_loc2, adjust = "holm")## contrast odds.ratio SE df null z.ratio p.value
## inside / edge 0.346 0.379 Inf 1 -0.970 0.9960
## inside / outside 0.502 0.563 Inf 1 -0.615 1.0000
## edge / outside 1.452 1.430 Inf 1 0.378 1.0000
##
## P value adjustment: holm method for 3 tests
## Tests are performed on the log odds ratio scaleemm_loc2## location prob SE df asymp.LCL asymp.UCL
## inside 0.148 0.135 Inf 0.0210 0.584
## edge 0.334 0.190 Inf 0.0862 0.728
## outside 0.257 0.169 Inf 0.0577 0.661
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scalecld(emm_loc2, adjust = "holm", Letters = letters)## location prob SE df asymp.LCL asymp.UCL .group
## inside 0.148 0.135 Inf 0.0133 0.691 a
## outside 0.257 0.169 Inf 0.0401 0.741 a
## edge 0.334 0.190 Inf 0.0612 0.795 a
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## Intervals are back-transformed from the logit scale
## P value adjustment: holm method for 3 tests
## Tests are performed on the log odds ratio scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.#sim_res2 <- simulateResiduals(model2)
#plot(sim_res2)
#testDispersion(sim_res2)
#testZeroInflation(sim_res2)
#testOutliers(sim_res2)
#check_model(model2)Model 2 figure
emm_loc2 <- emmeans(model2, ~ location, type = "response")
cld_loc2 <- cld(emm_loc2, adjust = "holm", Letters = letters)
cld_df2 <- as.data.frame(cld_loc2)
cld_df2$.group <- gsub(" ", "", cld_df2$.group)
cld_df2$prob <- cld_df2$prob * 100
cld_df2$asymp.LCL <- cld_df2$asymp.LCL * 100
cld_df2$asymp.UCL <- cld_df2$asymp.UCL * 100
my_colors <- viridis(100, option = "D")[c(10, 40, 70, 90)]
ggplot(cld_df2, aes(x = location, y = prob, fill = location)) +
geom_col(width = 0.6, color = "black") +
geom_errorbar(
aes(ymin = asymp.LCL, ymax = asymp.UCL),
width = 0.2,
linewidth = 1
) +
geom_text(
aes(y = asymp.UCL + 5, label = .group),
size = 6
) +
scale_fill_manual(values = my_colors) +
annotate("text",
x = -Inf,
y = Inf,
label = "χ² = 0.942 df = 2 p = 0.6243",
hjust = -0.1,
vjust = 1.1,
size = 5) +
labs(
x = "Patch Location",
y = "Probability of Mortality (%)"
) +
ylim(0, 100) +
theme_classic() +
theme(
legend.position = "none",
axis.title = element_text(size = 16),
axis.text = element_text(size = 14)
)
Model 3: immediate post fire seedling mortality
chpt2_r_data$mortality2_beta <- pmin(pmax(chpt2_r_data$mortality2, 0.001), 0.999)
model3 <- glmmTMB(
mortality2_beta ~ location + season + max.temp + ros + (1 | patch) + (1 | block),
data = chpt2_r_data,
family = beta_family(link = "logit")
)## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluationsummary(model3)## Family: beta ( logit )
## Formula:
## mortality2_beta ~ location + season + max.temp + ros + (1 | patch) +
## (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## -34.4 -20.1 26.2 -52.4 27
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 9.974e-01 9.987e-01
## block (Intercept) 6.269e-10 2.504e-05
## Number of obs: 36, groups: patch, 12; block, 4
##
## Dispersion parameter for beta family (): 10.1
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.706166 0.853068 0.828 0.40779
## locationedge -0.732088 0.343723 -2.130 0.03318 *
## locationoutside -0.859141 0.321892 -2.669 0.00761 **
## seasongrowing -1.151780 0.680672 -1.692 0.09062 .
## max.temp 0.002491 0.001147 2.172 0.02985 *
## ros -2.942300 2.333018 -1.261 0.20725
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(model3, type = 2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: mortality2_beta
## Chisq Df Pr(>Chisq)
## location 7.7052 2 0.02122 *
## season 2.8633 1 0.09062 .
## max.temp 4.7179 1 0.02985 *
## ros 1.5905 1 0.20725
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm_loc3 <- emmeans(model3, ~ location, type = "response")
pairs(emm_loc3, adjust = "holm")## contrast odds.ratio SE df null z.ratio p.value
## inside / edge 2.08 0.715 Inf 1 2.130 0.0664
## inside / outside 2.36 0.760 Inf 1 2.669 0.0228
## edge / outside 1.14 0.355 Inf 1 0.407 0.6841
##
## Results are averaged over the levels of: season
## P value adjustment: holm method for 3 tests
## Tests are performed on the log odds ratio scaleemm_loc3## location response SE df asymp.LCL asymp.UCL
## inside 0.805 0.0613 Inf 0.657 0.898
## edge 0.664 0.0860 Inf 0.482 0.808
## outside 0.636 0.0866 Inf 0.456 0.784
##
## Results are averaged over the levels of: season
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scalecld(emm_loc3, adjust = "holm", Letters = letters)## location response SE df asymp.LCL asymp.UCL .group
## outside 0.636 0.0866 Inf 0.416 0.810 a
## edge 0.664 0.0860 Inf 0.440 0.833 ab
## inside 0.805 0.0613 Inf 0.618 0.913 b
##
## Results are averaged over the levels of: season
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## Intervals are back-transformed from the logit scale
## P value adjustment: holm method for 3 tests
## Tests are performed on the log odds ratio scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.#sim_res3 <- simulateResiduals(model3)
#plot(sim_res3)
#testDispersion(sim_res3)
#testZeroInflation(sim_res3)
#testOutliers(sim_res3)
#check_model(model3)Model 3 figure- location
emm_loc3 <- emmeans(model3, ~ location, type = "response")
cld_loc3 <- cld(emm_loc3, adjust = "holm", Letters = letters)
cld_df3 <- as.data.frame(cld_loc3)
cld_df3$.group <- gsub(" ", "", cld_df3$.group)
cld_df3$response <- cld_df3$response * 100
cld_df3$asymp.LCL <- cld_df3$asymp.LCL * 100
cld_df3$asymp.UCL <- cld_df3$asymp.UCL * 100
my_colors <- viridis(100, option = "D")[c(10, 40, 70, 90)]
ggplot(cld_df3, aes(x = location, y = response, fill = location)) +
geom_col(width = 0.65, color = "black") +
geom_errorbar(
aes(ymin = asymp.LCL, ymax = asymp.UCL),
width = 0.2,
linewidth = 0.8
) +
geom_text(
aes(y = asymp.UCL + 5, label = .group),
size = 6
) +
scale_fill_manual(values = my_colors) +
labs(
x = "Patch Location",
y = "Probability of Mortality (%)"
) +
ylim(0, 100) +
annotate("text",
x = -Inf,
y = Inf,
label = "χ² = 7.705 df = 2 p = 0.0212",
hjust = -0.1,
vjust = 1.1,
size = 4) +
theme_classic() +
theme(
legend.position = "none",
axis.title = element_text(size = 16),
axis.text = element_text(size = 14)
)
model 3 figure- max temp
library(ggeffects)
library(ggplot2)
library(viridis)
# Predictions from beta regression model
pred3 <- ggpredict(model3, terms = "max.temp")## You are calculating adjusted predictions on the population-level (i.e.
## `type = "fixed"`) for a *generalized* linear mixed model.
## This may produce biased estimates due to Jensen's inequality. Consider
## setting `bias_correction = TRUE` to correct for this bias.
## See also the documentation of the `bias_correction` argument.pred_df <- as.data.frame(pred3)
pred_df$predicted <- pred_df$predicted * 100
pred_df$conf.low <- pred_df$conf.low * 100
pred_df$conf.high <- pred_df$conf.high * 100
ggplot(pred_df, aes(x = x, y = predicted)) +
# Solid ribbon
geom_ribbon(aes(ymin = conf.low,
ymax = conf.high),
fill = "steelblue",
alpha = 0.35,
colour = NA) +
# Solid line
geom_line(color = "steelblue",
linewidth = 1.3) +
# Labels
labs(
x = "Maximum Temperature (°C)",
y = "Predicted Seedling Mortality (%)"
) +
annotate("text",
x = max(pred_df$x),
y = min(pred_df$predicted)-10,
hjust = 1,
vjust = 0,
label = "χ² = 4.72\n df = 1\n p = 0.0298",
size = 5) +
theme_classic(base_size = 14) +
theme(
axis.title = element_text(face = "bold"),
axis.text = element_text(color = "black")
)
Model 5: severity (tree scorch), first order fire effect
###create the binary variable (1 for greater than or equal to 0.7 scorch, 0 for less than 0.7)
chpt2_r_data$tree.scorch_prop <- chpt2_r_data$tree.scorch / 100
chpt2_r_data$high_scorch <- ifelse(chpt2_r_data$tree.scorch_prop >= 0.7, 1, 0)
model5 <- glmmTMB(
high_scorch ~ season + location + ros + max.temp + (1|patch) + (1|block),
data = chpt2_r_data,
family = binomial(link = "logit")
)
summary(model5)## Family: binomial ( logit )
## Formula:
## high_scorch ~ season + location + ros + max.temp + (1 | patch) +
## (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## 51.5 64.2 -17.8 35.5 28
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 1.029 1.015
## block (Intercept) 1.098 1.048
## Number of obs: 36, groups: patch, 12; block, 4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.376955 2.890840 -0.130 0.8963
## seasongrowing 3.628621 2.168315 1.674 0.0942 .
## locationedge -2.143485 1.434037 -1.495 0.1350
## locationoutside -1.174080 1.248948 -0.940 0.3472
## ros 8.003390 8.952915 0.894 0.3714
## max.temp -0.000185 0.004011 -0.046 0.9632
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm5 <- emmeans(model5, ~ location, type = "response")
pairs(emm5, adjust = "holm")## contrast odds.ratio SE df null z.ratio p.value
## inside / edge 8.529 12.200 Inf 1 1.495 0.4050
## inside / outside 3.235 4.040 Inf 1 0.940 0.6944
## edge / outside 0.379 0.461 Inf 1 -0.797 0.6944
##
## Results are averaged over the levels of: season
## P value adjustment: holm method for 3 tests
## Tests are performed on the log odds ratio scalecld(emm5, adjust = "holm", Letters = letters)## location prob SE df asymp.LCL asymp.UCL .group
## edge 0.380 0.290 Inf 0.0314 0.921 a
## outside 0.618 0.277 Inf 0.0885 0.964 a
## inside 0.840 0.167 Inf 0.2131 0.990 a
##
## Results are averaged over the levels of: season
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## Intervals are back-transformed from the logit scale
## P value adjustment: holm method for 3 tests
## Tests are performed on the log odds ratio scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.#sim_res5 <- simulateResiduals(model5)
#plot(sim_res5)
#testDispersion(sim_res5)
#testZeroInflation(sim_res5)
#testOutliers(sim_res5)
#check_model(model5)model 5 figure
chpt2_r_data$tree.scorch_prop <- chpt2_r_data$tree.scorch / 100
chpt2_r_data$high_scorch <- ifelse(chpt2_r_data$tree.scorch_prop >= 0.7, 1, 0)
emm5 <- emmeans(model5, ~ location, type = "response")
cld5 <- cld(emm5, adjust = "holm", Letters = letters)
cld_df5 <- as.data.frame(cld5)
cld_df5$.group <- gsub(" ", "", cld_df5$.group)
cld_df5$prob <- cld_df5$prob * 100
cld_df5$asymp.LCL <- cld_df5$asymp.LCL * 100
cld_df5$asymp.UCL <- cld_df5$asymp.UCL * 100
my_colors <- viridis(100, option = "D")[c(10, 40, 70, 90)]
ggplot(cld_df5, aes(x = location, y = prob, fill = location)) +
geom_col(width = 0.65, color = "black") +
geom_errorbar(
aes(ymin = asymp.LCL, ymax = asymp.UCL),
width = 0.2,
linewidth = 0.8
) +
geom_text(
aes(y = asymp.UCL + 5, label = .group),
size = 6
) +
scale_fill_manual(values = my_colors) +
labs(
x = "Patch Location",
y = "Probability of High Tree Scorch (≥70%)"
) +
coord_cartesian(ylim = c(0, 100)) +
theme_classic() +
theme(
legend.position = "none",
axis.title = element_text(size = 16),
axis.text = element_text(size = 14)
)
Model 6: impacts on severity (char height), first order fire effect
model6 <- glmmTMB(
tree.char ~ season + location + ros + max.temp + (1 | patch ) + (1 | block),
data = chpt2_r_data,
family = Gamma(link = "log")
)## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluationsummary(model6)## Family: Gamma ( log )
## Formula:
## tree.char ~ season + location + ros + max.temp + (1 | patch) + (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## 97.8 112.0 -39.9 79.8 27
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 1.413e-01 3.759e-01
## block (Intercept) 1.930e-10 1.389e-05
## Number of obs: 36, groups: patch, 12; block, 4
##
## Dispersion estimate for Gamma family (sigma^2): 0.169
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3142579 0.4572089 -0.687 0.49187
## seasongrowing -0.6508829 0.2861865 -2.274 0.02295 *
## locationedge -0.0506114 0.1872238 -0.270 0.78691
## locationoutside -0.0698913 0.1744770 -0.401 0.68873
## ros -0.6598510 1.3904424 -0.475 0.63510
## max.temp 0.0018213 0.0006712 2.713 0.00666 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(model6, type = 2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: tree.char
## Chisq Df Pr(>Chisq)
## season 5.1726 1 0.02295 *
## location 0.1663 2 0.92023
## ros 0.2252 1 0.63510
## max.temp 7.3625 1 0.00666 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm6 <- emmeans(model6, ~ location, type = "response")
pairs(emm6, adjust = "holm")## contrast ratio SE df null z.ratio p.value
## inside / edge 1.05 0.197 Inf 1 0.270 1.0000
## inside / outside 1.07 0.187 Inf 1 0.401 1.0000
## edge / outside 1.02 0.183 Inf 1 0.108 1.0000
##
## Results are averaged over the levels of: season
## P value adjustment: holm method for 3 tests
## Tests are performed on the log scalecld(emm6, adjust = "holm", Letters = letters)## location response SE df asymp.LCL asymp.UCL .group
## outside 1.34 0.227 Inf 0.890 2.01 a
## edge 1.36 0.243 Inf 0.889 2.09 a
## inside 1.43 0.245 Inf 0.951 2.16 a
##
## Results are averaged over the levels of: season
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## Intervals are back-transformed from the log scale
## P value adjustment: holm method for 3 tests
## Tests are performed on the log scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.#sim_res6 <- simulateResiduals(model6)
#plot(sim_res6)
#testDispersion(sim_res6)
#testZeroInflation(sim_res6)
#testOutliers(sim_res6)
#check_model(model6)Model 6 figure
library(ggeffects)
library(ggplot2)
library(viridis)
# Generate predictions (holding other vars at typical values)
pred6 <- ggpredict(model6, terms = "max.temp")
# Convert to dataframe (ggeffects already is, but keeps things explicit)
pred_df <- as.data.frame(pred6)
# Plot
ggplot(pred_df, aes(x = x, y = predicted)) +
# Confidence ribbon
geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
fill = viridis(1, option = "D"),
alpha = 0.25) +
# Prediction line
geom_line(color = viridis(1, option = "D"),
linewidth = 1.3) +
# Labels
labs(
x = "Maximum Temperature (°C)",
y = "Predicted Tree Char Height (m)"
) +
# Annotation (replace with real values later)
annotate("text",
x = min(pred_df$x),
y = max(pred_df$predicted),
hjust = 0,
vjust = 1,
label = "χ² = 7.36\n df = 1\n p = 0.0067",
size = 5) +
# Theme
theme_classic(base_size = 14) +
theme(
axis.title = element_text(face = "bold"),
axis.text = element_text(color = "black")
)
model 6, season and char
library(emmeans)
library(ggplot2)
library(viridis)
library(multcomp)
# Estimated marginal means for season
emm_season6 <- emmeans(model6, ~ season, type = "response")
cld_season6 <- cld(emm_season6, adjust = "holm", Letters = letters)
cld_df6 <- as.data.frame(cld_season6)
cld_df6$.group <- gsub(" ", "", cld_df6$.group)
# No scaling needed (Gamma model → continuous response)
cld_df6$response <- cld_df6$response
cld_df6$asymp.LCL <- cld_df6$asymp.LCL
cld_df6$asymp.UCL <- cld_df6$asymp.UCL
my_colors <- viridis(100, option = "D")[c(25,60)]
# Plot
ggplot(cld_df6, aes(x = season, y = response, fill = season)) +
geom_col(width = 0.65, color = "black") +
geom_errorbar(
aes(ymin = asymp.LCL, ymax = asymp.UCL),
width = 0.2,
linewidth = 0.8
) +
geom_text(
aes(y = asymp.UCL + 0.05 * max(asymp.UCL), label = .group),
size = 6
) +
scale_fill_manual(values = my_colors) +
labs(
x = "Season",
y = "Tree Char Height (m)"
) +
annotate("text",
x = -Inf,
y = Inf,
label = "χ² = 5.173, df = 1, p = 0.0229",
hjust = -1.5,
vjust = 1.1,
size = 5) +
theme_classic() +
theme(
legend.position = "none",
axis.title = element_text(size = 16),
axis.text = element_text(size = 14)
)
combines first order fire effects figures
Model 7: impacts on delayed planted seedling mortality: 2nd order fire effect
model7 <- glmmTMB(
mortality3 ~ location * season + difn + ros + max.temp,
data = chpt2_r_data,
family = ordbeta(link = "logit")
)
summary(model7)## Family: ordbeta ( logit )
## Formula: mortality3 ~ location * season + difn + ros + max.temp
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## 58.3 75.5 -17.2 34.3 19
##
##
## Dispersion parameter for ordbeta family (): 10.1
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.315449 2.639746 0.877 0.3804
## locationedge 0.821250 0.795584 1.032 0.3020
## locationoutside 0.019961 0.832481 0.024 0.9809
## seasongrowing -0.177888 0.887174 -0.200 0.8411
## difn -0.022757 0.027428 -0.830 0.4067
## ros 1.733842 2.618028 0.662 0.5078
## max.temp -0.003492 0.002680 -1.303 0.1926
## locationedge:seasongrowing -2.751582 1.358321 -2.026 0.0428 *
## locationoutside:seasongrowing -0.115826 1.035156 -0.112 0.9109
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(model7, type = "III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: mortality3
## Chisq Df Pr(>Chisq)
## (Intercept) 0.7694 1 0.3804
## location 2.2292 2 0.3281
## season 0.0402 1 0.8411
## difn 0.6884 1 0.4067
## ros 0.4386 1 0.5078
## max.temp 1.6978 1 0.1926
## location:season 4.9973 2 0.0822 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm7 <- emmeans(model7, ~ season * location, type = "response")
pairs(emm7, adjust = "FDR")## contrast odds.ratio SE df null z.ratio p.value
## dormant inside / growing inside 1.195 1.060 Inf 1 0.201 0.9677
## dormant inside / dormant edge 0.440 0.350 Inf 1 -1.032 0.5033
## dormant inside / growing edge 8.234 12.500 Inf 1 1.392 0.3343
## dormant inside / dormant outside 0.980 0.816 Inf 1 -0.024 0.9809
## dormant inside / growing outside 1.315 1.210 Inf 1 0.298 0.9677
## growing inside / dormant edge 0.368 0.268 Inf 1 -1.373 0.3343
## growing inside / growing edge 6.892 9.160 Inf 1 1.453 0.3343
## growing inside / dormant outside 0.820 0.599 Inf 1 -0.271 0.9677
## growing inside / growing outside 1.101 0.867 Inf 1 0.122 0.9677
## dormant edge / growing edge 18.718 22.800 Inf 1 2.405 0.2423
## dormant edge / dormant outside 2.228 1.330 Inf 1 1.346 0.3343
## dormant edge / growing outside 2.989 1.960 Inf 1 1.671 0.3343
## growing edge / dormant outside 0.119 0.133 Inf 1 -1.904 0.3343
## growing edge / growing outside 0.160 0.167 Inf 1 -1.754 0.3343
## dormant outside / growing outside 1.341 0.792 Inf 1 0.497 0.9285
##
## P value adjustment: fdr method for 15 tests
## Tests are performed on the log odds ratio scalecld(emm7, adjust = "FDR", Letters = letters)## season location response SE df asymp.LCL asymp.UCL .group
## growing edge 0.044 0.0486 Inf 0.00217 0.493 a
## growing outside 0.224 0.0864 Inf 0.07191 0.517 a
## growing inside 0.241 0.1160 Inf 0.05643 0.627 a
## dormant inside 0.275 0.1450 Inf 0.05273 0.720 a
## dormant outside 0.279 0.0864 Inf 0.11058 0.546 a
## dormant edge 0.463 0.1200 Inf 0.19462 0.754 a
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 6 estimates
## Intervals are back-transformed from the logit scale
## P value adjustment: fdr method for 15 tests
## Tests are performed on the log odds ratio scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.#sim_res7 <- simulateResiduals(model7)
#plot(sim_res7)
#testDispersion(sim_res7)
#testZeroInflation(sim_res7)
#testOutliers(sim_res7)
#check_model(model7)model 7 fig
library(ggeffects)
library(ggplot2)
library(viridis)
library(grid) # for unit()
# ----------------------------
# Predictions
# ----------------------------
pred7 <- ggpredict(model7, terms = c("location", "season"))
pred_df <- as.data.frame(pred7)
# percent scale
pred_df$predicted <- pred_df$predicted * 100
pred_df$conf.low <- pred_df$conf.low * 100
pred_df$conf.high <- pred_df$conf.high * 100
pred_df$season <- pred_df$group
# ----------------------------
# VIRIDIS (EXPLICIT CONTROL)
# ----------------------------
base_cols <- viridis(10, option = "D")
loc_levels <- sort(unique(pred_df$x))
# pick specific viridis positions
cols <- c(
inside = base_cols[3],
edge = base_cols[6],
outside = base_cols[9]
)
# ensure names match your factor levels
cols <- cols[loc_levels]
# ----------------------------
# CHI-SQUARE TEXT (dormant only)
# ----------------------------
chi_df <- data.frame(
season = "growing",
x = 1.55, # centered between 3 bars (adjust if needed)
y = max(pred_df$predicted) * 1.45,
label = "χ² = 4.997, df = 2, p = 0.0822"
)
# ----------------------------
# PLOT
# ----------------------------
ggplot(pred_df,
aes(x = x,
y = predicted,
fill = x)) +
# bars
geom_col(width = 0.7, color = "black") +
# error bars
geom_errorbar(aes(ymin = conf.low,
ymax = conf.high),
width = 0.2,
linewidth = 0.8) +
# chi-square annotation (ONLY dormant facet)
geom_text(
data = chi_df,
aes(x = x, y = y, label = label),
inherit.aes = FALSE,
size = 4,
fontface = "italic"
) +
# facets
facet_wrap(
~season,
nrow = 1,
labeller = labeller(season = c(
growing = "Growing",
dormant = "Dormant"
))
) +
# viridis manual colors
scale_fill_manual(values = cols) +
labs(
x = "Location",
y = "Predicted Mortality (%)"
) +
theme_classic(base_size = 14) +
theme(
legend.position = "none",
panel.spacing = unit(1.5, "lines"),
panel.border = element_rect(color = "black", fill = NA),
strip.background = element_rect(color = "black", fill = "white"),
strip.text = element_text(size = 14, face = "bold"),
axis.title.y = element_text(size = 14),
axis.text.x = element_text(size = 12)
)
model 9: planted seedling growth post fire, 2nd order fire effect
model9_data <- chpt2_r_data[!is.na(chpt2_r_data$rcd20), ]
model9_data$rcd20_ig <- model9_data$rcd20 + 1e-6
model9 <- glmmTMB(
rcd20_ig ~ location + max.temp + ros + tree.scorch + tree.char +
difn + green.needles + season + (1 | patch) + (1 | block),
data = model9_data,
family = gaussian(link = "log"),
control = glmmTMBControl(optimizer = optim, optArgs = list(method = "BFGS")
))
summary(model9)## Family: gaussian ( log )
## Formula:
## rcd20_ig ~ location + max.temp + ros + tree.scorch + tree.char +
## difn + green.needles + season + (1 | patch) + (1 | block)
## Data: model9_data
##
## AIC BIC logLik -2*log(L) df.resid
## 149.3 167.1 -61.6 123.3 16
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 0.1422 0.377
## block (Intercept) 1.1316 1.064
## Residual 1.7531 1.324
## Number of obs: 29, groups: patch, 12; block, 4
##
## Dispersion estimate for gaussian family (sigma^2): 1.75
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1432882 1.2428518 -0.920 0.35763
## locationedge 0.0279611 0.2126652 0.132 0.89540
## locationoutside -0.0265196 0.2078471 -0.128 0.89847
## max.temp 0.0018267 0.0006062 3.013 0.00258 **
## ros -1.4837444 0.9358266 -1.585 0.11285
## tree.scorch -0.0001214 0.0081062 -0.015 0.98805
## tree.char 0.1510273 0.1039283 1.453 0.14617
## difn 0.0041250 0.0117000 0.353 0.72441
## green.needles 0.0024901 0.0217064 0.115 0.90867
## seasongrowing 0.8096581 1.3640937 0.594 0.55281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(model9, type = "II")## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: rcd20_ig
## Chisq Df Pr(>Chisq)
## location 0.0629 2 0.969034
## max.temp 9.0799 1 0.002584 **
## ros 2.5138 1 0.112855
## tree.scorch 0.0002 1 0.988047
## tree.char 2.1118 1 0.146172
## difn 0.1243 1 0.724414
## green.needles 0.0132 1 0.908669
## season 0.3523 1 0.552813
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#emm9 <- emmeans(model9, ~ location, type = "response")
#pairs(emm9, adjust = "FDR")
#cld(emm9, adjust = "FDR", Letters = letters)
#sim_res9 <- simulateResiduals(model9)
#plot(sim_res9)
#testDispersion(sim_res9)
#testZeroInflation(sim_res9)
#testOutliers(sim_res9)
#check_model(model9)model 9 figure
library(viridis)
# Predictions from model9
pred9 <- ggpredict(model9, terms = "max.temp")
pred_df <- as.data.frame(pred9)
# Pick a nice mid-range viridis color
col <- viridis(5, option = "D")[3]
ggplot(pred_df, aes(x = x, y = predicted)) +
# Solid ribbon (CI)
geom_ribbon(aes(ymin = conf.low,
ymax = conf.high),
fill = col,
alpha = 0.35,
colour = NA) +
# Solid line
geom_line(color = col,
linewidth = 1.3) +
# Labels
labs(
x = "Maximum Temperature (°C)",
y = "Predicted Change Root Collar Diameter (mm)"
) +
annotate("text",
x = min(pred_df$x),
y = max(pred_df$predicted),
hjust = 0,
vjust = -3,
label = "χ² = 9.079 \n df = 1 \n p = 0.0026",
size = 5) +
theme_classic(base_size = 12) +
theme(
axis.title = element_text(face = "bold"),
axis.text.x = element_text(color = "black", size = 14),
axis.text.y = element_text(color = "black", size = 12)
)
combined 7 and 9 figure