Reworked Biodiversity
Alan Ivory
2026-08-11
pacman::p_load(randomForest, Boruta, caret, tidyverse, vip, pdp,
sf, corrplot, gridExtra, scales, grid, 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 25 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" ...
## $ Field_LandCover : chr "Evergreen" "Mixed" "Evergreen" "Mixed" ...
## $ 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 Field_LandCover Region
## Length:206 Length:206 Length:206 Length:206
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## Status SoilType Soil_Group TreeCanopy_NLCD
## Length:206 Alfisols : 31 Min. : 18 Min. : 0.00
## Class :character Entisols : 31 1st Qu.:377 1st Qu.:51.50
## Mode :character Spodosols: 1 Median :381 Median :75.50
## Ultisols :143 Mean :330 Mean :66.22
## 3rd Qu.:389 3rd Qu.:86.75
## Max. :389 Max. :98.00
##
## aspect elevation ndvi pr
## Min. : 0.0 Min. : 17.00 Min. :0.2728 Min. :3.850
## 1st Qu.: 90.0 1st Qu.: 38.00 1st Qu.:0.4303 1st Qu.:4.345
## Median :180.0 Median : 59.00 Median :0.4624 Median :4.592
## Mean :171.7 Mean : 57.68 Mean :0.4503 Mean :4.566
## 3rd Qu.:270.0 3rd Qu.: 75.00 3rd Qu.:0.4846 3rd Qu.:4.796
## Max. :350.6 Max. :111.00 Max. :0.5357 Max. :5.294
##
## slope sph srad tmmn
## Min. : 0.000 Min. :0.01051 Min. :204.7 Min. :285.7
## 1st Qu.: 2.149 1st Qu.:0.01074 1st Qu.:206.8 1st Qu.:286.6
## Median : 2.984 Median :0.01091 Median :207.0 Median :287.2
## Mean : 3.601 Mean :0.01145 Mean :210.6 Mean :287.4
## 3rd Qu.: 4.688 3rd Qu.:0.01183 3rd Qu.:216.6 3rd Qu.:288.6
## Max. :18.429 Max. :0.01325 Max. :224.2 Max. :290.3
##
## agb avg_pop_den_1km road_density_1km mean_canopy_cover_plot
## Min. : 6.90 Min. : 0.0000 Min. :0.000 Min. : 0.00
## 1st Qu.: 21.21 1st Qu.: 0.5956 1st Qu.:1.771 1st Qu.:20.00
## Median : 35.82 Median : 2.1841 Median :3.103 Median :32.86
## Mean : 38.51 Mean : 8.3667 Mean :3.037 Mean :34.16
## 3rd Qu.: 50.89 3rd Qu.: 5.2578 3rd Qu.:3.810 3rd Qu.:48.57
## Max. :112.06 Max. :122.9795 Max. :6.932 Max. :85.71
## NA's :30Check 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$Field_LandCover <- as.numeric(as.factor(model_data$Field_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 28 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 ...
## $ Field_LandCover : num 2 3 2 3 2 2 2 3 2 2 ...
## $ 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")]
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 : 14cat(sprintf(" Total predictors : %d\n", ncol(all_predictors)))## Total predictors : 14cat(sprintf(" Observations (post-NA) : %d\n\n", nrow(model_df)))## Observations (post-NA) : 206Correlation 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 6 and column 8 with corr 0.916
## Means: 0.406 vs 0.208 so flagging column 6
## Compare row 8 and column 5 with corr 0.913
## Means: 0.358 vs 0.177 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 2 highly correlated predictors (|r| > 0.9) while keeping 'pr': sph, tmmnRename Variables for Readability
# Lookup table: raw name -> clear label (for display only)
var_labels <- c(
Cogongrass_Cover = "Cogongrass Cover",
elevation = "Elevation",
aspect = "Aspect",
slope = "Slope",
pr = "Precipitation",
sph = "Specific Humidity",
srad = "Surface Radiation",
tmmn = "Minimum Temperature",
TreeCanopy_NLCD = "Tree Canopy Cover",
NLCD_LandCover = "Land Cover Type",
SoilType = "Soil Order",
ndvi = "NDVI",
avg_pop_den_1km = "Population Density",
road_density_1km = "Road Density"
)
# Helper: map raw names -> clear labels, falling back to the raw name if unmapped
clean_label <- function(x) ifelse(x %in% names(var_labels), var_labels[x], x)
cat("Variable display labels defined:\n")## Variable display labels defined:print(var_labels)## Cogongrass_Cover elevation aspect
## "Cogongrass Cover" "Elevation" "Aspect"
## slope pr sph
## "Slope" "Precipitation" "Specific Humidity"
## srad tmmn TreeCanopy_NLCD
## "Surface Radiation" "Minimum Temperature" "Tree Canopy Cover"
## NLCD_LandCover SoilType ndvi
## "Land Cover Type" "Soil Order" "NDVI"
## avg_pop_den_1km road_density_1km
## "Population Density" "Road Density"Variable 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.93 secs:## confirmed 5 attributes: avg_pop_den_1km, Cogongrass_Cover, pr, srad, TreeCanopy_NLCD;## rejected 1 attribute: aspect;## still have 6 attributes left.## After 27 iterations, +2.1 secs:## rejected 1 attribute: road_density_1km;## still have 5 attributes left.## After 41 iterations, +3.2 secs:## confirmed 1 attribute: SoilType;## rejected 1 attribute: NLCD_LandCover;## still have 3 attributes left.## After 82 iterations, +6.1 secs:## confirmed 1 attribute: ndvi;## still have 2 attributes left.## After 86 iterations, +6.4 secs:## confirmed 1 attribute: elevation;## still have 1 attribute left.boruta_final <- TentativeRoughFix(boruta_out)
selected_vars <- getSelectedAttributes(boruta_final, withTentative = FALSE)
selected_env <- selected_vars
cat(sprintf("\n Variables confirmed : %d\n", length(selected_vars)))##
## Variables confirmed : 9cat(sprintf(" Selected predictors : %d (%s)\n",
length(selected_env), paste(clean_label(selected_env), collapse = ", ")))## Selected predictors : 9 (Cogongrass Cover, Elevation, Slope, Precipitation, Surface Radiation, Tree Canopy Cover, Soil Order, NDVI, Population Density)# The base-graphics Boruta plot is no longer written here. It is rebuilt in
# ggplot as panel A of the combined importance figure (RF_Variable_Importance
# chunk), so that both stages appear in one figure. Writing it here as well
# would leave a stale boruta_importance.png on disk that could be pasted into
# the manuscript by mistake.Train/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.2076294
## % Var explained: 23.58Evaluate 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.5305cat(sprintf(" MAE : %.4f\n", mae_val))## MAE : 0.3938cat(sprintf(" Bias : %.4f\n", bias_val))## Bias : 0.0106cat(sprintf(" Relative Bias : %.2f%%\n", rel_bias))## Relative Bias : 0.45%cat(sprintf(" R² : %.4f\n\n", r2_val))## R² : 0.2039# 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)Full-Data Model for Importance and Partial Dependence
## Importance and partial dependence come from a model fitted to ALL 206 plots.
## rf_final (80% training split) is kept only for the test-set metrics: across
## 100 random splits its ranking for cogongrass ranged from 1st to 9th, because
## cogongrass cover is high at relatively few plots. See Importance_Rank_Stability.
##
## ntree = 5000 with a seed: permutation importance carries Monte Carlo noise,
## and two unseeded 1000-tree fits to the same data gave surface radiation 20.4
## and 17.5 - enough to reorder predictors that sit close together.
##
## CAUTION: with scale = TRUE (the default) %IncMSE is divided by a standard
## error that depends on ntree, so the values scale roughly with sqrt(ntree).
## They are comparable among predictors here, but NOT against a model fitted
## with a different number of trees. Report ntree wherever these values appear.
## For values that do not move with ntree, use importance(rf_pdp, scale = FALSE).
set.seed(97)
rf_pdp <- randomForest(Shannon_Diversity ~ ., data = model_df_sel,
ntree = 5000, mtry = best_mtry, importance = TRUE)
pdp_train <- model_df_sel
imp_df <- as.data.frame(importance(rf_pdp)) %>%
rownames_to_column("Variable") %>%
mutate(Label = clean_label(Variable)) %>%
arrange(desc(`%IncMSE`))
print(imp_df[, c("Label", "%IncMSE")])## Label %IncMSE
## 1 Precipitation 46.53139
## 2 Surface Radiation 44.23703
## 3 Cogongrass Cover 40.88013
## 4 Population Density 31.03536
## 5 Tree Canopy Cover 26.28788
## 6 Elevation 23.15662
## 7 Soil Order 20.88005
## 8 NDVI 17.00884
## 9 Slope 13.99268## Unscaled equivalent, for reference: raw mean decrease in MSE, ntree-independent
imp_raw <- as.data.frame(importance(rf_pdp, scale = FALSE)) %>%
rownames_to_column("Variable") %>%
mutate(Label = clean_label(Variable)) %>%
arrange(desc(`%IncMSE`))
cat("\nUnscaled (raw mean decrease in MSE, independent of ntree):\n")##
## Unscaled (raw mean decrease in MSE, independent of ntree):print(imp_raw[, c("Label", "%IncMSE")])## Label %IncMSE
## 1 Precipitation 0.03702372
## 2 Surface Radiation 0.03698700
## 3 Cogongrass Cover 0.03123550
## 4 Population Density 0.02371575
## 5 Tree Canopy Cover 0.01909553
## 6 Elevation 0.01670811
## 7 NDVI 0.01277797
## 8 Slope 0.01034288
## 9 Soil Order 0.01032864Variable Importance Plot (combined Boruta + final RF)
## One figure, two stages of model construction:
## A) Boruta selection - which predictors were kept, and why
## B) Final random forest - how the kept predictors ranked in the model used
## for prediction
##
## The Boruta panel is rebuilt in ggplot rather than using plot.Boruta(), which
## is base graphics and cannot be combined with patchwork. Rebuilding also lets
## both panels run horizontally so the variable names read normally.
##
## NOTE the two panels are different metrics from different models: Boruta
## reports Z-scores of mean decrease in accuracy (ranger, 12 predictors + 12
## shadows, all 206 plots), the final forest reports %IncMSE (randomForest,
## 9 predictors, 5,000 trees). Axis labels and the caption both say so.
conf_green <- "#1B7837"
rej_red <- "#B2182B"
shadow_blue <- "#4393C3"
## ---- Panel A: Boruta ----------------------------------------------------
imp_long <- as.data.frame(boruta_final$ImpHistory) %>%
tidyr::pivot_longer(dplyr::everything(),
names_to = "Variable", values_to = "Importance") %>%
dplyr::filter(is.finite(Importance)) # -Inf once an attribute resolves
dec <- data.frame(Variable = names(boruta_final$finalDecision),
Decision = as.character(boruta_final$finalDecision),
stringsAsFactors = FALSE)
imp_long <- imp_long %>%
dplyr::mutate(
is_shadow = grepl("^shadow", Variable),
Decision = ifelse(is_shadow, "Shadow",
dec$Decision[match(Variable, dec$Variable)]),
Label = ifelse(is_shadow,
sub("^shadow", "Shadow ", Variable),
unname(clean_label(Variable)))
)
## order by median importance so the panel reads like the original Boruta plot
ord <- imp_long %>%
dplyr::group_by(Label) %>%
dplyr::summarise(med = median(Importance), .groups = "drop") %>%
dplyr::arrange(med) %>%
dplyr::pull(Label)
imp_long$Label <- factor(imp_long$Label, levels = ord)
p_boruta <- ggplot(imp_long, aes(Label, Importance, fill = Decision)) +
geom_hline(yintercept = 0, colour = "grey70", linewidth = 0.3) +
geom_boxplot(outlier.size = 0.5, linewidth = 0.3, width = 0.7) +
coord_flip() +
scale_fill_manual(values = c(Confirmed = conf_green,
Rejected = rej_red,
Shadow = shadow_blue,
Tentative = "#F4A582"),
breaks = c("Confirmed", "Rejected", "Shadow"),
name = NULL) +
labs(x = NULL, y = "Boruta importance",
title = "A) Boruta Feature Selection") +
theme_bw(base_size = 11) +
theme(panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold", size = 11),
legend.position = c(0.82, 0.16),
legend.background = element_rect(fill = scales::alpha("white", 0.8),
colour = NA),
legend.key.size = ggplot2::unit(0.4, "cm"),
legend.text = element_text(size = 8))
## ---- Panel B: final random forest ---------------------------------------
## imp_df is built in Full_Data_Model above, from rf_pdp (5,000 trees, n = 206).
p_rf <- ggplot(imp_df, aes(reorder(Label, `%IncMSE`), `%IncMSE`)) +
geom_col(fill = conf_green, alpha = 0.9, width = 0.7) +
coord_flip() +
labs(x = NULL, y = "Increase in MSE (%)",
title = "B) Final Random Forest") +
theme_bw(base_size = 11) +
theme(panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold", size = 11))
## ---- Combine -------------------------------------------------------------
## Panel B reuses the "Confirmed" green so the link is visible: every bar in B
## is one of the green boxes in A.
fig_importance <- p_boruta + p_rf + patchwork::plot_layout(widths = c(1.1, 1))
ggsave("variable_importance_combined.png", fig_importance,
width = 11, height = 5.5, units = "in", dpi = 300, bg = "white")
## individual versions retained in case either is wanted alone
ggsave("boruta_importance_gg.png", p_boruta, width = 6, height = 5.5,
units = "in", dpi = 300, bg = "white")
ggsave("variable_importance.png", p_rf, width = 6, height = 4.5,
units = "in", dpi = 300, bg = "white")
cat("Combined importance figure saved.\n")## Combined importance figure saved.print(imp_df[, c("Label", "%IncMSE")])## Label %IncMSE
## 1 Precipitation 46.53139
## 2 Surface Radiation 44.23703
## 3 Cogongrass Cover 40.88013
## 4 Population Density 31.03536
## 5 Tree Canopy Cover 26.28788
## 6 Elevation 23.15662
## 7 Soil Order 20.88005
## 8 NDVI 17.00884
## 9 Slope 13.99268Partial Dependence Plots (Top 4 Variables)
top4_env <- imp_df %>%
slice_head(n = 4) %>%
pull(Variable) # raw names, e.g. "Cogongrass_Cover"
pdp_plots <- lapply(top4_env, function(var) {
pd <- partial(rf_final, pred.var = var, train = train_data)
autoplot(pd) +
labs(title = clean_label(var), x = clean_label(var), y = "Shannon Diversity") +
theme_bw(base_size = 11)
})Partial Dependence Plots
gator_blue <- "#0021A5"
gator_orange <- "#0021A5" # was the blue hex, so ribbon and line matched
n_boot <- 50
## ---- 1. Panel order: focal variable first, then importance --------------
## Panels follow the order the Results paragraph discusses them, so the text
## reads Figure 5A, 5B, 5C... in sequence. Cogongrass leads because it is the
## subject of the paper; everything after it is ordered by %IncMSE, and each
## panel carries its importance value so the ranking is still visible.
## Importance and partial dependence are computed from a model fitted to ALL
## 206 plots. rf_final (80% training split) is kept only for the test-set
## performance metrics: across 100 random splits its importance ranking for
## cogongrass ranged from 1st to 9th, because cogongrass cover is high at
## relatively few plots. See the Importance_Rank_Stability chunk.
## rf_pdp / imp_df are built in the Full_Data_Model chunk above, so Figure 4,
## Table 1 and these panels all report the same model.
focus_var <- "Cogongrass_Cover"
selected_env <- c(focus_var, setdiff(imp_df$Variable, focus_var))
imp_lookup <- setNames(imp_df$`%IncMSE`, imp_df$Variable)
print(imp_df[, c("Variable", "%IncMSE")])## Variable %IncMSE
## 1 pr 46.53139
## 2 srad 44.23703
## 3 Cogongrass_Cover 40.88013
## 4 avg_pop_den_1km 31.03536
## 5 TreeCanopy_NLCD 26.28788
## 6 elevation 23.15662
## 7 SoilType 20.88005
## 8 ndvi 17.00884
## 9 slope 13.99268## ---- 2. Declare which predictors are actually categorical ---------------
categorical_vars <- c("SoilType")
# as.numeric(as.factor()) assigns codes alphabetically
soil_levels <- c("Alfisols", "Entisols", "Spodosols", "Ultisols")
var_labels <- list(
Cogongrass_Cover = "Cogongrass Cover (%)",
pr = "Precipitation (mm)",
srad = expression(paste("Surface Radiation (W/", m^2, ")")),
TreeCanopy_NLCD = "Tree Canopy (%)",
elevation = "Elevation (m)",
slope = "Slope (degrees)",
SoilType = "Soil Order",
avg_pop_den_1km = expression(paste("Population Density (people/", km^2, ")")),
ndvi = "NDVI",
sph = "Specific Humidity (g/kg)" # kept in case it returns
)
## ---- 3. Fixed prediction grid per variable ------------------------------
## Every bootstrap PDP must be evaluated at the SAME x values. Otherwise the
## 2.5/97.5 quantiles are taken over a single value at most grid points and
## the ribbon collapses onto the mean line.
make_grid <- function(var) {
x <- pdp_train[[var]]
g <- if (var %in% categorical_vars) {
sort(unique(x))
} else {
seq(min(x, na.rm = TRUE), max(x, na.rm = TRUE), length.out = 30)
}
setNames(data.frame(g), var)
}
pred_grids <- setNames(lapply(selected_env, make_grid), selected_env)
## ---- 4. Bootstrap the forest once, read every PDP off each fit ----------
## The original refit a forest inside the per-variable loop: 9 x 50 = 450 fits
## where 50 will do.
## Seeded separately so the bootstrap resamples, and therefore every partial
## dependence curve and span, stay reproducible even if rf_pdp above is changed.
set.seed(97)
boot_fits <- lapply(seq_len(n_boot), function(j) {
idx <- sample(nrow(pdp_train), replace = TRUE)
tb <- pdp_train[idx, ]
list(fit = randomForest(Shannon_Diversity ~ ., data = tb, ntree = 500),
dat = tb)
})
pdp_data_list <- lapply(selected_env, function(var) {
pd_all <- map_dfr(seq_along(boot_fits), function(j) {
pd <- partial(boot_fits[[j]]$fit,
pred.var = var,
pred.grid = pred_grids[[var]],
train = boot_fits[[j]]$dat)
pd$iter <- j
pd
})
pd_summary <- pd_all %>%
group_by(x = .data[[var]]) %>%
summarise(y = mean(yhat),
ymin = quantile(yhat, 0.025),
ymax = quantile(yhat, 0.975),
.groups = "drop") %>%
arrange(x) %>%
mutate(Variable = var)
## smooth only where there are enough distinct x values to justify it
if (!(var %in% categorical_vars) && nrow(pd_summary) >= 10) {
pd_summary %>%
mutate(y_smooth = loess(y ~ x, span = 0.6)$fitted,
ymin_smooth = loess(ymin ~ x, span = 0.6)$fitted,
ymax_smooth = loess(ymax ~ x, span = 0.6)$fitted)
} else {
pd_summary %>%
mutate(y_smooth = y, ymin_smooth = ymin, ymax_smooth = ymax)
}
})
names(pdp_data_list) <- selected_env
## ---- 5. Shared y-axis across panels -------------------------------------
all_pdp <- bind_rows(pdp_data_list)
y_rng <- range(c(all_pdp$ymin_smooth, all_pdp$ymax_smooth), na.rm = TRUE)
y_pad <- 0.05 * diff(y_rng)
y_limits <- c(y_rng[1] - y_pad, y_rng[2] + y_pad)
## ---- 6. Build panels ----------------------------------------------------
panel_labels <- LETTERS[seq_along(selected_env)]
pdp_boot_list <- lapply(seq_along(selected_env), function(i) {
var <- selected_env[i]
d <- pdp_data_list[[var]]
lab <- sprintf("%s %%IncMSE = %.1f", panel_labels[i], imp_lookup[[var]])
if (var %in% categorical_vars) {
d$x_lab <- factor(soil_levels[d$x], levels = soil_levels)
p <- ggplot(d, aes(x = x_lab)) +
geom_errorbar(aes(ymin = ymin, ymax = ymax),
width = 0.15, colour = gator_blue) +
geom_point(aes(y = y), colour = gator_blue, size = 2.5)
} else {
p <- ggplot(d, aes(x = x)) +
geom_ribbon(aes(ymin = ymin_smooth, ymax = ymax_smooth),
fill = gator_orange, alpha = 0.25) +
geom_line(aes(y = y_smooth), colour = gator_blue, linewidth = 1.2)
}
p +
coord_cartesian(ylim = y_limits) +
labs(x = var_labels[[var]], y = NULL) +
annotate("text", x = -Inf, y = Inf, label = lab,
hjust = -0.08, vjust = 1.5, fontface = "bold", size = 3.4) +
theme_classic(base_size = 12) +
theme(
axis.title.x = element_text(size = 10),
axis.text = element_text(size = 9),
axis.text.x = if (var %in% categorical_vars)
element_text(angle = 45, hjust = 1, size = 8) else element_text(size = 9),
plot.margin = ggplot2::margin(5, 5, 5, 5)
)
})
## ---- 7. Save ------------------------------------------------------------
## 3 columns keeps each panel narrow, which is what Victoria asked for in
## comment 106 - narrower panels make the trends read steeper.
png("partial_dependence_bootstrap.png",
width = 1050, height = 1150, res = 150)
grid.arrange(
textGrob("Shannon Diversity", rot = 90,
gp = gpar(fontsize = 14, fontface = "bold")),
arrangeGrob(grobs = pdp_boot_list, ncol = 3),
ncol = 2, widths = c(0.05, 0.95)
)
dev.off()## png
## 2cat("\nPanel order (A onward):\n")##
## Panel order (A onward):print(data.frame(Panel = panel_labels,
Variable = selected_env,
IncMSE = round(unname(imp_lookup[selected_env]), 2)))## Panel Variable IncMSE
## 1 A Cogongrass_Cover 40.88
## 2 B pr 46.53
## 3 C srad 44.24
## 4 D avg_pop_den_1km 31.04
## 5 E TreeCanopy_NLCD 26.29
## 6 F elevation 23.16
## 7 G SoilType 20.88
## 8 H ndvi 17.01
## 9 I slope 13.99Effect Magnitude
## Span of the bootstrapped partial dependence mean curve for each predictor.
## The Abstract and Results now rest on this: cogongrass ranks 4th by %IncMSE
## but should show the largest span here. Requires pdp_data_list and imp_df
## from the PDP_pub chunk above.
## Plain-text labels. Do NOT use clean_label() here: the PDP_pub chunk above
## overwrites `var_labels` with a list of plot expressions, so clean_label()
## returns 'expression(paste("Surface Radiation (W/", m^2, ")"))' as a string.
plain_labels <- c(
Cogongrass_Cover = "Cogongrass Cover",
pr = "Precipitation",
srad = "Surface Radiation",
TreeCanopy_NLCD = "Tree Canopy Cover",
elevation = "Elevation",
slope = "Slope",
SoilType = "Soil Order",
avg_pop_den_1km = "Population Density",
ndvi = "NDVI",
sph = "Specific Humidity"
)
pdp_spans <- bind_rows(pdp_data_list) %>%
group_by(Variable) %>%
summarise(
pdp_min = min(y_smooth, na.rm = TRUE),
pdp_max = max(y_smooth, na.rm = TRUE),
span = pdp_max - pdp_min,
.groups = "drop"
) %>%
left_join(imp_df[, c("Variable", "%IncMSE")], by = "Variable") %>%
mutate(Label = unname(plain_labels[Variable])) %>%
arrange(desc(span)) %>%
mutate(across(where(is.numeric), ~round(.x, 3)))
print(as.data.frame(pdp_spans[, c("Label", "span", "pdp_min", "pdp_max", "%IncMSE")]))## Label span pdp_min pdp_max %IncMSE
## 1 Cogongrass Cover 0.517 1.847 2.364 40.880
## 2 Precipitation 0.360 2.125 2.484 46.531
## 3 Elevation 0.291 2.069 2.360 23.157
## 4 Population Density 0.268 2.082 2.349 31.035
## 5 NDVI 0.197 2.161 2.358 17.009
## 6 Surface Radiation 0.166 2.286 2.452 44.237
## 7 Slope 0.156 2.198 2.353 13.993
## 8 Tree Canopy Cover 0.135 2.228 2.363 26.288
## 9 Soil Order 0.075 2.321 2.396 20.880cog <- pdp_spans %>% filter(Variable == "Cogongrass_Cover")
top_env <- pdp_spans %>% filter(Variable != "Cogongrass_Cover") %>% slice_max(span, n = 1)
cat(sprintf(
"\n Cogongrass span : %.3f units\n Largest environmental span : %.3f units (%s)\n Ratio : %.1fx\n",
cog$span, top_env$span, top_env$Label, cog$span / top_env$span))##
## Cogongrass span : 0.517 units
## Largest environmental span : 0.360 units (Precipitation)
## Ratio : 1.4xSoil Order Contrast
## The Results claim soil order "showed no meaningful differentiation, with
## overlapping confidence intervals across all four orders". This quantifies it.
##
## The right statistic is a SIGNED pairwise contrast from the bootstrap, whose
## interval can include zero. A max-minus-min range cannot include zero by
## construction, so it would look significant no matter what.
##
## Requires boot_fits, pred_grids and soil_levels from the PDP_pub chunk.
soil_boot <- map_dfr(seq_along(boot_fits), function(j) {
pd <- partial(boot_fits[[j]]$fit,
pred.var = "SoilType",
pred.grid = pred_grids[["SoilType"]],
train = boot_fits[[j]]$dat)
pd$iter <- j
pd
}) %>%
mutate(Order = soil_levels[SoilType])
## Every pairwise contrast, with bootstrap 95% intervals
soil_wide <- soil_boot %>%
dplyr::select(iter, Order, yhat) %>%
pivot_wider(names_from = Order, values_from = yhat)
pairs_tbl <- combn(soil_levels, 2, simplify = FALSE) %>%
map_dfr(function(pr) {
d <- soil_wide[[pr[1]]] - soil_wide[[pr[2]]]
tibble(contrast = paste(pr[1], "-", pr[2]),
median = median(d),
lo = quantile(d, 0.025),
hi = quantile(d, 0.975),
excl_0 = !(quantile(d, 0.025) < 0 & quantile(d, 0.975) > 0))
}) %>%
mutate(across(where(is.numeric), ~round(.x, 3)))
cat("===== Soil order: pairwise partial dependence contrasts =====\n")## ===== Soil order: pairwise partial dependence contrasts =====print(as.data.frame(pairs_tbl))## contrast median lo hi excl_0
## 1 Alfisols - Entisols 0.054 0.009 0.138 TRUE
## 2 Alfisols - Spodosols 0.064 0.003 0.143 TRUE
## 3 Alfisols - Ultisols 0.068 0.004 0.154 TRUE
## 4 Entisols - Spodosols 0.005 -0.006 0.017 FALSE
## 5 Entisols - Ultisols 0.012 -0.011 0.060 FALSE
## 6 Spodosols - Ultisols 0.009 -0.015 0.045 FALSEcat(sprintf("\nContrasts whose 95%% interval excludes zero: %d of %d\n",
sum(pairs_tbl$excl_0), nrow(pairs_tbl)))##
## Contrasts whose 95% interval excludes zero: 3 of 6## Same style of contrast for cogongrass, for comparison: predicted diversity
## at 0%% cover vs at the highest observed cover.
cog_boot <- map_dfr(seq_along(boot_fits), function(j) {
pd <- partial(boot_fits[[j]]$fit,
pred.var = "Cogongrass_Cover",
pred.grid = pred_grids[["Cogongrass_Cover"]],
train = boot_fits[[j]]$dat)
pd$iter <- j
pd
})
cog_ends <- cog_boot %>%
group_by(iter) %>%
summarise(d = yhat[which.max(Cogongrass_Cover)] - yhat[which.min(Cogongrass_Cover)],
.groups = "drop")
cat(sprintf(
"\nCogongrass: highest observed cover vs 0%%\n median %.3f units, 95%% interval %.3f to %.3f\n",
median(cog_ends$d), quantile(cog_ends$d, 0.025), quantile(cog_ends$d, 0.975)))##
## Cogongrass: highest observed cover vs 0%
## median -0.541 units, 95% interval -0.810 to -0.112## Data-side check, independent of the model. Spodosols has n = 1, so it is
## excluded; report the test both ways if you prefer.
soil_obs <- model_df_sel %>%
mutate(Order = soil_levels[SoilType]) %>%
filter(Order != "Spodosols")
kw <- kruskal.test(Shannon_Diversity ~ factor(Order), data = soil_obs)
cat(sprintf(
"\nObserved Shannon diversity by soil order (Spodosols excluded, n = %d):\n Kruskal-Wallis chi-squared = %.3f, df = %d, P = %.3f\n",
nrow(soil_obs), kw$statistic, kw$parameter, kw$p.value))##
## Observed Shannon diversity by soil order (Spodosols excluded, n = 205):
## Kruskal-Wallis chi-squared = 8.998, df = 2, P = 0.011print(soil_obs %>% group_by(Order) %>%
summarise(n = n(), median_SHDI = round(median(Shannon_Diversity), 3),
.groups = "drop") %>% as.data.frame())## Order n median_SHDI
## 1 Alfisols 31 2.645
## 2 Entisols 31 2.264
## 3 Ultisols 143 2.368Importance Rank Stability
## The manuscript states cogongrass ranks 4th. That comes from a single
## createDataPartition draw, and Boruta (all 206 plots) put it 2nd. This
## repeats the split to see how stable the rank actually is.
set.seed(97)
n_rep <- 100
rank_reps <- map_dfr(seq_len(n_rep), function(i) {
idx <- createDataPartition(model_df_sel$Shannon_Diversity, p = 0.8, list = FALSE)
rf <- randomForest(Shannon_Diversity ~ .,
data = model_df_sel[idx, ],
ntree = 1000, mtry = best_mtry, importance = TRUE)
as.data.frame(importance(rf)) %>%
rownames_to_column("Variable") %>%
mutate(rank = rank(-`%IncMSE`), iter = i) %>%
dplyr::select(Variable, IncMSE = `%IncMSE`, rank, iter)
})
rank_summary <- rank_reps %>%
group_by(Variable) %>%
summarise(median_rank = median(rank),
rank_min = min(rank),
rank_max = max(rank),
pct_top1 = 100 * mean(rank == 1),
median_imp = median(IncMSE),
imp_lo = quantile(IncMSE, 0.025),
imp_hi = quantile(IncMSE, 0.975),
.groups = "drop") %>%
arrange(median_rank) %>%
mutate(across(where(is.numeric), ~round(.x, 2)))
cat("===== Importance rank across", n_rep, "train/test splits =====\n")## ===== Importance rank across 100 train/test splits =====print(as.data.frame(rank_summary))## Variable median_rank rank_min rank_max pct_top1 median_imp imp_lo
## 1 pr 1 1 4 56 18.20 13.28
## 2 srad 2 1 4 26 16.48 11.90
## 3 Cogongrass_Cover 3 1 9 18 14.49 2.36
## 4 avg_pop_den_1km 4 2 9 0 12.23 6.10
## 5 TreeCanopy_NLCD 5 3 9 0 10.83 7.34
## 6 elevation 6 4 9 0 8.73 4.90
## 7 SoilType 7 4 9 0 7.90 3.95
## 8 ndvi 8 2 9 0 6.45 2.73
## 9 slope 9 4 9 0 5.08 0.66
## imp_hi
## 1 22.62
## 2 21.80
## 3 20.86
## 4 16.77
## 5 14.38
## 6 11.12
## 7 12.71
## 8 11.33
## 9 9.64cog_r <- rank_reps %>% filter(Variable == "Cogongrass_Cover")
cat(sprintf(
"\nCogongrass cover: median rank %.0f (range %d-%d), ranked 1st in %.0f%% of splits\n",
median(cog_r$rank), min(cog_r$rank), max(cog_r$rank), 100 * mean(cog_r$rank == 1)))##
## Cogongrass cover: median rank 3 (range 1-9), ranked 1st in 18% of splits## Split-free ranking: the all-206-plot model already fitted in PDP_pub.
## Permutation importance is computed out-of-bag, so this is legitimate.
full_imp <- as.data.frame(importance(rf_pdp)) %>%
rownames_to_column("Variable") %>%
arrange(desc(`%IncMSE`)) %>%
mutate(across(where(is.numeric), ~round(.x, 2)))
cat("\n===== Importance from a model fit to all 206 plots =====\n")##
## ===== Importance from a model fit to all 206 plots =====print(as.data.frame(full_imp[, c("Variable", "%IncMSE")]))## Variable %IncMSE
## 1 pr 46.53
## 2 srad 44.24
## 3 Cogongrass_Cover 40.88
## 4 avg_pop_den_1km 31.04
## 5 TreeCanopy_NLCD 26.29
## 6 elevation 23.16
## 7 SoilType 20.88
## 8 ndvi 17.01
## 9 slope 13.99Summary 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 : 9cat("------------------------------------------------------------\n")## ------------------------------------------------------------cat(sprintf(" Test RMSE : %.4f\n", rmse_val))## Test RMSE : 0.5305cat(sprintf(" Test MAE : %.4f\n", mae_val))## Test MAE : 0.3938cat(sprintf(" Test Bias : %.4f\n", bias_val))## Test Bias : 0.0106cat(sprintf(" Relative Bias : %.2f%%\n", rel_bias))## Relative Bias : 0.45%cat(sprintf(" Test R² : %.4f\n", r2_val))## Test R² : 0.2039cat("============================================================\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'