Reworked Biodiversity
Alan Ivory
2026-07-09
pacman::p_load(randomForest, Boruta, caret, tidyverse, vip, pdp,
sf, corrplot, gridExtra, spdep, adespatial, tibble,
readxl, terra, patchwork)
EnvironmentalOutputs <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/03_Biodiversity/06_Processing/Environmental_Outputs.xlsx")
social_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/03_Biodiversity/06_Processing/Social_Outputs.xlsx")
# Merge Environmental and Social Data
EnvironmentalOutputs <- EnvironmentalOutputs %>%
left_join(social_data, by = "Plot")
# Tree PCQ Data
tree_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
sheet = "Tree_PCQ")
# Veg Data
Veg_Cover <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
sheet = "Veg_Cover")
# Shrub Cover Data
shrub_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
sheet = "Shrub_Cover")
# Site Data
CogonSites <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/CogonSites_FL_AL_MS.xlsx")Tree Canopy
canopy_plot_avg <- Veg_Cover %>%
distinct(Plot, Quadrat, Canopy_Cover) %>%
group_by(Plot) %>%
summarise(
mean_canopy_cover_plot = mean(Canopy_Cover, na.rm = TRUE),
.groups = "drop"
)
EnvironmentalOutputs <- EnvironmentalOutputs %>%
left_join(canopy_plot_avg, by = "Plot")
cor(EnvironmentalOutputs$TreeCanopy_NLCD,
EnvironmentalOutputs$mean_canopy_cover_plot,
use = "complete.obs")## [1] 0.1742495Extract data from tif
agb_raster <- rast("I:/Cogongrass_Project/data/TreeMap2022_CONUS_DRYBIO_L.tif")
points <- vect(EnvironmentalOutputs,
geom = c("Longitude", "Latitude"),
crs = "EPSG:4326")
agb_vals <- extract(agb_raster, points)## Warning: [extract] transforming vector data to the CRS of the rasterEnvironmentalOutputs$agb <- agb_vals[,2]Extract Soil Data
# Load raster
r_out <- rast("I:/Cogongrass_Project/data/Soils/taxorder_CONUS_30m.tif")
# Reproject points to match raster CRS
points_proj <- project(points, crs(r_out))
# Soil value extract
soil_vals <- extract(r_out, points_proj)
# Add soil type to EnvironmentalOutputs
EnvironmentalOutputs$SoilType <- soil_vals$taxorder
table(EnvironmentalOutputs$SoilType)##
## 1 4 9 10
## 31 31 1 143# Convert numeric IDs to labeled factor
soil_lookup <- data.frame(
ID = 1:11,
SoilOrder = c("Alfisols","Andisols","Aridisols","Entisols","Gelisols",
"Histosols","Inceptisols","Mollisols","Spodosols","Ultisols","Vertisols")
)
soil_vals$SoilType <- soil_lookup$SoilOrder[soil_vals$taxorder]
EnvironmentalOutputs$SoilType <- as.factor(soil_vals$SoilType)NLCD
nlcd <- rast("I:/Cogongrass_Project/data/nlcd/Annual_NLCD_LndCov_2024_CU_C1V1.tif")
# Lookup table
nlcd_table <- data.frame(
value = c(11,12,21,22,23,24,31,41,42,43,51,52,71,72,73,74,81,82,90,95),
class = c(
"Open Water", "Perennial Ice/Snow",
"Developed Open", "Developed Low", "Developed Medium", "Developed High",
"Barren",
"Deciduous Forest", "Evergreen Forest", "Mixed Forest",
"Dwarf Scrub", "Shrub/Scrub",
"Grassland", "Sedge", "Lichens", "Moss",
"Pasture/Hay", "Cultivated Crops",
"Woody Wetlands", "Emergent Wetlands"
)
)
levels(nlcd) <- nlcd_table
plot(nlcd)
points(points_proj, col = "red", pch = 16)
nlcd_vals <- extract(nlcd, points_proj, bind = TRUE)## Warning: [extract] transforming vector data to the CRS of the rasterstr(nlcd_vals)## S4 class 'SpatVector' [package "terra"]nrow(nlcd_vals)## [1] 206nlcd_df <- as.data.frame(nlcd_vals)
str(nlcd_df)## 'data.frame': 206 obs. of 24 variables:
## $ Plot : chr "BI200" "BI201" "BI202" "BI97" ...
## $ Site : chr "BRSF" "BRSF" "BRSF" "BRSF" ...
## $ BurnYear : num 2022 2022 2022 2022 2022 ...
## $ Camera : chr NA "S-4K" NA NA ...
## $ Muname : chr "TroupLoamySand_0_to_5_percent_slopes" "TroupLoamySand_0_to_5_percent_slopes" "TroupLoamySand_0_to_5_percent_slopes" "FuquayLoamySand_0_to_5_percent_slopes" ...
## $ NLCD_LandCover : chr "Evergreen" "Evergreen" "Evergreen" "Evergreen" ...
## $ Region : chr "FA" "FA" "FA" "FA" ...
## $ Status : chr "Invaded" "Invaded" "Invaded" "Invaded" ...
## $ SoilType : chr "Ultisols" "Ultisols" "Ultisols" "Ultisols" ...
## $ Soil_Group : num 381 381 381 389 381 381 381 381 381 389 ...
## $ TreeCanopy_NLCD : num 66 94 70 59 83 63 66 71 82 54 ...
## $ aspect : num 0 246.8 180 246.8 13.1 ...
## $ elevation : num 72 70 71 48 58 71 65 55 61 43 ...
## $ ndvi : num 0.46 0.475 0.461 0.468 0.466 ...
## $ pr : num 4.76 4.76 4.81 4.8 4.82 ...
## $ slope : num 0.927 2.349 1.854 2.351 4.751 ...
## $ sph : num 0.0109 0.0109 0.0108 0.0107 0.0107 ...
## $ srad : num 206 206 207 207 207 ...
## $ tmmn : num 287 287 287 287 287 ...
## $ agb : num 22.49 22.49 13.73 9.6 9.23 ...
## $ avg_pop_den_1km : num 0.00479 0.00305 0.57717 2.00298 1.05855 ...
## $ road_density_1km : num 6.07 6.44 6.29 2.34 5.35 ...
## $ mean_canopy_cover_plot: num 15.7 37.1 32.9 37.1 27.1 ...
## $ class : Factor w/ 20 levels "Open Water","Perennial Ice/Snow",..: 9 9 9 9 9 9 9 9 9 3 ...nlcd_column <- nlcd_df[, ncol(nlcd_df)]
nlcd_vals <- extract(nlcd, points_proj)## Warning: [extract] transforming vector data to the CRS of the rastercomparison <- cbind(EnvironmentalOutputs,
NLCD_extracted = nlcd_vals[,2])
comparison <- comparison %>%
dplyr::select(Plot, NLCD_LandCover, NLCD_extracted)
comparison$NLCD_LandCover <- as.character(comparison$NLCD_extracted)Number of quadrats sampled per plot
# Count the total number of quadrats per plot
quadrat_count <- Veg_Cover %>%
group_by(Plot) %>%
summarize(total_quadrats = n_distinct(Quadrat), .groups = "drop")Filter All data to only include specified species (Per PLANTS database)
#Filter tree data to only include trees with "tree" in the growth column
tree_data <- dplyr::filter(tree_data, Growth == "Tree")
#Filter Veg Cover to exclude Shrubs and Trees
Veg_Cover <- dplyr::filter(Veg_Cover, Growth != "Shrub" & Growth != "Tree")
#Filter Shrub Cover to only include Shrubs and Trees
shrub_data <- dplyr::filter(shrub_data, Growth == "Shrub" | Growth == "Tree")Shrub Cover Conversion
# Total length of Shrub cover at a site
shrub_cover <- shrub_data %>%
mutate(Cover = Line_End - Line_Start) %>%
group_by(Species_Name, Plot, InvStatus) %>%
summarise(Total_Cover = sum(Cover, na.rm = TRUE), .groups = "drop") %>%
mutate(Percent_Cover = Total_Cover / 3000 * 100)Herbacous Cover Conversion
# Combine Plot and Quadrat columns
Veg_Cover <- Veg_Cover %>%
mutate(Plot_Quadrat = paste(Plot, Quadrat, sep = '_'))
# Join with CogonSites to get site information
Veg_Cover <- Veg_Cover %>%
left_join(CogonSites, by = "Plot")
# Sum species cover across quadrats for each species at each plot
veg_cover_summed <- Veg_Cover %>%
group_by(Plot, Species_Name, InvStatus) %>%
summarize(total_cover = sum(Cover_Per, na.rm = TRUE), .groups = "drop")
# Calculate average herbaceous species cover
avg_species_cover <- veg_cover_summed %>%
left_join(quadrat_count, by = "Plot") %>%
mutate(avg_cover = total_cover / total_quadrats)Merging Herb cover with Shrub
This species matrix includes herbaceous and shrub species
# Merge shrub cover with herbaceous average cover
combined_cover <- avg_species_cover %>%
full_join(
shrub_cover %>%
dplyr::select(Plot, Species_Name, Percent_Cover, InvStatus),
by = c("Plot", "Species_Name")
) %>%
mutate(
overlap_flag = ifelse(!is.na(avg_cover) & !is.na(Percent_Cover), TRUE, FALSE),
final_cover = case_when(
!is.na(avg_cover) & is.na(Percent_Cover) ~ avg_cover, # Use herbaceous cover if no shrub data
is.na(avg_cover) & !is.na(Percent_Cover) ~ Percent_Cover, # Use shrub cover if no herbaceous data
TRUE ~ NA_real_ # Leave as NA where overlaps exist
)
)Extract Cogongrass Cover
# Extract cogongrass cover
cogongrass_cover <- combined_cover %>%
filter(Species_Name == "Imperata_cylindrica") %>%
dplyr::select(Plot, final_cover) %>%
rename(Cogongrass_Cover = final_cover)Species Matrix
combined_cover <- combined_cover %>%
filter(Species_Name != "Imperata_cylindrica" , Species_Name != "Imperata_cylindrica_Live") # Remove cogongrass from species matrix
## Remove any non_native species
combined_cover <- combined_cover %>%
filter(InvStatus.x != "Non_Native") # Remove non-native species from species matrix
species_matrix <- combined_cover %>%
dplyr::select(Plot, Species_Name, final_cover) %>%
pivot_wider(
names_from = Species_Name,
values_from = final_cover,
values_fill = 0
)Shannon Diversity
# Calculate Shannon diversity index for each site
shannon_diversity <- species_matrix %>%
dplyr::select(-Plot) %>%
vegan::diversity(index = "shannon") %>%
as.data.frame() %>%
setNames("Shannon_Diversity") %>%
mutate(Plot = species_matrix$Plot)Merge Shannon Diversity with Cogongrass Cover and Environmental Data
# Merge Shannon diversity with cogongrass cover, tree canopy cover, and environmental outputs
model_data <- shannon_diversity %>%
left_join(cogongrass_cover, by = "Plot") %>%
left_join(EnvironmentalOutputs, by = "Plot") %>%
mutate(Cogongrass_Cover = ifelse(is.na(Cogongrass_Cover), 0, Cogongrass_Cover))
summary(model_data)## Shannon_Diversity Plot Cogongrass_Cover Site
## Min. :0.000 Length:206 Min. : 0.00 Length:206
## 1st Qu.:2.065 Class :character 1st Qu.: 0.00 Class :character
## Median :2.407 Mode :character Median : 0.00 Mode :character
## Mean :2.333 Mean :17.08
## 3rd Qu.:2.720 3rd Qu.:29.82
## Max. :3.552 Max. :92.14
##
## BurnYear Camera Latitude Longitude
## Min. :2000 Length:206 Min. :28.48 Min. :-89.80
## 1st Qu.:2014 Class :character 1st Qu.:29.56 1st Qu.:-88.91
## Median :2022 Mode :character Median :30.77 Median :-87.14
## Mean :2018 Mean :30.38 Mean :-86.72
## 3rd Qu.:2022 3rd Qu.:30.98 3rd Qu.:-83.58
## Max. :2024 Max. :31.59 Max. :-81.94
## NA's :42
## Muname NLCD_LandCover Region Status
## Length:206 Length:206 Length:206 Length:206
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## SoilType Soil_Group TreeCanopy_NLCD aspect elevation
## Alfisols : 31 Min. : 18 Min. : 0.00 Min. : 0.0 Min. : 17.00
## Entisols : 31 1st Qu.:377 1st Qu.:51.50 1st Qu.: 90.0 1st Qu.: 38.00
## Spodosols: 1 Median :381 Median :75.50 Median :180.0 Median : 59.00
## Ultisols :143 Mean :330 Mean :66.22 Mean :171.7 Mean : 57.68
## 3rd Qu.:389 3rd Qu.:86.75 3rd Qu.:270.0 3rd Qu.: 75.00
## Max. :389 Max. :98.00 Max. :350.6 Max. :111.00
##
## ndvi pr slope sph
## Min. :0.2728 Min. :3.850 Min. : 0.000 Min. :0.01051
## 1st Qu.:0.4303 1st Qu.:4.345 1st Qu.: 2.149 1st Qu.:0.01074
## Median :0.4624 Median :4.592 Median : 2.984 Median :0.01091
## Mean :0.4503 Mean :4.566 Mean : 3.601 Mean :0.01145
## 3rd Qu.:0.4846 3rd Qu.:4.796 3rd Qu.: 4.688 3rd Qu.:0.01183
## Max. :0.5357 Max. :5.294 Max. :18.429 Max. :0.01325
##
## srad tmmn agb avg_pop_den_1km
## Min. :204.7 Min. :285.7 Min. : 6.90 Min. : 0.0000
## 1st Qu.:206.8 1st Qu.:286.6 1st Qu.: 21.21 1st Qu.: 0.5956
## Median :207.0 Median :287.2 Median : 35.82 Median : 2.1841
## Mean :210.6 Mean :287.4 Mean : 38.51 Mean : 8.3667
## 3rd Qu.:216.6 3rd Qu.:288.6 3rd Qu.: 50.89 3rd Qu.: 5.2578
## Max. :224.2 Max. :290.3 Max. :112.06 Max. :122.9795
## NA's :30
## road_density_1km mean_canopy_cover_plot
## Min. :0.000 Min. : 0.00
## 1st Qu.:1.771 1st Qu.:20.00
## Median :3.103 Median :32.86
## Mean :3.037 Mean :34.16
## 3rd Qu.:3.810 3rd Qu.:48.57
## Max. :6.932 Max. :85.71
## Check Distribution of Response Variable
hist(model_data$Shannon_Diversity, breaks=30, main="Histogram of Shannon Diversity", xlab="Shannon Diversity")
summary(model_data$Shannon_Diversity)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2.065 2.407 2.333 2.720 3.552Model Data Preparation
set.seed(97)
# convert categorical variables to numeric
model_data$NLCD_LandCover <- as.numeric(as.factor(model_data$NLCD_LandCover))
model_data$SoilType <- as.numeric(as.factor(model_data$SoilType))
model_data$Site <- as.numeric(as.factor(model_data$Site))
str(model_data)## 'data.frame': 206 obs. of 27 variables:
## $ Shannon_Diversity : num 2.69 2.65 2.54 1.61 2.97 ...
## $ Plot : chr "BI200" "BI201" "BI202" "BI97" ...
## $ Cogongrass_Cover : num 22.1 10.7 39.6 21.1 2.5 ...
## $ Site : num 2 2 2 2 2 2 2 2 2 2 ...
## $ BurnYear : num 2022 2022 2022 2022 2022 ...
## $ Camera : chr "NA" "S-4K" "NA" "NA" ...
## $ Latitude : num 30.8 30.8 30.8 30.9 30.8 ...
## $ Longitude : num -86.9 -86.9 -86.8 -86.9 -86.9 ...
## $ Muname : chr "TroupLoamySand_0_to_5_percent_slopes" "TroupLoamySand_0_to_5_percent_slopes" "TroupLoamySand_0_to_5_percent_slopes" "FuquayLoamySand_0_to_5_percent_slopes" ...
## $ NLCD_LandCover : num 4 4 4 4 4 4 4 4 4 4 ...
## $ Region : chr "FA" "FA" "FA" "FA" ...
## $ Status : chr "Invaded" "Invaded" "Invaded" "Invaded" ...
## $ SoilType : num 4 4 4 4 4 4 4 4 4 4 ...
## $ Soil_Group : num 381 381 381 389 381 381 381 381 381 389 ...
## $ TreeCanopy_NLCD : num 66 94 70 59 83 63 66 71 82 54 ...
## $ aspect : num 0 246.8 180 246.8 13.1 ...
## $ elevation : num 72 70 71 48 58 71 65 55 61 43 ...
## $ ndvi : num 0.46 0.475 0.461 0.468 0.466 ...
## $ pr : num 4.76 4.76 4.81 4.8 4.82 ...
## $ slope : num 0.927 2.349 1.854 2.351 4.751 ...
## $ sph : num 0.0109 0.0109 0.0108 0.0107 0.0107 ...
## $ srad : num 206 206 207 207 207 ...
## $ tmmn : num 287 287 287 287 287 ...
## $ agb : num 45 45 27.5 19.2 18.5 ...
## $ avg_pop_den_1km : num 0.00479 0.00305 0.57717 2.00298 1.05855 ...
## $ road_density_1km : num 6.07 6.44 6.29 2.34 5.35 ...
## $ mean_canopy_cover_plot: num 15.7 37.1 32.9 37.1 27.1 ...env_predictors <- model_data[, c("Cogongrass_Cover", "elevation", "aspect",
"slope", "pr", "sph", "srad", "tmmn",
"TreeCanopy_NLCD", "NLCD_LandCover",
"SoilType", "ndvi", "avg_pop_den_1km",
"road_density_1km", "agb")]
all_predictors <- cbind(env_predictors)
model_df <- cbind(Shannon_Diversity = model_data$Shannon_Diversity,
all_predictors) %>%
na.omit()
cat(sprintf(" Environmental predictors : %d\n", ncol(env_predictors)))## Environmental predictors : 15cat(sprintf(" Total predictors : %d\n", ncol(all_predictors)))## Total predictors : 15cat(sprintf(" Observations (post-NA) : %d\n\n", nrow(model_df)))## Observations (post-NA) : 176Correlation Check
cor_matrix <- cor(model_df, use = "complete.obs")
png("correlation_matrix.png", width = 900, height = 800, res = 120)
corrplot(cor_matrix, method = "color", type = "upper", tl.cex = 0.75,
addCoef.col = "black", number.cex = 0.55,
title = "Predictor Correlation Matrix", mar = c(0,0,1,0))
dev.off()## png
## 2cat("Correlation matrix saved.\n")## Correlation matrix saved.# Compute correlation matrix for predictors only (exclude response)
pred_cor <- cor(model_df[, -1], use = "complete.obs") # remove Shannon_Diversity
# Find pairs with |r| > 0.9
high_cor_pairs <- findCorrelation(pred_cor, cutoff = 0.9, names = TRUE, verbose = TRUE)## Compare row 8 and column 6 with corr 0.913
## Means: 0.382 vs 0.196 so flagging column 8
## All correlations <= 0.9if ("pr" %in% high_cor_pairs) {
# Get all variables highly correlated with 'pr'
pr_cor_vars <- names(which(abs(pred_cor["pr", ]) > 0.9))
pr_cor_vars <- setdiff(pr_cor_vars, "pr")
# Remove these variables instead of pr
high_cor_pairs <- setdiff(high_cor_pairs, "pr")
high_cor_pairs <- unique(c(high_cor_pairs, pr_cor_vars))
}
# Remove highly correlated predictors
if (length(high_cor_pairs) > 0) {
cat(sprintf(" Removing %d highly correlated predictors (|r| > 0.9) while keeping 'pr': %s\n",
length(high_cor_pairs),
paste(high_cor_pairs, collapse = ", ")))
model_df <- model_df[, !(colnames(model_df) %in% high_cor_pairs)]
} else {
cat(" No highly correlated predictors (|r| > 0.9) found — no removals necessary.\n")
}## Removing 1 highly correlated predictors (|r| > 0.9) while keeping 'pr': tmmnVariable Selection- Boruta
cat("========== STEP 3: Boruta Variable Selection ==========\n")## ========== STEP 3: Boruta Variable Selection ==========cat(" Running Boruta (~ 1 minute)...\n")## Running Boruta (~ 1 minute)...boruta_out <- Boruta(Shannon_Diversity ~ ., data = model_df,
doTrace = 1, maxRuns = 150, num.trees = 500)## After 11 iterations, +0.99 secs:## confirmed 3 attributes: avg_pop_den_1km, Cogongrass_Cover, pr;## rejected 3 attributes: agb, NLCD_LandCover, slope;## still have 8 attributes left.## After 15 iterations, +1.2 secs:## confirmed 1 attribute: srad;## rejected 1 attribute: aspect;## still have 6 attributes left.## After 18 iterations, +1.4 secs:## confirmed 1 attribute: TreeCanopy_NLCD;## still have 5 attributes left.## After 36 iterations, +2.6 secs:## rejected 1 attribute: SoilType;## still have 4 attributes left.## After 49 iterations, +3.4 secs:## confirmed 1 attribute: sph;## still have 3 attributes left.## After 52 iterations, +3.6 secs:## rejected 1 attribute: road_density_1km;## still have 2 attributes left.## After 104 iterations, +6.6 secs:## confirmed 2 attributes: elevation, ndvi;## no more attributes left.boruta_final <- TentativeRoughFix(boruta_out)## Warning in TentativeRoughFix(boruta_out): There are no Tentative attributes!
## Returning original object.selected_vars <- getSelectedAttributes(boruta_final, withTentative = FALSE)
# Selected environmental vars
selected_env <- selected_vars
cat(sprintf("\n Variables confirmed : %d\n", length(selected_vars)))##
## Variables confirmed : 8cat(sprintf(" Selected predictors : %d (%s)\n",
length(selected_env), paste(selected_env, collapse = ", ")))## Selected predictors : 8 (Cogongrass_Cover, elevation, pr, sph, srad, TreeCanopy_NLCD, ndvi, avg_pop_den_1km)# Boruta importance plot
png("boruta_importance.png", width = 1000, height = 600, res = 120)
par(mar = c(11, 4, 3, 1))
plot(boruta_final, las = 2, xlab = "",
main = "")
dev.off()## png
## 2Train/Test Split
model_df_sel <- model_df[, c("Shannon_Diversity", selected_vars)]
train_idx <- createDataPartition(model_df_sel$Shannon_Diversity, p = 0.8, list = FALSE)
train_data <- model_df_sel[ train_idx, ]
test_data <- model_df_sel[-train_idx, ]Tune Random Forest via Cross-Validation
cat("========== STEP 4: Tuning Random Forest (10-fold CV) ==========\n")## ========== STEP 4: Tuning Random Forest (10-fold CV) ==========ctrl <- trainControl(method = "cv", number = 10, verboseIter = FALSE)
tune_grid <- expand.grid(mtry = seq(2, max(2, floor(sqrt(length(selected_vars))) + 3), by = 1))
rf_tuned <- train(Shannon_Diversity ~ .,
data = train_data,
method = "rf",
trControl = ctrl,
tuneGrid = tune_grid,
ntree = 500,
importance = TRUE)
best_mtry <- rf_tuned$bestTune$mtry
cat(sprintf(" Best mtry: %d\n\n", best_mtry))## Best mtry: 2rf_final <- randomForest(Shannon_Diversity ~ .,
data = train_data,
ntree = 1000,
mtry = best_mtry,
importance = TRUE)
print(rf_final)##
## Call:
## randomForest(formula = Shannon_Diversity ~ ., data = train_data, ntree = 1000, mtry = best_mtry, importance = TRUE)
## Type of random forest: regression
## Number of trees: 1000
## No. of variables tried at each split: 2
##
## Mean of squared residuals: 0.2463338
## % Var explained: 18.56Evaluate on Test Set
cat("========== STEP 5: Model Evaluation ==========\n")## ========== STEP 5: Model Evaluation ==========pred_test <- predict(rf_final, newdata = test_data)
bias_val <- mean(pred_test - test_data$Shannon_Diversity)
rel_bias <- bias_val / mean(test_data$Shannon_Diversity) * 100
rmse_val <- sqrt(mean((pred_test - test_data$Shannon_Diversity)^2))
mae_val <- mean(abs(pred_test - test_data$Shannon_Diversity))
ss_res <- sum((pred_test - test_data$Shannon_Diversity)^2)
ss_tot <- sum((test_data$Shannon_Diversity - mean(test_data$Shannon_Diversity))^2)
r2_val <- 1 - ss_res / ss_tot
cat(sprintf(" RMSE : %.4f\n", rmse_val))## RMSE : 0.4418cat(sprintf(" MAE : %.4f\n", mae_val))## MAE : 0.3510cat(sprintf(" Bias : %.4f\n", bias_val))## Bias : 0.0290cat(sprintf(" Relative Bias : %.2f%%\n", rel_bias))## Relative Bias : 1.27%cat(sprintf(" R² : %.4f\n\n", r2_val))## R² : 0.2661# Predicted vs observed plot
pred_obs_df <- data.frame(Observed = test_data$Shannon_Diversity,
Predicted = pred_test)
p_pred <- ggplot(pred_obs_df, aes(x = Observed, y = Predicted)) +
geom_point(alpha = 0.7, colour = "#2E8B57", size = 2.5) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "grey30") +
annotate("text", x = min(pred_obs_df$Observed) + 0.05,
y = max(pred_obs_df$Predicted) - 0.1,
label = sprintf("R² = %.3f\nRMSE = %.3f", r2_val, rmse_val),
hjust = 0, size = 4.5) +
labs(title = "",
x = "Observed", y = "Predicted") +
theme_bw(base_size = 13)
ggsave("predicted_vs_observed.png", p_pred, width = 6, height = 5, dpi = 150)Variable Importance Plot
imp_df <- as.data.frame(importance(rf_final)) %>%
rownames_to_column("Variable") %>%
arrange(desc(`%IncMSE`))
p_imp <- ggplot(imp_df, aes(x = reorder(Variable, `%IncMSE`),
y = `%IncMSE`)) +
geom_col(fill = "#2E8B57", alpha = 0.85) +
coord_flip() +
labs(title = "",
x = NULL,
y = "% Increase in MSE") +
theme_bw(base_size = 13)
ggsave("variable_importance.png", p_imp, width = 8, height = 5, dpi = 150)
cat("Variable importance plot saved.\n")## Variable importance plot saved.Partial Dependence Plots
top4_env <- imp_df %>%
slice_head(n = 4) %>%
pull(Variable)
pdp_plots <- lapply(top4_env, function(var) {
pd <- partial(rf_final, pred.var = var, train = train_data)
autoplot(pd) +
labs(title = var, x = var, y = "Shannon Diversity") +
theme_bw(base_size = 11)
})
png("partial_dependence_spatial_top4.png", width = 1100, height = 900, res = 130)
do.call(grid.arrange, c(pdp_plots, ncol = 2,
top = "Partial Dependence — Top 4 Environmental Predictors"))
dev.off()## png
## 2cat("Partial dependence plots saved.\n\n")## Partial dependence plots saved.Partial Dependence Plots
gator_blue <- "#0021A5"
gator_orange <- "#0021A5"
selected_env <- c(
"pr", "srad", "TreeCanopy_NLCD", "Cogongrass_Cover",
"elevation", "avg_pop_den_1km", "ndvi", "sph"
)
n_boot <- 50
var_labels <- list(
pr = "Precipitation (mm)",
srad = expression(paste("Surface Radiation (W/", m^2, ")")),
TreeCanopy_NLCD = "Tree Canopy (%)",
Cogongrass_Cover = "Cogongrass Cover (%)",
elevation = "Elevation (m)",
avg_pop_den_1km = expression(paste("Population Density (people/", km^2, ")")),
ndvi = "NDVI",
sph = "Specific Humidity (g/kg)"
)
#-----------------------------
# Generate ALL PDP summaries
#-----------------------------
pdp_data_list <- lapply(seq_along(selected_env), function(i) {
var <- selected_env[i]
boot_res <- lapply(1:n_boot, function(j) {
idx <- sample(1:nrow(train_data), replace = TRUE)
train_boot <- train_data[idx, ]
rf_boot <- randomForest(
Shannon_Diversity ~ .,
data = train_boot,
ntree = 500
)
pd <- partial(rf_boot, pred.var = var, train = train_boot)
pd$iter <- j
pd
})
pd_all <- bind_rows(boot_res)
pd_summary <- pd_all %>%
group_by(!!sym(var)) %>%
summarise(
y = mean(yhat),
ymin = quantile(yhat, 0.025),
ymax = quantile(yhat, 0.975),
.groups = "drop"
) %>%
arrange(!!sym(var)) %>%
mutate(
y_smooth = loess(y ~ get(var), span = 0.6)$fitted,
ymin_smooth = loess(ymin ~ get(var), span = 0.6)$fitted,
ymax_smooth = loess(ymax ~ get(var), span = 0.6)$fitted,
Variable = var
)
return(pd_summary)
})
#-----------------------------
# GLOBAL y-axis limits
#-----------------------------
all_pdp <- bind_rows(pdp_data_list)
y_min <- min(all_pdp$ymin_smooth, na.rm = TRUE)
y_max <- max(all_pdp$ymax_smooth, na.rm = TRUE)
y_pad <- 0.05 * (y_max - y_min)
y_limits <- c(y_min - y_pad, y_max + y_pad)
#-----------------------------
# Build plots
#-----------------------------
panel_labels <- LETTERS[1:length(pdp_data_list)]
pdp_boot_list <- lapply(seq_along(pdp_data_list), function(i) {
pd_summary <- pdp_data_list[[i]]
var <- selected_env[i]
ggplot(pd_summary, aes(x = .data[[var]])) +
geom_ribbon(
aes(ymin = ymin_smooth, ymax = ymax_smooth),
fill = gator_orange,
alpha = 0.25
) +
geom_line(
aes(y = y_smooth),
color = gator_blue,
linewidth = 1.2
) +
coord_cartesian(ylim = y_limits) +
labs(
x = var_labels[[var]],
y = NULL
) +
annotate("text",
x = Inf, y = Inf,
label = panel_labels[i],
hjust = 1.2, vjust = 1.5,
fontface = "bold", size = 5) +
theme_classic(base_size = 12) +
theme(
axis.title.x = element_text(size = 11),
axis.text = element_text(size = 10),
plot.margin = margin(5, 5, 5, 5)
)
})
#-----------------------------
# Save figure
#-----------------------------
png("partial_dependence_bootstrap.png",
width = 1300, height = 900, res = 150)
gridExtra::grid.arrange(
grid::textGrob("Shannon Diversity", rot = 90,
gp = grid::gpar(fontsize = 14, fontface = "bold")),
gridExtra::arrangeGrob(
grobs = pdp_boot_list,
ncol = 2
),
ncol = 2,
widths = c(0.05, 0.95)
)
dev.off()## png
## 2Summary Output
cat("============================================================\n")## ============================================================cat(" RANDOM FOREST MODEL SUMMARY\n")## RANDOM FOREST MODEL SUMMARYcat("============================================================\n")## ============================================================cat(sprintf(" Selected predictors : %d\n", length(selected_vars)))## Selected predictors : 8cat("------------------------------------------------------------\n")## ------------------------------------------------------------cat(sprintf(" Test RMSE : %.4f\n", rmse_val))## Test RMSE : 0.4418cat(sprintf(" Test MAE : %.4f\n", mae_val))## Test MAE : 0.3510cat(sprintf(" Test Bias : %.4f\n", bias_val))## Test Bias : 0.0290cat(sprintf(" Relative Bias : %.2f%%\n", rel_bias))## Relative Bias : 1.27%cat(sprintf(" Test R² : %.4f\n", r2_val))## Test R² : 0.2661cat("============================================================\n")## ============================================================saveRDS(rf_final, "shannon_diversity_rf_model.rds")
cat("\nModel saved as 'shannon_diversity_rf_model.rds'\n")##
## Model saved as 'shannon_diversity_rf_model.rds'