Refugia Landscape Modeling
Alan Ivory
2026-04-08
Selected Variables
# 1. Define the selected variable sets
selected_models <- list(
Clumpiness = c("ag_1985", "wood_cov", "road_dens_km_km2",
"nearest_road_m", "NDVI_pre_90d", "NDVI_pre_365d", "mean_hydro_dist_m",
"mean_slope", "soil_order_1", "vs_max_3b_2a", "vs_mea_60b_2a",
"tmmn_max_3b_2a"),
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")
)Explore Distributions of Landscape Metrics
summary(landscape_data$clumpiness)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.2619 0.8533 0.8750 0.8705 0.8972 1.0000 5summary(landscape_data$edge_dens)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 141.6 489.8 786.0 1000.8 1164.6 10354.5summary(landscape_data$avg_patch_area)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.08912 0.17099 0.28256 0.43624 0.53432 4.91433hist(landscape_data$avg_patch_area)
hist(landscape_data$edge_dens)
hist(landscape_data$clumpiness)
Prepare Data for Modeling
# Shrink clumpiness away from 0/1 for quasibinomial
eps <- 0.001
landscape_data$clumpiness_glm <- pmin(pmax(landscape_data$clumpiness, eps), 1 - eps)Model Function
make_formula <- function(response, predictors) {
as.formula(
paste(response, "~", paste(predictors, collapse = " + "))
)
}GLM Models
# CLUMPINESS ------------------------------------------------
clump_form <- make_formula(
"clumpiness",
selected_models$Clumpiness
)
clump_glm <- glm(
clump_form,
data = landscape_data,
family = gaussian(),
na.action = na.omit
)
# EDGE DENSITY ----------------------------------------------
edge_form <- make_formula(
"edge_dens",
selected_models$Edge_Density
)
edge_glm <- glm(
edge_form,
data = landscape_data,
family = Gamma(link = "log"),
na.action = na.omit
)
# AVERAGE PATCH AREA ----------------------------------------
area_form <- make_formula(
"avg_patch_area",
selected_models$Avg_Patch_Area
)
area_glm <- glm(
area_form,
data = landscape_data,
family = Gamma(link = "log"),
na.action = na.omit
)Model Diagnostics
# Model summaries
summary(clump_glm)##
## Call:
## glm(formula = clump_form, family = gaussian(), data = landscape_data,
## na.action = na.omit)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.102e-01 1.477e-01 4.131 4.2e-05 ***
## ag_1985 -1.132e-01 8.620e-02 -1.314 0.18953
## wood_cov -4.643e-04 1.984e-04 -2.340 0.01964 *
## road_dens_km_km2 1.159e-04 3.631e-03 0.032 0.97456
## nearest_road_m -1.045e-05 8.186e-06 -1.276 0.20250
## NDVI_pre_90d 2.102e-02 1.818e-02 1.156 0.24802
## NDVI_pre_365d 1.143e-02 3.803e-02 0.300 0.76396
## mean_hydro_dist_m -2.687e-06 7.400e-06 -0.363 0.71663
## mean_slope -3.981e-03 1.473e-03 -2.703 0.00711 **
## soil_order_1 -1.496e-03 7.795e-04 -1.920 0.05542 .
## vs_max_3b_2a -1.415e-03 2.043e-03 -0.693 0.48885
## vs_mea_60b_2a -1.809e-03 6.006e-03 -0.301 0.76333
## tmmn_max_3b_2a 1.142e-03 4.950e-04 2.307 0.02145 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.003005342)
##
## Null deviance: 1.7044 on 532 degrees of freedom
## Residual deviance: 1.5628 on 520 degrees of freedom
## (5 observations deleted due to missingness)
## AIC: -1567.9
##
## Number of Fisher Scoring iterations: 2summary(edge_glm)##
## Call:
## glm(formula = edge_form, family = Gamma(link = "log"), data = landscape_data,
## na.action = na.omit)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.353e+00 2.300e+00 1.458 0.1455
## ag_1985 3.221e+00 1.248e+00 2.580 0.0101 *
## wood_cov 1.457e-02 3.000e-03 4.856 1.59e-06 ***
## nearest_road_m 1.979e-04 1.210e-04 1.635 0.1027
## NDVI_pre_90d 8.931e-02 2.917e-01 0.306 0.7596
## NDVI_pre_365d 2.876e-02 6.204e-01 0.046 0.9630
## mean_hydro_dist_m -4.415e-05 1.130e-04 -0.391 0.6962
## mean_slope 1.641e-02 2.689e-02 0.610 0.5420
## mean_elevation 8.082e-04 6.614e-04 1.222 0.2223
## soil_order_1 1.147e-02 1.173e-02 0.978 0.3286
## vs_mea_3b_2a -7.248e-02 5.832e-02 -1.243 0.2145
## vs_mea_30b_2a 1.873e-02 8.950e-02 0.209 0.8343
## tmmn_max_3b_2a 1.175e-02 7.929e-03 1.482 0.1389
## pr_max_3b_2a -1.517e-04 2.057e-03 -0.074 0.9412
## pr_max_10b_2a -5.897e-04 1.711e-03 -0.345 0.7306
## pr_mea_60b_2a -3.805e-02 2.798e-02 -1.360 0.1744
## rmin_mea_60b_2a -1.976e-02 1.020e-02 -1.938 0.0532 .
## spi1y 8.733e-02 4.839e-02 1.805 0.0717 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.671667)
##
## Null deviance: 267.17 on 537 degrees of freedom
## Residual deviance: 223.67 on 520 degrees of freedom
## AIC: 8296.1
##
## Number of Fisher Scoring iterations: 12summary(area_glm)##
## Call:
## glm(formula = area_form, family = Gamma(link = "log"), data = landscape_data,
## na.action = na.omit)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.383e+00 4.047e+00 0.342 0.73278
## ag_1985 -1.233e+00 1.439e+00 -0.857 0.39196
## wood_cov -8.113e-03 3.622e-03 -2.240 0.02553 *
## road_dens_km_km2 1.105e-01 6.167e-02 1.792 0.07368 .
## nearest_road_m -1.142e-04 1.428e-04 -0.800 0.42426
## NDVI_pre_90d -1.568e-01 3.296e-01 -0.476 0.63438
## NDVI_pre_365d 4.425e-01 6.908e-01 0.641 0.52204
## mean_hydro_dist_m -6.585e-05 1.262e-04 -0.522 0.60198
## mean_elevation -1.314e-03 6.290e-04 -2.090 0.03714 *
## soil_order_1 -3.231e-02 1.325e-02 -2.438 0.01510 *
## vs_min_3b_2a -1.049e-01 5.741e-02 -1.827 0.06820 .
## tmmn_mea_60b_2a -4.200e-02 1.320e-02 -3.181 0.00155 **
## tmmx_max_3b_2a 3.550e-02 1.384e-02 2.566 0.01058 *
## pr_mea_10b_2a -1.728e-02 1.930e-02 -0.896 0.37089
## pr_mea_3b_2a 1.015e-02 1.152e-02 0.881 0.37853
## vpd_min_10b_2a -2.597e-01 2.276e-01 -1.141 0.25426
## spei1y 1.029e-01 4.712e-02 2.184 0.02944 *
## spei14d 9.167e-02 6.604e-02 1.388 0.16570
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.8585633)
##
## Null deviance: 376.58 on 537 degrees of freedom
## Residual deviance: 328.05 on 520 degrees of freedom
## AIC: 83.51
##
## Number of Fisher Scoring iterations: 7# Check multicollinearity
vif(clump_glm)## ag_1985 wood_cov road_dens_km_km2 nearest_road_m
## 1.249362 1.812614 1.128367 1.070279
## NDVI_pre_90d NDVI_pre_365d mean_hydro_dist_m mean_slope
## 1.244923 1.554912 1.291469 1.198413
## soil_order_1 vs_max_3b_2a vs_mea_60b_2a tmmn_max_3b_2a
## 1.142051 1.558555 1.712144 1.197812vif(edge_glm)## ag_1985 wood_cov nearest_road_m NDVI_pre_90d
## 1.215773 1.872931 1.046878 1.457930
## NDVI_pre_365d mean_hydro_dist_m mean_slope mean_elevation
## 1.857524 1.349505 1.813655 2.019399
## soil_order_1 vs_mea_3b_2a vs_mea_30b_2a tmmn_max_3b_2a
## 1.160489 2.170767 2.092714 1.388618
## pr_max_3b_2a pr_max_10b_2a pr_mea_60b_2a rmin_mea_60b_2a
## 1.365825 1.677048 2.148358 1.515802
## spi1y
## 1.993432vif(area_glm)## ag_1985 wood_cov road_dens_km_km2 nearest_road_m
## 1.264064 2.135892 1.169023 1.140318
## NDVI_pre_90d NDVI_pre_365d mean_hydro_dist_m mean_elevation
## 1.456628 1.801586 1.315881 1.428774
## soil_order_1 vs_min_3b_2a tmmn_mea_60b_2a tmmx_max_3b_2a
## 1.158571 1.196786 3.083414 2.242347
## pr_mea_10b_2a pr_mea_3b_2a vpd_min_10b_2a spei1y
## 2.800097 1.632142 2.444759 1.705995
## spei14d
## 2.351243# Residual checks
par(mfrow = c(2, 2))
plot(clump_glm)
plot(edge_glm)
plot(area_glm)
par(mfrow = c(1, 1))Residual Diagnostics
shapiro.test(residuals(clump_glm))##
## Shapiro-Wilk normality test
##
## data: residuals(clump_glm)
## W = 0.79991, p-value < 2.2e-16shapiro.test(residuals(edge_glm))##
## Shapiro-Wilk normality test
##
## data: residuals(edge_glm)
## W = 0.92828, p-value = 2.315e-15shapiro.test(residuals(area_glm))##
## Shapiro-Wilk normality test
##
## data: residuals(area_glm)
## W = 0.95681, p-value = 1.878e-11# Gaussian
sim_clump <- simulateResiduals(clump_glm)
plot(sim_clump)
# Gamma
sim_edge <- simulateResiduals(edge_glm)
plot(sim_edge)## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
sim_area <- simulateResiduals(area_glm)
plot(sim_area)
plot(simulateResiduals(clump_glm))
plot(simulateResiduals(edge_glm))## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
plot(simulateResiduals(area_glm))
range(predict(clump_glm, type = "response"))## [1] 0.8005707 0.9143867range(predict(edge_glm, type = "response"))## [1] 363.1029 3284.7507range(predict(area_glm, type = "response"))## [1] 0.1486563 0.8872436r2(clump_glm)## R2: 0.083r2(edge_glm)## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.198r2(area_glm)## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.171GAM Versions
clump_gam <- gam(
clumpiness ~
ag_1985 +
wood_cov +
road_dens_km_km2 +
nearest_road_m +
NDVI_pre_90d +
NDVI_pre_365d +
mean_hydro_dist_m +
soil_order_1 +
vs_max_3b_2a +
vs_mea_60b_2a +
s(tmmn_max_3b_2a, k = 6) +
s(mean_slope, k = 5),
family = betar(link = "logit"),
data = landscape_data,
method = "REML"
)## Warning in family$saturated.ll(y, prior.weights, theta): saturated likelihood
## may be inaccurateedge_gam <- gam(
edge_dens ~
s(ag_1985) +
s(wood_cov) +
nearest_road_m +
s(NDVI_pre_90d) +
s(NDVI_pre_365d) +
s(mean_hydro_dist_m) +
s(mean_slope) +
s(mean_elevation) +
soil_order_1 +
s(vs_mea_3b_2a) +
s(vs_mea_30b_2a) +
s(tmmn_max_3b_2a) +
s(pr_max_3b_2a) +
s(pr_max_10b_2a) +
s(pr_mea_60b_2a) +
s(rmin_mea_60b_2a) +
s(spi1y),
family = Gamma(link = "log"),
data = landscape_data,
method = "REML"
)
area_gam <- gam(
avg_patch_area ~
s(ag_1985) +
s(wood_cov) +
s(road_dens_km_km2) +
nearest_road_m +
s(NDVI_pre_90d) +
s(NDVI_pre_365d) +
s(mean_hydro_dist_m) +
s(mean_elevation) +
soil_order_1 +
s(vs_min_3b_2a) +
s(tmmn_mea_60b_2a) +
s(tmmx_max_3b_2a) +
s(pr_mea_10b_2a) +
s(vpd_min_10b_2a) +
s(spei1y) +
s(spei14d),
family = Gamma(link = "log"),
data = landscape_data,
method = "REML"
)GAM Model Diagnostics
# Model summaries
summary(clump_gam)##
## Family: Beta regression(9.683)
## Link function: logit
##
## Formula:
## clumpiness ~ ag_1985 + wood_cov + road_dens_km_km2 + nearest_road_m +
## NDVI_pre_90d + NDVI_pre_365d + mean_hydro_dist_m + soil_order_1 +
## vs_max_3b_2a + vs_mea_60b_2a + s(tmmn_max_3b_2a, k = 6) +
## s(mean_slope, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.7144125 0.3709431 10.013 < 2e-16 ***
## ag_1985 -1.7453154 1.2069159 -1.446 0.14815
## wood_cov 0.0088246 0.0030050 2.937 0.00332 **
## road_dens_km_km2 -0.1078618 0.0530364 -2.034 0.04198 *
## nearest_road_m -0.0003408 0.0001075 -3.170 0.00153 **
## NDVI_pre_90d 0.8422525 0.2706934 3.111 0.00186 **
## NDVI_pre_365d -1.2889113 0.5713494 -2.256 0.02408 *
## mean_hydro_dist_m -0.0001893 0.0001049 -1.805 0.07113 .
## soil_order_1 -0.0471512 0.0121005 -3.897 9.75e-05 ***
## vs_max_3b_2a 0.0044153 0.0300483 0.147 0.88318
## vs_mea_60b_2a -0.3840900 0.0888266 -4.324 1.53e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(tmmn_max_3b_2a) 4.254 4.754 22.32 0.000263 ***
## s(mean_slope) 3.846 3.986 57.05 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = -0.413 Deviance explained = 49.7%
## -REML = -649.56 Scale est. = 1 n = 533gam.check(clump_gam)## Warning in object$family$saturated.ll(y, wts, object$family$getTheta(TRUE)):
## saturated likelihood may be inaccurate
## Warning in object$family$saturated.ll(y, wts, object$family$getTheta(TRUE)):
## saturated likelihood may be inaccurate
##
## Method: REML Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-1.677398e-05,0.0001089247]
## (score -649.5582 & scale 1).
## Hessian positive definite, eigenvalue range [1.013798,243.3512].
## Model rank = 20 / 20## Warning in object$family$saturated.ll(y, wts, object$family$getTheta(TRUE)):
## saturated likelihood may be inaccurate## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(tmmn_max_3b_2a) 5.00 4.25 0.86 0.255
## s(mean_slope) 4.00 3.85 0.80 0.015 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sim_clump <- simulateResiduals(clump_gam, n = 1000)
plot(sim_clump)
testUniformity(sim_clump)
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.22026, p-value < 2.2e-16
## alternative hypothesis: two-sidedtestDispersion(sim_clump)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.434, p-value < 2.2e-16
## alternative hypothesis: two.sidedtestZeroInflation(sim_clump)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sidedconcurvity(clump_gam, full = TRUE)## para s(tmmn_max_3b_2a) s(mean_slope)
## worst 0.991084 0.2306100 0.2253097
## observed 0.991084 0.1606834 0.1414294
## estimate 0.991084 0.1493649 0.1871834plot(clump_gam, pages = 1, shade = TRUE, rug = TRUE)
pred <- predict(clump_gam, type = "response")
range(pred)## [1] 0.6602795 0.9852276summary(edge_gam)##
## Family: Gamma
## Link function: log
##
## Formula:
## edge_dens ~ s(ag_1985) + s(wood_cov) + nearest_road_m + s(NDVI_pre_90d) +
## s(NDVI_pre_365d) + s(mean_hydro_dist_m) + s(mean_slope) +
## s(mean_elevation) + soil_order_1 + s(vs_mea_3b_2a) + s(vs_mea_30b_2a) +
## s(tmmn_max_3b_2a) + s(pr_max_3b_2a) + s(pr_max_10b_2a) +
## s(pr_mea_60b_2a) + s(rmin_mea_60b_2a) + s(spi1y)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.7298504 0.0848894 79.278 <2e-16 ***
## nearest_road_m 0.0001349 0.0001034 1.305 0.193
## soil_order_1 0.0102162 0.0100008 1.022 0.307
## ---
## 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.150 2.661 4.315 0.00752 **
## s(wood_cov) 2.338 2.950 10.626 1.66e-06 ***
## s(NDVI_pre_90d) 1.002 1.004 0.953 0.32927
## s(NDVI_pre_365d) 1.447 1.776 0.170 0.84333
## s(mean_hydro_dist_m) 2.116 2.558 1.994 0.14150
## s(mean_slope) 1.001 1.003 0.064 0.80311
## s(mean_elevation) 7.259 8.235 6.439 < 2e-16 ***
## s(vs_mea_3b_2a) 1.002 1.003 0.426 0.51516
## s(vs_mea_30b_2a) 1.001 1.002 0.023 0.88162
## s(tmmn_max_3b_2a) 1.841 2.346 2.155 0.10671
## s(pr_max_3b_2a) 1.154 1.286 0.851 0.31578
## s(pr_max_10b_2a) 4.410 5.362 1.840 0.08757 .
## s(pr_mea_60b_2a) 2.215 2.816 1.132 0.37140
## s(rmin_mea_60b_2a) 6.362 7.569 2.283 0.02076 *
## s(spi1y) 1.260 1.470 1.789 0.12646
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.194 Deviance explained = 37.9%
## -REML = 4137 Scale est. = 0.45837 n = 538summary(area_gam)##
## Family: Gamma
## Link function: log
##
## Formula:
## avg_patch_area ~ s(ag_1985) + s(wood_cov) + s(road_dens_km_km2) +
## nearest_road_m + s(NDVI_pre_90d) + s(NDVI_pre_365d) + s(mean_hydro_dist_m) +
## s(mean_elevation) + soil_order_1 + s(vs_min_3b_2a) + s(tmmn_mea_60b_2a) +
## s(tmmx_max_3b_2a) + s(pr_mea_10b_2a) + s(vpd_min_10b_2a) +
## s(spei1y) + s(spei14d)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.584e-01 1.056e-01 -6.233 9.65e-10 ***
## nearest_road_m -7.038e-05 1.452e-04 -0.485 0.6280
## soil_order_1 -3.090e-02 1.245e-02 -2.482 0.0134 *
## ---
## 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.000 1.000 0.860 0.3543
## s(wood_cov) 2.409 3.028 3.234 0.0225 *
## s(road_dens_km_km2) 2.601 3.314 2.353 0.0650 .
## s(NDVI_pre_90d) 2.068 2.640 0.631 0.4951
## s(NDVI_pre_365d) 2.269 2.885 1.114 0.3311
## s(mean_hydro_dist_m) 1.814 2.213 0.788 0.4392
## s(mean_elevation) 1.457 1.760 4.454 0.0127 *
## s(vs_min_3b_2a) 1.433 1.756 3.062 0.0760 .
## s(tmmn_mea_60b_2a) 1.000 1.000 10.068 0.0016 **
## s(tmmx_max_3b_2a) 2.863 3.668 1.616 0.1607
## s(pr_mea_10b_2a) 1.000 1.000 0.211 0.6465
## s(vpd_min_10b_2a) 2.522 3.146 2.216 0.0836 .
## s(spei1y) 5.496 6.577 2.268 0.0291 *
## s(spei14d) 1.893 2.397 0.744 0.4469
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0994 Deviance explained = 22%
## -REML = 60.451 Scale est. = 0.71422 n = 538Quasibinomial GLM
clump_quasi <- glm(
clump_form,
data = landscape_data,
family = quasibinomial(),
na.action = na.omit
)
summary(clump_quasi)##
## Call:
## glm(formula = clump_form, family = quasibinomial(), data = landscape_data,
## na.action = na.omit)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.998e-01 1.275e+00 -0.235 0.81421
## ag_1985 -1.023e+00 7.357e-01 -1.391 0.16488
## wood_cov -4.381e-03 1.789e-03 -2.449 0.01466 *
## road_dens_km_km2 4.905e-04 3.193e-02 0.015 0.98775
## nearest_road_m -8.614e-05 6.423e-05 -1.341 0.18046
## NDVI_pre_90d 1.896e-01 1.580e-01 1.201 0.23044
## NDVI_pre_365d 9.795e-02 3.305e-01 0.296 0.76704
## mean_hydro_dist_m -2.444e-05 6.765e-05 -0.361 0.71807
## mean_slope -3.164e-02 1.223e-02 -2.587 0.00995 **
## soil_order_1 -1.387e-02 6.957e-03 -1.993 0.04676 *
## vs_max_3b_2a -1.310e-02 1.781e-02 -0.735 0.46245
## vs_mea_60b_2a -1.444e-02 5.219e-02 -0.277 0.78222
## tmmn_max_3b_2a 9.819e-03 4.273e-03 2.298 0.02197 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 0.02560982)
##
## Null deviance: 14.527 on 532 degrees of freedom
## Residual deviance: 13.285 on 520 degrees of freedom
## (5 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5# 1. Check for remaining patterns
par(mfrow = c(2,2))
plot(clump_glm)
# 2. Check for influential points
plot(cooks.distance(clump_glm), type = "h")
abline(h = 4/nrow(landscape_data), col = "red")
# 3. Examine quantile-quantile plot more carefully
car::qqPlot(residuals(clump_glm, type = "deviance"))## 92 142
## 92 139# 4. Look for patterns by predictor
residualPlots(clump_glm)
## Test stat Pr(>|Test stat|)
## ag_1985 0.0064 0.9362
## wood_cov 0.0055 0.9409
## road_dens_km_km2 0.0017 0.9669
## nearest_road_m 0.0008 0.9774
## NDVI_pre_90d 0.0004 0.9838
## NDVI_pre_365d 0.0060 0.9384
## mean_hydro_dist_m 0.0006 0.9802
## mean_slope 0.0123 0.9116
## soil_order_1 0.0084 0.9268
## vs_max_3b_2a 0.0002 0.9885
## vs_mea_60b_2a 0.0001 0.9928
## tmmn_max_3b_2a 0.0026 0.9596Gam Options
Helper Functions
run_diagnostics <- function(mod, name) {
gam.check(mod)
res <- simulateResiduals(mod, n = 1000)
plot(res)
print(testDispersion(res))
print(testOutliers(res))
print(check_concurvity(mod))
message(glue("{name} R² = {round(r2(mod)$R2, 3)}"))
png(glue("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/08_Refugia/Refugia/04_Figure/{name}_diag.png"),
width = 1600, height = 1200, res = 200)
plot(res)
dev.off()
}
plot_smooths_and_preds <- function(mod, name) {
png(glue("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/08_Refugia/Refugia/04_Figure/{name}_smooths.png"),
width = 2000, height = 1500, res = 200)
plot(mod, pages = 1, scheme = 2, shade = TRUE)
dev.off()
terms <- attr(terms(mod), "term.labels")
preds <- map(terms, ~ ggpredict(mod, terms = .x))
p <- map(preds, ~ plot(.x) +
theme_classic() +
labs(title = NULL)) |>
wrap_plots() +
plot_layout(axis_titles = "collect")
ggsave(
glue("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/08_Refugia/Refugia/04_Figure/{name}_preds.png"),
p, width = 14, height = 10
)
}Clumpiness Model
I believe there were only 8 non-zero values for nearest_road_m, thus I did not include a smoothing term.
Beta regression fits the data (between 0 and 1), but we were left with large residuals at the extremes. For now, we have a quasibinomial model, but this does not seem to fit much better.
# Moves values slightly away from 0 and 1
eps <- 1e-4
landscape_data <- landscape_data %>%
mutate(clumpiness_beta = clumpiness * (1 - 2 * eps) + eps)
mod_clumpiness_qb <- gam(
clumpiness_beta ~
s(ag_1985) +
s(wood_cov) +
road_dens_km_km2 +
nearest_road_m +
s(NDVI_pre_90d) +
s(NDVI_pre_365d) +
s(mean_hydro_dist_m) +
s(mean_slope) +
soil_order_1 +
s(vs_max_3b_2a) +
s(vs_mea_60b_2a) +
s(tmmn_max_3b_2a),
data = landscape_data,
family = quasibinomial(link = "logit"),
method = "REML",
na.action = na.omit
)
plot(mod_clumpiness_qb, residuals = TRUE, pch = 16)
qq.gam(mod_clumpiness_qb)
plot_smooths_and_preds(mod_clumpiness_qb, "Clumpiness_beta")Edge Density
edge_formula <- edge_dens ~
s(ag_1985) +
s(wood_cov) +
nearest_road_m +
s(NDVI_pre_90d) +
s(NDVI_pre_365d) +
s(mean_hydro_dist_m) +
s(mean_slope) +
s(mean_elevation) +
soil_order_1 +
s(vs_mea_3b_2a) +
s(vs_mea_30b_2a) +
s(tmmn_max_3b_2a) +
s(pr_max_3b_2a) +
s(pr_max_10b_2a) +
s(pr_mea_60b_2a) +
s(rmin_mea_60b_2a) +
s(spi1y)
mod_edge <- gam(
edge_formula,
data = landscape_data,
family = Gamma(link = "log"),
method = "REML",
na.action = na.omit,
optimizer = c("outer", "bfgs")
)
run_diagnostics(mod_edge, "Edge_Density")
##
## Method: REML Optimizer: outer bfgs
## full convergence after 45 iterations.
## Gradient range [-0.002485165,0.002674089]
## (score 4136.951 & scale 0.4584247).
## Hessian positive definite, eigenvalue range [0.0006969184,460.3965].
## Model rank = 138 / 138
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(ag_1985) 9.00 2.15 0.91 0.080 .
## s(wood_cov) 9.00 2.34 1.01 0.845
## s(NDVI_pre_90d) 9.00 1.00 0.93 0.220
## s(NDVI_pre_365d) 9.00 1.46 0.90 0.035 *
## s(mean_hydro_dist_m) 9.00 2.11 0.99 0.755
## s(mean_slope) 9.00 1.00 0.95 0.430
## s(mean_elevation) 9.00 7.26 1.04 0.975
## s(vs_mea_3b_2a) 9.00 1.00 0.97 0.515
## s(vs_mea_30b_2a) 9.00 1.00 0.97 0.510
## s(tmmn_max_3b_2a) 9.00 1.84 0.97 0.485
## s(pr_max_3b_2a) 9.00 1.15 0.83 <2e-16 ***
## s(pr_max_10b_2a) 9.00 4.41 0.97 0.575
## s(pr_mea_60b_2a) 9.00 2.21 0.86 <2e-16 ***
## s(rmin_mea_60b_2a) 9.00 6.36 0.93 0.265
## s(spi1y) 9.00 1.26 0.97 0.575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.2783, p-value = 0.052
## alternative hypothesis: two.sided
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: res
## outliers at both margin(s) = 3, observations = 538, p-value = 0.0944
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.001151428 0.016208905
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.005576208
##
##
## Low Concurvity
##
## Term VIF VIF %
## s(ag_1985) 2.28 0.46
## s(wood_cov) 3.57 0.67
## s(NDVI_pre_90d) 2.28 0.50
## s(NDVI_pre_365d) 3.05 0.59
## s(mean_slope) 3.39 0.60
## s(mean_elevation) 3.76 0.64
## s(vs_mea_3b_2a) 3.47 0.61
## s(vs_mea_30b_2a) 4.06 0.64
## s(tmmn_max_3b_2a) 2.54 0.54
## s(pr_max_3b_2a) 2.44 0.48
## s(pr_mea_60b_2a) 3.74 0.64
## s(rmin_mea_60b_2a) 2.41 0.46
## s(spi1y) 3.84 0.65
##
## High Concurvity
##
## Term VIF VIF %
## Parametric 10.76 0.91
## s(mean_hydro_dist_m) 81.21 0.65
## s(pr_max_10b_2a) 81.72 0.69## Edge_Density R² = 0.194## png
## 2plot_smooths_and_preds(mod_edge, "Edge_Density")Average Patch Area
patch_formula <- avg_patch_area ~
s(ag_1985) +
s(wood_cov) +
s(road_dens_km_km2) +
nearest_road_m +
s(NDVI_pre_90d) +
s(NDVI_pre_365d) +
s(mean_hydro_dist_m) +
s(mean_elevation) +
soil_order_1 +
s(vs_min_3b_2a) +
s(tmmn_mea_60b_2a) +
s(tmmx_max_3b_2a) +
s(pr_mea_10b_2a) +
s(pr_mea_3b_2a) +
s(vpd_min_10b_2a) +
s(spei1y) +
s(spei14d)
mod_patch <- gam(
patch_formula,
data = landscape_data,
family = Gamma(link = "log"),
method = "REML",
na.action = na.omit,
optimizer = c("outer", "bfgs")
)
run_diagnostics(mod_patch, "Avg_Patch_Area")
##
## Method: REML Optimizer: outer bfgs
## full convergence after 85 iterations.
## Gradient range [-4.202996e-06,2.457007e-06]
## (score 59.37161 & scale 0.7006064).
## Hessian positive definite, eigenvalue range [2.815699e-06,341.8756].
## Model rank = 138 / 138
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(ag_1985) 9.00 1.00 0.88 0.070 .
## s(wood_cov) 9.00 1.00 0.93 0.305
## s(road_dens_km_km2) 9.00 2.55 0.95 0.620
## s(NDVI_pre_90d) 9.00 2.20 1.00 0.925
## s(NDVI_pre_365d) 9.00 2.25 0.95 0.580
## s(mean_hydro_dist_m) 9.00 1.84 0.98 0.775
## s(mean_elevation) 9.00 1.00 0.88 0.035 *
## s(vs_min_3b_2a) 9.00 1.97 0.91 0.180
## s(tmmn_mea_60b_2a) 9.00 1.00 0.91 0.285
## s(tmmx_max_3b_2a) 9.00 2.56 0.94 0.470
## s(pr_mea_10b_2a) 9.00 2.03 0.90 0.115
## s(pr_mea_3b_2a) 9.00 2.90 0.93 0.330
## s(vpd_min_10b_2a) 9.00 2.59 0.93 0.350
## s(spei1y) 9.00 6.39 0.94 0.490
## s(spei14d) 9.00 2.19 0.95 0.590
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1405, p-value = 0.328
## alternative hypothesis: two.sided
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: res
## outliers at both margin(s) = 2, observations = 538, p-value = 0.2918
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.0004505206 0.0133637398
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.003717472
##
##
## Low Concurvity
##
## Term VIF VIF %
## s(ag_1985) 2.24 0.44
## s(wood_cov) 3.54 0.67
## s(road_dens_km_km2) 2.05 0.32
## s(NDVI_pre_90d) 2.11 0.47
## s(NDVI_pre_365d) 2.91 0.57
## s(mean_hydro_dist_m) 2.98 0.53
## s(mean_elevation) 3.44 0.55
## s(vs_min_3b_2a) 1.86 0.41
## s(pr_mea_3b_2a) 2.66 0.55
## s(vpd_min_10b_2a) 4.40 0.70
## s(spei1y) 2.71 0.57
## s(spei14d) 4.28 0.66
##
## Moderate Concurvity
##
## Term VIF VIF %
## s(tmmn_mea_60b_2a) 6.34 0.77
## s(tmmx_max_3b_2a) 5.55 0.69
## s(pr_mea_10b_2a) 5.33 0.72
##
## High Concurvity
##
## Term VIF VIF %
## Parametric 10.90 0.91## Avg_Patch_Area R² = 0.12## png
## 2plot_smooths_and_preds(mod_patch, "Avg_Patch_Area")
summary(mod_patch)##
## Family: Gamma
## Link function: log
##
## Formula:
## avg_patch_area ~ s(ag_1985) + s(wood_cov) + s(road_dens_km_km2) +
## nearest_road_m + s(NDVI_pre_90d) + s(NDVI_pre_365d) + s(mean_hydro_dist_m) +
## s(mean_elevation) + soil_order_1 + s(vs_min_3b_2a) + s(tmmn_mea_60b_2a) +
## s(tmmx_max_3b_2a) + s(pr_mea_10b_2a) + s(pr_mea_3b_2a) +
## s(vpd_min_10b_2a) + s(spei1y) + s(spei14d)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.457e-01 1.051e-01 -6.146 1.62e-09 ***
## nearest_road_m -5.485e-05 1.441e-04 -0.381 0.704
## soil_order_1 -3.356e-02 1.239e-02 -2.708 0.007 **
## ---
## 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.000 1.000 0.829 0.36304
## s(wood_cov) 1.000 1.000 5.570 0.01865 *
## s(road_dens_km_km2) 2.547 3.247 2.451 0.05833 .
## s(NDVI_pre_90d) 2.203 2.811 0.757 0.43430
## s(NDVI_pre_365d) 2.251 2.863 0.893 0.42901
## s(mean_hydro_dist_m) 1.840 2.241 0.811 0.43952
## s(mean_elevation) 1.000 1.000 7.857 0.00526 **
## s(vs_min_3b_2a) 1.972 2.511 1.818 0.18123
## s(tmmn_mea_60b_2a) 1.000 1.000 10.834 0.00107 **
## s(tmmx_max_3b_2a) 2.558 3.297 1.423 0.22965
## s(pr_mea_10b_2a) 2.025 2.588 0.929 0.38185
## s(pr_mea_3b_2a) 2.903 3.581 2.275 0.08251 .
## s(vpd_min_10b_2a) 2.588 3.229 2.368 0.06822 .
## s(spei1y) 6.386 7.490 2.528 0.01625 *
## s(spei14d) 2.192 2.781 1.112 0.30245
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.12 Deviance explained = 24.3%
## -REML = 59.372 Scale est. = 0.70061 n = 538Clumpiness log-normal
mod_clum_ln <- gam(
log(clumpiness) ~
s(ag_1985) +
s(wood_cov) +
s(road_dens_km_km2) +
nearest_road_m +
s(NDVI_pre_90d) +
s(NDVI_pre_365d) +
s(mean_hydro_dist_m) +
soil_order_1 +
s(vs_max_3b_2a) +
s(vs_mea_60b_2a) +
s(tmmn_max_3b_2a) +
s(mean_slope),
data = landscape_data,
family = gaussian(),
method = "REML"
)
summary(mod_clum_ln)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(clumpiness) ~ s(ag_1985) + s(wood_cov) + s(road_dens_km_km2) +
## nearest_road_m + s(NDVI_pre_90d) + s(NDVI_pre_365d) + s(mean_hydro_dist_m) +
## soil_order_1 + s(vs_max_3b_2a) + s(vs_mea_60b_2a) + s(tmmn_max_3b_2a) +
## s(mean_slope)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.272e-01 9.433e-03 -13.487 <2e-16 ***
## nearest_road_m -1.141e-05 1.169e-05 -0.976 0.329
## soil_order_1 -1.742e-03 1.110e-03 -1.569 0.117
## ---
## 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.113 2.625 1.689 0.2214
## s(wood_cov) 1.000 1.001 4.066 0.0443 *
## s(road_dens_km_km2) 1.000 1.001 0.057 0.8122
## s(NDVI_pre_90d) 1.000 1.000 0.774 0.3792
## s(NDVI_pre_365d) 1.344 1.617 0.151 0.7554
## s(mean_hydro_dist_m) 1.000 1.000 0.011 0.9180
## s(vs_max_3b_2a) 1.000 1.000 0.629 0.4279
## s(vs_mea_60b_2a) 1.000 1.000 0.285 0.5939
## s(tmmn_max_3b_2a) 3.201 4.026 2.111 0.0804 .
## s(mean_slope) 2.100 2.683 3.172 0.0221 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0654 Deviance explained = 9.48%
## -REML = -542.04 Scale est. = 0.00599 n = 533# Model Fit
run_diagnostics(mod_clum_ln, "Clumpiness_LogNormal")
##
## Method: REML Optimizer: outer newton
## full convergence after 12 iterations.
## Gradient range [-4.450082e-05,0.0002723412]
## (score -542.0401 & scale 0.005989958).
## Hessian positive definite, eigenvalue range [5.127599e-06,260.0069].
## Model rank = 93 / 93
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(ag_1985) 9.00 2.11 1.05 0.90
## s(wood_cov) 9.00 1.00 1.00 0.46
## s(road_dens_km_km2) 9.00 1.00 0.90 0.01 **
## s(NDVI_pre_90d) 9.00 1.00 0.97 0.20
## s(NDVI_pre_365d) 9.00 1.34 1.02 0.64
## s(mean_hydro_dist_m) 9.00 1.00 0.99 0.34
## s(vs_max_3b_2a) 9.00 1.00 1.01 0.60
## s(vs_mea_60b_2a) 9.00 1.00 1.02 0.64
## s(tmmn_max_3b_2a) 9.00 3.20 1.02 0.60
## s(mean_slope) 9.00 2.10 0.97 0.19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.96645, p-value = 0.622
## alternative hypothesis: two.sided
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: res
## outliers at both margin(s) = 3, observations = 533, p-value = 0.09244
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.001162243 0.016360135
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.005628518
##
##
## Low Concurvity
##
## Term VIF VIF %
## s(ag_1985) 2.24 0.41
## s(wood_cov) 2.89 0.59
## s(road_dens_km_km2) 1.97 0.24
## s(NDVI_pre_90d) 1.69 0.36
## s(NDVI_pre_365d) 2.24 0.45
## s(mean_hydro_dist_m) 2.34 0.37
## s(vs_max_3b_2a) 1.89 0.40
## s(vs_mea_60b_2a) 2.81 0.52
## s(tmmn_max_3b_2a) 1.80 0.33
## s(mean_slope) 2.09 0.32
##
## Moderate Concurvity
##
## Term VIF VIF %
## Parametric 9.39 0.89## Clumpiness_LogNormal R² = 0.065## png
## 2Log-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.046e+00 1.684e+00 2.403 0.01661 *
## wood_cov 1.151e-02 2.251e-03 5.115 4.47e-07 ***
## nearest_road_m 2.395e-04 8.503e-05 2.817 0.00504 **
## NDVI_pre_90d -1.210e-01 2.068e-01 -0.585 0.55869
## soil_order_1 6.192e-03 8.230e-03 0.752 0.45217
## vs_mea_3b_2a -9.245e-03 4.090e-02 -0.226 0.82128
## vs_mea_30b_2a -1.420e-02 6.487e-02 -0.219 0.82679
## tmmn_max_3b_2a 7.153e-03 5.651e-03 1.266 0.20615
## pr_max_3b_2a 1.858e-03 1.462e-03 1.271 0.20442
## rmin_mea_60b_2a -5.934e-03 7.264e-03 -0.817 0.41438
## ---
## 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.053 2.545 6.806 0.000529 ***
## s(NDVI_pre_365d) 2.103 2.675 1.082 0.313253
## s(mean_hydro_dist_m) 2.508 3.036 4.503 0.003764 **
## s(mean_slope) 1.651 2.085 0.293 0.712259
## s(mean_elevation) 6.308 7.414 7.942 < 2e-16 ***
## s(pr_max_10b_2a) 2.813 3.546 3.173 0.018005 *
## s(pr_mea_60b_2a) 1.639 2.058 0.691 0.504395
## s(spi1y) 2.519 3.161 1.746 0.156460
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.287 Deviance explained = 32.8%
## -REML = 523.58 Scale est. = 0.3195 n = 538# Model Fit
run_diagnostics(mod_edge_ln, "Edge_Density_LogNormal")
##
## Method: REML Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [-2.95986e-08,2.881336e-08]
## (score 523.5755 & scale 0.3194992).
## Hessian positive definite, eigenvalue range [0.08372431,260.0378].
## Model rank = 82 / 82
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(ag_1985) 9.00 2.05 0.95 0.10 .
## s(NDVI_pre_365d) 9.00 2.10 0.96 0.14
## s(mean_hydro_dist_m) 9.00 2.51 1.01 0.59
## s(mean_slope) 9.00 1.65 0.95 0.10
## s(mean_elevation) 9.00 6.31 1.06 0.95
## s(pr_max_10b_2a) 9.00 2.81 1.01 0.57
## s(pr_mea_60b_2a) 9.00 1.64 0.89 <2e-16 ***
## s(spi1y) 9.00 2.52 0.99 0.32
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.94094, p-value = 0.334
## alternative hypothesis: two.sided
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: res
## outliers at both margin(s) = 4, observations = 538, p-value = 0.02379
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.002029379 0.018926228
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.007434944
##
##
## Low Concurvity
##
## Term VIF VIF %
## s(ag_1985) 1.80 0.26
## s(NDVI_pre_365d) 2.55 0.49
## s(mean_slope) 2.73 0.51
## s(mean_elevation) 3.00 0.57
## s(pr_mea_60b_2a) 2.89 0.54
## s(spi1y) 2.94 0.55
##
## High Concurvity
##
## Term VIF VIF %
## Parametric 5796.49 1.00
## s(mean_hydro_dist_m) 64.31 0.56
## s(pr_max_10b_2a) 64.17 0.57## Edge_Density_LogNormal R² = 0.287## png
## 2performance::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)
landscape_data_clean <- landscape_data[-outlier_idx, ]
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) 4.2328068 1.5486239 2.733 0.00649 **
## wood_cov 0.0095938 0.0021062 4.555 6.59e-06 ***
## nearest_road_m 0.0002476 0.0000791 3.130 0.00185 **
## NDVI_pre_90d -0.1397748 0.1917718 -0.729 0.46643
## soil_order_1 0.0041987 0.0076844 0.546 0.58504
## vs_mea_3b_2a -0.0074816 0.0380738 -0.197 0.84430
## vs_mea_30b_2a 0.0105044 0.0604795 0.174 0.86218
## tmmn_max_3b_2a 0.0059790 0.0052639 1.136 0.25657
## pr_max_3b_2a 0.0017911 0.0013595 1.318 0.18827
## pr_mea_60b_2a -0.0096153 0.0187928 -0.512 0.60912
## rmin_mea_60b_2a -0.0004783 0.0068544 -0.070 0.94439
## ---
## 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.165 2.677 7.639 0.000156 ***
## s(NDVI_pre_365d) 2.559 3.247 2.167 0.085175 .
## s(mean_hydro_dist_m) 2.897 3.525 6.363 0.000143 ***
## s(mean_slope) 2.154 2.759 0.859 0.374510
## s(mean_elevation) 5.975 7.089 9.048 < 2e-16 ***
## s(pr_max_10b_2a) 2.532 3.206 4.635 0.002679 **
## s(spi1y) 2.820 3.529 2.277 0.066347 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.319 Deviance explained = 35.9%
## -REML = 479.82 Scale est. = 0.27464 n = 531# Model Fit
run_diagnostics(mod_edge_clean, "Edge_Density_LogNormal")
##
## Method: REML Optimizer: outer newton
## full convergence after 6 iterations.
## Gradient range [-0.0001555174,2.595653e-05]
## (score 479.8249 & scale 0.2746411).
## Hessian positive definite, eigenvalue range [0.2768978,256.5383].
## Model rank = 74 / 74
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(ag_1985) 9.00 2.17 0.96 0.23
## s(NDVI_pre_365d) 9.00 2.56 0.97 0.26
## s(mean_hydro_dist_m) 9.00 2.90 0.99 0.38
## s(mean_slope) 9.00 2.15 0.96 0.18
## s(mean_elevation) 9.00 5.97 1.10 0.97
## s(pr_max_10b_2a) 9.00 2.53 1.00 0.45
## s(spi1y) 9.00 2.82 0.98 0.32## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9393, p-value = 0.312
## alternative hypothesis: two.sided
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: res
## outliers at both margin(s) = 1, observations = 531, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 4.767835e-05 1.044768e-02
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.001883239
##
##
## Low Concurvity
##
## Term VIF VIF %
## s(ag_1985) 1.75 0.23
## s(NDVI_pre_365d) 2.48 0.47
## s(mean_slope) 2.57 0.50
## s(mean_elevation) 2.91 0.54
## s(spi1y) 2.90 0.54
##
## High Concurvity
##
## Term VIF VIF %
## Parametric 5295.64 1.00
## s(mean_hydro_dist_m) 70.71 0.55
## s(pr_max_10b_2a) 70.40 0.54## Edge_Density_LogNormal R² = 0.319## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully## png
## 2performance::check_outliers(mod_edge_clean)## OK: No outliers detected.
## - Based on the following method and threshold: cook (1).
## - For variable: (Whole model)## Plot
plot(mod_edge_clean, 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 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 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 Backtransformed - linear relationships
# Function to create baseline data
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])
}))
}
# Function to get effect with back-transformation
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)
# Back-transform
newdat$fit <- exp(pred$fit + 0.5 * pred$se.fit^2)
newdat$upper <- exp(pred$fit + 2 * pred$se.fit)
newdat$lower <- exp(pred$fit - 2 * pred$se.fit)
data.frame(
x = x_seq,
fit = newdat$fit,
upper = newdat$upper,
lower = newdat$lower,
variable = var
)
}
# Variables to plot
focus_vars <- c("wood_cov", "nearest_road_m")
effects_list <- lapply(focus_vars, function(v) {
get_effect(mod_edge_clean, landscape_data_clean, v)
})
effects_df <- do.call(rbind, effects_list)
# Custom x-axis labels with units
x_labels <- c("wood_cov" = "Wood Cover (%)", "nearest_road_m" = "Distance to Nearest Road (m)")
# Plot
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") +
theme_classic(base_size = 14) + # base font size
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"
)
})
# Combine with labels A and B
combined_plot <- plots[[1]] + plots[[2]] + plot_layout(ncol = 2) + plot_annotation(tag_levels = "A")
combined_plot
# Save the combined plot
ggsave(
"C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/08_Refugia/Refugia/04_Figure/Edge_Density_Effects.png",
combined_plot, width = 12, height = 6
)Edge Density Backtransformed - non-linear relationships
get_smooth_effect <- function(model, data, var, n = 100) {
# Generate sequence along the smooth variable
x <- data[[var]]
x <- x[is.finite(x)]
if(length(x) < 2) return(NULL)
x_seq <- seq(min(x), max(x), length.out = n)
# Create baseline data
newdat <- make_baseline(data)
newdat <- newdat[rep(1, n), ]
# Set focal variable to sequence
newdat[[var]] <- x_seq
# Predict for smooth term only
pred <- predict(model, newdata = newdat, se.fit = TRUE, type = "link")
# Back-transform
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)
return(data.frame(
x = x_seq,
fit = fit,
upper = upper,
lower = lower,
variable = var
))
}
smooth_vars <- c("ag_1985", "mean_hydro_dist_m", "mean_elevation")
smooth_effects <- lapply(smooth_vars, function(v) get_smooth_effect(mod_edge_clean, landscape_data_clean, v))
smooth_effects <- do.call(rbind, smooth_effects)
x_labels_smooth <- c(
"ag_1985" = "Agricultural Land Use (1985, %)",
"mean_hydro_dist_m" = "Distance to Hydrologic Feature (m)",
"mean_elevation" = "Elevation (m)"
)
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") +
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"
)
})
library(patchwork)
combined_smooth_plot <- smooth_plots[[1]] + smooth_plots[[2]] + smooth_plots[[3]] +
plot_layout(ncol = 3) +
plot_annotation(tag_levels = "A") # Labels A, B, C
combined_smooth_plot
ggsave(
"C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/08_Refugia/Refugia/04_Figure/Edge_Density_Smooth_Effects.png",
combined_smooth_plot, width = 18, height = 6
)Edge Density back-transformed table
# Extract parametric table
param_table <- summary(mod_edge_clean)$p.table
# Convert to data frame
param_df <- as.data.frame(param_table)
# Back-transform estimates
param_df$exp_estimate <- exp(param_df$Estimate)
# Percent change for interpretation
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)
)
param_df## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.2328067906 1.548624e+00 2.73326967 6.493588e-03
## wood_cov 0.0095938154 2.106166e-03 4.55510851 6.589517e-06
## nearest_road_m 0.0002475648 7.909992e-05 3.12977286 1.851875e-03
## NDVI_pre_90d -0.1397748026 1.917718e-01 -0.72886017 4.664292e-01
## soil_order_1 0.0041986829 7.684370e-03 0.54639258 5.850404e-01
## vs_mea_3b_2a -0.0074816008 3.807376e-02 -0.19650283 8.442966e-01
## vs_mea_30b_2a 0.0105043590 6.047949e-02 0.17368464 8.621838e-01
## tmmn_max_3b_2a 0.0059789952 5.263930e-03 1.13584242 2.565679e-01
## pr_max_3b_2a 0.0017911418 1.359473e-03 1.31752661 1.882669e-01
## pr_mea_60b_2a -0.0096152856 1.879282e-02 -0.51164689 6.091245e-01
## rmin_mea_60b_2a -0.0004783424 6.854393e-03 -0.06978626 9.443917e-01
## exp_estimate percent_change
## (Intercept) 68.910 6791.0
## wood_cov 1.010 1.0
## nearest_road_m 1.000 0.0
## NDVI_pre_90d 0.870 -13.0
## soil_order_1 1.004 0.4
## vs_mea_3b_2a 0.993 -0.7
## vs_mea_30b_2a 1.011 1.1
## tmmn_max_3b_2a 1.006 0.6
## pr_max_3b_2a 1.002 0.2
## pr_mea_60b_2a 0.990 -1.0
## rmin_mea_60b_2a 1.000 0.0smooth_table <- summary(mod_edge_clean)$s.table
smooth_df <- as.data.frame(smooth_table)
smooth_df <- smooth_df %>%
dplyr::mutate(
edf = round(edf, 2),
F = round(F, 2),
p_value = round(`p-value`, 4)
)
smooth_df## edf Ref.df F p-value p_value
## s(ag_1985) 2.17 2.676754 7.64 0.0001562843 0.0002
## s(NDVI_pre_365d) 2.56 3.247215 2.17 0.0851751803 0.0852
## s(mean_hydro_dist_m) 2.90 3.524709 6.36 0.0001432534 0.0001
## s(mean_slope) 2.15 2.758910 0.86 0.3745100109 0.3745
## s(mean_elevation) 5.97 7.088744 9.05 0.0000000000 0.0000
## s(pr_max_10b_2a) 2.53 3.205610 4.63 0.0026789036 0.0027
## s(spi1y) 2.82 3.529157 2.28 0.0663467090 0.0663write.csv(param_df, "backtransformed_param_summary.csv", row.names = TRUE)
write.csv(smooth_df, "smooth_summary.csv", row.names = TRUE)Patch Area Log-normal
Model fit seems appropriate.
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.253e+00 3.457e+00 1.230 0.21912
## ag_1985 -1.143e+00 1.199e+00 -0.954 0.34052
## wood_cov -8.444e-03 3.120e-03 -2.706 0.00703 **
## nearest_road_m -7.284e-05 1.285e-04 -0.567 0.57117
## NDVI_pre_90d -2.895e-01 2.745e-01 -1.054 0.29215
## NDVI_pre_365d 4.057e-01 5.782e-01 0.702 0.48316
## mean_elevation -1.287e-03 5.281e-04 -2.436 0.01517 *
## soil_order_1 -1.914e-02 1.105e-02 -1.732 0.08388 .
## vs_min_3b_2a -7.859e-02 4.776e-02 -1.646 0.10045
## tmmn_mea_60b_2a -3.399e-02 1.104e-02 -3.079 0.00219 **
## tmmx_max_3b_2a 1.717e-02 1.189e-02 1.444 0.14930
## pr_mea_10b_2a 2.543e-03 1.328e-02 0.191 0.84821
## spei1y 7.349e-02 3.911e-02 1.879 0.06081 .
## spei14d 4.654e-02 5.201e-02 0.895 0.37130
## ---
## 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.645 3.374 3.238 0.0179 *
## s(mean_hydro_dist_m) 2.016 2.462 1.523 0.2187
## s(vpd_min_10b_2a) 2.382 2.974 2.546 0.0674 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0951 Deviance explained = 12.9%
## -REML = 666.42 Scale est. = 0.58951 n = 538# Model Fit
run_diagnostics(mod_pa_ln, "Avg_Patch_Area_LogNormal")
##
## Method: REML Optimizer: outer newton
## full convergence after 6 iterations.
## Gradient range [-2.527887e-07,3.435896e-10]
## (score 666.4209 & scale 0.5895059).
## Hessian positive definite, eigenvalue range [0.4351183,260.5054].
## Model rank = 41 / 41
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(road_dens_km_km2) 9.00 2.65 1.01 0.64
## s(mean_hydro_dist_m) 9.00 2.02 1.05 0.88
## s(vpd_min_10b_2a) 9.00 2.38 0.96 0.13
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.96055, p-value = 0.51
## alternative hypothesis: two.sided
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: res
## outliers at both margin(s) = 0, observations = 538, p-value = 0.6328
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.0000000 0.0068332
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0
##
##
## Low Concurvity
##
## Term VIF VIF %
## s(road_dens_km_km2) 1.56 0.16
## s(mean_hydro_dist_m) 1.56 0.17
## s(vpd_min_10b_2a) 2.80 0.53
##
## High Concurvity
##
## Term VIF VIF %
## Parametric 11671.84 1.00## Avg_Patch_Area_LogNormal R² = 0.095## png
## 2performance::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 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 Backtransformed
# Extract parametric table
param_table_pa <- summary(mod_pa_ln)$p.table
# Convert to data frame
param_df_pa <- as.data.frame(param_table_pa)
# Back-transform estimates
param_df_pa$exp_estimate <- exp(param_df_pa$Estimate)
# Percent change for interpretation
param_df_pa$percent_change <- (param_df_pa$exp_estimate - 1) * 100
# Round nicely for reporting
param_df_pa <- param_df_pa %>%
mutate(
exp_estimate = round(exp_estimate, 3),
percent_change = round(percent_change, 1)
)
param_df_pa## Estimate Std. Error t value Pr(>|t|) exp_estimate
## (Intercept) 4.253067e+00 3.4567609522 1.2303620 0.219121135 70.321
## ag_1985 -1.143421e+00 1.1985343960 -0.9540162 0.340521301 0.319
## wood_cov -8.443526e-03 0.0031202025 -2.7060826 0.007033159 0.992
## nearest_road_m -7.284403e-05 0.0001285439 -0.5666861 0.571173402 1.000
## NDVI_pre_90d -2.894978e-01 0.2745391105 -1.0544866 0.292152776 0.749
## NDVI_pre_365d 4.057233e-01 0.5781743269 0.7017318 0.483162233 1.500
## mean_elevation -1.286655e-03 0.0005281116 -2.4363309 0.015173673 0.999
## soil_order_1 -1.914329e-02 0.0110531620 -1.7319288 0.083882746 0.981
## vs_min_3b_2a -7.858946e-02 0.0477573080 -1.6456008 0.100453793 0.924
## tmmn_mea_60b_2a -3.398817e-02 0.0110396156 -3.0787454 0.002189165 0.967
## tmmx_max_3b_2a 1.716742e-02 0.0118874550 1.4441626 0.149299026 1.017
## pr_mea_10b_2a 2.543449e-03 0.0132819138 0.1914972 0.848211260 1.003
## spei1y 7.349321e-02 0.0391137640 1.8789604 0.060812285 1.076
## spei14d 4.654231e-02 0.0520125628 0.8948282 0.371295368 1.048
## percent_change
## (Intercept) 6932.1
## ag_1985 -68.1
## wood_cov -0.8
## nearest_road_m 0.0
## NDVI_pre_90d -25.1
## NDVI_pre_365d 50.0
## mean_elevation -0.1
## soil_order_1 -1.9
## vs_min_3b_2a -7.6
## tmmn_mea_60b_2a -3.3
## tmmx_max_3b_2a 1.7
## pr_mea_10b_2a 0.3
## spei1y 7.6
## spei14d 4.8# Smooth term table
smooth_table_pa <- summary(mod_pa_ln)$s.table
# Convert to data frame
smooth_df_pa <- as.data.frame(smooth_table_pa)
# Optional: round for presentation
smooth_df_pa <- smooth_df_pa %>%
mutate(
edf = round(edf, 2),
F = round(F, 2),
p_value = round(`p-value`, 4)
)
smooth_df_pa## edf Ref.df F p-value p_value
## s(road_dens_km_km2) 2.65 3.374415 3.24 0.01786407 0.0179
## s(mean_hydro_dist_m) 2.02 2.461627 1.52 0.21868769 0.2187
## s(vpd_min_10b_2a) 2.38 2.974471 2.55 0.06738564 0.0674## Save csv
write.csv(param_df_pa, "patch_area_backtransformed_param_summary.csv", row.names = TRUE)
write.csv(smooth_df_pa, "patch_area_smooth_summary.csv", row.names = TRUE)Patch Area Backtransformed- linear effects
get_param_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)
# Baseline data
newdat <- make_baseline(data)
newdat <- newdat[rep(1, n), ]
# Vary focal variable
newdat[[var]] <- x_seq
# Predict full model
pred <- predict(model, newdata = newdat, se.fit = TRUE, type = "link")
# Back-transform
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)
data.frame(
x = x_seq,
fit = fit,
upper = upper,
lower = lower,
variable = var
)
}
param_vars <- c("wood_cov", "mean_elevation", "tmmn_mea_60b_2a")
param_effects <- lapply(param_vars, function(v) get_param_effect(mod_pa_ln, landscape_data, v))
param_effects <- do.call(rbind, param_effects)
x_labels_pa <- c(
"wood_cov" = "Wood Cover (%)",
"mean_elevation" = "Elevation (m)",
"tmmn_mea_60b_2a" = "Min Temperature (K)"
)
pa_plots <- lapply(param_vars, function(var) {
df <- param_effects %>% filter(variable == var)
ggplot(df, aes(x = x, y = fit)) +
geom_line(linewidth = 1.2, color = "steelblue") +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "lightblue") +
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_pa[var],
y = "Predicted Patch Area",
tag = ""
)
})
# Combine with patchwork, labeled A, B, C
combined_pa_plot <- pa_plots[[1]] + pa_plots[[2]] + pa_plots[[3]] +
plot_layout(ncol = 3) +
plot_annotation(tag_levels = "A") # A, B, C
combined_pa_plot
ggsave(
"C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/08_Refugia/Refugia/04_Figure/Patch_Area_Effects.png",
combined_pa_plot, width = 18, height = 6
)Patch Area Backtransformed- non-linear effects
get_smooth_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)
# Baseline data
newdat <- make_baseline(data)
newdat <- newdat[rep(1, n), ]
# Vary focal smooth variable
newdat[[var]] <- x_seq
# Predict full model (type = "link" to back-transform)
pred <- predict(model, newdata = newdat, se.fit = TRUE, type = "link")
# Back-transform
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)
data.frame(
x = x_seq,
fit = fit,
upper = upper,
lower = lower,
variable = var
)
}
road_effect <- get_smooth_effect(mod_pa_ln, landscape_data, "road_dens_km_km2")
road_plot <- ggplot(road_effect, aes(x = x, y = fit)) +
geom_line(linewidth = 1.2, color = "#0072B2") +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "#0072B2") +
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 = "Road Density (km/km²)",
y = "Predicted Patch Area"
)
road_plot
ggsave(
"C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/08_Refugia/Refugia/04_Figure/Patch_Area_Road_Effect.png",
road_plot, width = 8, height = 6
)