Cogongrass & Fire Arrival Time: Cover x Distribution and Climate Scenarios (Objectives 1-2)
Alan Ivory
2026-07-16
1 Purpose
This analysis asks how cogongrass cover, its spatial arrangement, and their interaction shape predicted fire arrival time, and how those effects depend on background climate and fuel-load conditions. We do this with six scenario models (three fixing fuel moisture at Low/Mean/High while fuel load varies; three fixing fuel load at Low/Mean/High while moisture varies). Cover, arrangement, and the cover x arrangement interaction are estimated within each scenario, so no separate pooled cover x arrangement model is needed.
Reclassification: the 0% and 100% cover landscapes have no defined spatial arrangement; we treat complete dominance as the limiting case of clumping and fold them into the Clumped arrangement.
Reporting: because the simulations are a deterministic full factorial, every term is “significant”; results are reported with the estimate (or χ²), its standard error, and the p-value, and interpreted through effect magnitudes rather than significance.
2 Load & Aggregate to the Scenario Level
con <- dbConnect(duckdb())
scenario_dt <- dbGetQuery(con, sprintf("
SELECT cover_percent, distribution,
inv_moist_live, inv_moist_dead, ni_moist_live, ni_moist_dead,
inv_fuel_live, inv_fuel_dead, ni_fuel_live, ni_fuel_dead,
AVG(arrival_mean) AS arrival_avg,
STDDEV(arrival_mean) AS arrival_sd,
COUNT(*) AS n_ign
FROM read_csv_auto('%s')
GROUP BY ALL
", input_path))
dbDisconnect(con, shutdown = TRUE)
setDT(scenario_dt)
cat("Scenarios:", nrow(scenario_dt), "\n")## Scenarios: 190269print(table(scenario_dt$n_ign))##
## 50
## 1902692.1 Sanity checks
stopifnot(all(scenario_dt$n_ign == 50))
stopifnot(all(scenario_dt$arrival_avg > 0))
stopifnot(sum(is.na(scenario_dt$arrival_avg)) == 0)
cover_for_na <- scenario_dt[is.na(distribution), unique(cover_percent)]
cat("cover_percent with distribution == NA (folded into Clump):",
paste(cover_for_na, collapse = ", "), "\n")## cover_percent with distribution == NA (folded into Clump): 100, 0stopifnot(setequal(cover_for_na, range(scenario_dt$cover_percent)))2.2 Factor & contrast set-up
# Fold the 0%/100% (arrangement-free) landscapes into "Clump"
scenario_dt[, distribution := fifelse(is.na(distribution), "Clump", distribution)]
scenario_dt[, distribution := factor(distribution,
levels = c("Random", "ModClump", "Clump"))]
mf_cols <- c("inv_moist_live","inv_moist_dead","ni_moist_live","ni_moist_dead",
"inv_fuel_live","inv_fuel_dead","ni_fuel_live","ni_fuel_dead")
inv_terms <- c("inv_moist_live","inv_moist_dead","inv_fuel_live","inv_fuel_dead")
ni_terms <- c("ni_moist_live","ni_moist_dead","ni_fuel_live","ni_fuel_dead")
scenario_dt[, (mf_cols) := lapply(.SD, factor, levels = c("L","H","M")), .SDcols = mf_cols]
scenario_dt[, cover_c := cover_percent - 50]
options(contrasts = c("contr.sum", "contr.poly"))
named_contr_sum <- function(lvls) {
n <- length(lvls); cm <- contr.sum(n)
colnames(cm) <- lvls[1:(n - 1)]; cm
}
contrasts(scenario_dt$distribution) <- named_contr_sum(levels(scenario_dt$distribution))
for (col in mf_cols) contrasts(scenario_dt[[col]]) <- named_contr_sum(levels(scenario_dt[[col]]))
cat("Rows per arrangement (Clump includes folded 0/100):\n")## Rows per arrangement (Clump includes folded 0/100):print(table(scenario_dt$distribution))##
## Random ModClump Clump
## 59049 59049 72171fuel_labels <- c(
inv_moist_live = "Invaded Live Fuel Moisture", inv_moist_dead = "Invaded Dead Fuel Moisture",
inv_fuel_live = "Invaded Live Fuel Load", inv_fuel_dead = "Invaded Dead Fuel Load",
ni_moist_live = "Non-Invaded Live Fuel Moisture", ni_moist_dead = "Non-Invaded Dead Fuel Moisture",
ni_fuel_live = "Non-Invaded Live Fuel Load", ni_fuel_dead = "Non-Invaded Dead Fuel Load")
arr_labels <- c(Random = "Random", ModClump = "Moderately Clumped", Clump = "Clumped")3 Cover, Arrangement, and Fuel Conditions Across Six Scenarios
Three scenarios fix all four moisture terms at one quantile (drought = Low, moderate = Mean, high = High) while fuel load varies; three fix all four fuel-load terms (Low/Mean/High) while moisture varies. Each subset also contains its 0%/100% Clumped rows. Within each scenario we model arrival time as cover x arrangement plus the four still-varying fuel factors, and read the result on its own terms.
moist_terms <- c("inv_moist_live","inv_moist_dead","ni_moist_live","ni_moist_dead")
load_terms <- c("inv_fuel_live","inv_fuel_dead","ni_fuel_live","ni_fuel_dead")
fit_scenario <- function(fixed_terms, fixed_level, varying_terms, label) {
keep <- Reduce(`&`, lapply(fixed_terms, function(cn) scenario_dt[[cn]] == fixed_level))
d <- scenario_dt[keep]
d[, distribution := droplevels(distribution)]
contrasts(d$distribution) <- named_contr_sum(levels(d$distribution))
form <- as.formula(paste0(
"arrival_avg ~ cover_c * distribution + ", paste(varying_terms, collapse = " + ")))
m <- glm(form, data = d, family = Gamma(link = "log"))
cs <- summary(m)$coefficients["cover_c", ] # Estimate, SE, t, p
beta <- cs["Estimate"]; se <- cs["Std. Error"]
pct10 <- (exp(10 * beta) - 1) * 100
pct10_se <- abs(10 * exp(10 * beta) * se) * 100 # delta-method SE on % scale
list(label = label, model = m, data = d, varying = varying_terms,
n = nrow(d), n_params = length(coef(m)),
cover_beta = unname(beta), cover_se = unname(se), cover_p = unname(cs["Pr(>|t|)"]),
pct10 = unname(pct10), pct10_se = unname(pct10_se))
}
scenarios <- list(
fit_scenario(moist_terms, "L", load_terms, "Low Moisture"),
fit_scenario(moist_terms, "M", load_terms, "Moderate Moisture"),
fit_scenario(moist_terms, "H", load_terms, "High Moisture"),
fit_scenario(load_terms, "L", moist_terms, "Low Fuel Load"),
fit_scenario(load_terms, "M", moist_terms, "Mean Fuel Load"),
fit_scenario(load_terms, "H", moist_terms, "High Fuel Load"))3.1 Dispersion diagnostics (DHARMa)
disp_tbl <- rbindlist(lapply(scenarios, function(s) {
sr <- simulateResiduals(s$model, n = 250, seed = 1)
td <- testDispersion(sr, plot = FALSE)
data.table(scenario = s$label, dispersion = round(unname(td$statistic), 3),
disp_p = signif(td$p.value, 3))
}))
print(disp_tbl)## scenario dispersion disp_p
## <char> <num> <num>
## 1: Low Moisture 1.191 0.000
## 2: Moderate Moisture 1.168 0.000
## 3: High Moisture 0.994 0.928
## 4: Low Fuel Load 1.264 0.000
## 5: Mean Fuel Load 1.733 0.000
## 6: High Fuel Load 1.377 0.0003.2 Cover effect (per scenario)
o_summary <- rbindlist(lapply(scenarios, function(s) data.table(
scenario = s$label,
n = s$n,
cover_pct_per10 = round(s$pct10, 2),
cover_pct_se = round(s$pct10_se, 2),
cover_p = signif(s$cover_p, 3)))) # Wald p for the cover slope
print(o_summary)## scenario n cover_pct_per10 cover_pct_se cover_p
## <char> <int> <num> <num> <num>
## 1: Low Moisture 2349 -7.46 0.07 0.00e+00
## 2: Moderate Moisture 2349 -5.73 0.08 0.00e+00
## 3: High Moisture 2349 -8.63 0.17 0.00e+00
## 4: Low Fuel Load 2349 -7.45 0.20 2.11e-230
## 5: Mean Fuel Load 2349 -8.16 0.15 0.00e+00
## 6: High Fuel Load 2349 -7.85 0.14 0.00e+00fwrite(o_summary, file.path(out_dir, "scenario_summary.csv"))3.3 Figure 1 — predicted arrival time by cover across scenarios
curves <- rbindlist(lapply(scenarios, function(s) {
ref <- setNames(as.list(rep("M", length(s$varying))), s$varying)
emm <- emmeans(s$model, ~ cover_c,
at = c(list(cover_c = seq(-40, 40, 5), distribution = "ModClump"), ref),
type = "response")
d <- as.data.table(as.data.frame(emm))
lcl <- grep("LCL$|lower\\.CL$", names(d), value = TRUE)[1]
ucl <- grep("UCL$|upper\\.CL$", names(d), value = TRUE)[1]
data.table(scenario = s$label, cover = d$cover_c + 50, response = d$response,
LCL = d[[lcl]], UCL = d[[ucl]],
family = fifelse(grepl("moisture", s$label, ignore.case = TRUE),
"Moisture scenarios (fuel load varying)",
"Fuel-load scenarios (moisture varying)"))
}))
ggplot(curves, aes(cover, response, color = scenario, fill = scenario)) +
geom_line(linewidth = 1) +
geom_ribbon(aes(ymin = LCL, ymax = UCL), alpha = 0.12, color = NA) +
facet_wrap(~ family, scales = "free_y") +
labs(x = "Cogongrass cover (%)", y = "Predicted arrival time (min)",
color = NULL, fill = NULL) +
theme_minimal(base_size = 12) + theme(legend.position = "bottom")
ggsave(file.path(out_dir, "Fig1_scenario_cover_curves.png"), width = 10, height = 6, dpi = 300)3.4 Figure 2 — how cover x arrangement interact within each scenario
Spatial-arrangement differences are the differences among the arrangement-specific cover slopes. We plot each arrangement’s cover slope (% change in arrival time per +10% cover, 95% CI) within each scenario and annotate Tukey-adjusted compact-letter groupings: arrangements that share a letter do not differ (p >= 0.05).
# Compact-letter display for 3 arrangements, from their pairwise Tukey p-values
cld3 <- function(p_RM, p_RC, p_MC) {
items <- c("Random", "ModClump", "Clump")
nd <- c(RM = p_RM >= 0.05, RC = p_RC >= 0.05, MC = p_MC >= 0.05)
if (all(nd)) return(setNames(rep("a", 3), items))
cliques <- list()
if (nd["RM"]) cliques <- c(cliques, list(c("Random", "ModClump")))
if (nd["RC"]) cliques <- c(cliques, list(c("Random", "Clump")))
if (nd["MC"]) cliques <- c(cliques, list(c("ModClump", "Clump")))
covered <- unique(unlist(cliques))
for (it in items) if (!(it %in% covered)) cliques <- c(cliques, list(it))
cliques <- cliques[order(sapply(cliques, function(cl) min(match(cl, items))))]
lab <- setNames(rep("", 3), items)
for (i in seq_along(cliques)) for (it in cliques[[i]]) lab[it] <- paste0(lab[it], letters[i])
lab
}
slope_by_scn <- rbindlist(lapply(scenarios, function(s) {
emt <- emtrends(s$model, ~ distribution, var = "cover_c")
d <- as.data.table(as.data.frame(emt))
pr <- as.data.frame(pairs(emt)) # Tukey-adjusted
# exact contrast match (avoid "Clump" matching the substring in "ModClump")
getp <- function(a, b) pr$p.value[pr$contrast %in% c(paste(a, "-", b), paste(b, "-", a))][1]
lab <- cld3(getp("Random", "ModClump"), getp("Random", "Clump"), getp("ModClump", "Clump"))
d[, `:=`(scenario = s$label,
letter = lab[as.character(distribution)],
pct = (exp(10 * cover_c.trend) - 1) * 100,
pct_lo = (exp(10 * (cover_c.trend - 1.96 * SE)) - 1) * 100,
pct_hi = (exp(10 * (cover_c.trend + 1.96 * SE)) - 1) * 100)]
d
}))
scn_levels <- c("Low Moisture","Moderate Moisture","High Moisture",
"Low Fuel Load","Mean Fuel Load","High Fuel Load")
slope_by_scn[, scenario := factor(scenario, levels = rev(scn_levels))]
slope_by_scn[, distribution := factor(arr_labels[as.character(distribution)], levels = arr_labels)]
arr_cols <- c(Random = "#F2A900", `Moderately Clumped` = "#22884C", Clumped = "#002657")
ggplot(slope_by_scn, aes(x = pct, y = scenario, color = distribution)) +
geom_pointrange(aes(xmin = pct_lo, xmax = pct_hi),
position = position_dodge(width = 0.6), size = 0.5) +
geom_text(aes(x = pct_hi, label = letter), position = position_dodge(width = 0.6),
hjust = -0.6, size = 3.2, fontface = "bold", show.legend = FALSE) +
scale_color_manual(values = arr_cols, name = "Arrangement") +
labs(x = "Change in Predicted Arrival Time Per +10% Cover (%)",
y = NULL) +
theme_minimal(base_size = 14) +
theme(panel.grid.major.y = element_blank(), legend.position = c(0.85, 0.22))
ggsave(file.path(out_dir, "Fig2_interaction_by_scenario.png"), width = 9, height = 6, dpi = 300)3.5 Pairwise arrangement contrasts (which arrangements differ)
These Tukey-adjusted pairwise contrasts of the arrangement-specific cover slopes are the sole test of which arrangements differ within each scenario (no omnibus interaction test is used); the compact-letter groupings in Figure 2 are derived from them. Tukey adjustment controls the family-wise error rate across the three comparisons.
pairwise_arr <- rbindlist(lapply(scenarios, function(s) {
pr <- as.data.table(as.data.frame(
pairs(emtrends(s$model, ~ distribution, var = "cover_c")))) # Tukey by default
pr[, scenario := s$label]
pr
}), fill = TRUE)
ratio_col <- grep("ratio", names(pairwise_arr), value = TRUE)[1]
print(pairwise_arr[, .(scenario, contrast,
slope_diff = round(estimate, 6), SE = round(SE, 6),
stat = round(get(ratio_col), 2),
p_tukey = signif(p.value, 3))])## scenario contrast slope_diff SE stat p_tukey
## <char> <char> <num> <num> <num> <num>
## 1: Low Moisture Random - ModClump -0.001383 0.000197 -7.01 1.10e-11
## 2: Low Moisture Random - Clump -0.001072 0.000173 -6.18 2.29e-09
## 3: Low Moisture ModClump - Clump 0.000311 0.000173 1.79 1.72e-01
## 4: Moderate Moisture Random - ModClump -0.001061 0.000220 -4.81 4.74e-06
## 5: Moderate Moisture Random - Clump -0.000835 0.000194 -4.31 5.04e-05
## 6: Moderate Moisture ModClump - Clump 0.000226 0.000194 1.16 4.75e-01
## 7: High Moisture Random - ModClump -0.001051 0.000483 -2.17 7.59e-02
## 8: High Moisture Random - Clump -0.000106 0.000425 -0.25 9.66e-01
## 9: High Moisture ModClump - Clump 0.000945 0.000425 2.22 6.75e-02
## 10: Low Fuel Load Random - ModClump -0.000783 0.000565 -1.39 3.49e-01
## 11: Low Fuel Load Random - Clump -0.000096 0.000497 -0.19 9.80e-01
## 12: Low Fuel Load ModClump - Clump 0.000687 0.000497 1.38 3.50e-01
## 13: Mean Fuel Load Random - ModClump -0.001267 0.000446 -2.84 1.26e-02
## 14: Mean Fuel Load Random - Clump -0.000505 0.000392 -1.29 4.01e-01
## 15: Mean Fuel Load ModClump - Clump 0.000761 0.000392 1.94 1.27e-01
## 16: High Fuel Load Random - ModClump -0.001411 0.000405 -3.48 1.46e-03
## 17: High Fuel Load Random - Clump -0.000662 0.000356 -1.86 1.51e-01
## 18: High Fuel Load ModClump - Clump 0.000749 0.000356 2.10 8.95e-02fwrite(pairwise_arr, file.path(out_dir, "pairwise_arrangement_slopes.csv"))3.6 Invaded vs. non-invaded fuel impact within each scenario
inv_eff <- rbindlist(lapply(scenarios, function(s) {
vars_LH <- intersect(s$varying, c(inv_terms, ni_terms))
cmp <- avg_comparisons(s$model,
variables = setNames(lapply(vars_LH, function(v) c("L","H")), vars_LH),
type = "response")
d <- as.data.table(cmp)
d[, `:=`(scenario = s$label,
type = fifelse(term %in% inv_terms, "Invaded", "Non-invaded"),
label = fuel_labels[term])]
d[, .(scenario, label, type,
estimate = round(estimate, 2), std.error = round(std.error, 2),
conf.low = round(conf.low, 2), conf.high = round(conf.high, 2),
p.value = signif(p.value, 3))]
}))
print(inv_eff) # signed per-factor effects (detail behind Table 3)## scenario label type estimate
## <char> <char> <char> <num>
## 1: Low Moisture Invaded Dead Fuel Load Invaded -11.01
## 2: Low Moisture Invaded Live Fuel Load Invaded -1.46
## 3: Low Moisture Non-Invaded Dead Fuel Load Non-invaded -6.67
## 4: Low Moisture Non-Invaded Live Fuel Load Non-invaded -0.66
## 5: Moderate Moisture Invaded Dead Fuel Load Invaded -21.77
## 6: Moderate Moisture Invaded Live Fuel Load Invaded 9.67
## 7: Moderate Moisture Non-Invaded Dead Fuel Load Non-invaded -8.63
## 8: Moderate Moisture Non-Invaded Live Fuel Load Non-invaded 0.51
## 9: High Moisture Invaded Dead Fuel Load Invaded -123.01
## 10: High Moisture Invaded Live Fuel Load Invaded 45.69
## 11: High Moisture Non-Invaded Dead Fuel Load Non-invaded -59.72
## 12: High Moisture Non-Invaded Live Fuel Load Non-invaded 16.02
## 13: Low Fuel Load Invaded Dead Fuel Moisture Invaded 50.02
## 14: Low Fuel Load Invaded Live Fuel Moisture Invaded 27.83
## 15: Low Fuel Load Non-Invaded Dead Fuel Moisture Non-invaded 75.01
## 16: Low Fuel Load Non-Invaded Live Fuel Moisture Non-invaded 3.39
## 17: Mean Fuel Load Invaded Dead Fuel Moisture Invaded 24.91
## 18: Mean Fuel Load Invaded Live Fuel Moisture Invaded 16.18
## 19: Mean Fuel Load Non-Invaded Dead Fuel Moisture Non-invaded 40.59
## 20: Mean Fuel Load Non-Invaded Live Fuel Moisture Non-invaded 11.31
## 21: High Fuel Load Invaded Dead Fuel Moisture Invaded 22.13
## 22: High Fuel Load Invaded Live Fuel Moisture Invaded 12.01
## 23: High Fuel Load Non-Invaded Dead Fuel Moisture Non-invaded 30.15
## 24: High Fuel Load Non-Invaded Live Fuel Moisture Non-invaded 4.22
## scenario label type estimate
## std.error conf.low conf.high p.value
## <num> <num> <num> <num>
## 1: 0.18 -11.37 -10.65 0.000000e+00
## 2: 0.18 -1.81 -1.10 3.900000e-16
## 3: 0.18 -7.02 -6.32 2.630000e-302
## 4: 0.18 -1.01 -0.31 2.310000e-04
## 5: 0.31 -22.38 -21.15 0.000000e+00
## 6: 0.31 9.07 10.27 1.890000e-216
## 7: 0.31 -9.23 -8.03 1.590000e-175
## 8: 0.30 -0.09 1.10 9.630000e-02
## 9: 2.26 -127.44 -118.57 0.000000e+00
## 10: 2.08 41.62 49.76 2.350000e-107
## 11: 2.03 -63.71 -55.74 1.350000e-189
## 12: 2.01 12.08 19.96 1.650000e-15
## 13: 1.50 47.08 52.96 6.900000e-244
## 14: 1.37 25.15 30.51 2.140000e-92
## 15: 1.60 71.87 78.16 0.000000e+00
## 16: 1.36 0.72 6.06 1.290000e-02
## 17: 0.80 23.35 26.47 2.360000e-215
## 18: 0.76 14.70 17.66 8.510000e-102
## 19: 0.85 38.93 42.25 0.000000e+00
## 20: 0.77 9.80 12.82 1.300000e-48
## 21: 0.58 20.99 23.26 5.160022e-319
## 22: 0.55 10.94 13.08 3.030000e-107
## 23: 0.60 28.97 31.33 0.000000e+00
## 24: 0.55 3.13 5.31 2.670000e-14
## std.error conf.low conf.high p.valuefwrite(inv_eff, file.path(out_dir, "invaded_effects_by_scenario.csv"))3.7 Invaded vs. non-invaded contrasts (Table 3)
We test whether each invaded fuel factor’s effect differs from its non-invaded counterpart by contrasting their signed average marginal effects (matched dead and live components) with the delta method applied to the comparisons’ joint covariance. These contrasts are Table 3 of the main text.
inv_ni_contrasts <- rbindlist(lapply(scenarios, function(s) {
rbindlist(lapply(c("dead", "live"), function(ld) {
fi <- s$varying[grepl(paste0("_", ld, "$"), s$varying) & grepl("^inv_", s$varying)]
fn <- s$varying[grepl(paste0("_", ld, "$"), s$varying) & grepl("^ni_", s$varying)]
cmp <- avg_comparisons(s$model,
variables = setNames(list(c("L","H"), c("L","H")), c(fi, fn)),
type = "response")
V <- vcov(cmp) # joint covariance of the 2 effects
Lc <- ifelse(cmp$term == fi, 1, ifelse(cmp$term == fn, -1, 0))
diff <- sum(Lc * cmp$estimate) # invaded - non-invaded (delta method)
se <- sqrt(as.numeric(t(Lc) %*% V %*% Lc))
data.table(scenario = s$label, component = ld,
invaded = round(cmp$estimate[cmp$term == fi], 2),
non_invaded = round(cmp$estimate[cmp$term == fn], 2),
difference = round(diff, 2), SE = round(se, 2),
p = signif(2 * pnorm(-abs(diff / se)), 3))
}))
}))
print(inv_ni_contrasts)## scenario component invaded non_invaded difference SE p
## <char> <char> <num> <num> <num> <num> <num>
## 1: Low Moisture dead -11.01 -6.67 -4.34 0.25 3.12e-66
## 2: Low Moisture live -1.46 -0.66 -0.80 0.25 1.59e-03
## 3: Moderate Moisture dead -21.77 -8.63 -13.14 0.43 1.87e-203
## 4: Moderate Moisture live 9.67 0.51 9.16 0.43 2.15e-99
## 5: High Moisture dead -123.01 -59.72 -63.29 2.85 3.83e-109
## 6: High Moisture live 45.69 16.02 29.68 2.87 4.74e-25
## 7: Low Fuel Load dead 50.02 75.01 -24.99 1.99 4.85e-36
## 8: Low Fuel Load live 27.83 3.39 24.44 1.92 5.50e-37
## 9: Mean Fuel Load dead 24.91 40.59 -15.68 1.10 4.16e-46
## 10: Mean Fuel Load live 16.18 11.31 4.87 1.07 5.20e-06
## 11: High Fuel Load dead 22.13 30.15 -8.02 0.79 3.62e-24
## 12: High Fuel Load live 12.01 4.22 7.79 0.77 8.33e-24fwrite(inv_ni_contrasts, file.path(out_dir, "invaded_vs_noninvaded_contrasts.csv"))4 Verification & Saved Outputs
stopifnot(nrow(scenario_dt) == 190269)
cat("Arrangement counts (Clump includes folded 0/100):\n")## Arrangement counts (Clump includes folded 0/100):print(table(scenario_dt$distribution))##
## Random ModClump Clump
## 59049 59049 72171cat("\nCover shortens arrival in every scenario (all cover_pct_per10 < 0):\n")##
## Cover shortens arrival in every scenario (all cover_pct_per10 < 0):print(o_summary[, .(scenario, cover_pct_per10, cover_pct_se, cover_p)])## scenario cover_pct_per10 cover_pct_se cover_p
## <char> <num> <num> <num>
## 1: Low Moisture -7.46 0.07 0.00e+00
## 2: Moderate Moisture -5.73 0.08 0.00e+00
## 3: High Moisture -8.63 0.17 0.00e+00
## 4: Low Fuel Load -7.45 0.20 2.11e-230
## 5: Mean Fuel Load -8.16 0.15 0.00e+00
## 6: High Fuel Load -7.85 0.14 0.00e+00cat("\nArrangement differences are tested pairwise (Tukey); significant contrasts:\n")##
## Arrangement differences are tested pairwise (Tukey); significant contrasts:print(pairwise_arr[p.value < 0.05, .(scenario, contrast, p_tukey = signif(p.value, 3))])## scenario contrast p_tukey
## <char> <char> <num>
## 1: Low Moisture Random - ModClump 1.10e-11
## 2: Low Moisture Random - Clump 2.29e-09
## 3: Moderate Moisture Random - ModClump 4.74e-06
## 4: Moderate Moisture Random - Clump 5.04e-05
## 5: Mean Fuel Load Random - ModClump 1.26e-02
## 6: High Fuel Load Random - ModClump 1.46e-03saveRDS(setNames(lapply(scenarios, `[[`, "model"),
sapply(scenarios, `[[`, "label")),
file.path(out_dir, "scenario_models.rds"))
cat("Saved models, tables, and figures to:\n", out_dir, "\n")## Saved models, tables, and figures to:
## C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/02_FuelSpatial/08_ArrivalTime_Analysis/Outputs5 Session Info
sessionInfo()## R version 4.5.3 (2026-03-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] marginaleffects_0.32.0 moments_0.14.1 patchwork_1.3.2
## [4] emmeans_1.11.2-8 DHARMa_0.4.7 lubridate_1.9.4
## [7] forcats_1.0.0 stringr_1.5.2 dplyr_1.1.4
## [10] purrr_1.1.0 readr_2.1.5 tidyr_1.3.1
## [13] tibble_3.3.0 ggplot2_4.0.0 tidyverse_2.0.0
## [16] data.table_1.17.8 duckdb_1.5.4.2 DBI_1.2.3
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.53 bslib_0.9.0 insight_1.4.4
## [5] lattice_0.22-9 tzdb_0.5.0 vctrs_0.6.5 tools_4.5.3
## [9] Rdpack_2.6.4 generics_0.1.4 sandwich_3.1-1 pacman_0.5.1
## [13] pkgconfig_2.0.3 Matrix_1.7-4 checkmate_2.3.3 RColorBrewer_1.1-3
## [17] S7_0.2.0 lifecycle_1.0.5 compiler_4.5.3 farver_2.1.2
## [21] textshaping_1.0.3 codetools_0.2-20 htmltools_0.5.8.1 sass_0.4.10
## [25] yaml_2.3.10 pillar_1.11.1 nloptr_2.2.1 jquerylib_0.1.4
## [29] MASS_7.3-65 cachem_1.1.0 reformulas_0.4.1 boot_1.3-32
## [33] multcomp_1.4-30 nlme_3.1-168 tidyselect_1.2.1 digest_0.6.37
## [37] mvtnorm_1.3-3 stringi_1.8.7 labeling_0.4.3 splines_4.5.3
## [41] fastmap_1.2.0 grid_4.5.3 cli_3.6.5 magrittr_2.0.4
## [45] survival_3.8-6 TH.data_1.1-5 withr_3.0.2 backports_1.5.0
## [49] scales_1.4.0 timechange_0.3.0 estimability_1.5.1 rmarkdown_2.29
## [53] lme4_1.1-37 ragg_1.5.0 zoo_1.8-14 hms_1.1.3
## [57] coda_0.19-4.1 evaluate_1.0.5 knitr_1.50 rbibutils_2.3
## [61] rlang_1.1.6 Rcpp_1.1.0 xtable_1.8-4 glue_1.8.0
## [65] rstudioapi_0.17.1 minqa_1.2.8 jsonlite_2.0.0 R6_2.6.1
## [69] systemfonts_1.3.1