1 Data Background

arrival_results_first50.csv contains fire arrival-time results (minutes to reach the target cells) for 50 ignitions per landscape scenario, crossed over the full design: 8 fuel/moisture factors (invaded/non-invaded x live/dead x fuel load/moisture, each Low/Medium/High -> 3^8 = 6,561 combinations), 11 cogongrass cover levels (0-100% by 10s), and up to 3 spatial distributions (Random, ModClump, Clump; undefined at 0% and 100% cover). That is 190,269 unique landscape scenarios x 50 ignitions = 9,513,450 rows.

2 Load & Aggregate to the Scenario Level

Each landscape scenario has 50 ignition points, but ignition location doesn’t change any of the design variables (cover, distribution, fuel/moisture). Fitting a regression on all 9.5M ignition-level rows would treat 50 pseudo-replicates as independent observations for every landscape-level predictor, inflating degrees of freedom. The correct unit of analysis is the scenario (n = 190,269), with the 50 ignitions summarized into a mean (and SD) per scenario.

To speed up the aggregation step, this uses DuckDB to read the CSV and compute the scenario-level summary statistics in SQL. If you don’t have DuckDB installed, the commented-out fread alternative below does the same thing.

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: 190269
print(table(scenario_dt$n_ign))   # should all be 50
## 
##     50 
## 190269
# ---- Alternative if you don't have the duckdb package / prefer fread ----
# dt <- fread(input_path)
# group_cols <- c("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")
# scenario_dt <- dt[, .(arrival_avg = mean(arrival_mean), arrival_sd = sd(arrival_mean),
#                        n_ign = .N), by = group_cols]

2.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 values with distribution == NA:", paste(cover_for_na, collapse = ", "), "\n")
## cover_percent values with distribution == NA: 100, 0
stopifnot(setequal(cover_for_na, range(scenario_dt$cover_percent)))

2.2 Factor set-up

scenario_dt[, distribution := fifelse(is.na(distribution), "CompleteDominance", distribution)]

scenario_dt[, distribution := factor(distribution,
              levels = c("Random", "ModClump", "Clump", "CompleteDominance"))]

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")

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]]))
}

3 Exploratory

p1 <- ggplot(scenario_dt, aes(arrival_avg)) +
  geom_histogram(bins = 60, fill = "steelblue") +
  labs(title = "Distribution of scenario-mean arrival time", x = "Arrival time (min)") +
  theme_minimal()

p2 <- ggplot(scenario_dt, aes(sample = arrival_avg)) +
  stat_qq() + stat_qq_line(color = "red") +
  labs(title = "QQ plot (raw arrival time)") +
  theme_minimal()

p1 + p2

cat("Skewness:", moments::skewness(scenario_dt$arrival_avg), "\n")
## Skewness: 2.124421
cat("Mean:", mean(scenario_dt$arrival_avg), " Median:", median(scenario_dt$arrival_avg), "\n")
## Mean: 75.386  Median: 58.46665

Arrival time is strictly positive and right-skewed (mean > median) — the motivation for a Gamma GLM rather than a Gaussian-errors model.

ggplot(scenario_dt, aes(x = factor(cover_percent), y = arrival_avg, fill = distribution)) +
  geom_boxplot(outlier.size = 0.3) +
  labs(title = "Arrival time by cogongrass cover and distribution",
       x = "Cogongrass cover (%)", y = "Arrival time (min)", fill = "Distribution") +
  theme_minimal()

5 Effects of Interest

5.1 Type III tests

car::Anova(m_gamma, type = "III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: arrival_avg
##                             LR Chisq Df Pr(>Chisq)    
## cover_c                       339457  1  < 2.2e-16 ***
## distribution                   10886  3  < 2.2e-16 ***
## inv_moist_live                 68629  2  < 2.2e-16 ***
## inv_moist_dead                173802  2  < 2.2e-16 ***
## ni_moist_live                   7820  2  < 2.2e-16 ***
## ni_moist_dead                 387130  2  < 2.2e-16 ***
## inv_fuel_live                  18122  2  < 2.2e-16 ***
## inv_fuel_dead                 136092  2  < 2.2e-16 ***
## ni_fuel_live                     624  2  < 2.2e-16 ***
## ni_fuel_dead                   60751  2  < 2.2e-16 ***
## cover_c:distribution            1016  3  < 2.2e-16 ***
## cover_c:inv_moist_live         33479  2  < 2.2e-16 ***
## cover_c:inv_moist_dead         79305  2  < 2.2e-16 ***
## cover_c:ni_moist_live           3941  2  < 2.2e-16 ***
## cover_c:ni_moist_dead         199867  2  < 2.2e-16 ***
## cover_c:inv_fuel_live           9307  2  < 2.2e-16 ***
## cover_c:inv_fuel_dead          65669  2  < 2.2e-16 ***
## cover_c:ni_fuel_live             329  2  < 2.2e-16 ***
## cover_c:ni_fuel_dead           30565  2  < 2.2e-16 ***
## distribution:inv_moist_live      327  6  < 2.2e-16 ***
## distribution:inv_moist_dead      788  6  < 2.2e-16 ***
## distribution:ni_moist_live       198  6  < 2.2e-16 ***
## distribution:ni_moist_dead     10293  6  < 2.2e-16 ***
## distribution:inv_fuel_live        50  6  5.828e-09 ***
## distribution:inv_fuel_dead       632  6  < 2.2e-16 ***
## distribution:ni_fuel_live         17  6   0.008125 ** 
## distribution:ni_fuel_dead       1539  6  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.2 Cover x Distribution interaction

ref_at <- setNames(as.list(rep("M", length(mf_cols))), mf_cols)

emm_grid <- emmeans(m_gamma, ~ cover_c * distribution,
                     at = c(list(cover_c = seq(-40, 40, 10)), ref_at),
                     type = "response")
emm_df <- as.data.frame(emm_grid)
emm_df <- emm_df[emm_df$distribution != "CompleteDominance", ]
emm_df$distribution <- droplevels(emm_df$distribution)

lcl_col <- grep("LCL$|lower\\.CL$", names(emm_df), value = TRUE)[1]
ucl_col <- grep("UCL$|upper\\.CL$", names(emm_df), value = TRUE)[1]
emm_df$LCL <- emm_df[[lcl_col]]
emm_df$UCL <- emm_df[[ucl_col]]

ggplot(emm_df, aes(x = cover_c + 50, y = response, color = distribution, fill = distribution)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = LCL, ymax = UCL), alpha = 0.15, color = NA) +
  labs(title = "Predicted fire arrival time by cogongrass cover and distribution",
       x = "Cogongrass cover (%)", y = "Predicted arrival time (min)",
       color = "Distribution", fill = "Distribution") +
  theme_minimal()

ggsave(file.path(out_dir, "cover_by_distribution_effect.png"), width = 8, height = 5, dpi = 300)


emm_extremes <- emmeans(m_gamma, ~ cover_c,
                         at = c(list(cover_c = c(-50, 50), distribution = "CompleteDominance"), ref_at),
                         type = "response")
print(as.data.frame(emm_extremes))
##  cover_c response        SE     df lower.CL upper.CL
##      -50 75.50552 0.4824440 190181 74.56584 76.45704
##       50 41.75765 0.2668113 190181 41.23797 42.28389
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

5.3 Cover

gamma_summary <- summary(m_gamma)
coef_tbl <- coef(gamma_summary)
cover_row <- coef_tbl["cover_c", ]
print(cover_row)
##      Estimate    Std. Error       t value      Pr(>|t|) 
## -7.511584e-03  1.289652e-05 -5.824504e+02  0.000000e+00
beta_cover <- cover_row["Estimate"]
se_cover   <- cover_row["Std. Error"]
pct_per_10 <- (exp(10 * beta_cover) - 1) * 100
pct_lo     <- (exp(10 * (beta_cover - 1.96 * se_cover)) - 1) * 100
pct_hi     <- (exp(10 * (beta_cover + 1.96 * se_cover)) - 1) * 100
cat(sprintf("Cover coefficient: %.5f (SE = %.5f)\n", beta_cover, se_cover))
## Cover coefficient: -0.00751 (SE = 0.00001)
cat(sprintf("=> %.2f%% change in arrival time per 10-point increase in cover [%.2f%%, %.2f%%]\n",
            pct_per_10, pct_lo, pct_hi))
## => -7.24% change in arrival time per 10-point increase in cover [-7.26%, -7.21%]

5.4 Fuel moisture / fuel load main effects

fuel_terms <- c("inv_moist_live","ni_moist_live","inv_moist_dead","ni_moist_dead",
                 "inv_fuel_live","ni_fuel_live","inv_fuel_dead","ni_fuel_dead")

fuel_labels <- c(
  inv_moist_live = "Invaded Live Fuel Moisture",
  ni_moist_live  = "Non-Invaded Live Fuel Moisture",
  inv_moist_dead = "Invaded Dead Fuel Moisture",
  ni_moist_dead  = "Non-Invaded Dead Fuel Moisture",
  inv_fuel_live  = "Invaded Live Fuel Load",
  ni_fuel_live   = "Non-Invaded Live Fuel Load",
  inv_fuel_dead  = "Invaded Dead Fuel Load",
  ni_fuel_dead   = "Non-Invaded Dead Fuel Load"
)

gator_colors <- c("Invaded" = "#FA4616", "Non-invaded" = "#0021A5")
fuel_pred_list <- lapply(fuel_terms, function(v) {
  at_list <- ref_at
  at_list[[v]] <- NULL
  at_list$distribution <- "Random"
  emm <- emmeans(m_gamma, as.formula(paste("~", v)), at = at_list, type = "response")
  d <- as.data.frame(emm)
  lcl_col <- grep("LCL$|lower\\.CL$", names(d), value = TRUE)[1]
  ucl_col <- grep("UCL$|upper\\.CL$", names(d), value = TRUE)[1]
  data.table(factor = v,
             type   = fifelse(grepl("^inv_", v), "Invaded", "Non-invaded"),
             level  = d[[v]],
             response = d$response,
             LCL = d[[lcl_col]],
             UCL = d[[ucl_col]])
})
fuel_pred_dt <- rbindlist(fuel_pred_list)
fuel_pred_dt[, level := factor(level, levels = c("L", "M", "H"))]
fuel_pred_dt[, factor_label := factor(fuel_labels[factor], levels = fuel_labels[fuel_terms])]

ggplot(fuel_pred_dt, aes(x = level, y = response, color = type)) +
  geom_pointrange(aes(ymin = LCL, ymax = UCL)) +
  facet_wrap(~ factor_label, ncol = 4) +
  scale_color_manual(values = gator_colors) +
  labs(title = "Fuel moisture / fuel load main effects",
       x = "Quantile level (Low / Medium / High)", y = "Predicted Arrival Time (min)",
       color = "Fuel type") +
  theme_minimal(base_size = 10)

ggsave(file.path(out_dir, "fuel_moisture_main_effects.png"), width = 12, height = 6, dpi = 300)

5.5 Sensitivity: how much does each fuel factor move arrival time?

sensitivity_tbl <- dcast(fuel_pred_dt[level %in% c("L", "H")],
                          factor + type ~ level, value.var = "response")
setnames(sensitivity_tbl, c("L", "H"), c("arrival_L", "arrival_H"))
sensitivity_tbl[, delta_min   := arrival_H - arrival_L]     # minutes, H - L
sensitivity_tbl[, pct_change  := (arrival_H / arrival_L - 1) * 100]
sensitivity_tbl[, abs_delta   := abs(delta_min)]

print(sensitivity_tbl[order(-abs_delta), .(factor, type, arrival_L, arrival_H, delta_min, pct_change)])
##            factor        type arrival_L arrival_H  delta_min pct_change
##            <char>      <char>     <num>     <num>      <num>      <num>
## 1:  ni_moist_dead Non-invaded  46.24802  78.80974  32.561724  70.406748
## 2: inv_moist_dead     Invaded  45.73579  72.42467  26.688875  58.354464
## 3:  inv_fuel_dead     Invaded  66.21608  43.26832 -22.947759 -34.655872
## 4: inv_moist_live     Invaded  41.75481  56.60361  14.848799  35.561889
## 5:   ni_fuel_dead Non-invaded  58.28785  46.28092 -12.006932 -20.599373
## 6:  inv_fuel_live     Invaded  48.72589  56.68413   7.958230  16.332651
## 7:  ni_moist_live Non-invaded  50.27682  54.53087   4.254045   8.461244
## 8:   ni_fuel_live Non-invaded  50.92869  52.18897   1.260289   2.474615
fwrite(sensitivity_tbl, file.path(out_dir, "fuel_factor_sensitivity.csv"))


sensitivity_tbl[, factor_label := factor(fuel_labels[factor], levels = rev(fuel_labels[fuel_terms]))]

ggplot(sensitivity_tbl, aes(x = factor_label, y = delta_min, fill = type)) +
  geom_col() +
  coord_flip() +
  scale_fill_manual(values = gator_colors) +
  labs(title = "Sensitivity: change in predicted arrival time, low -> high fuel quantile",
       x = NULL, y = "Change in arrival time (minutes)", fill = "Fuel type") +
  theme_minimal()

ggsave(file.path(out_dir, "fuel_factor_sensitivity_tornado.png"), width = 8, height = 5, dpi = 300)

6 Invaded vs. Non-Invaded Fuel Factors: Relative Impact

Are cogongrass (invaded) fuel conditions more important for fire arrival time than the non-invaded fuel’s conditions?

6.1 Overall impact: nested likelihood-ratio test

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")

build_formula <- function(terms) {
  as.formula(paste0(
    "arrival_avg ~ cover_c*distribution + cover_c*(", paste(terms, collapse = "+"), ") + ",
    "distribution*(", paste(terms, collapse = "+"), ")"
  ))
}

m_no_invaded    <- glm(build_formula(ni_terms),  data = scenario_dt, family = Gamma(link = "log"))
m_no_noninvaded <- glm(build_formula(inv_terms), data = scenario_dt, family = Gamma(link = "log"))

cat("Impact of INVADED fuel terms (full vs. without them):\n")
## Impact of INVADED fuel terms (full vs. without them):
print(anova(m_no_invaded, m_gamma, test = "LRT"))
## Analysis of Deviance Table
## 
## Model 1: arrival_avg ~ cover_c * distribution + cover_c * (ni_moist_live + 
##     ni_moist_dead + ni_fuel_live + ni_fuel_dead) + distribution * 
##     (ni_moist_live + ni_moist_dead + ni_fuel_live + ni_fuel_dead)
## Model 2: arrival_avg ~ cover_c * distribution + cover_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) + 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)
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1    190221    26736.4                          
## 2    190181     4412.6 40    22324 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nImpact of NON-INVADED fuel terms (full vs. without them):\n")
## 
## Impact of NON-INVADED fuel terms (full vs. without them):
print(anova(m_no_noninvaded, m_gamma, test = "LRT"))
## Analysis of Deviance Table
## 
## Model 1: arrival_avg ~ cover_c * distribution + cover_c * (inv_moist_live + 
##     inv_moist_dead + inv_fuel_live + inv_fuel_dead) + distribution * 
##     (inv_moist_live + inv_moist_dead + inv_fuel_live + inv_fuel_dead)
## Model 2: arrival_avg ~ cover_c * distribution + cover_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) + 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)
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1    190221    29132.9                          
## 2    190181     4412.6 40    24720 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.2 Does the impact shift with cover level?

factors <- c(inv_terms, ni_terms)
cover_levels <- seq(10, 90, 20)

get_pct_change <- function(fac, cov) {
  at_list <- setNames(as.list(rep("M", length(mf_cols))), mf_cols)
  at_list[[fac]] <- c("L", "H")
  at_list$cover_c <- cov - 50
  at_list$distribution <- "Random"
  emm <- emmeans(m_gamma, as.formula(paste("~", fac)), at = at_list, type = "response")
  d <- as.data.frame(emm)
  lv <- d$response[d[[fac]] == "L"]
  hv <- d$response[d[[fac]] == "H"]
  (hv / lv - 1) * 100
}

effect_grid <- expand.grid(factor = factors, cover = cover_levels, stringsAsFactors = FALSE)
effect_grid$pct_change <- mapply(get_pct_change, effect_grid$factor, effect_grid$cover)
effect_grid$type <- ifelse(grepl("^inv_", effect_grid$factor), "Invaded", "Non-invaded")
effect_grid$abs_pct <- abs(effect_grid$pct_change)

summary_df <- aggregate(abs_pct ~ type + cover, data = effect_grid, FUN = mean)
print(summary_df)
##           type cover    abs_pct
## 1      Invaded    10  8.5842178
## 2  Non-invaded    10 61.9422758
## 3      Invaded    30 21.9251482
## 4  Non-invaded    30 41.8335170
## 5      Invaded    50 36.2262192
## 6  Non-invaded    50 25.4854950
## 7      Invaded    70 51.7899497
## 8  Non-invaded    70 11.9566429
## 9      Invaded    90 68.9471743
## 10 Non-invaded    90  0.5158219
ggplot(summary_df, aes(x = cover, y = abs_pct, color = type)) +
  geom_line(linewidth = 1.2) + geom_point(size = 2) +
  labs(title = "Relative impact of invaded vs. non-invaded fuel conditions by cover level",
       x = "Cogongrass cover (%)",
       y = "Mean |% change in arrival time| (lowest vs. highest fuel quantile)",
       color = "Fuel type") +
  theme_minimal()

ggsave(file.path(out_dir, "invaded_vs_noninvaded_impact_by_cover.png"), width = 7, height = 5, dpi = 300)

7 Sensitivity Check

p_map <- c(Random = 0, ModClump = 0.3, Clump = 0.55)
sens_dt <- scenario_dt[distribution != "CompleteDominance"]
sens_dt[, clustering_p := p_map[as.character(distribution)]]

m_gamma_p <- glm(arrival_avg ~ cover_c * clustering_p +
                    inv_moist_live + inv_moist_dead + ni_moist_live + ni_moist_dead +
                    inv_fuel_live + inv_fuel_dead + ni_fuel_live + ni_fuel_dead,
                  data = sens_dt, family = Gamma(link = "log"))
summary(m_gamma_p)
## 
## Call:
## glm(formula = arrival_avg ~ cover_c * clustering_p + inv_moist_live + 
##     inv_moist_dead + ni_moist_live + ni_moist_dead + inv_fuel_live + 
##     inv_fuel_dead + ni_fuel_live + ni_fuel_dead, family = Gamma(link = "log"), 
##     data = sens_dt)
## 
## Coefficients:
##                        Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)           4.155e+00  1.156e-03 3593.796   <2e-16 ***
## cover_c              -7.627e-03  4.478e-05 -170.311   <2e-16 ***
## clustering_p          1.335e-01  3.197e-03   41.766   <2e-16 ***
## inv_moist_liveL      -1.665e-01  1.016e-03 -163.811   <2e-16 ***
## inv_moist_liveH       1.299e-01  1.016e-03  127.775   <2e-16 ***
## inv_moist_deadL      -1.833e-01  1.016e-03 -180.374   <2e-16 ***
## inv_moist_deadH       2.592e-01  1.016e-03  255.013   <2e-16 ***
## ni_moist_liveL       -4.284e-02  1.016e-03  -42.145   <2e-16 ***
## ni_moist_liveH        5.502e-02  1.016e-03   54.134   <2e-16 ***
## ni_moist_deadL       -2.469e-01  1.016e-03 -242.897   <2e-16 ***
## ni_moist_deadH        3.892e-01  1.016e-03  382.946   <2e-16 ***
## inv_fuel_liveL       -6.896e-02  1.016e-03  -67.840   <2e-16 ***
## inv_fuel_liveH        8.279e-02  1.016e-03   81.450   <2e-16 ***
## inv_fuel_deadL        2.208e-01  1.016e-03  217.243   <2e-16 ***
## inv_fuel_deadH       -1.907e-01  1.016e-03 -187.660   <2e-16 ***
## ni_fuel_liveL        -1.495e-02  1.016e-03  -14.704   <2e-16 ***
## ni_fuel_liveH         1.401e-02  1.016e-03   13.787   <2e-16 ***
## ni_fuel_deadL         1.449e-01  1.016e-03  142.542   <2e-16 ***
## ni_fuel_deadH        -1.317e-01  1.016e-03 -129.583   <2e-16 ***
## cover_c:clustering_p  2.461e-04  1.238e-04    1.987   0.0469 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.09151222)
## 
##     Null deviance: 53985  on 177146  degrees of freedom
## Residual deviance: 13333  on 177127  degrees of freedom
## AIC: 1519017
## 
## Number of Fisher Scoring iterations: 6

7.1 Marginal Effects

cmp <- avg_comparisons(
  m_gamma,
  variables = list(
    inv_moist_live = c("L","H"),
    ni_moist_live  = c("L","H"),
    inv_moist_dead = c("L","H"),
    ni_moist_dead  = c("L","H"),
    inv_fuel_live  = c("L","H"),
    ni_fuel_live   = c("L","H"),
    inv_fuel_dead  = c("L","H"),
    ni_fuel_dead   = c("L","H")
  ),
  type = "response"
)

cmp
## 
##            Term Estimate Std. Error      z Pr(>|z|)     S  2.5 % 97.5 %
##  inv_fuel_dead    -26.29     0.0704 -373.2   <0.001   Inf -26.43 -26.15
##  inv_fuel_live      9.45     0.0687  137.6   <0.001   Inf   9.32   9.59
##  inv_moist_dead    29.08     0.0719  404.5   <0.001   Inf  28.94  29.22
##  inv_moist_live    18.25     0.0680  268.2   <0.001   Inf  18.12  18.38
##  ni_fuel_dead     -22.82     0.0703 -324.6   <0.001   Inf -22.96 -22.69
##  ni_fuel_live       2.37     0.0682   34.7   <0.001 874.6   2.23   2.50
##  ni_moist_dead     55.21     0.0818  674.6   <0.001   Inf  55.05  55.37
##  ni_moist_live      8.09     0.0689  117.4   <0.001   Inf   7.96   8.23
## 
## Type: response
## Comparison: H - L
hypotheses(
  cmp,
  hypothesis =
"(b1 + b2 + b3 + b4)/4 = (b5 + b6 + b7 + b8)/4"
)
## 
##                       Hypothesis Estimate Std. Error     z Pr(>|z|)   S 2.5 %
##  (b1+b2+b3+b4)/4=(b5+b6+b7+b8)/4    -3.09     0.0487 -63.3   <0.001 Inf -3.18
##  97.5 %
##   -2.99

8 Appendix — Mixed-Effects Cross-Check on Ignition-Level Data

con <- dbConnect(duckdb())
ign_dt <- setDT(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, ignition, arrival_mean
  FROM read_csv_auto('%s')
", input_path)))
dbDisconnect(con, shutdown = TRUE)

ign_dt[, scenario_id := .GRP, by = .(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)]
ign_dt[, distribution := fifelse(is.na(distribution), "CompleteDominance", distribution)]
ign_dt[, distribution := factor(distribution, levels = c("CompleteDominance","Random","ModClump","Clump"))]
ign_dt[, (mf_cols) := lapply(.SD, factor, levels = c("L","M","H")), .SDcols = mf_cols]
ign_dt[, cover_c := cover_percent - 50]

m_glmm <- glmmTMB(arrival_mean ~ cover_c * 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 +
                     (1 | scenario_id),
                   data = ign_dt, family = Gamma(link = "log"))
summary(m_glmm)

9 Save Model Objects

saveRDS(m_gamma,          file.path(out_dir, "gamma_glm_primary_model.rds"))
saveRDS(m_gamma_p,        file.path(out_dir, "gamma_glm_continuous_p_sensitivity.rds"))
saveRDS(m_no_invaded,     file.path(out_dir, "gamma_glm_no_invaded_terms.rds"))
saveRDS(m_no_noninvaded,  file.path(out_dir, "gamma_glm_no_noninvaded_terms.rds"))
fwrite(scenario_dt, file.path(out_dir, "scenario_level_data.csv"))
fwrite(effect_grid, file.path(out_dir, "invaded_vs_noninvaded_effect_grid.csv"))