Set-up

Install necessary packages and import appropriate data

knitr::opts_chunk$set(echo = TRUE)

pacman::p_load(tidyverse, readxl, raster, vegan, tigris, sf, sjPlot, sp, spOccupancy, ggrepel, lme4, lmerTest, MuMIn, brms, MCMCvis, bayesplot, patchwork)

set.seed(97)

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

# CWD Data
cwd_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
                       sheet = "CWD")

# Soil Data
fuel_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Fuel_Sampling")

# 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
CameraData <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/04_Wildlife/02_Data/CameraData.xlsx")

CameraLoc <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/04_Wildlife/02_Data/CameraLoc.xlsx",
                  sheet = "CameraLocations")

# Add effort data
effort_matrix <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/04_Wildlife/02_Data/CameraLoc.xlsx",
                            sheet = "Effort_Matrix_Full") %>%
  pivot_longer(cols = matches("^202[4-5]-"), names_to = "week", values_to = "days") %>%
  filter(days == "7") %>%
  dplyr::select(Plot, week)

Coordinate Transformation

coords_sf <- CameraLoc %>%
  dplyr::select(Plot, Long, Lat) %>%
  st_as_sf(coords = c("Long", "Lat"), crs = 4326) %>%
  st_transform(crs = 32616)  # UTM zone 16N

coords <- st_coordinates(coords_sf)
rownames(coords) <- CameraLoc$Plot

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

Provides the option to rarify vegetative species

# Calculate the total number of sites
total_sites <- nrow(CameraLoc)

# Function to filter data by frequency
filter_by_frequency <- function(df) {
  # Group data by species and calculate the frequency
  freq <- df %>%
    group_by(Species) %>%
    summarise(Frequency = n_distinct(Plot) / nrow(CameraLoc) * 100) %>%
    filter(Frequency >= 0)
  
  # Filter the original data to include only species with frequency >= XX%
  filtered_df <- df %>%
    filter(Species %in% freq$Species)
  
  return(filtered_df)
}

# Filter tree data by frequency
tree_data <- filter_by_frequency(tree_data)

# Filter Veg Cover data by frequency
Veg_Cover <- filter_by_frequency(Veg_Cover)

# Filter Shrub Cover data by frequency
shrub_data <- filter_by_frequency(shrub_data)

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) %>%
  summarise(Shrub_Total_Cover = sum(Cover, na.rm = TRUE), .groups = "drop") %>%
  mutate(Shrub_Percent_Cover = Shrub_Total_Cover / 3000 * 100)

# Summed length of shrub over at a site
shrub_cover_summed <- shrub_cover %>%
  group_by(Plot) %>%
  summarize(total_shrub_cover = sum(Shrub_Total_Cover, na.rm = TRUE), .groups = "drop")

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(CameraLoc, by = "Plot")

# Sum species cover across quadrats for each species at each plot
veg_cover_summed <- Veg_Cover %>%
  group_by(Plot, Species_Name) %>%
  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, Shrub_Percent_Cover),
    by = c("Plot", "Species_Name")
  ) %>%
  mutate(
    overlap_flag = ifelse(!is.na(avg_cover) & !is.na(Shrub_Percent_Cover), TRUE, FALSE), # Flag overlaps
    final_cover = case_when(
      !is.na(avg_cover) & is.na(Shrub_Percent_Cover) ~ avg_cover,  # Use herbaceous cover if no shrub data
      is.na(avg_cover) & !is.na(Shrub_Percent_Cover) ~ Shrub_Percent_Cover, # Use shrub cover if no herbaceous data
      TRUE ~ NA_real_ # NA where overlaps exist
    )
  )

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

Summarize Cogongrass Cover

avg_cogongrass_cover <- species_matrix %>%
  group_by(Plot) %>%
  summarize(Avg_Cogongrass_Cover = sum(Imperata_cylindrica, na.rm = TRUE) / n(), .groups = "drop")

Herbacous Shannon Diversity Index

# Summarize species cover by site
site_species_cover <- Veg_Cover %>%
  group_by(Plot, Species_Name) %>%
  summarize(total_cover = sum(Cover_Per, na.rm = TRUE)) %>%
  ungroup()
## `summarise()` has grouped output by 'Plot'. You can override using the
## `.groups` argument.
## Remove all Imperata_cylindrica_Live and Imperata_cylindrica from species
site_species_cover <- site_species_cover %>%
  filter(Species_Name != "Imperata_cylindrica_Live" & Species_Name != "Imperata_cylindrica")

# Calculate Shannon diversity per site
Veg_shannon_diversity <- site_species_cover %>%
  group_by(Plot) %>%
  mutate(proportion = total_cover / sum(total_cover)) %>%
  summarize(Veg_shannon_index = -sum(proportion * log(proportion), na.rm = TRUE))

print(Veg_shannon_diversity)
## # A tibble: 206 × 2
##    Plot  Veg_shannon_index
##    <chr>             <dbl>
##  1 BI200              2.75
##  2 BI201              2.70
##  3 BI202              2.59
##  4 BI97               1.61
##  5 BI99               2.97
##  6 BN210              2.97
##  7 BN211              2.43
##  8 BN212              2.22
##  9 BN96               3.05
## 10 BN98               2.79
## # ℹ 196 more rows

Vegetation Height

if (!is.numeric(fuel_data$Height)) {
  fuel_data$Height <- as.numeric(as.character(fuel_data$Height))
}
## Warning: NAs introduced by coercion
# Calculate average vegetation height per plot
veg_height <- fuel_data %>%
  group_by(Plot) %>%
  summarize(avg_veg_height = mean(Height, na.rm = TRUE), .groups = "drop")

Tree Metrics

# Tree density from point-centered quarter data
if (!is.numeric(tree_data$Distance)) {
  tree_data$Distance <- as.numeric(as.character(tree_data$Distance))
}

tree_density_data <- tree_data %>%
  group_by(Plot) %>%
  summarize(Average_Distance = mean(Distance) / 100,  # Convert to meters
            Tree_Density = 10000 / (Average_Distance^2))  # Convert to trees per hectare

# Average canopy cover from vegetation quadrats
tree_canopy_data <- Veg_Cover %>%
  distinct(Plot, Quadrat, .keep_all = TRUE) %>%  # Ensure each quadrat counts once per plot
  group_by(Plot) %>%
  summarize(Avg_Canopy_Cover = mean(Canopy_Cover, na.rm = TRUE), .groups = "drop") # Calculate the average canopy cover per plot

cor(tree_density_data$Tree_Density, tree_canopy_data$Avg_Canopy_Cover)
## [1] 0.2742307

CWD Clean

cwd_pa <- cwd_data %>%
  group_by(Plot) %>%
  summarise(
    CWD_presence = as.integer(any(Line_Start > 0 & Line_End > 0)),
    .groups = "drop"
  )

Plot Species level relationship with cogongrass cover

# Collecting scaled atrributes
cogon.mean <- attr(covariates_matrix, "scaled:center")["Avg_Cogongrass_Cover"]
cogon.sd <- attr(covariates_matrix, "scaled:scale")["Avg_Cogongrass_Cover"]

cogon.mean <- mean(CameraLoc$Avg_Cogongrass_Cover)
cogon.sd <- sd(CameraLoc$Avg_Cogongrass_Cover)

# Unscaled prediction values
cogon.pred.vals <- seq(0, 100, length.out = 100)

# Scale prediction values
cogon.pred.scale <- (cogon.pred.vals - cogon.mean) / cogon.sd
cogon.pred.scale2 <- cogon.pred.scale^2

canopy.mean <- attr(covariates_matrix, "scaled:center")["Avg_Canopy_Cover"]
canopy.sd <- attr(covariates_matrix, "scaled:scale")["Avg_Canopy_Cover"]

pred.df <- as.matrix(data.frame(
  intercept = 1,
  Avg_Cogongrass_Cover = cogon.pred.scale,
  I.Avg_Cogongrass_Cover.2. = cogon.pred.scale2,
  Avg_Canopy_Cover = 0  
))

out.pred <- predict(
  ms_cover_canopyQ,
  X.0 = pred.df,
  ignore.RE = TRUE
)


# Extract the species-specific betas for Coyote
coyote_index <- which(ms_cover_canopyQ$sp.names == "Canis_latrans")

# Get column indices for the parameters we need
intercept_col <- paste0("(Intercept)-", "Canis_latrans")
cogon_col <- paste0("Avg_Cogongrass_Cover-", "Canis_latrans")
cogon2_col <- paste0("I(Avg_Cogongrass_Cover^2)-", "Canis_latrans")
canopy_col <- paste0("Avg_Canopy_Cover-", "Canis_latrans")

# Extract posterior samples for each coefficient
beta_intercept <- ms_cover_canopyQ$beta.samples[, intercept_col]
beta_cogon <- ms_cover_canopyQ$beta.samples[, cogon_col]
beta_cogon2 <- ms_cover_canopyQ$beta.samples[, cogon2_col]
beta_canopy <- ms_cover_canopyQ$beta.samples[, canopy_col]

# Create prediction gradient (0 to 100%)
cogon.pred.vals <- seq(0, 100, length.out = 1000)

# Scale using same transformation
cogon.pred.scale <- (cogon.pred.vals - cogon.mean) / cogon.sd
cogon.pred.scale2 <- cogon.pred.scale^2

n_samples <- nrow(ms_cover_canopyQ$beta.samples)
n_pred <- length(cogon.pred.vals)

logit_psi <- matrix(NA, nrow = n_samples, ncol = n_pred)

for(i in 1:n_samples) {
  logit_psi[i, ] <- beta_intercept[i] + 
                    beta_cogon[i] * cogon.pred.scale + 
                    beta_cogon2[i] * cogon.pred.scale2 + 
                    beta_canopy[i] * 0  # Setting to mean (0 on scaled)
}

# Convert to probability scale
psi_samples <- plogis(logit_psi)

# Calculate quantiles across MCMC samples for each prediction point
psi.coyote.quants <- apply(psi_samples, 2, quantile, 
                            probs = c(0.025, 0.5, 0.975))

psi.coyote.dat <- data.frame(
  Avg_Cogongrass_Cover = cogon.pred.vals,
  psi.low  = psi.coyote.quants[1, ],
  psi.med  = psi.coyote.quants[2, ],
  psi.high = psi.coyote.quants[3, ]
)

ggplot(psi.coyote.dat, aes(x = Avg_Cogongrass_Cover, y = psi.med)) +
  geom_ribbon(aes(ymin = psi.low, ymax = psi.high), fill = "grey70", alpha = 0.5) +
  geom_line(size = 1.2) +
  theme_bw() +
  labs(x = "Average Cogongrass Cover (%)",
       y = "Occupancy Probability (Canis latrans)") +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_continuous(limits = c(0, 100))

summary(psi.coyote.dat)
##  Avg_Cogongrass_Cover    psi.low          psi.med          psi.high     
##  Min.   :  0          Min.   :0.1597   Min.   :0.3878   Min.   :0.6769  
##  1st Qu.: 25          1st Qu.:0.2400   1st Qu.:0.5267   1st Qu.:0.8324  
##  Median : 50          Median :0.5635   Median :0.9855   Median :1.0000  
##  Mean   : 50          Mean   :0.5581   Mean   :0.8087   Mean   :0.9222  
##  3rd Qu.: 75          3rd Qu.:0.8687   3rd Qu.:1.0000   3rd Qu.:1.0000  
##  Max.   :100          Max.   :0.9733   Max.   :1.0000   Max.   :1.0000

Community Occupancy

# Community-level (average across species)
comm_intercept <- ms_cover_canopyQ$beta.comm.samples[, "(Intercept)"]
comm_cogon <- ms_cover_canopyQ$beta.comm.samples[, "Avg_Cogongrass_Cover"]
comm_cogon2 <- ms_cover_canopyQ$beta.comm.samples[, "I(Avg_Cogongrass_Cover^2)"]
comm_canopy <- ms_cover_canopyQ$beta.comm.samples[, "Avg_Canopy_Cover"]

# Same prediction approach
logit_psi_comm <- matrix(NA, nrow = nrow(ms_cover_canopyQ$beta.comm.samples), 
                         ncol = n_pred)

for(i in 1:nrow(ms_cover_canopyQ$beta.comm.samples)) {
  logit_psi_comm[i, ] <- comm_intercept[i] + 
                         comm_cogon[i] * cogon.pred.scale + 
                         comm_cogon2[i] * cogon.pred.scale2 + 
                         comm_canopy[i] * 0
}

psi_comm_samples <- plogis(logit_psi_comm)

psi.comm.quants <- apply(psi_comm_samples, 2, quantile, 
                         probs = c(0.025, 0.5, 0.975))

psi.comm.dat <- data.frame(
  Avg_Cogongrass_Cover = cogon.pred.vals,
  psi.low  = psi.comm.quants[1, ],
  psi.med  = psi.comm.quants[2, ],
  psi.high = psi.comm.quants[3, ]
)

ggplot(psi.comm.dat, aes(x = Avg_Cogongrass_Cover, y = psi.med)) +
  geom_ribbon(aes(ymin = psi.low, ymax = psi.high), fill = "grey70", alpha = 0.5) +
  geom_line(size = 1.2) +
  theme_bw() +
  labs(x = "Average Cogongrass Cover (%)",
       y = "Community Occupancy Probability") +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_continuous(limits = c(0, 100))

Loop through all species in the community

# Prediction gradient
cogon.pred.vals <- seq(0, 100, length.out = 1000)
cogon.pred.scale <- (cogon.pred.vals - cogon.mean) / cogon.sd
cogon.pred.scale2 <- cogon.pred.scale^2

n_samples <- nrow(ms_cover_canopyQ$beta.samples)
n_pred <- length(cogon.pred.vals)

# List to store results
species_predictions <- list()
species_plots <- list()

# Loop through all species
for(sp in ms_cover_canopyQ$sp.names) {
  
  # Get column names for this species
  intercept_col <- paste0("(Intercept)-", sp)
  cogon_col <- paste0("Avg_Cogongrass_Cover-", sp)
  cogon2_col <- paste0("I(Avg_Cogongrass_Cover^2)-", sp)
  canopy_col <- paste0("Avg_Canopy_Cover-", sp)
  
  # Extract posterior samples
  beta_intercept <- ms_cover_canopyQ$beta.samples[, intercept_col]
  beta_cogon <- ms_cover_canopyQ$beta.samples[, cogon_col]
  beta_cogon2 <- ms_cover_canopyQ$beta.samples[, cogon2_col]
  beta_canopy <- ms_cover_canopyQ$beta.samples[, canopy_col]
  
  # Calculate linear predictor
  logit_psi <- matrix(NA, nrow = n_samples, ncol = n_pred)
  
  for(i in 1:n_samples) {
    logit_psi[i, ] <- beta_intercept[i] + 
                      beta_cogon[i] * cogon.pred.scale + 
                      beta_cogon2[i] * cogon.pred.scale2 + 
                      beta_canopy[i] * 0
  }
  
  # Convert to probability
  psi_samples <- plogis(logit_psi)
  
  # Calculate quantiles
  psi.quants <- apply(psi_samples, 2, quantile, probs = c(0.025, 0.5, 0.975))
  
  # Store results
  species_predictions[[sp]] <- data.frame(
    species = sp,
    Avg_Cogongrass_Cover = cogon.pred.vals,
    psi.low  = psi.quants[1, ],
    psi.med  = psi.quants[2, ],
    psi.high = psi.quants[3, ]
  )
  
  # Create plot
  species_plots[[sp]] <- ggplot(species_predictions[[sp]], 
                                 aes(x = Avg_Cogongrass_Cover, y = psi.med)) +
    geom_ribbon(aes(ymin = psi.low, ymax = psi.high), fill = "grey70", alpha = 0.5) +
    geom_line(size = 1.2) +
    theme_bw() +
    labs(x = "Average Cogongrass Cover (%)",
         y = "Occupancy Probability",
         title = gsub("_", " ", sp)) +  # Replace underscores with spaces
    scale_y_continuous(limits = c(0, 1)) +
    scale_x_continuous(limits = c(0, 100))
}

# View individual plots
species_plots[["Canis_latrans"]]

species_plots[["Procyon_lotor"]]

# Combine all plots into a grid
all_species_plot <- wrap_plots(species_plots, ncol = 4)
all_species_plot

# Save combined plot
ggsave("all_species_cogongrass_response.png", 
       all_species_plot, 
       width = 12, height = 16, dpi = 300)

# Combine all predictions into one dataframe for faceted plot
all_predictions <- do.call(rbind, species_predictions)

# Create faceted plot
faceted_plot <- ggplot(all_predictions, aes(x = Avg_Cogongrass_Cover, y = psi.med)) +
  geom_ribbon(aes(ymin = psi.low, ymax = psi.high), fill = "grey70", alpha = 0.5) +
  geom_line(size = 1.2) +
  facet_wrap(~ species, ncol = 2, labeller = labeller(species = function(x) gsub("_", " ", x))) +
  theme_bw() +
  labs(x = "Average Cogongrass Cover (%)",
       y = "Occupancy Probability") +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_continuous(limits = c(0, 100)) +
  theme(strip.text = element_text(face = "italic"))  # Italicize species names

faceted_plot

ggsave("all_species_faceted.png", 
       faceted_plot, 
       width = 10, height = 12, dpi = 300)

Species Lookup

SpeciesLookup <- tibble::tribble(
  ~Name, ~Scientific_Name, ~Common_Name,
  "Canis_latrans", "Canis latrans", "Coyote",
  "Dasypus_novemcinctus", "Dasypus novemcinctus", "Nine-banded Armadillo",
  "Didelphis_virginiana", "Didelphis virginiana", "Virginia Opossum",
  "Lynx_rufus", "Lynx rufus", "Bobcat",
  "Meleagris_gallopavo", "Meleagris gallopavo", "Wild Turkey",
  "Procyon_lotor", "Procyon lotor", "Raccoon",
  "Sciurus_carolinensis", "Sciurus carolinensis", "Eastern Gray Squirrel",
  "Sylvilagus_floridanus", "Sylvilagus floridanus", "Eastern Cottontail Rabbit",
)

Varying cogongrass and canopy cover

# Prediction gradient for cogongrass
cogon.pred.vals <- seq(0, 100, length.out = 1000)
cogon.pred.scale <- (cogon.pred.vals - cogon.mean) / cogon.sd
cogon.pred.scale2 <- cogon.pred.scale^2

# Define canopy cover scenarios with factor levels
canopy_scenarios <- data.frame(
  scenario = factor(c("Low", "Average", "High"),
                   levels = c("Low", "Average", "High")),
  canopy_scaled = c(-1, 0, 1),
  color = c("#D55E00", "#009E73", "#0072B2")
)

n_samples <- nrow(ms_cover_canopyQ$beta.samples)
n_pred <- length(cogon.pred.vals)

# Lists for results
species_predictions <- list()
species_plots <- list()

# Loop through all species
for(sp in ms_cover_canopyQ$sp.names) {
  
  # Get column names for this species
  intercept_col <- paste0("(Intercept)-", sp)
  cogon_col <- paste0("Avg_Cogongrass_Cover-", sp)
  cogon2_col <- paste0("I(Avg_Cogongrass_Cover^2)-", sp)
  canopy_col <- paste0("Avg_Canopy_Cover-", sp)
  
  # Extract posterior samples
  beta_intercept <- ms_cover_canopyQ$beta.samples[, intercept_col]
  beta_cogon <- ms_cover_canopyQ$beta.samples[, cogon_col]
  beta_cogon2 <- ms_cover_canopyQ$beta.samples[, cogon2_col]
  beta_canopy <- ms_cover_canopyQ$beta.samples[, canopy_col]
  
  # Store predictions for each canopy scenario
  scenario_predictions <- list()
  
  for(j in 1:nrow(canopy_scenarios)) {
    
    canopy_val <- canopy_scenarios$canopy_scaled[j]
    
    # Calculate linear predictor
    logit_psi <- matrix(NA, nrow = n_samples, ncol = n_pred)
    
    for(i in 1:n_samples) {
      logit_psi[i, ] <- beta_intercept[i] + 
                        beta_cogon[i] * cogon.pred.scale + 
                        beta_cogon2[i] * cogon.pred.scale2 + 
                        beta_canopy[i] * canopy_val
    }
    
    # Convert to probability
    psi_samples <- plogis(logit_psi)
    
    # Calculate quantiles
    psi.quants <- apply(psi_samples, 2, quantile, probs = c(0.025, 0.5, 0.975))
    
    # Store results
    scenario_predictions[[j]] <- data.frame(
      species = sp,
      scenario = canopy_scenarios$scenario[j],
      Avg_Cogongrass_Cover = cogon.pred.vals,
      psi.low  = psi.quants[1, ],
      psi.med  = psi.quants[2, ],
      psi.high = psi.quants[3, ]
    )
  }
  
  # Combine all scenarios for this species
  species_predictions[[sp]] <- do.call(rbind, scenario_predictions)
  
  all_predictions <- all_predictions %>%
  left_join(SpeciesLookup, by = c("species" = "Name"))
  # Create plot for this species
  species_plots[[sp]] <- ggplot(species_predictions[[sp]], 
                                 aes(x = Avg_Cogongrass_Cover, y = psi.med, 
                                     color = scenario, fill = scenario)) +
    geom_ribbon(aes(ymin = psi.low, ymax = psi.high), alpha = 0.2, color = NA) +
    geom_line(size = 1.2) +
    scale_color_manual(values = canopy_scenarios$color) +
    scale_fill_manual(values = canopy_scenarios$color) +
    theme_bw() +
    labs(x = "Average Cogongrass Cover (%)",
         y = "Occupancy Probability",
         title = gsub("_", " ", sp),
         color = "Canopy Cover",
         fill = "Canopy Cover") +
    scale_y_continuous(limits = c(0, 1)) +
    scale_x_continuous(limits = c(0, 100)) +
    theme(legend.position = "bottom")
  
}

# View individual plots
species_plots[["Canis_latrans"]]

species_plots[["Procyon_lotor"]]

# Combine all plots into a grid
all_species_plot <- wrap_plots(species_plots, ncol = 2, guides = "collect") & 
  theme(legend.position = "bottom")

all_species_plot

ggsave("all_species_cogongrass_canopy.png", 
       all_species_plot, 
       width = 14, height = 16, dpi = 300)

# Combine all predictions into one dataframe
all_predictions <- do.call(rbind, species_predictions)

all_predictions <- all_predictions %>%
  left_join(
    SpeciesLookup %>% dplyr::select(Name, Common_Name),
    by = c("species" = "Name")
  )

all_predictions$scenario <- factor(all_predictions$scenario,
                                   levels = c("Low", "Average", "High"))

# Create faceted plot
faceted_plot <- ggplot(all_predictions, 
                       aes(x = Avg_Cogongrass_Cover, y = psi.med, 
                           color = scenario, fill = scenario)) +
  geom_ribbon(aes(ymin = psi.low, ymax = psi.high), alpha = 0.2, color = NA) +
  geom_line(size = 1.2) +
  facet_wrap(~ Common_Name, ncol = 2) +
  scale_color_manual(values = canopy_scenarios$color) +
  scale_fill_manual(values = canopy_scenarios$color) +
  theme_bw() +
  labs(x = "Average Cogongrass Cover (%)",
       y = "Occupancy Probability",
       color = "Canopy Cover",
       fill = "Canopy Cover") +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_continuous(limits = c(0, 100)) +
  theme(strip.text = element_text(face = "bold"),
      legend.position = "bottom")


faceted_plot

ggsave("all_species_faceted_canopy.png", 
       faceted_plot, 
       width = 12, height = 14, dpi = 600)

Community Occupancy Plot- Cogongrass and Canopy Cover

# Prediction gradient for cogongrass
cogon.pred.vals <- seq(0, 100, length.out = 1000)
cogon.pred.scale <- (cogon.pred.vals - cogon.mean) / cogon.sd
cogon.pred.scale2 <- cogon.pred.scale^2

# Define canopy cover scenarios
canopy_scenarios <- data.frame(
  scenario = c("Low", "Average", "High"),
  canopy_scaled = c(-1, 0, 1),
  color = c("#D55E00", "#009E73", "#0072B2")
)

n_pred <- length(cogon.pred.vals)

# Get community-level parameters
comm_intercept <- ms_cover_canopyQ$beta.comm.samples[, "(Intercept)"]
comm_cogon <- ms_cover_canopyQ$beta.comm.samples[, "Avg_Cogongrass_Cover"]
comm_cogon2 <- ms_cover_canopyQ$beta.comm.samples[, "I(Avg_Cogongrass_Cover^2)"]
comm_canopy <- ms_cover_canopyQ$beta.comm.samples[, "Avg_Canopy_Cover"]

# Store predictions for each canopy scenario
community_predictions <- list()

for(j in 1:nrow(canopy_scenarios)) {
  
  canopy_val <- canopy_scenarios$canopy_scaled[j]
  
  # Calculate linear predictor
  logit_psi_comm <- matrix(NA, 
                           nrow = nrow(ms_cover_canopyQ$beta.comm.samples), 
                           ncol = n_pred)
  
  for(i in 1:nrow(ms_cover_canopyQ$beta.comm.samples)) {
    logit_psi_comm[i, ] <- comm_intercept[i] + 
                           comm_cogon[i] * cogon.pred.scale + 
                           comm_cogon2[i] * cogon.pred.scale2 + 
                           comm_canopy[i] * canopy_val
  }
  
  # Convert to probability
  psi_comm_samples <- plogis(logit_psi_comm)
  
  # Calculate quantiles
  psi.comm.quants <- apply(psi_comm_samples, 2, quantile, 
                           probs = c(0.025, 0.5, 0.975))
  
  # Store results
  community_predictions[[j]] <- data.frame(
    scenario = canopy_scenarios$scenario[j],
    Avg_Cogongrass_Cover = cogon.pred.vals,
    psi.low  = psi.comm.quants[1, ],
    psi.med  = psi.comm.quants[2, ],
    psi.high = psi.comm.quants[3, ]
  )
}

# Combine all scenarios
psi.comm.dat <- do.call(rbind, community_predictions)

# Reorder the scenario factor levels
psi.comm.dat$scenario <- factor(psi.comm.dat$scenario, 
                                levels = c("Low", "Average", "High"))

community_plot <- ggplot(psi.comm.dat, 
                         aes(x = Avg_Cogongrass_Cover, y = psi.med, 
                             color = scenario, fill = scenario)) +
  geom_ribbon(aes(ymin = psi.low, ymax = psi.high), alpha = 0.2, color = NA) +
  geom_line(size = 1.5) +
  scale_color_manual(values = canopy_scenarios$color) +
  scale_fill_manual(values = canopy_scenarios$color) +
  theme_bw(base_size = 14) +
  labs(x = "Average Cogongrass Cover (%)",
       y = "Community Occupancy Probability",
       color = "Canopy Cover",
       fill = "Canopy Cover") +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_continuous(limits = c(0, 100)) +
  theme(legend.position = "bottom",
        plot.title = element_text(hjust = 0.5))

community_plot

ggsave("community_cogongrass_canopy.png", 
       community_plot, 
       width = 10, height = 7, dpi = 300)

write.csv(psi.comm.dat, "community_predictions_by_canopy.csv", row.names = FALSE)