Refugia Edge Density & Patch Area Modeling
Alan Ivory
2026-07-30
Selected Variables
selected_models <- list(
Edge_Density = c("ag_1985", "wood_cov", "nearest_road_m",
"NDVI_pre_90d", "NDVI_pre_365d", "mean_hydro_dist_m", "mean_slope",
"mean_elevation", "soil_order_1", "vs_mea_3b_2a",
"vs_mea_30b_2a", "tmmn_max_3b_2a", "pr_max_3b_2a", "pr_max_10b_2a",
"pr_mea_60b_2a", "rmin_mea_60b_2a", "spi1y"),
Avg_Patch_Area = c("ag_1985", "wood_cov", "road_dens_km_km2", "nearest_road_m",
"NDVI_pre_90d", "NDVI_pre_365d",
"mean_hydro_dist_m", "mean_elevation", "soil_order_1", "vs_min_3b_2a",
"tmmn_mea_60b_2a", "tmmx_max_3b_2a", "pr_mea_10b_2a",
"pr_mea_3b_2a", "vpd_min_10b_2a", "spei1y", "spei14d")
)
selected_models## $Edge_Density
## [1] "ag_1985" "wood_cov" "nearest_road_m"
## [4] "NDVI_pre_90d" "NDVI_pre_365d" "mean_hydro_dist_m"
## [7] "mean_slope" "mean_elevation" "soil_order_1"
## [10] "vs_mea_3b_2a" "vs_mea_30b_2a" "tmmn_max_3b_2a"
## [13] "pr_max_3b_2a" "pr_max_10b_2a" "pr_mea_60b_2a"
## [16] "rmin_mea_60b_2a" "spi1y"
##
## $Avg_Patch_Area
## [1] "ag_1985" "wood_cov" "road_dens_km_km2"
## [4] "nearest_road_m" "NDVI_pre_90d" "NDVI_pre_365d"
## [7] "mean_hydro_dist_m" "mean_elevation" "soil_order_1"
## [10] "vs_min_3b_2a" "tmmn_mea_60b_2a" "tmmx_max_3b_2a"
## [13] "pr_mea_10b_2a" "pr_mea_3b_2a" "vpd_min_10b_2a"
## [16] "spei1y" "spei14d"Prediction Helpers
Uses the lognormal mean correction
exp(fit + 0.5 * se^2).
make_baseline <- function(data) {
as.data.frame(lapply(data, function(x) {
if (is.numeric(x)) return(mean(x, na.rm = TRUE))
if (is.factor(x)) return(levels(x)[1])
if (is.character(x)) return(na.omit(x)[1])
return(x[1])
}))
}
# Marginal effect
get_effect <- function(model, data, var, n = 100) {
x <- data[[var]]
x <- x[is.finite(x)]
if (length(x) < 2) return(NULL)
x_seq <- seq(min(x), max(x), length.out = n)
newdat <- make_baseline(data)
newdat <- newdat[rep(1, n), ]
newdat[[var]] <- x_seq
pred <- predict(model, newdata = newdat, se.fit = TRUE, type = "link")
data.frame(
x = x_seq,
fit = exp(pred$fit + 0.5 * pred$se.fit^2),
upper = exp(pred$fit + 2 * pred$se.fit),
lower = exp(pred$fit - 2 * pred$se.fit),
variable = var
)
}
# Interaction
get_interaction_lines <- function(model, data, var1, var2,
probs = c(0.1, 0.25, 0.5, 0.75, 0.9), n = 100) {
x_seq <- seq(min(data[[var1]], na.rm = TRUE),
max(data[[var1]], na.rm = TRUE), length.out = n)
y_vals <- quantile(data[[var2]], probs = probs, na.rm = TRUE)
df_list <- lapply(seq_along(y_vals), function(i) {
newdat <- make_baseline(data)
newdat <- newdat[rep(1, n), ]
newdat[[var1]] <- x_seq
newdat[[var2]] <- y_vals[i]
pred <- predict(model, newdata = newdat, se.fit = TRUE, type = "link")
data.frame(
x = x_seq,
fit = exp(pred$fit + 0.5 * pred$se.fit^2),
upper = exp(pred$fit + 2 * pred$se.fit),
lower = exp(pred$fit - 2 * pred$se.fit),
group = paste0(names(y_vals)[i], " (", round(y_vals[i], 2), ")")
)
})
do.call(rbind, df_list)
}
# Back-transformed parametric / smooth summary tables
backtransform_tables <- function(model) {
param_df <- as.data.frame(summary(model)$p.table)
param_df$exp_estimate <- exp(param_df$Estimate)
param_df$percent_change <- (param_df$exp_estimate - 1) * 100
param_df <- param_df %>%
dplyr::mutate(
exp_estimate = round(exp_estimate, 3),
percent_change = round(percent_change, 1)
)
smooth_df <- as.data.frame(summary(model)$s.table) %>%
dplyr::mutate(
edf = round(edf, 2),
F = round(F, 2),
p_value = round(`p-value`, 4)
)
list(param = param_df, smooth = smooth_df)
}Supplemental Figure Helpers
Build a single faceted figure showing the fitted relationship for every one-dimensional term in a model, with significance stars taken from the model summary. Two-dimensional tensor interactions are excluded here; they appear in their own figures.
# Pretty axis labels, shared by both models
var_labels <- c(
ag_1985 = "Agricultural Land Use (1985, %)",
wood_cov = "Wood Cover (%)",
road_dens_km_km2 = "Road Density (km/km²)",
nearest_road_m = "Distance to Nearest Road (m)",
NDVI_pre_90d = "NDVI (90-day)",
NDVI_pre_365d = "NDVI (365-day)",
mean_hydro_dist_m = "Distance to Hydrologic Feature (m)",
mean_slope = "Slope",
mean_elevation = "Elevation (m)",
soil_order_1 = "Dominant Soil Order (code)",
vs_mea_3b_2a = "Mean Wind Speed (3-day, m/s)",
vs_mea_30b_2a = "Mean Wind Speed (30-day, m/s)",
vs_min_3b_2a = "Minimum Wind Speed (3-day, m/s)",
tmmn_max_3b_2a = "Max of Minimum Temperature (3-day, K)",
tmmn_mea_60b_2a = "Mean Minimum Temperature (60-day, K)",
tmmx_max_3b_2a = "Maximum Temperature (3-day, K)",
pr_max_3b_2a = "Maximum Precipitation (3-day, mm)",
pr_max_10b_2a = "Maximum Precipitation (10-day, mm)",
pr_mea_10b_2a = "Mean Precipitation (10-day, mm)",
pr_mea_60b_2a = "Mean Precipitation (60-day, mm)",
rmin_mea_60b_2a = "Mean Minimum Relative Humidity (60-day, %)",
vpd_min_10b_2a = "Minimum VPD (10-day, kPa)",
spi1y = "SPI (1 year)",
spei1y = "SPEI (1 year)",
spei14d = "SPEI (14 day)"
)
# p-value for every 1-D term, named by the underlying variable.
# Parametric terms use pTerms.pv, which gives one joint p-value per term - so a
# factor like soil_order_1 gets a single test rather than one row per level.
term_pvalues <- function(model) {
s <- summary(model)
p_par <- s$pTerms.pv
if (is.null(p_par)) { # fallback if no parametric terms
p_par <- numeric(0)
}
p_par <- p_par[!grepl(":", names(p_par))] # parametric interactions excluded
p_sm <- s$s.table[, "p-value"]
nm <- gsub("^(s|ti|te|t2)\\(|\\)$", "", names(p_sm))
keep <- !grepl(",", nm) # 2-D tensor terms excluded
p_sm <- p_sm[keep]
names(p_sm) <- nm[keep]
c(p_par, p_sm)
}
sig_stars <- function(p) {
ifelse(is.na(p), "",
ifelse(p < 0.001, "***",
ifelse(p < 0.01, "**",
ifelse(p < 0.05, "*",
ifelse(p < 0.10, ".", "")))))
}
# Panel title with significance stars
panel_title <- function(v, p, labels = var_labels) {
lab <- if (v %in% names(labels)) labels[[v]] else v
trimws(paste(lab, sig_stars(p)))
}
# Marginal effects for every continuous 1-D term, with starred panel titles
supplemental_effects <- function(model, data, labels = var_labels) {
pvals <- term_pvalues(model)
vars <- names(pvals)
vars <- vars[vars %in% names(data)]
vars <- vars[vapply(data[vars], is.numeric, logical(1))]
eff <- do.call(rbind, lapply(vars, function(v) get_effect(model, data, v)))
panel <- setNames(vapply(vars, function(v) panel_title(v, pvals[[v]], labels),
character(1)), vars)
eff$panel <- factor(panel[eff$variable], levels = panel[vars])
eff
}
# Fitted value at each level of a categorical predictor, other terms at baseline
get_factor_effect <- function(model, data, var) {
lv <- levels(data[[var]])
newdat <- make_baseline(data)
newdat <- newdat[rep(1, length(lv)), ]
newdat[[var]] <- factor(lv, levels = lv)
pred <- predict(model, newdata = newdat, se.fit = TRUE, type = "link")
data.frame(
level = factor(lv, levels = lv),
fit = exp(pred$fit + 0.5 * pred$se.fit^2),
upper = exp(pred$fit + 2 * pred$se.fit),
lower = exp(pred$fit - 2 * pred$se.fit),
variable = var
)
}
# All categorical terms in a model
supplemental_factor_effects <- function(model, data, labels = var_labels) {
pvals <- term_pvalues(model)
vars <- names(pvals)
vars <- vars[vars %in% names(data)]
vars <- vars[vapply(data[vars], is.factor, logical(1))]
if (length(vars) == 0) return(NULL)
eff <- do.call(rbind, lapply(vars, function(v) get_factor_effect(model, data, v)))
panel <- setNames(vapply(vars, function(v) panel_title(v, pvals[[v]], labels),
character(1)), vars)
eff$panel <- factor(panel[eff$variable], levels = panel[vars])
eff
}
plot_factor_effects <- function(eff, y_lab, colour = "#0072B2") {
ggplot(eff, aes(x = level, y = fit)) +
geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.15, colour = colour) +
geom_point(size = 3, colour = colour) +
facet_wrap(~ panel, scales = "free") +
theme_classic(base_size = 12) +
theme(
strip.background = element_blank(),
strip.text = element_text(size = 11, face = "bold"),
axis.title = element_text(size = 13, face = "bold")
) +
labs(x = NULL, y = y_lab)
}
plot_supplemental <- function(eff, y_lab, ncol = 4, colour = "#0072B2") {
ggplot(eff, aes(x = x, y = fit)) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = colour) +
geom_line(linewidth = 1, colour = colour) +
facet_wrap(~ panel, scales = "free", ncol = ncol) +
theme_classic(base_size = 12) +
theme(
strip.background = element_blank(),
strip.text = element_text(size = 11, face = "bold"),
axis.title = element_text(size = 13, face = "bold")
) +
labs(x = NULL, y = y_lab)
}Cross-Validation Helpers
cv_gam <- function(formula, data, K = 10, reps = 5, seed = 1) {
vars <- all.vars(formula)
d <- data[complete.cases(data[, vars, drop = FALSE]), , drop = FALSE]
y <- eval(formula[[2]], envir = d) # observed log response
set.seed(seed)
rep_stats <- lapply(seq_len(reps), function(r) {
folds <- sample(rep(seq_len(K), length.out = nrow(d)))
pred <- rep(NA_real_, nrow(d))
for (k in seq_len(K)) {
train <- d[folds != k, ]
test <- d[folds == k, ]
# A held-out row whose factor level is absent from the training fold cannot
# be predicted; drop those rows rather than losing the whole fold.
fac_vars <- vars[vapply(d[vars], is.factor, logical(1))]
ok_rows <- rep(TRUE, nrow(test))
for (fv in fac_vars) {
ok_rows <- ok_rows & test[[fv]] %in% unique(train[[fv]])
}
fit <- tryCatch(gam(formula, data = train, method = "REML"),
error = function(e) NULL)
if (is.null(fit) || !any(ok_rows)) next
idx <- which(folds == k)[ok_rows]
pred[idx] <- tryCatch(
as.numeric(predict(fit, newdata = test[ok_rows, , drop = FALSE], type = "link")),
error = function(e) NA_real_
)
}
ok <- is.finite(pred) & is.finite(y)
res <- y[ok] - pred[ok]
data.frame(
rep = r,
RMSE = sqrt(mean(res^2)),
MAE = mean(abs(res)),
R2 = 1 - sum(res^2) / sum((y[ok] - mean(y[ok]))^2),
n_pred = sum(ok)
)
})
do.call(rbind, rep_stats)
}
cv_compare <- function(model_list, data, K = 10, reps = 5, seed = 1) {
res <- lapply(names(model_list), function(nm) {
s <- cv_gam(formula(model_list[[nm]]), data, K = K, reps = reps, seed = seed)
data.frame(
model = nm,
CV_RMSE = mean(s$RMSE),
RMSE_SD = sd(s$RMSE),
CV_MAE = mean(s$MAE),
CV_R2 = mean(s$R2),
R2_SD = sd(s$R2)
)
})
res <- do.call(rbind, res)
res <- res[order(res$CV_RMSE), ]
res$delta_RMSE <- res$CV_RMSE - min(res$CV_RMSE)
rownames(res) <- NULL
res
}EDGE DENSITY
Log-normal edge density
mod_edge_ln <- gam(
log(edge_dens) ~
s(ag_1985) + wood_cov + nearest_road_m +
NDVI_pre_90d + s(NDVI_pre_365d) + s(mean_hydro_dist_m) +
s(mean_slope) + s(mean_elevation) + soil_order_1 +
vs_mea_3b_2a + vs_mea_30b_2a + tmmn_max_3b_2a +
pr_max_3b_2a + s(pr_max_10b_2a) + s(pr_mea_60b_2a) +
rmin_mea_60b_2a + s(spi1y),
data = landscape_data,
family = gaussian(),
method = "REML"
)
summary(mod_edge_ln)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(edge_dens) ~ s(ag_1985) + wood_cov + nearest_road_m + NDVI_pre_90d +
## s(NDVI_pre_365d) + s(mean_hydro_dist_m) + s(mean_slope) +
## s(mean_elevation) + soil_order_1 + vs_mea_3b_2a + vs_mea_30b_2a +
## tmmn_max_3b_2a + pr_max_3b_2a + s(pr_max_10b_2a) + s(pr_mea_60b_2a) +
## rmin_mea_60b_2a + s(spi1y)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.058e+00 1.678e+00 2.418 0.015942 *
## wood_cov 9.226e-03 2.495e-03 3.698 0.000241 ***
## nearest_road_m 2.007e-04 8.627e-05 2.327 0.020372 *
## NDVI_pre_90d -6.439e-02 2.071e-01 -0.311 0.756004
## soil_order_11 3.891e-02 8.344e-02 0.466 0.641189
## soil_order_14 -2.332e-01 1.155e-01 -2.019 0.043994 *
## soil_order_19 1.694e-01 1.071e-01 1.582 0.114237
## soil_order_1Other 6.434e-03 1.174e-01 0.055 0.956325
## vs_mea_3b_2a -7.905e-03 4.073e-02 -0.194 0.846178
## vs_mea_30b_2a -1.251e-02 6.569e-02 -0.190 0.849015
## tmmn_max_3b_2a 7.596e-03 5.648e-03 1.345 0.179230
## pr_max_3b_2a 1.955e-03 1.459e-03 1.340 0.180980
## rmin_mea_60b_2a -4.650e-03 7.255e-03 -0.641 0.521797
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(ag_1985) 2.056 2.549 5.465 0.00233 **
## s(NDVI_pre_365d) 2.232 2.839 1.162 0.27034
## s(mean_hydro_dist_m) 2.268 2.742 2.697 0.06220 .
## s(mean_slope) 1.800 2.295 0.394 0.64004
## s(mean_elevation) 6.371 7.481 8.277 < 2e-16 ***
## s(pr_max_10b_2a) 2.859 3.600 2.987 0.02205 *
## s(pr_mea_60b_2a) 1.170 1.318 0.134 0.73958
## s(spi1y) 2.756 3.451 2.055 0.09755 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.294 Deviance explained = 33.8%
## -REML = 521.33 Scale est. = 0.31647 n = 538performance::check_outliers(mod_edge_ln)## OK: No outliers detected.
## - Based on the following method and threshold: cook (1).
## - For variable: (Whole model)Cleaned up log-normal edge density
sim_res <- simulateResiduals(mod_edge_ln)
outlier_idx <- outliers(sim_res)
# Drops NA rows
rows_used_edge <- as.numeric(rownames(mod_edge_ln$model))
if (length(outlier_idx) > 0) {
landscape_data_clean <- landscape_data[-rows_used_edge[outlier_idx], ]
} else {
landscape_data_clean <- landscape_data
}
nrow(landscape_data) - nrow(landscape_data_clean) # n outliers dropped## [1] 8## Base model, outliers removed
mod_edge_clean <- gam(log(edge_dens) ~
s(ag_1985) + wood_cov + nearest_road_m +
NDVI_pre_90d + s(NDVI_pre_365d) + s(mean_hydro_dist_m) +
s(mean_slope) + s(mean_elevation) + soil_order_1 +
vs_mea_3b_2a + vs_mea_30b_2a + tmmn_max_3b_2a +
pr_max_3b_2a + s(pr_max_10b_2a) + pr_mea_60b_2a +
rmin_mea_60b_2a + s(spi1y),
data = landscape_data_clean,
method = "REML")
summary(mod_edge_clean)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(edge_dens) ~ s(ag_1985) + wood_cov + nearest_road_m + NDVI_pre_90d +
## s(NDVI_pre_365d) + s(mean_hydro_dist_m) + s(mean_slope) +
## s(mean_elevation) + soil_order_1 + vs_mea_3b_2a + vs_mea_30b_2a +
## tmmn_max_3b_2a + pr_max_3b_2a + s(pr_max_10b_2a) + pr_mea_60b_2a +
## rmin_mea_60b_2a + s(spi1y)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.990e+00 1.536e+00 2.598 0.00967 **
## wood_cov 7.439e-03 2.325e-03 3.199 0.00147 **
## nearest_road_m 2.025e-04 7.981e-05 2.537 0.01147 *
## NDVI_pre_90d -8.212e-02 1.918e-01 -0.428 0.66872
## soil_order_11 6.812e-02 7.711e-02 0.883 0.37744
## soil_order_14 -2.872e-01 1.082e-01 -2.655 0.00818 **
## soil_order_19 1.785e-01 1.003e-01 1.780 0.07570 .
## soil_order_1Other 5.086e-02 1.087e-01 0.468 0.64015
## vs_mea_3b_2a -2.641e-03 3.777e-02 -0.070 0.94429
## vs_mea_30b_2a -7.508e-03 6.137e-02 -0.122 0.90268
## tmmn_max_3b_2a 7.114e-03 5.233e-03 1.360 0.17459
## pr_max_3b_2a 2.504e-03 1.348e-03 1.858 0.06372 .
## pr_mea_60b_2a -1.068e-02 1.869e-02 -0.572 0.56792
## rmin_mea_60b_2a 2.853e-03 6.793e-03 0.420 0.67468
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(ag_1985) 2.110 2.612 5.803 0.00137 **
## s(NDVI_pre_365d) 2.590 3.284 2.225 0.07666 .
## s(mean_hydro_dist_m) 2.388 2.888 2.907 0.04126 *
## s(mean_slope) 2.085 2.672 0.877 0.41404
## s(mean_elevation) 6.244 7.360 9.853 < 2e-16 ***
## s(pr_max_10b_2a) 2.643 3.341 3.679 0.00846 **
## s(spi1y) 3.571 4.435 2.825 0.03417 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.335 Deviance explained = 37.8%
## -REML = 474.65 Scale est. = 0.26859 n = 530## NDVI x drought interaction
mod_edge_NDVI_Climate <- gam(log(edge_dens) ~
s(ag_1985) + wood_cov + nearest_road_m +
NDVI_pre_90d + s(NDVI_pre_365d) + s(mean_hydro_dist_m) +
s(mean_slope) + s(mean_elevation) + soil_order_1 +
vs_mea_3b_2a + vs_mea_30b_2a + tmmn_max_3b_2a +
pr_max_3b_2a + s(pr_max_10b_2a) + pr_mea_60b_2a +
rmin_mea_60b_2a + s(spi1y) +
ti(NDVI_pre_365d, spi1y),
data = landscape_data_clean,
method = "REML")
summary(mod_edge_NDVI_Climate)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(edge_dens) ~ s(ag_1985) + wood_cov + nearest_road_m + NDVI_pre_90d +
## s(NDVI_pre_365d) + s(mean_hydro_dist_m) + s(mean_slope) +
## s(mean_elevation) + soil_order_1 + vs_mea_3b_2a + vs_mea_30b_2a +
## tmmn_max_3b_2a + pr_max_3b_2a + s(pr_max_10b_2a) + pr_mea_60b_2a +
## rmin_mea_60b_2a + s(spi1y) + ti(NDVI_pre_365d, spi1y)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.296e+00 1.528e+00 2.812 0.00513 **
## wood_cov 8.072e-03 2.327e-03 3.468 0.00057 ***
## nearest_road_m 2.045e-04 7.942e-05 2.575 0.01032 *
## NDVI_pre_90d -1.280e-01 1.907e-01 -0.671 0.50242
## soil_order_11 7.953e-02 7.673e-02 1.037 0.30044
## soil_order_14 -2.515e-01 1.072e-01 -2.347 0.01932 *
## soil_order_19 2.123e-01 9.967e-02 2.130 0.03365 *
## soil_order_1Other 6.301e-02 1.077e-01 0.585 0.55866
## vs_mea_3b_2a -1.571e-02 3.745e-02 -0.419 0.67507
## vs_mea_30b_2a 6.891e-03 6.178e-02 0.112 0.91124
## tmmn_max_3b_2a 6.056e-03 5.226e-03 1.159 0.24702
## pr_max_3b_2a 2.538e-03 1.331e-03 1.907 0.05708 .
## pr_mea_60b_2a -1.809e-02 1.947e-02 -0.929 0.35324
## rmin_mea_60b_2a 2.353e-03 6.772e-03 0.347 0.72839
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(ag_1985) 2.002 2.482 6.355 0.000937 ***
## s(NDVI_pre_365d) 2.364 2.994 1.817 0.144227
## s(mean_hydro_dist_m) 2.167 2.617 2.004 0.139193
## s(mean_slope) 1.922 2.454 0.710 0.500594
## s(mean_elevation) 6.173 7.290 11.083 < 2e-16 ***
## s(pr_max_10b_2a) 2.546 3.227 3.105 0.021598 *
## s(spi1y) 2.384 2.984 1.047 0.353945
## ti(NDVI_pre_365d,spi1y) 7.621 9.640 2.325 0.011105 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.356 Deviance explained = 40.5%
## -REML = 471.01 Scale est. = 0.25981 n = 530## Land cover x topography interactions
mod_edge_LandCov_Slope <- gam(log(edge_dens) ~
s(ag_1985) + wood_cov + nearest_road_m +
NDVI_pre_90d + s(NDVI_pre_365d) + s(mean_hydro_dist_m) +
s(mean_slope) + s(mean_elevation) + soil_order_1 +
vs_mea_3b_2a + vs_mea_30b_2a + tmmn_max_3b_2a +
pr_max_3b_2a + s(pr_max_10b_2a) + pr_mea_60b_2a +
rmin_mea_60b_2a + s(spi1y) +
ti(wood_cov, mean_slope) +
ti(ag_1985, mean_slope),
data = landscape_data_clean,
method = "REML")
summary(mod_edge_LandCov_Slope)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(edge_dens) ~ s(ag_1985) + wood_cov + nearest_road_m + NDVI_pre_90d +
## s(NDVI_pre_365d) + s(mean_hydro_dist_m) + s(mean_slope) +
## s(mean_elevation) + soil_order_1 + vs_mea_3b_2a + vs_mea_30b_2a +
## tmmn_max_3b_2a + pr_max_3b_2a + s(pr_max_10b_2a) + pr_mea_60b_2a +
## rmin_mea_60b_2a + s(spi1y) + ti(wood_cov, mean_slope) + ti(ag_1985,
## mean_slope)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.182e+00 1.528e+00 2.736 0.00644 **
## wood_cov 7.068e-03 2.341e-03 3.019 0.00266 **
## nearest_road_m 2.121e-04 7.916e-05 2.680 0.00761 **
## NDVI_pre_90d -2.795e-02 1.916e-01 -0.146 0.88406
## soil_order_11 5.130e-02 7.673e-02 0.669 0.50407
## soil_order_14 -2.596e-01 1.088e-01 -2.386 0.01739 *
## soil_order_19 1.559e-01 9.995e-02 1.560 0.11936
## soil_order_1Other 2.005e-02 1.085e-01 0.185 0.85351
## vs_mea_3b_2a -4.959e-03 3.745e-02 -0.132 0.89471
## vs_mea_30b_2a 3.100e-03 6.115e-02 0.051 0.95958
## tmmn_max_3b_2a 6.404e-03 5.200e-03 1.232 0.21871
## pr_max_3b_2a 2.697e-03 1.340e-03 2.013 0.04465 *
## pr_mea_60b_2a -8.444e-03 1.852e-02 -0.456 0.64865
## rmin_mea_60b_2a 2.789e-03 6.725e-03 0.415 0.67852
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(ag_1985) 1.868 2.314 6.569 0.000949 ***
## s(NDVI_pre_365d) 2.471 3.133 1.982 0.112081
## s(mean_hydro_dist_m) 2.356 2.849 2.669 0.059373 .
## s(mean_slope) 1.152 1.281 0.103 0.766239
## s(mean_elevation) 6.215 7.329 11.243 < 2e-16 ***
## s(pr_max_10b_2a) 2.602 3.293 4.344 0.003523 **
## s(spi1y) 3.586 4.452 3.246 0.015804 *
## ti(wood_cov,mean_slope) 3.170 4.306 2.482 0.038955 *
## ti(ag_1985,mean_slope) 1.001 1.002 5.677 0.017537 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.348 Deviance explained = 39.4%
## -REML = 467.08 Scale est. = 0.26317 n = 530## Combined interaction model
mod_edge_Int <- gam(log(edge_dens) ~
s(ag_1985) + wood_cov + nearest_road_m +
NDVI_pre_90d + s(NDVI_pre_365d) + s(mean_hydro_dist_m) +
s(mean_slope) + s(mean_elevation) + soil_order_1 +
vs_mea_3b_2a + vs_mea_30b_2a + tmmn_max_3b_2a +
pr_max_3b_2a + s(pr_max_10b_2a) + pr_mea_60b_2a +
rmin_mea_60b_2a + s(spi1y) +
ti(NDVI_pre_365d, spi1y) +
ti(wood_cov, mean_slope) +
ag_1985:mean_slope,
data = landscape_data_clean,
method = "REML")
summary(mod_edge_Int)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(edge_dens) ~ s(ag_1985) + wood_cov + nearest_road_m + NDVI_pre_90d +
## s(NDVI_pre_365d) + s(mean_hydro_dist_m) + s(mean_slope) +
## s(mean_elevation) + soil_order_1 + vs_mea_3b_2a + vs_mea_30b_2a +
## tmmn_max_3b_2a + pr_max_3b_2a + s(pr_max_10b_2a) + pr_mea_60b_2a +
## rmin_mea_60b_2a + s(spi1y) + ti(NDVI_pre_365d, spi1y) + ti(wood_cov,
## mean_slope) + ag_1985:mean_slope
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.402e+00 1.523e+00 2.891 0.00402 **
## wood_cov 7.629e-03 2.347e-03 3.251 0.00123 **
## nearest_road_m 2.150e-04 7.889e-05 2.726 0.00664 **
## NDVI_pre_90d -6.404e-02 1.910e-01 -0.335 0.73761
## soil_order_11 6.686e-02 7.654e-02 0.873 0.38283
## soil_order_14 -2.426e-01 1.076e-01 -2.256 0.02454 *
## soil_order_19 1.813e-01 9.967e-02 1.819 0.06948 .
## soil_order_1Other 2.705e-02 1.077e-01 0.251 0.80188
## vs_mea_3b_2a -1.507e-02 3.725e-02 -0.405 0.68600
## vs_mea_30b_2a 1.344e-02 6.160e-02 0.218 0.82738
## tmmn_max_3b_2a 5.911e-03 5.198e-03 1.137 0.25608
## pr_max_3b_2a 2.743e-03 1.327e-03 2.067 0.03928 *
## pr_mea_60b_2a -1.958e-02 1.928e-02 -1.016 0.31023
## rmin_mea_60b_2a 2.223e-03 6.723e-03 0.331 0.74098
## ag_1985:mean_slope -1.338e+00 6.093e-01 -2.196 0.02855 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(ag_1985) 1.765 2.183 5.963 0.00211 **
## s(NDVI_pre_365d) 2.289 2.892 1.485 0.20687
## s(mean_hydro_dist_m) 2.131 2.574 1.910 0.15955
## s(mean_slope) 1.005 1.010 0.904 0.33943
## s(mean_elevation) 6.097 7.212 12.595 < 2e-16 ***
## s(pr_max_10b_2a) 2.647 3.343 3.539 0.01112 *
## s(spi1y) 2.416 3.024 1.188 0.31117
## ti(NDVI_pre_365d,spi1y) 7.547 9.743 1.998 0.03592 *
## ti(wood_cov,mean_slope) 3.500 4.739 2.535 0.03026 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.367 Deviance explained = 41.9%
## -REML = 464.79 Scale est. = 0.25545 n = 530## Model Selection
mod_selection_edge <- model.sel(mod_edge_clean, mod_edge_NDVI_Climate,
mod_edge_LandCov_Slope, mod_edge_Int)
mod_selection_edge## Model selection table
## (Int) NDV_pre_90d nrs_rod_m pr_max_3b_2a pr_mea_60b_2a
## mod_edge_Int 4.402 -0.06404 0.0002150 0.002743 -0.019580
## mod_edge_NDVI_Climate 4.296 -0.12800 0.0002045 0.002538 -0.018090
## mod_edge_LandCov_Slope 4.182 -0.02795 0.0002121 0.002697 -0.008444
## mod_edge_clean 3.990 -0.08212 0.0002025 0.002504 -0.010680
## rmn_mea_60b_2a s(ag_1985) s(men_elv) s(men_hyd_dst_m)
## mod_edge_Int 0.002223 + + +
## mod_edge_NDVI_Climate 0.002353 + + +
## mod_edge_LandCov_Slope 0.002789 + + +
## mod_edge_clean 0.002853 + + +
## s(men_slp) s(NDV_pre_365) s(pr_max_10b_2a) s(sp1)
## mod_edge_Int + + + +
## mod_edge_NDVI_Climate + + + +
## mod_edge_LandCov_Slope + + + +
## mod_edge_clean + + + +
## sol_ord_1 tmm_max_3b_2a vs_mea_30b_2a vs_mea_3b_2a
## mod_edge_Int + 0.005911 0.013440 -0.015070
## mod_edge_NDVI_Climate + 0.006056 0.006891 -0.015710
## mod_edge_LandCov_Slope + 0.006404 0.003100 -0.004959
## mod_edge_clean + 0.007114 -0.007508 -0.002641
## wod_cov ti(NDV_pre_365,sp1) ti(ag_1985,men_slp)
## mod_edge_Int 0.007629 +
## mod_edge_NDVI_Climate 0.008072 +
## mod_edge_LandCov_Slope 0.007068 +
## mod_edge_clean 0.007439
## ti(wod_cov,men_slp) ag_1985:men_slp df logLik AICc
## mod_edge_Int + -1.338 52 -367.205 851.7
## mod_edge_NDVI_Climate 48 -373.431 854.3
## mod_edge_LandCov_Slope + 44 -378.331 855.1
## mod_edge_clean 41 -385.239 860.9
## delta weight
## mod_edge_Int 0.00 0.679
## mod_edge_NDVI_Climate 2.57 0.188
## mod_edge_LandCov_Slope 3.37 0.126
## mod_edge_clean 9.19 0.007
## Models ranked by AICc(x)## Selected model
mod_edge_best <- mod_edge_Int
summary(mod_edge_best)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(edge_dens) ~ s(ag_1985) + wood_cov + nearest_road_m + NDVI_pre_90d +
## s(NDVI_pre_365d) + s(mean_hydro_dist_m) + s(mean_slope) +
## s(mean_elevation) + soil_order_1 + vs_mea_3b_2a + vs_mea_30b_2a +
## tmmn_max_3b_2a + pr_max_3b_2a + s(pr_max_10b_2a) + pr_mea_60b_2a +
## rmin_mea_60b_2a + s(spi1y) + ti(NDVI_pre_365d, spi1y) + ti(wood_cov,
## mean_slope) + ag_1985:mean_slope
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.402e+00 1.523e+00 2.891 0.00402 **
## wood_cov 7.629e-03 2.347e-03 3.251 0.00123 **
## nearest_road_m 2.150e-04 7.889e-05 2.726 0.00664 **
## NDVI_pre_90d -6.404e-02 1.910e-01 -0.335 0.73761
## soil_order_11 6.686e-02 7.654e-02 0.873 0.38283
## soil_order_14 -2.426e-01 1.076e-01 -2.256 0.02454 *
## soil_order_19 1.813e-01 9.967e-02 1.819 0.06948 .
## soil_order_1Other 2.705e-02 1.077e-01 0.251 0.80188
## vs_mea_3b_2a -1.507e-02 3.725e-02 -0.405 0.68600
## vs_mea_30b_2a 1.344e-02 6.160e-02 0.218 0.82738
## tmmn_max_3b_2a 5.911e-03 5.198e-03 1.137 0.25608
## pr_max_3b_2a 2.743e-03 1.327e-03 2.067 0.03928 *
## pr_mea_60b_2a -1.958e-02 1.928e-02 -1.016 0.31023
## rmin_mea_60b_2a 2.223e-03 6.723e-03 0.331 0.74098
## ag_1985:mean_slope -1.338e+00 6.093e-01 -2.196 0.02855 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(ag_1985) 1.765 2.183 5.963 0.00211 **
## s(NDVI_pre_365d) 2.289 2.892 1.485 0.20687
## s(mean_hydro_dist_m) 2.131 2.574 1.910 0.15955
## s(mean_slope) 1.005 1.010 0.904 0.33943
## s(mean_elevation) 6.097 7.212 12.595 < 2e-16 ***
## s(pr_max_10b_2a) 2.647 3.343 3.539 0.01112 *
## s(spi1y) 2.416 3.024 1.188 0.31117
## ti(NDVI_pre_365d,spi1y) 7.547 9.743 1.998 0.03592 *
## ti(wood_cov,mean_slope) 3.500 4.739 2.535 0.03026 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.367 Deviance explained = 41.9%
## -REML = 464.79 Scale est. = 0.25545 n = 530# Model Fit
performance::check_collinearity(mod_edge_best)## Model has interaction terms. VIFs might be inflated.
## Try to center the variables used for the interaction, or check
## multicollinearity among predictors of a model without interaction terms.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
## nearest_road_m 1.17 [ 1.09, 1.32] 1.08 0.86 [0.76, 0.92]
## NDVI_pre_90d 1.63 [ 1.47, 1.84] 1.28 0.61 [0.54, 0.68]
## vs_mea_3b_2a 2.30 [ 2.04, 2.62] 1.52 0.44 [0.38, 0.49]
## vs_mea_30b_2a 2.57 [ 2.27, 2.94] 1.60 0.39 [0.34, 0.44]
## tmmn_max_3b_2a 1.56 [ 1.41, 1.76] 1.25 0.64 [0.57, 0.71]
## pr_max_3b_2a 1.48 [ 1.35, 1.67] 1.22 0.67 [0.60, 0.74]
## pr_mea_60b_2a 2.64 [ 2.33, 3.02] 1.62 0.38 [0.33, 0.43]
## rmin_mea_60b_2a 1.69 [ 1.52, 1.91] 1.30 0.59 [0.52, 0.66]
##
## Moderate Correlation
##
## Term VIF VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
## s(ag_1985) 8.34 [ 7.18, 9.73] 2.89 0.12 [0.10, 0.14]
## s(mean_hydro_dist_m) 5.41 [ 4.69, 6.29] 2.33 0.18 [0.16, 0.21]
## soil_order_1 6.01 [ 5.19, 6.99] 1.25 0.17 [0.14, 0.19]
##
## High Correlation
##
## Term VIF VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
## wood_cov 57.47 [ 48.97, 67.48] 7.58 0.02 [0.01, 0.02]
## s(NDVI_pre_365d) 190.45 [162.08, 223.81] 13.80 5.25e-03 [0.00, 0.01]
## s(mean_slope) 14.06 [ 12.04, 16.45] 3.75 0.07 [0.06, 0.08]
## s(mean_elevation) 34.63 [ 29.54, 40.63] 5.88 0.03 [0.02, 0.03]
## s(pr_max_10b_2a) 30.39 [ 25.93, 35.65] 5.51 0.03 [0.03, 0.04]
## s(spi1y) 149.23 [127.02, 175.36] 12.22 6.70e-03 [0.01, 0.01]performance::check_outliers(mod_edge_best)## OK: No outliers detected.
## - Based on the following method and threshold: cook (1).
## - For variable: (Whole model)## Plot
plot(mod_edge_best, pages = 4, shade = TRUE, shade.col = "lightblue", all.terms = TRUE)## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter
## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in axis(1, at = seq_along(ll), labels = ll, ...): "shade" is not a
## graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter
## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter
Edge Density Cross-Validation
edge_candidates <- list(
mod_edge_clean = mod_edge_clean,
mod_edge_NDVI_Climate = mod_edge_NDVI_Climate,
mod_edge_LandCov_Slope = mod_edge_LandCov_Slope,
mod_edge_Int = mod_edge_Int
)
cv_edge <- cv_compare(edge_candidates, landscape_data_clean,
K = 10, reps = 5, seed = 1)
cv_edge## model CV_RMSE RMSE_SD CV_MAE CV_R2 R2_SD
## 1 mod_edge_NDVI_Climate 0.5428980 0.001507820 0.4323507 0.2684994 0.004064240
## 2 mod_edge_clean 0.5447582 0.002380710 0.4333763 0.2634712 0.006436699
## 3 mod_edge_Int 0.5461756 0.005574657 0.4358400 0.2595832 0.015179926
## 4 mod_edge_LandCov_Slope 0.5470742 0.003787106 0.4352183 0.2571783 0.010314839
## delta_RMSE
## 1 0.000000000
## 2 0.001860201
## 3 0.003277582
## 4 0.004176179write.csv(cv_edge, "edge_density_cv_results.csv", row.names = FALSE)Edge Density back-transformed table
edge_tables <- backtransform_tables(mod_edge_best)
edge_tables$param## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.4015533243 1.522628e+00 2.8907614 0.004015262
## wood_cov 0.0076291534 2.347053e-03 3.2505242 0.001231910
## nearest_road_m 0.0002150471 7.888781e-05 2.7259864 0.006642655
## NDVI_pre_90d -0.0640384530 1.910360e-01 -0.3352167 0.737606360
## soil_order_11 0.0668603091 7.654498e-02 0.8734774 0.382834645
## soil_order_14 -0.2426018234 1.075574e-01 -2.2555571 0.024542683
## soil_order_19 0.1813229694 9.966722e-02 1.8192839 0.069483796
## soil_order_1Other 0.0270504639 1.077487e-01 0.2510515 0.801880470
## vs_mea_3b_2a -0.0150668574 3.724572e-02 -0.4045258 0.686004189
## vs_mea_30b_2a 0.0134406572 6.160483e-02 0.2181754 0.827384120
## tmmn_max_3b_2a 0.0059107148 5.198286e-03 1.1370507 0.256077858
## pr_max_3b_2a 0.0027433480 1.327347e-03 2.0667903 0.039283392
## pr_mea_60b_2a -0.0195803530 1.927562e-02 -1.0158094 0.310225993
## rmin_mea_60b_2a 0.0022234833 6.722750e-03 0.3307401 0.740983436
## ag_1985:mean_slope -1.3380551876 6.092582e-01 -2.1962040 0.028549415
## exp_estimate percent_change
## (Intercept) 81.577 8057.7
## wood_cov 1.008 0.8
## nearest_road_m 1.000 0.0
## NDVI_pre_90d 0.938 -6.2
## soil_order_11 1.069 6.9
## soil_order_14 0.785 -21.5
## soil_order_19 1.199 19.9
## soil_order_1Other 1.027 2.7
## vs_mea_3b_2a 0.985 -1.5
## vs_mea_30b_2a 1.014 1.4
## tmmn_max_3b_2a 1.006 0.6
## pr_max_3b_2a 1.003 0.3
## pr_mea_60b_2a 0.981 -1.9
## rmin_mea_60b_2a 1.002 0.2
## ag_1985:mean_slope 0.262 -73.8edge_tables$smooth## edf Ref.df F p-value p_value
## s(ag_1985) 1.76 2.183396 5.96 0.002106984 0.0021
## s(NDVI_pre_365d) 2.29 2.892089 1.48 0.206872322 0.2069
## s(mean_hydro_dist_m) 2.13 2.573509 1.91 0.159546537 0.1595
## s(mean_slope) 1.01 1.010253 0.90 0.339428772 0.3394
## s(mean_elevation) 6.10 7.212241 12.60 0.000000000 0.0000
## s(pr_max_10b_2a) 2.65 3.343162 3.54 0.011118777 0.0111
## s(spi1y) 2.42 3.024152 1.19 0.311166937 0.3112
## ti(NDVI_pre_365d,spi1y) 7.55 9.743193 2.00 0.035921540 0.0359
## ti(wood_cov,mean_slope) 3.50 4.739282 2.54 0.030263742 0.0303write.csv(edge_tables$param, "backtransformed_param_summary.csv", row.names = TRUE)
write.csv(edge_tables$smooth, "smooth_summary.csv", row.names = TRUE)Edge Density Backtransformed - linear relationships
focus_vars <- c("wood_cov", "nearest_road_m")
effects_df <- do.call(rbind, lapply(focus_vars, function(v) {
get_effect(mod_edge_best, landscape_data_clean, v)
}))
x_labels <- c("wood_cov" = "Wood Cover (%)",
"nearest_road_m" = "Distance to Nearest Road (m)")
# Shared y-axis across panels
y_range_linear <- range(c(effects_df$lower, effects_df$upper), na.rm = TRUE)
plots <- lapply(focus_vars, function(var) {
df <- effects_df %>% filter(variable == var)
ggplot(df, aes(x = x, y = fit)) +
geom_line(linewidth = 1.2, color = "#D55E00") +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "#D55E00") +
coord_cartesian(ylim = y_range_linear) +
theme_classic(base_size = 14) +
theme(
axis.title = element_text(size = 16, face = "bold"),
axis.text = element_text(size = 14),
plot.tag = element_text(size = 18, face = "bold")
) +
labs(x = x_labels[var], y = "Predicted Edge Density")
})
combined_plot <- plots[[1]] + plots[[2]] +
plot_layout(ncol = 2) +
plot_annotation(tag_levels = "A")
combined_plot
ggsave(file.path(fig_dir, "Edge_Density_Effects.png"),
combined_plot, width = 12, height = 6)Edge Density Backtransformed - non-linear relationships
smooth_vars <- c("ag_1985", "mean_hydro_dist_m", "mean_elevation")
smooth_effects <- do.call(rbind, lapply(smooth_vars, function(v) {
get_effect(mod_edge_best, landscape_data_clean, v)
}))
x_labels_smooth <- c(
"ag_1985" = "Agricultural Land Use (1985, %)",
"mean_hydro_dist_m" = "Distance to Hydrologic Feature (m)",
"mean_elevation" = "Elevation (m)"
)
# Shared y-axis across panels
y_range_smooth_edge <- range(c(smooth_effects$lower, smooth_effects$upper), na.rm = TRUE)
smooth_plots <- lapply(smooth_vars, function(var) {
df <- smooth_effects %>% filter(variable == var)
ggplot(df, aes(x = x, y = fit)) +
geom_line(linewidth = 1.2, color = "forestgreen") +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "lightgreen") +
coord_cartesian(ylim = y_range_smooth_edge) +
theme_classic(base_size = 14) +
theme(
axis.title = element_text(size = 16, face = "bold"),
axis.text = element_text(size = 14),
plot.tag = element_text(size = 18, face = "bold")
) +
labs(x = x_labels_smooth[var], y = "Predicted Edge Density")
})
combined_smooth_plot <- smooth_plots[[1]] + smooth_plots[[2]] + smooth_plots[[3]] +
plot_layout(ncol = 3) +
plot_annotation(tag_levels = "A")
combined_smooth_plot
ggsave(file.path(fig_dir, "Edge_Density_Smooth_Effects.png"),
combined_smooth_plot, width = 18, height = 6)Edge Density Significant Main Effects
Significant terms that do NOT appear in the interaction figure. Distance to the nearest hydrologic feature is excluded because it is not significant in the retained model (p = 0.160); both precipitation terms are included because they are (p = 0.011 and 0.039). Woody cover and historic agriculture are significant but are shown in the interaction figure instead, so they are not repeated here.
smooth_vars_all <- c("mean_elevation", "nearest_road_m",
"pr_max_10b_2a", "pr_max_3b_2a")
smooth_effects_all <- do.call(rbind, lapply(smooth_vars_all, function(v) {
get_effect(mod_edge_best, landscape_data_clean, v)
}))
x_labels_all <- c(
"mean_elevation" = "Elevation (m)",
"nearest_road_m" = "Distance to Nearest Road (m)",
"pr_max_10b_2a" = "Maximum Precipitation (10-day, mm)",
"pr_max_3b_2a" = "Maximum Precipitation (3-day, mm)"
)
# Shared y-axis
y_range_smooth <- range(c(smooth_effects_all$lower, smooth_effects_all$upper), na.rm = TRUE)
all_plots <- lapply(smooth_vars_all, function(var) {
df <- smooth_effects_all %>% filter(variable == var)
ggplot(df, aes(x = x, y = fit)) +
geom_line(linewidth = 1.2, color = "#D55E00") +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "#D55E00") +
coord_cartesian(ylim = y_range_smooth) +
theme_classic(base_size = 14) +
theme(
axis.title = element_text(size = 16, face = "bold"),
axis.text = element_text(size = 14),
plot.tag = element_text(size = 18, face = "bold")
) +
labs(x = x_labels_all[var], y = "Edge Density")
})
combined_patchwork <- all_plots[[1]] + all_plots[[2]] +
all_plots[[3]] + all_plots[[4]] +
plot_layout(ncol = 2) +
plot_annotation(tag_levels = "A")
combined_patchwork
ggsave(file.path(fig_dir, "Edge_Density_Significant_Effects.png"),
combined_patchwork, width = 12, height = 10)Edge Density Interaction Plots
## Pastel palette for the five quantile groups, running cool (low) to warm (high)
## so the ordering of the quantiles reads intuitively. Change here to restyle all
## three panels at once.
pastel_pal <- c("#92C5DE", "#C5E0B4", "#F6D98A", "#F4A582", "#C994C7")
## 1. NDVI x SPI
ndvi_spi_lines <- get_interaction_lines(
mod_edge_best, landscape_data_clean, "NDVI_pre_365d", "spi1y"
)
ndvi_spi_plot <- ggplot(ndvi_spi_lines, aes(x = x, y = fit, color = group, fill = group)) +
geom_line(linewidth = 1.2) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.15, color = NA) +
theme_classic(base_size = 14) +
theme(
axis.title = element_text(size = 16, face = "bold"),
axis.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12)
) +
scale_color_manual(values = pastel_pal) +
scale_fill_manual(values = pastel_pal) +
labs(x = "NDVI (365-day)", y = "Edge Density",
color = "SPI (quantiles)", fill = "SPI (quantiles)")
## 2. Wood Cover x Slope
wood_slope_df <- get_interaction_lines(
mod_edge_best, landscape_data_clean, "wood_cov", "mean_slope"
)
wood_slope_plot <- ggplot(wood_slope_df, aes(x = x, y = fit, color = group, fill = group)) +
geom_line(linewidth = 1.2) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.15, color = NA) +
theme_classic(base_size = 14) +
theme(
axis.title = element_text(size = 16, face = "bold"),
axis.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12)
) +
scale_color_manual(values = pastel_pal) +
scale_fill_manual(values = pastel_pal) +
labs(x = "Wood Cover (%)", y = "Edge Density",
color = "Slope (quantiles)", fill = "Slope (quantiles)")
## 3. Agriculture x Slope
ag_slope_df <- get_interaction_lines(
mod_edge_best, landscape_data_clean, "ag_1985", "mean_slope"
)
ag_slope_plot <- ggplot(ag_slope_df, aes(x = x, y = fit, color = group, fill = group)) +
geom_line(linewidth = 1.2) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.15, color = NA) +
theme_classic(base_size = 14) +
theme(
axis.title = element_text(size = 16, face = "bold"),
axis.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12)
) +
scale_color_manual(values = pastel_pal) +
scale_fill_manual(values = pastel_pal) +
labs(x = "Agricultural Land Use (1985, %)", y = "Edge Density",
color = "Slope (quantiles)", fill = "Slope (quantiles)")
## Shared y-axis scale
combined_int_data <- rbind(ndvi_spi_lines, wood_slope_df, ag_slope_df)
y_range_int <- range(c(combined_int_data$lower, combined_int_data$upper), na.rm = TRUE)
ndvi_spi_plot <- ndvi_spi_plot + coord_cartesian(ylim = y_range_int)
wood_slope_plot <- wood_slope_plot + coord_cartesian(ylim = y_range_int)
ag_slope_plot <- ag_slope_plot + coord_cartesian(ylim = y_range_int)
final_plot <- ndvi_spi_plot + wood_slope_plot + ag_slope_plot +
plot_layout(ncol = 1) +
plot_annotation(tag_levels = "A")
final_plot
ggsave(file.path(fig_dir, "Edge_Density_Interactions.png"),
final_plot, width = 10, height = 18)PATCH AREA
Patch Area Log-normal
Baseline model, all observations.
mod_pa_ln <- gam(
log(avg_patch_area) ~
ag_1985 +
wood_cov +
s(road_dens_km_km2) +
nearest_road_m +
NDVI_pre_90d +
NDVI_pre_365d +
s(mean_hydro_dist_m) +
mean_elevation +
soil_order_1 +
vs_min_3b_2a +
tmmn_mea_60b_2a +
tmmx_max_3b_2a +
pr_mea_10b_2a +
s(vpd_min_10b_2a) +
spei1y +
spei14d,
data = landscape_data,
method = "REML"
)
summary(mod_pa_ln)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(avg_patch_area) ~ ag_1985 + wood_cov + s(road_dens_km_km2) +
## nearest_road_m + NDVI_pre_90d + NDVI_pre_365d + s(mean_hydro_dist_m) +
## mean_elevation + soil_order_1 + vs_min_3b_2a + tmmn_mea_60b_2a +
## tmmx_max_3b_2a + pr_mea_10b_2a + s(vpd_min_10b_2a) + spei1y +
## spei14d
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.3506115 3.4469229 1.262 0.20746
## ag_1985 -1.2137005 1.2104532 -1.003 0.31649
## wood_cov -0.0097382 0.0033543 -2.903 0.00385 **
## nearest_road_m -0.0000426 0.0001298 -0.328 0.74288
## NDVI_pre_90d -0.2726924 0.2778005 -0.982 0.32675
## NDVI_pre_365d 0.4372282 0.5831493 0.750 0.45374
## mean_elevation -0.0013383 0.0005324 -2.514 0.01226 *
## soil_order_11 0.1766769 0.1108523 1.594 0.11159
## soil_order_14 0.0351500 0.1443108 0.244 0.80766
## soil_order_19 -0.1429263 0.1271527 -1.124 0.26152
## soil_order_1Other 0.3192304 0.1537520 2.076 0.03837 *
## vs_min_3b_2a -0.0885299 0.0479380 -1.847 0.06536 .
## tmmn_mea_60b_2a -0.0338649 0.0110683 -3.060 0.00233 **
## tmmx_max_3b_2a 0.0164805 0.0118824 1.387 0.16606
## pr_mea_10b_2a 0.0022377 0.0134478 0.166 0.86791
## spei1y 0.0710355 0.0391043 1.817 0.06987 .
## spei14d 0.0494428 0.0525899 0.940 0.34758
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(road_dens_km_km2) 2.601 3.321 3.184 0.0198 *
## s(mean_hydro_dist_m) 2.011 2.461 1.460 0.2272
## s(vpd_min_10b_2a) 2.409 3.007 2.532 0.0562 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.1 Deviance explained = 13.9%
## -REML = 664.27 Scale est. = 0.58606 n = 538performance::check_outliers(mod_pa_ln)## OK: No outliers detected.
## - Based on the following method and threshold: cook (1).
## - For variable: (Whole model)## Plot
plot(mod_pa_ln, pages = 4, shade = TRUE, shade.col = "lightblue", all.terms = TRUE)## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter
## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter
## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in axis(1, at = seq_along(ll), labels = ll, ...): "shade" is not a
## graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter
## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter
Cleaned up log-normal patch area
sim_res_pa <- simulateResiduals(mod_pa_ln)
outlier_idx_pa <- outliers(sim_res_pa)
rows_used_pa <- as.numeric(rownames(mod_pa_ln$model))
if (length(outlier_idx_pa) > 0) {
landscape_data_pa_clean <- landscape_data[-rows_used_pa[outlier_idx_pa], ]
} else {
landscape_data_pa_clean <- landscape_data
}
nrow(landscape_data) - nrow(landscape_data_pa_clean) # n outliers dropped## [1] 1## Base model, outliers removed
mod_pa_clean <- gam(
log(avg_patch_area) ~
ag_1985 + wood_cov + s(road_dens_km_km2) + nearest_road_m +
NDVI_pre_90d + NDVI_pre_365d + s(mean_hydro_dist_m) +
mean_elevation + soil_order_1 +
vs_min_3b_2a + tmmn_mea_60b_2a + tmmx_max_3b_2a +
pr_mea_10b_2a + s(vpd_min_10b_2a) + spei1y + spei14d,
data = landscape_data_pa_clean,
method = "REML"
)
summary(mod_pa_clean)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(avg_patch_area) ~ ag_1985 + wood_cov + s(road_dens_km_km2) +
## nearest_road_m + NDVI_pre_90d + NDVI_pre_365d + s(mean_hydro_dist_m) +
## mean_elevation + soil_order_1 + vs_min_3b_2a + tmmn_mea_60b_2a +
## tmmx_max_3b_2a + pr_mea_10b_2a + s(vpd_min_10b_2a) + spei1y +
## spei14d
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.001e+00 3.432e+00 1.166 0.24422
## ag_1985 -1.205e+00 1.203e+00 -1.001 0.31715
## wood_cov -9.154e-03 3.342e-03 -2.739 0.00638 **
## nearest_road_m -5.139e-05 1.293e-04 -0.398 0.69113
## NDVI_pre_90d -3.066e-01 2.765e-01 -1.109 0.26799
## NDVI_pre_365d 4.259e-01 5.799e-01 0.734 0.46307
## mean_elevation -1.359e-03 5.295e-04 -2.567 0.01055 *
## soil_order_11 1.709e-01 1.102e-01 1.550 0.12167
## soil_order_14 -5.375e-03 1.443e-01 -0.037 0.97029
## soil_order_19 -1.444e-01 1.265e-01 -1.142 0.25419
## soil_order_1Other 3.177e-01 1.529e-01 2.078 0.03824 *
## vs_min_3b_2a -9.795e-02 4.778e-02 -2.050 0.04088 *
## tmmn_mea_60b_2a -3.348e-02 1.101e-02 -3.042 0.00247 **
## tmmx_max_3b_2a 1.728e-02 1.183e-02 1.461 0.14469
## pr_mea_10b_2a 3.443e-03 1.338e-02 0.257 0.79701
## spei1y 7.300e-02 3.888e-02 1.878 0.06101 .
## spei14d 4.183e-02 5.236e-02 0.799 0.42473
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(road_dens_km_km2) 2.534 3.238 2.940 0.0292 *
## s(mean_hydro_dist_m) 2.062 2.523 1.629 0.1830
## s(vpd_min_10b_2a) 2.477 3.091 2.671 0.0452 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.101 Deviance explained = 13.9%
## -REML = 660.1 Scale est. = 0.57918 n = 537## NDVI x drought interaction
mod_pa_NDVI_Climate <- gam(
log(avg_patch_area) ~
ag_1985 + wood_cov + s(road_dens_km_km2) + nearest_road_m +
NDVI_pre_90d + NDVI_pre_365d + s(mean_hydro_dist_m) +
mean_elevation + soil_order_1 +
vs_min_3b_2a + tmmn_mea_60b_2a + tmmx_max_3b_2a +
pr_mea_10b_2a + s(vpd_min_10b_2a) + spei1y + spei14d +
ti(NDVI_pre_365d, spei1y),
data = landscape_data_pa_clean,
method = "REML"
)
summary(mod_pa_NDVI_Climate)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(avg_patch_area) ~ ag_1985 + wood_cov + s(road_dens_km_km2) +
## nearest_road_m + NDVI_pre_90d + NDVI_pre_365d + s(mean_hydro_dist_m) +
## mean_elevation + soil_order_1 + vs_min_3b_2a + tmmn_mea_60b_2a +
## tmmx_max_3b_2a + pr_mea_10b_2a + s(vpd_min_10b_2a) + spei1y +
## spei14d + ti(NDVI_pre_365d, spei1y)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.394e+00 3.417e+00 1.286 0.19900
## ag_1985 -1.300e+00 1.197e+00 -1.087 0.27773
## wood_cov -8.983e-03 3.336e-03 -2.693 0.00732 **
## nearest_road_m -5.274e-05 1.301e-04 -0.405 0.68534
## NDVI_pre_90d -3.162e-01 2.769e-01 -1.142 0.25398
## NDVI_pre_365d 1.477e-01 6.135e-01 0.241 0.80986
## mean_elevation -1.366e-03 5.292e-04 -2.581 0.01013 *
## soil_order_11 1.631e-01 1.099e-01 1.484 0.13845
## soil_order_14 1.641e-03 1.432e-01 0.011 0.99086
## soil_order_19 -1.418e-01 1.274e-01 -1.113 0.26625
## soil_order_1Other 2.986e-01 1.521e-01 1.963 0.05020 .
## vs_min_3b_2a -9.924e-02 4.750e-02 -2.089 0.03717 *
## tmmn_mea_60b_2a -3.307e-02 1.097e-02 -3.014 0.00271 **
## tmmx_max_3b_2a 1.605e-02 1.181e-02 1.360 0.17456
## pr_mea_10b_2a 2.558e-03 1.331e-02 0.192 0.84768
## spei1y 7.543e-02 3.873e-02 1.948 0.05201 .
## spei14d 4.173e-02 5.215e-02 0.800 0.42402
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(road_dens_km_km2) 2.633 3.360 3.206 0.0190 *
## s(mean_hydro_dist_m) 1.890 2.315 1.116 0.3254
## s(vpd_min_10b_2a) 2.651 3.309 3.112 0.0219 *
## ti(NDVI_pre_365d,spei1y) 2.772 3.294 2.145 0.0717 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.112 Deviance explained = 15.5%
## -REML = 658.14 Scale est. = 0.57157 n = 537## Land cover x topography interactions
mod_pa_LandCov_Elev <- gam(
log(avg_patch_area) ~
ag_1985 + wood_cov + s(road_dens_km_km2) + nearest_road_m +
NDVI_pre_90d + NDVI_pre_365d + s(mean_hydro_dist_m) +
mean_elevation + soil_order_1 +
vs_min_3b_2a + tmmn_mea_60b_2a + tmmx_max_3b_2a +
pr_mea_10b_2a + s(vpd_min_10b_2a) + spei1y + spei14d +
ti(wood_cov, mean_elevation) +
ti(ag_1985, mean_elevation),
data = landscape_data_pa_clean,
method = "REML"
)
summary(mod_pa_LandCov_Elev)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(avg_patch_area) ~ ag_1985 + wood_cov + s(road_dens_km_km2) +
## nearest_road_m + NDVI_pre_90d + NDVI_pre_365d + s(mean_hydro_dist_m) +
## mean_elevation + soil_order_1 + vs_min_3b_2a + tmmn_mea_60b_2a +
## tmmx_max_3b_2a + pr_mea_10b_2a + s(vpd_min_10b_2a) + spei1y +
## spei14d + ti(wood_cov, mean_elevation) + ti(ag_1985, mean_elevation)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.449e+00 3.451e+00 1.000 0.31794
## ag_1985 -1.351e+00 1.288e+00 -1.049 0.29461
## wood_cov -9.905e-03 3.512e-03 -2.821 0.00498 **
## nearest_road_m -5.756e-05 1.294e-04 -0.445 0.65658
## NDVI_pre_90d -3.122e-01 2.764e-01 -1.130 0.25912
## NDVI_pre_365d 4.661e-01 5.848e-01 0.797 0.42579
## mean_elevation -8.599e-04 8.923e-04 -0.964 0.33571
## soil_order_11 1.759e-01 1.107e-01 1.590 0.11252
## soil_order_14 5.794e-02 1.504e-01 0.385 0.70021
## soil_order_19 -1.326e-01 1.274e-01 -1.041 0.29837
## soil_order_1Other 3.144e-01 1.527e-01 2.059 0.04002 *
## vs_min_3b_2a -1.039e-01 4.785e-02 -2.170 0.03044 *
## tmmn_mea_60b_2a -3.208e-02 1.111e-02 -2.888 0.00404 **
## tmmx_max_3b_2a 1.788e-02 1.188e-02 1.505 0.13303
## pr_mea_10b_2a 4.583e-03 1.341e-02 0.342 0.73261
## spei1y 7.849e-02 3.902e-02 2.011 0.04482 *
## spei14d 2.926e-02 5.290e-02 0.553 0.58042
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(road_dens_km_km2) 2.450 3.131 2.361 0.0670 .
## s(mean_hydro_dist_m) 1.976 2.416 1.323 0.2512
## s(vpd_min_10b_2a) 2.513 3.134 2.818 0.0366 *
## ti(wood_cov,mean_elevation) 1.901 2.565 0.971 0.4085
## ti(ag_1985,mean_elevation) 1.386 1.669 0.358 0.5457
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.103 Deviance explained = 14.7%
## -REML = 656.33 Scale est. = 0.57743 n = 537## Combined interaction model
mod_pa_Int <- gam(
log(avg_patch_area) ~
ag_1985 + wood_cov + s(road_dens_km_km2) + nearest_road_m +
NDVI_pre_90d + NDVI_pre_365d + s(mean_hydro_dist_m) +
mean_elevation + soil_order_1 +
vs_min_3b_2a + tmmn_mea_60b_2a + tmmx_max_3b_2a +
pr_mea_10b_2a + s(vpd_min_10b_2a) + spei1y + spei14d +
ti(NDVI_pre_365d, spei1y) +
ti(wood_cov, mean_elevation) +
ag_1985:mean_elevation,
data = landscape_data_pa_clean,
method = "REML"
)
summary(mod_pa_Int)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(avg_patch_area) ~ ag_1985 + wood_cov + s(road_dens_km_km2) +
## nearest_road_m + NDVI_pre_90d + NDVI_pre_365d + s(mean_hydro_dist_m) +
## mean_elevation + soil_order_1 + vs_min_3b_2a + tmmn_mea_60b_2a +
## tmmx_max_3b_2a + pr_mea_10b_2a + s(vpd_min_10b_2a) + spei1y +
## spei14d + ti(NDVI_pre_365d, spei1y) + ti(wood_cov, mean_elevation) +
## ag_1985:mean_elevation
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.954e+00 3.439e+00 1.150 0.25076
## ag_1985 -3.943e-01 1.627e+00 -0.242 0.80865
## wood_cov -1.007e-02 3.488e-03 -2.887 0.00406 **
## nearest_road_m -6.027e-05 1.301e-04 -0.463 0.64344
## NDVI_pre_90d -3.134e-01 2.772e-01 -1.131 0.25879
## NDVI_pre_365d 1.360e-01 6.151e-01 0.221 0.82516
## mean_elevation -5.297e-04 8.860e-04 -0.598 0.55020
## soil_order_11 1.601e-01 1.100e-01 1.455 0.14631
## soil_order_14 4.973e-02 1.491e-01 0.333 0.73893
## soil_order_19 -1.285e-01 1.282e-01 -1.002 0.31660
## soil_order_1Other 2.998e-01 1.522e-01 1.970 0.04936 *
## vs_min_3b_2a -1.012e-01 4.756e-02 -2.127 0.03392 *
## tmmn_mea_60b_2a -3.190e-02 1.108e-02 -2.879 0.00416 **
## tmmx_max_3b_2a 1.655e-02 1.187e-02 1.395 0.16367
## pr_mea_10b_2a 2.813e-03 1.332e-02 0.211 0.83284
## spei1y 7.950e-02 3.890e-02 2.044 0.04149 *
## spei14d 3.593e-02 5.252e-02 0.684 0.49416
## ag_1985:mean_elevation -1.289e-02 1.846e-02 -0.699 0.48517
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(road_dens_km_km2) 2.553 3.260 2.725 0.0387 *
## s(mean_hydro_dist_m) 1.825 2.237 0.954 0.3790
## s(vpd_min_10b_2a) 2.640 3.295 3.104 0.0226 *
## ti(NDVI_pre_365d,spei1y) 2.737 3.260 1.995 0.0894 .
## ti(wood_cov,mean_elevation) 1.009 1.017 1.252 0.2674
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.111 Deviance explained = 15.7%
## -REML = 659.34 Scale est. = 0.57262 n = 537## Model Selection
mod_selection_pa <- model.sel(mod_pa_clean, mod_pa_NDVI_Climate,
mod_pa_LandCov_Elev, mod_pa_Int)
mod_selection_pa## Model selection table
## (Int) ag_1985 men_elv NDV_pre_365 NDV_pre_90d nrs_rod_m
## mod_pa_NDVI_Climate 4.394 -1.3000 -0.0013660 0.1477 -0.3162 -5.274e-05
## mod_pa_clean 4.001 -1.2050 -0.0013590 0.4259 -0.3066 -5.139e-05
## mod_pa_Int 3.954 -0.3943 -0.0005297 0.1360 -0.3134 -6.027e-05
## mod_pa_LandCov_Elev 3.449 -1.3510 -0.0008598 0.4661 -0.3122 -5.756e-05
## pr_mea_10b_2a s(men_hyd_dst_m) s(rod_dns_km_km2)
## mod_pa_NDVI_Climate 0.002558 + +
## mod_pa_clean 0.003443 + +
## mod_pa_Int 0.002813 + +
## mod_pa_LandCov_Elev 0.004583 + +
## s(vpd_min_10b_2a) sol_ord_1 s14 sp1 tmmn_mea_60b_2a
## mod_pa_NDVI_Climate + + 0.04173 0.07543 -0.03307
## mod_pa_clean + + 0.04183 0.07300 -0.03348
## mod_pa_Int + + 0.03593 0.07950 -0.03190
## mod_pa_LandCov_Elev + + 0.02926 0.07849 -0.03208
## tmmx_max_3b_2a vs_min_3b_2a wod_cov ti(NDV_pre_365,sp1)
## mod_pa_NDVI_Climate 0.01605 -0.09924 -0.008983 +
## mod_pa_clean 0.01728 -0.09795 -0.009154
## mod_pa_Int 0.01655 -0.10120 -0.010070 +
## mod_pa_LandCov_Elev 0.01788 -0.10390 -0.009905
## ti(ag_1985,men_elv) ti(wod_cov,men_elv) ag_1985:men_elv df
## mod_pa_NDVI_Climate 30
## mod_pa_clean 26
## mod_pa_Int + -0.01289 32
## mod_pa_LandCov_Elev + + 30
## logLik AICc delta weight
## mod_pa_NDVI_Climate -597.959 1260.2 0.00 0.649
## mod_pa_clean -603.016 1262.7 2.46 0.190
## mod_pa_Int -597.492 1263.3 3.11 0.137
## mod_pa_LandCov_Elev -600.548 1266.8 6.62 0.024
## Models ranked by AICc(x)## Selected model: mod_pa_NDVI_Climate (lowest AICc, weight = 0.66;
## mod_pa_Int delta = 2.42, mod_pa_clean delta = 3.32)
mod_pa_best <- mod_pa_NDVI_Climate
summary(mod_pa_best)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(avg_patch_area) ~ ag_1985 + wood_cov + s(road_dens_km_km2) +
## nearest_road_m + NDVI_pre_90d + NDVI_pre_365d + s(mean_hydro_dist_m) +
## mean_elevation + soil_order_1 + vs_min_3b_2a + tmmn_mea_60b_2a +
## tmmx_max_3b_2a + pr_mea_10b_2a + s(vpd_min_10b_2a) + spei1y +
## spei14d + ti(NDVI_pre_365d, spei1y)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.394e+00 3.417e+00 1.286 0.19900
## ag_1985 -1.300e+00 1.197e+00 -1.087 0.27773
## wood_cov -8.983e-03 3.336e-03 -2.693 0.00732 **
## nearest_road_m -5.274e-05 1.301e-04 -0.405 0.68534
## NDVI_pre_90d -3.162e-01 2.769e-01 -1.142 0.25398
## NDVI_pre_365d 1.477e-01 6.135e-01 0.241 0.80986
## mean_elevation -1.366e-03 5.292e-04 -2.581 0.01013 *
## soil_order_11 1.631e-01 1.099e-01 1.484 0.13845
## soil_order_14 1.641e-03 1.432e-01 0.011 0.99086
## soil_order_19 -1.418e-01 1.274e-01 -1.113 0.26625
## soil_order_1Other 2.986e-01 1.521e-01 1.963 0.05020 .
## vs_min_3b_2a -9.924e-02 4.750e-02 -2.089 0.03717 *
## tmmn_mea_60b_2a -3.307e-02 1.097e-02 -3.014 0.00271 **
## tmmx_max_3b_2a 1.605e-02 1.181e-02 1.360 0.17456
## pr_mea_10b_2a 2.558e-03 1.331e-02 0.192 0.84768
## spei1y 7.543e-02 3.873e-02 1.948 0.05201 .
## spei14d 4.173e-02 5.215e-02 0.800 0.42402
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(road_dens_km_km2) 2.633 3.360 3.206 0.0190 *
## s(mean_hydro_dist_m) 1.890 2.315 1.116 0.3254
## s(vpd_min_10b_2a) 2.651 3.309 3.112 0.0219 *
## ti(NDVI_pre_365d,spei1y) 2.772 3.294 2.145 0.0717 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.112 Deviance explained = 15.5%
## -REML = 658.14 Scale est. = 0.57157 n = 537# Model Fit
performance::check_collinearity(mod_pa_best)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
## ag_1985 1.31 [ 1.21, 1.48] 1.15 0.76 [0.68, 0.83]
## wood_cov 2.70 [ 2.38, 3.10] 1.64 0.37 [0.32, 0.42]
## nearest_road_m 1.42 [ 1.30, 1.60] 1.19 0.70 [0.63, 0.77]
## NDVI_pre_90d 1.54 [ 1.40, 1.74] 1.24 0.65 [0.58, 0.72]
## mean_elevation 1.52 [ 1.38, 1.71] 1.23 0.66 [0.58, 0.73]
## soil_order_1 3.03 [ 2.66, 3.48] 1.15 0.33 [0.29, 0.38]
## vs_min_3b_2a 1.22 [ 1.13, 1.38] 1.11 0.82 [0.73, 0.88]
## tmmn_mea_60b_2a 3.20 [ 2.80, 3.68] 1.79 0.31 [0.27, 0.36]
## tmmx_max_3b_2a 2.45 [ 2.17, 2.80] 1.56 0.41 [0.36, 0.46]
## pr_mea_10b_2a 2.00 [ 1.79, 2.28] 1.41 0.50 [0.44, 0.56]
## s(vpd_min_10b_2a) 2.10 [ 1.87, 2.40] 1.45 0.48 [0.42, 0.53]
## spei1y 1.73 [ 1.56, 1.96] 1.32 0.58 [0.51, 0.64]
## spei14d 2.20 [ 1.96, 2.51] 1.48 0.45 [0.40, 0.51]
##
## Moderate Correlation
##
## Term VIF VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
## s(road_dens_km_km2) 9.68 [ 8.32, 11.30] 3.11 0.10 [0.09, 0.12]
##
## High Correlation
##
## Term VIF VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
## NDVI_pre_365d 31.44 [26.82, 36.88] 5.61 0.03 [0.03, 0.04]
## s(mean_hydro_dist_m) 19.83 [16.95, 23.23] 4.45 0.05 [0.04, 0.06]performance::check_outliers(mod_pa_best)## OK: No outliers detected.
## - Based on the following method and threshold: cook (1).
## - For variable: (Whole model)## Plot
plot(mod_pa_best, pages = 4, shade = TRUE, shade.col = "lightblue", all.terms = TRUE)## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter
## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in axis(1, at = seq_along(ll), labels = ll, ...): "shade" is not a
## graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "shade" is not a
## graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter
## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter## Warning in plot.window(...): "shade" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "shade" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "shade" is not a
## graphical parameter## Warning in box(...): "shade" is not a graphical parameter## Warning in title(...): "shade" is not a graphical parameter
Patch Area Cross-Validation
pa_candidates <- list(
mod_pa_clean = mod_pa_clean,
mod_pa_NDVI_Climate = mod_pa_NDVI_Climate,
mod_pa_LandCov_Elev = mod_pa_LandCov_Elev,
mod_pa_Int = mod_pa_Int
)
cv_pa <- cv_compare(pa_candidates, landscape_data_pa_clean,
K = 10, reps = 5, seed = 1)
cv_pa## model CV_RMSE RMSE_SD CV_MAE CV_R2 R2_SD
## 1 mod_pa_NDVI_Climate 0.7784013 0.003916184 0.6332445 0.05732665 0.009489843
## 2 mod_pa_clean 0.7802679 0.003579629 0.6362259 0.05280350 0.008689951
## 3 mod_pa_Int 0.7837982 0.002965821 0.6385766 0.04421812 0.007224407
## 4 mod_pa_LandCov_Elev 0.7904865 0.008415075 0.6432290 0.02775978 0.020794498
## delta_RMSE
## 1 0.000000000
## 2 0.001866563
## 3 0.005396868
## 4 0.012085168write.csv(cv_pa, "patch_area_cv_results.csv", row.names = FALSE)Patch Area back-transformed table
pa_tables <- backtransform_tables(mod_pa_best)
pa_tables$param## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.394206e+00 3.4167789673 1.28606676 0.199003266
## ag_1985 -1.300224e+00 1.1966041358 -1.08659520 0.277728985
## wood_cov -8.982941e-03 0.0033359428 -2.69277423 0.007319141
## nearest_road_m -5.273982e-05 0.0001300853 -0.40542501 0.685335128
## NDVI_pre_90d -3.162450e-01 0.2769142817 -1.14203206 0.253976529
## NDVI_pre_365d 1.476800e-01 0.6134721547 0.24072815 0.809862554
## mean_elevation -1.365839e-03 0.0005291711 -2.58109108 0.010126993
## soil_order_11 1.630897e-01 0.1099042766 1.48392462 0.138446595
## soil_order_14 1.641478e-03 0.1431732303 0.01146498 0.990856952
## soil_order_19 -1.417560e-01 0.1273673081 -1.11297043 0.266245379
## soil_order_1Other 2.985893e-01 0.1521133339 1.96293978 0.050196203
## vs_min_3b_2a -9.923651e-02 0.0474969275 -2.08932483 0.037173434
## tmmn_mea_60b_2a -3.307292e-02 0.0109731143 -3.01399580 0.002706580
## tmmx_max_3b_2a 1.605444e-02 0.0118082028 1.35960021 0.174557202
## pr_mea_10b_2a 2.558378e-03 0.0133123791 0.19218036 0.847677421
## spei1y 7.543448e-02 0.0387322631 1.94758773 0.052012612
## spei14d 4.172972e-02 0.0521546480 0.80011503 0.424016643
## exp_estimate percent_change
## (Intercept) 80.980 7998.0
## ag_1985 0.272 -72.8
## wood_cov 0.991 -0.9
## nearest_road_m 1.000 0.0
## NDVI_pre_90d 0.729 -27.1
## NDVI_pre_365d 1.159 15.9
## mean_elevation 0.999 -0.1
## soil_order_11 1.177 17.7
## soil_order_14 1.002 0.2
## soil_order_19 0.868 -13.2
## soil_order_1Other 1.348 34.8
## vs_min_3b_2a 0.906 -9.4
## tmmn_mea_60b_2a 0.967 -3.3
## tmmx_max_3b_2a 1.016 1.6
## pr_mea_10b_2a 1.003 0.3
## spei1y 1.078 7.8
## spei14d 1.043 4.3pa_tables$smooth## edf Ref.df F p-value p_value
## s(road_dens_km_km2) 2.63 3.359672 3.21 0.01895920 0.0190
## s(mean_hydro_dist_m) 1.89 2.314991 1.12 0.32541785 0.3254
## s(vpd_min_10b_2a) 2.65 3.308568 3.11 0.02194575 0.0219
## ti(NDVI_pre_365d,spei1y) 2.77 3.293882 2.15 0.07173445 0.0717write.csv(pa_tables$param, "patch_area_backtransformed_param_summary.csv", row.names = TRUE)
write.csv(pa_tables$smooth, "patch_area_smooth_summary.csv", row.names = TRUE)Patch Area Significant Effects
All six significant terms in the retained patch area model on one figure, sharing a common y-axis so effect magnitudes are directly comparable. Distance to the nearest hydrologic feature is excluded (p = 0.325), as is the NDVI x SPEI interaction (p = 0.072).
pa_sig_vars <- c("tmmn_mea_60b_2a", "vs_min_3b_2a", "vpd_min_10b_2a",
"wood_cov", "mean_elevation", "road_dens_km_km2")
pa_sig_effects <- do.call(rbind, lapply(pa_sig_vars, function(v) {
get_effect(mod_pa_best, landscape_data_pa_clean, v)
}))
x_labels_pa <- c(
"tmmn_mea_60b_2a" = "Mean Minimum Temperature (60-day, K)",
"vs_min_3b_2a" = "Minimum Wind Speed (3-day, m/s)",
"vpd_min_10b_2a" = "Minimum VPD (10-day, kPa)",
"wood_cov" = "Wood Cover (%)",
"mean_elevation" = "Elevation (m)",
"road_dens_km_km2" = "Road Density (km/km\u00b2)"
)
# Shared y-axis across all six panels
y_range_pa <- range(c(pa_sig_effects$lower, pa_sig_effects$upper), na.rm = TRUE)
pa_plots <- lapply(pa_sig_vars, function(var) {
df <- pa_sig_effects %>% filter(variable == var)
ggplot(df, aes(x = x, y = fit)) +
geom_line(linewidth = 1.2, color = "#0072B2") +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "#0072B2") +
coord_cartesian(ylim = y_range_pa) +
theme_classic(base_size = 14) +
theme(
axis.title = element_text(size = 15, face = "bold"),
axis.text = element_text(size = 13),
plot.tag = element_text(size = 18, face = "bold")
) +
labs(x = x_labels_pa[var], y = "Predicted Patch Area (ha)")
})
combined_pa_plot <- pa_plots[[1]] + pa_plots[[2]] + pa_plots[[3]] +
pa_plots[[4]] + pa_plots[[5]] + pa_plots[[6]] +
plot_layout(ncol = 3) +
plot_annotation(tag_levels = "A")
combined_pa_plot
ggsave(file.path(fig_dir, "Patch_Area_Significant_Effects.png"),
combined_pa_plot, width = 15, height = 9, dpi = 300)SUPPLEMENTAL FIGURES
Every one-dimensional term in each retained model, on the back-transformed response scale, with all other predictors held at their mean (numeric) or first level (categorical). Shaded bands are ±2 SE. Stars in the panel titles are the term’s significance in the model summary: *** p < 0.001, ** p < 0.01, * p < 0.05, . p < 0.10.
Supplemental: all edge density effects
supp_edge <- supplemental_effects(mod_edge_best, landscape_data_clean)
supp_edge_plot <- plot_supplemental(supp_edge,
y_lab = "Predicted Edge Density (m/ha)",
ncol = 4, colour = "#D55E00")
supp_edge_plot
ggsave(file.path(fig_dir, "FigS1_Edge_Density_All_Effects.png"),
supp_edge_plot, width = 16, height = 16, dpi = 300)
## Categorical predictors (soil order)
supp_edge_fac <- supplemental_factor_effects(mod_edge_best, landscape_data_clean)
if (!is.null(supp_edge_fac)) {
supp_edge_fac_plot <- plot_factor_effects(supp_edge_fac,
y_lab = "Predicted Edge Density (m/ha)",
colour = "#D55E00")
print(supp_edge_fac_plot)
ggsave(file.path(fig_dir, "FigS1b_Edge_Density_Factor_Effects.png"),
supp_edge_fac_plot, width = 7, height = 5, dpi = 300)
}
Supplemental: all patch area effects
supp_pa <- supplemental_effects(mod_pa_best, landscape_data_pa_clean)
supp_pa_plot <- plot_supplemental(supp_pa,
y_lab = "Predicted Patch Area (ha)",
ncol = 4, colour = "#0072B2")
supp_pa_plot
ggsave(file.path(fig_dir, "FigS2_Patch_Area_All_Effects.png"),
supp_pa_plot, width = 16, height = 16, dpi = 300)
## Categorical predictors (soil order)
supp_pa_fac <- supplemental_factor_effects(mod_pa_best, landscape_data_pa_clean)
if (!is.null(supp_pa_fac)) {
supp_pa_fac_plot <- plot_factor_effects(supp_pa_fac,
y_lab = "Predicted Patch Area (ha)",
colour = "#0072B2")
print(supp_pa_fac_plot)
ggsave(file.path(fig_dir, "FigS2b_Patch_Area_Factor_Effects.png"),
supp_pa_fac_plot, width = 7, height = 5, dpi = 300)
}
Supplemental: direction of each fitted relationship
A compact numeric summary of each panel — the fitted value at the low and high end of the observed range, and whether the relationship is monotonic — so the direction of every effect can be stated in the text without reading it off the figure by eye.
describe_direction <- function(eff) {
eff %>%
group_by(variable, panel) %>%
arrange(x, .by_group = TRUE) %>%
summarise(
fit_low = round(dplyr::first(fit), 3),
fit_high = round(dplyr::last(fit), 3),
fit_min = round(min(fit), 3),
fit_max = round(max(fit), 3),
x_at_max = round(x[which.max(fit)], 3),
x_at_min = round(x[which.min(fit)], 3),
shape = dplyr::case_when(
which.max(fit) == dplyr::n() & which.min(fit) == 1 ~ "increasing",
which.max(fit) == 1 & which.min(fit) == dplyr::n() ~ "decreasing",
which.max(fit) > 1 & which.max(fit) < dplyr::n() ~ "peak at intermediate",
TRUE ~ "trough at intermediate"
),
.groups = "drop"
)
}
direction_edge <- describe_direction(supp_edge)
direction_pa <- describe_direction(supp_pa)
direction_edge## # A tibble: 16 × 9
## variable panel fit_low fit_high fit_min fit_max x_at_max x_at_min shape
## <chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 NDVI_pre_365d NDVI… 792. 715. 497. 792. 0.227 0.462 trou…
## 2 NDVI_pre_90d NDVI… 525. 495. 495. 525. -0.147 0.712 decr…
## 3 ag_1985 Agri… 483. 1028. 483. 1028. 0.294 0 incr…
## 4 mean_elevation Elev… 879. 1479. 498. 1810. 259. 74.8 peak…
## 5 mean_hydro_di… Dist… 542. 892. 439. 892. 4072. 1316. trou…
## 6 mean_slope Slope 471. 562. 471. 572. 9.59 1.61 peak…
## 7 nearest_road_m Dist… 501. 2056. 501. 2056. 6033. 0 incr…
## 8 pr_max_10b_2a Maxi… 589. 774. 421. 774. 239. 86.9 trou…
## 9 pr_max_3b_2a Maxi… 483. 678. 483. 678. 120. 0 incr…
## 10 pr_mea_60b_2a Mean… 550. 450. 450. 550. 0.325 10.7 decr…
## 11 rmin_mea_60b_… Mean… 491. 519. 491. 519. 55.7 29.2 incr…
## 12 spi1y SPI … 651. 530. 497. 669. -1.41 0.19 peak…
## 13 tmmn_max_3b_2a Max … 449. 540. 449. 540. 298. 267. incr…
## 14 vs_mea_30b_2a Mean… 496. 522. 496. 522. 6.11 2.73 incr…
## 15 vs_mea_3b_2a Mean… 522. 483. 483. 522. 1.95 7.39 decr…
## 16 wood_cov Wood… 304. 591. 304. 591. 97.4 20.1 incr…direction_pa## # A tibble: 15 × 9
## variable panel fit_low fit_high fit_min fit_max x_at_max x_at_min shape
## <chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 NDVI_pre_365d NDVI… 0.286 0.304 0.286 0.304 0.651 0.227 incr…
## 2 NDVI_pre_90d NDVI… 0.353 0.267 0.267 0.353 -0.147 0.712 decr…
## 3 ag_1985 Agri… 0.296 0.214 0.214 0.296 0 0.294 decr…
## 4 mean_elevation Elev… 0.324 0.191 0.191 0.324 6.55 404. decr…
## 5 mean_hydro_di… Dist… 0.273 0.235 0.235 0.368 1645. 4072. peak…
## 6 nearest_road_m Dist… 0.293 0.291 0.27 0.293 0 3108. trou…
## 7 pr_mea_10b_2a Mean… 0.29 0.313 0.29 0.313 21.3 0 incr…
## 8 road_dens_km_… Road… 0.2 0.207 0.2 0.306 1.87 0 peak…
## 9 spei14d SPEI… 0.268 0.318 0.268 0.318 2.09 -2.09 incr…
## 10 spei1y SPEI… 0.229 0.396 0.229 0.396 2.09 -2.09 incr…
## 11 tmmn_mea_60b_… Mean… 0.435 0.18 0.18 0.435 269. 295. decr…
## 12 tmmx_max_3b_2a Maxi… 0.235 0.353 0.235 0.353 311. 285. incr…
## 13 vpd_min_10b_2a Mini… 0.375 0.439 0.284 0.439 1.18 0.371 trou…
## 14 vs_min_3b_2a Mini… 0.344 0.212 0.212 0.344 0.984 5.97 decr…
## 15 wood_cov Wood… 0.496 0.243 0.243 0.496 20.1 97.4 decr…write.csv(direction_edge, "edge_density_effect_directions.csv", row.names = FALSE)
write.csv(direction_pa, "patch_area_effect_directions.csv", row.names = FALSE)