C4_O3_Diversity

Set-up

Install necessary packages and import appropriate data

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

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

Number of quadrats sampled per plot

I moved this from a later section because the filtering process removed quadrats that did not capture any species. Rows labeled as “None” were removed, suggesting that the number of quadrats sampled per plot is not consistent across all plots.

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

Filter all data to only include species found at 3% of all sites

# 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 >= 3)
  
  # Filter the original data to include only species with frequency >= 3%
  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_ # Leave as 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.

# 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.28
##  2 BI201              2.20
##  3 BI202              1.50
##  4 BI97               1.82
##  5 BI99               3.06
##  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"
  )

Objective 1: Assessing if cogongrass invasion-related factors influence wildlife presence

Add Covariates

CameraLoc <- CameraLoc %>%
  left_join(cwd_pa, by = "Plot") %>%
  left_join(Veg_shannon_diversity, by = "Plot") %>%
  left_join(avg_cogongrass_cover, by = "Plot") %>%
  left_join(shrub_cover_summed %>% dplyr::select(Plot, total_shrub_cover), by = "Plot") %>%
  left_join(veg_height, by = "Plot") %>%
  left_join(tree_density_data %>% dplyr::select(Plot, Tree_Density), by = "Plot") %>%
  left_join(tree_canopy_data %>% dplyr::select(Plot, Avg_Canopy_Cover), by = "Plot") %>%
  dplyr::select(-Authority)

Subset to mammals with >1 observations

# Group by Name and count the number of observations
species_counts <- CameraData %>%
  filter(Class == "Mammalia" |
  Name == "Meleagris_gallopavo") %>%
  group_by(Name) %>%
  summarize(count = n(), .groups = "drop")

# Filter for species with count greater than 1
species_subset <- species_counts %>%
  filter(count > 1) %>%
  pull(Name)

# Filter CameraData to only include species with count > 1
CameraData <- CameraData %>%
  filter(Name %in% species_subset)

Create Detection Matrix

# Format Data Weekly
observations_weekly <- CameraData %>%
  group_by(Plot, week = format(as.Date(Date), "%Y-%U"), Name) %>%
  summarise(observations = n(), .groups = 'drop')

# Merge with Effort Matrix to include only valid weeks
observations_weekly <- effort_matrix %>%
  left_join(observations_weekly, by = c("Plot" = "Plot", "week")) %>%
  replace_na(list(observations = 0))

# Convert to wide format
observations_species <- observations_weekly %>%
  pivot_wider(names_from = Name, values_from = observations, values_fill = list(observations = 0)) %>%
  dplyr::select(-"NA")

# Create detection array
site_names <- sort(unique(observations_species$Plot))
species_names <- setdiff(colnames(observations_species), c("Plot", "week"))

num_sites <- length(site_names)
num_weeks <- length(unique(observations_species$week))
num_species <- length(species_names)

detection_array <- array(0, dim = c(num_sites, num_weeks, num_species))
dimnames(detection_array) <- list(site_names, unique(observations_species$week), species_names)

for (s in seq_along(species_names)) {
  species_col <- species_names[s]
  for (i in seq_along(site_names)) {
    site <- site_names[i]
    for (j in seq_along(unique(observations_species$week))) {
      week <- unique(observations_species$week)[j]
      detection_array[i, j, s] <- ifelse(
        any(observations_species$Plot == site & observations_species$week == week & observations_species[[species_col]] > 0),
        1, 0
      )
    }
  }
}

dim(detection_array) # Should be num_sites x num_weeks x num_species

## [1] 32 36 15

Prepare Site Covariates

# Duplicate CameraLoc to be used in Objective 2
CameraLoc_O2 <- CameraLoc

# Standardize the covariates
CameraLoc <- CameraLoc %>%
  dplyr::select(-Plot, -Camera, -Lat, -Long, -Status, - Start_Date)
covariates_matrix <- as.matrix(CameraLoc)
rownames(covariates_matrix) <- site_names

# Standardizing covariates
covariates_matrix <- scale(covariates_matrix)

Combine data into a list

This requires combining the presence data and the site covariate data into a single list. This also means that the presence data is in a 3-d format.

# Combine into a named list
data_list <- list(
  y = detection_array,
  occ.covs = covariates_matrix,
  det.covs = covariates_matrix
)

Convert occupancy and detection covariates to a dataframe (from matrix)

I am unsure why I only had an issue with total shrub cover, but this should fix the “cannot find” issue.

# Convert occupancy and detection covariates to a dataframe
data_list[["occ.covs"]] <- as.data.frame(data_list[["occ.covs"]])
data_list[["occ.covs"]]$total_shrub_cover <- as.numeric(data_list[["occ.covs"]]$total_shrub_cover)

data_list[["det.covs"]] <- as.data.frame(data_list[["det.covs"]])
data_list[["det.covs"]]$total_shrub_cover <- as.numeric(data_list[["det.covs"]]$total_shrub_cover)

# Make species the first dimension
data_list$y <- aperm(data_list$y, c(3, 1, 2))

dimnames(data_list$y) <- list(species = dimnames(data_list$y)[[1]],
                              site = dimnames(data_list$y)[[2]],
                              week = dimnames(data_list$y)[[3]])

str(data_list)

## List of 3
##  $ y       : num [1:15, 1:32, 1:36] 1 0 0 0 0 0 0 0 0 0 ...
##   ..- attr(*, "dimnames")=List of 3
##   .. ..$ species: chr [1:15] "Odocoileus_virginianus" "Canis_latrans" "Sciurus_niger" "Procyon_lotor" ...
##   .. ..$ site   : chr [1:32] "BI201" "BN211" "EI100" "EI102" ...
##   .. ..$ week   : chr [1:36] "2024-29" "2024-30" "2024-31" "2024-32" ...
##  $ occ.covs:'data.frame':    32 obs. of  12 variables:
##   ..$ Cogon_Patch_Size      : num [1:32] -0.237 -0.446 -0.337 -0.308 0.039 ...
##   ..$ VegetationDiversity   : num [1:32] -0.273 0.455 1.619 -0.273 2.929 ...
##   ..$ PostTreatmentDensities: num [1:32] 0.432 -0.729 0.432 2.169 1.13 ...
##   ..$ Auth                  : num [1:32] -2.28 -2.28 -1.04 -1.04 -1.04 ...
##   ..$ UTM_Zone              : num [1:32] -0.473 -0.473 -0.473 -0.473 -0.473 ...
##   ..$ CWD_presence          : num [1:32] 1.86 -0.521 -0.521 -0.521 -0.521 ...
##   ..$ Veg_shannon_index     : num [1:32] -0.143 0.359 0.604 -1.683 0.328 ...
##   ..$ Avg_Cogongrass_Cover  : num [1:32] -0.154 -0.708 0.308 2.045 1.121 ...
##   ..$ total_shrub_cover     : num [1:32] 1.526 0.5 -0.153 -1.003 -0.811 ...
##   ..$ avg_veg_height        : num [1:32] 0.466 -1.044 1.181 0.879 1.684 ...
##   ..$ Tree_Density          : num [1:32] -0.3629 -0.3564 -0.5111 3.5896 0.0958 ...
##   ..$ Avg_Canopy_Cover      : num [1:32] 0.1362 -0.0252 -0.9132 0.782 -1.9627 ...
##  $ det.covs:'data.frame':    32 obs. of  12 variables:
##   ..$ Cogon_Patch_Size      : num [1:32] -0.237 -0.446 -0.337 -0.308 0.039 ...
##   ..$ VegetationDiversity   : num [1:32] -0.273 0.455 1.619 -0.273 2.929 ...
##   ..$ PostTreatmentDensities: num [1:32] 0.432 -0.729 0.432 2.169 1.13 ...
##   ..$ Auth                  : num [1:32] -2.28 -2.28 -1.04 -1.04 -1.04 ...
##   ..$ UTM_Zone              : num [1:32] -0.473 -0.473 -0.473 -0.473 -0.473 ...
##   ..$ CWD_presence          : num [1:32] 1.86 -0.521 -0.521 -0.521 -0.521 ...
##   ..$ Veg_shannon_index     : num [1:32] -0.143 0.359 0.604 -1.683 0.328 ...
##   ..$ Avg_Cogongrass_Cover  : num [1:32] -0.154 -0.708 0.308 2.045 1.121 ...
##   ..$ total_shrub_cover     : num [1:32] 1.526 0.5 -0.153 -1.003 -0.811 ...
##   ..$ avg_veg_height        : num [1:32] 0.466 -1.044 1.181 0.879 1.684 ...
##   ..$ Tree_Density          : num [1:32] -0.3629 -0.3564 -0.5111 3.5896 0.0958 ...
##   ..$ Avg_Canopy_Cover      : num [1:32] 0.1362 -0.0252 -0.9132 0.782 -1.9627 ...

Objective 3: Assessing wildlife diversity and how it relates to site conditions

Wildlife Shannon Diversity Index

Wildlife Shannon Diversity Index

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

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

## take CWD presence data and merge to Camera Loc
cwd_pa <- cwd_data %>%
  group_by(Plot) %>%
  summarise(
    CWD_presence = as.integer(any(Line_Start > 0 & Line_End > 0)),
    .groups = "drop"
  )

CameraLoc <- CameraLoc %>%
  left_join(cwd_pa, by = "Plot")

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

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

# Observation per species, per site
site_species_counts <- CameraData %>%
  group_by(Plot, Name) %>%
  summarize(count = n(), .groups = "drop")

# Shannon Diversity per site

shannon_diversity_wildlife <- site_species_counts %>%
  group_by(Plot) %>%
  mutate(proportion = count / sum(count)) %>%
  summarize(Wild_shannon_index = -sum(proportion * log(proportion), na.rm = TRUE))

# View results
print(shannon_diversity_wildlife)

## # A tibble: 32 × 2
##    Plot  Wild_shannon_index
##    <chr>              <dbl>
##  1 BI201              0.305
##  2 BN211              0.303
##  3 EI100              0.292
##  4 EI102              0.412
##  5 EI104              0.310
##  6 EI106              0.446
##  7 EN101              0.117
##  8 EN103              0.194
##  9 EN107              0.861
## 10 JI07               0.764
## # ℹ 22 more rows

Merge Data

# Merge wildlife and vegetation diversity indices with invasion data
diversity_data <- shannon_diversity_wildlife %>%
  left_join(Veg_shannon_diversity, by = "Plot") %>%
  left_join(CameraLoc, by = "Plot")

# View merged data
print(diversity_data)

## # A tibble: 32 × 16
##    Plot  Wild_shannon_index Veg_shannon_index   Lat  Long Status     
##    <chr>              <dbl>             <dbl> <dbl> <dbl> <chr>      
##  1 BI201              0.305              2.20  30.8 -86.9 Invaded    
##  2 BN211              0.303              2.43  30.8 -86.9 Non_Invaded
##  3 EI100              0.292              2.54  31.0 -87.1 Invaded    
##  4 EI102              0.412              1.49  31.0 -87.1 Invaded    
##  5 EI104              0.310              2.42  31.0 -87.1 Invaded    
##  6 EI106              0.446              1.73  31.0 -87.1 Invaded    
##  7 EN101              0.117              2.07  31.0 -87.1 Non_Invaded
##  8 EN103              0.194              2.43  31.0 -87.1 Non_Invaded
##  9 EN107              0.861              2.69  31.0 -87.1 Non_Invaded
## 10 JI07               0.764              1.69  30.8 -87.1 Invaded    
## # ℹ 22 more rows
## # ℹ 10 more variables: Start_Date <dttm>, Camera <chr>, Cogon_Patch_Size <dbl>,
## #   VegetationDiversity <dbl>, PostTreatmentDensities <dbl>, Authority <chr>,
## #   Auth <dbl>, UTM_Zone <dbl>, CWD_presence <int>, Avg_Cogongrass_Cover <dbl>

summary(diversity_data)

##      Plot           Wild_shannon_index Veg_shannon_index      Lat       
##  Length:32          Min.   :0.04539    Min.   :1.188     Min.   :28.68  
##  Class :character   1st Qu.:0.29022    1st Qu.:2.024     1st Qu.:30.76  
##  Mode  :character   Median :0.44221    Median :2.369     Median :30.77  
##                     Mean   :0.56021    Mean   :2.264     Mean   :30.45  
##                     3rd Qu.:0.75653    3rd Qu.:2.522     3rd Qu.:30.90  
##                     Max.   :1.72721    Max.   :3.160     Max.   :31.01  
##       Long           Status            Start_Date                 
##  Min.   :-87.15   Length:32          Min.   :2024-05-09 00:00:00  
##  1st Qu.:-87.14   Class :character   1st Qu.:2024-06-03 18:00:00  
##  Median :-87.09   Mode  :character   Median :2024-06-06 12:00:00  
##  Mean   :-86.21                      Mean   :2024-06-26 04:30:00  
##  3rd Qu.:-86.89                      3rd Qu.:2024-06-27 00:00:00  
##  Max.   :-82.41                      Max.   :2024-11-04 00:00:00  
##     Camera          Cogon_Patch_Size  VegetationDiversity
##  Length:32          Min.   :   0.00   Min.   :11.00      
##  Class :character   1st Qu.:   0.00   1st Qu.:18.00      
##  Mode  :character   Median :  30.07   Median :20.50      
##                     Mean   : 458.39   Mean   :21.88      
##                     3rd Qu.: 233.84   3rd Qu.:25.00      
##                     Max.   :4168.92   Max.   :42.00      
##  PostTreatmentDensities  Authority              Auth          UTM_Zone    
##  Min.   :0.0000         Length:32          Min.   :1.000   Min.   :16.00  
##  1st Qu.:0.0000         Class :character   1st Qu.:2.000   1st Qu.:16.00  
##  Median :0.0000         Mode  :character   Median :3.000   Median :16.00  
##  Mean   :0.8978                            Mean   :2.844   Mean   :16.19  
##  3rd Qu.:1.4650                            3rd Qu.:3.000   3rd Qu.:16.00  
##  Max.   :3.7100                            Max.   :4.000   Max.   :17.00  
##   CWD_presence    Avg_Cogongrass_Cover
##  Min.   :0.0000   Min.   : 0.000      
##  1st Qu.:0.0000   1st Qu.: 0.000      
##  Median :0.0000   Median : 2.143      
##  Mean   :0.2188   Mean   :13.683      
##  3rd Qu.:0.0000   3rd Qu.:20.089      
##  Max.   :1.0000   Max.   :63.571

diversity_data <- diversity_data %>%
  mutate(Status = as.factor(Status))

# Standardizing the continuous variables
diversity_data <- diversity_data %>%
  mutate(
    Avg_Cogongrass_Cover = scale(Avg_Cogongrass_Cover),
    Veg_shannon_index = scale(Veg_shannon_index),
    Cogon_Patch_Size = scale(Cogon_Patch_Size) 
  )

diversity_data <- diversity_data %>%
  mutate(
    CWD_presence = as.factor(CWD_presence)
  )

Diversity Model

glmm_model <- lmer(
  Wild_shannon_index ~ Avg_Cogongrass_Cover * Veg_shannon_index + Cogon_Patch_Size + CWD_presence +
    (1 | Authority) + (1 | Camera),
  data = diversity_data
)

## boundary (singular) fit: see help('isSingular')

# Summarize the model
summary(glmm_model)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Wild_shannon_index ~ Avg_Cogongrass_Cover * Veg_shannon_index +  
##     Cogon_Patch_Size + CWD_presence + (1 | Authority) + (1 |      Camera)
##    Data: diversity_data
## 
## REML criterion at convergence: 40.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7111 -0.5786 -0.1424  0.3584  2.6998 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Authority (Intercept) 0.0000   0.0000  
##  Camera    (Intercept) 0.0229   0.1513  
##  Residual              0.1335   0.3654  
## Number of obs: 32, groups:  Authority, 5; Camera, 4
## 
## Fixed effects:
##                                         Estimate Std. Error        df t value
## (Intercept)                             0.458789   0.123390  5.039200   3.718
## Avg_Cogongrass_Cover                   -0.171326   0.159880 25.767254  -1.072
## Veg_shannon_index                       0.009569   0.111461 24.469771   0.086
## Cogon_Patch_Size                        0.002448   0.085673 25.069391   0.029
## CWD_presence1                          -0.064466   0.184902 25.963759  -0.349
## Avg_Cogongrass_Cover:Veg_shannon_index -0.218582   0.081616 25.988707  -2.678
##                                        Pr(>|t|)  
## (Intercept)                              0.0135 *
## Avg_Cogongrass_Cover                     0.2938  
## Veg_shannon_index                        0.9323  
## Cogon_Patch_Size                         0.9774  
## CWD_presence1                            0.7302  
## Avg_Cogongrass_Cover:Veg_shannon_index   0.0127 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Av_C_C Vg_sh_ Cg_P_S CWD_p1
## Avg_Cgngr_C  0.231                            
## Vg_shnnn_nd -0.036  0.613                     
## Cgn_Ptch_Sz -0.055 -0.494 -0.077              
## CWD_presnc1 -0.248  0.160  0.071 -0.349       
## Avg_C_C:V__  0.455  0.662  0.067 -0.352 -0.063
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

Plot the model results

# Extract the fixed effects
fixed_effects <- fixef(glmm_model)

# Extract the random effects
random_effects <- ranef(glmm_model)

# Plot the fixed effects
plot_model(glmm_model, type = "pred", terms = c("Avg_Cogongrass_Cover", "Veg_shannon_index"))

# Plot the random effects
plot_model(glmm_model, type = "re")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the sjPlot package.
##   Please report the issue at <https://github.com/strengejacke/sjPlot/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## [[1]]

## 
## [[2]]

Bayesian Diversity Model

bayesian_model <- brm(
  Wild_shannon_index ~ scale(Avg_Cogongrass_Cover) * scale(Veg_shannon_index) + scale(Cogon_Patch_Size) +
    CWD_presence +  (1 | Authority) + (1 | Camera),
  data = diversity_data,
  family = gaussian(),
  chains = 4,          
  iter = 5000,          
  warmup = 3000,       
  cores = parallel::detectCores(),  
  control = list(adapt_delta = 0.99, max_treedepth = 15),  # Further tuned parameters
  prior = c(
    prior(normal(0, 1), class = "b"),       # Regularizing prior for coefficients
    prior(normal(0, 1), class = "Intercept"),  # Weakly informative cauchy prior for intercept
    prior(normal(0, 1), class = "sd")       # For random effect standard deviations
  )
)

## Compiling Stan program...

## Start sampling

summary(bayesian_model)

##  Family: gaussian 
##   Links: mu = identity 
## Formula: Wild_shannon_index ~ scale(Avg_Cogongrass_Cover) * scale(Veg_shannon_index) + scale(Cogon_Patch_Size) + CWD_presence + (1 | Authority) + (1 | Camera) 
##    Data: diversity_data (Number of observations: 32) 
##   Draws: 4 chains, each with iter = 5000; warmup = 3000; thin = 1;
##          total post-warmup draws = 8000
## 
## Multilevel Hyperparameters:
## ~Authority (Number of levels: 5) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.16      0.16     0.01     0.59 1.00     3574     3497
## 
## ~Camera (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.27      0.23     0.01     0.90 1.00     2547     3487
## 
## Regression Coefficients:
##                                                  Estimate Est.Error l-95% CI
## Intercept                                            0.42      0.22    -0.04
## scaleAvg_Cogongrass_Cover                           -0.14      0.18    -0.49
## scaleVeg_shannon_index                               0.02      0.13    -0.23
## scaleCogon_Patch_Size                               -0.02      0.10    -0.21
## CWD_presence1                                       -0.05      0.20    -0.46
## scaleAvg_Cogongrass_Cover:scaleVeg_shannon_index    -0.21      0.09    -0.40
##                                                  u-95% CI Rhat Bulk_ESS
## Intercept                                            0.86 1.00     4050
## scaleAvg_Cogongrass_Cover                            0.22 1.00     4862
## scaleVeg_shannon_index                               0.27 1.00     6874
## scaleCogon_Patch_Size                                0.18 1.00     5879
## CWD_presence1                                        0.35 1.00     5690
## scaleAvg_Cogongrass_Cover:scaleVeg_shannon_index    -0.02 1.00     5241
##                                                  Tail_ESS
## Intercept                                            3582
## scaleAvg_Cogongrass_Cover                            5317
## scaleVeg_shannon_index                               5334
## scaleCogon_Patch_Size                                5260
## CWD_presence1                                        5192
## scaleAvg_Cogongrass_Cover:scaleVeg_shannon_index     5484
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.39      0.06     0.29     0.53 1.00     5808     5330
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

plot(bayesian_model)

# Plot the fixed effects
plot(bayesian_model, "Avg_Cogongrass_Cover")

## Warning: Argument 'pars' is deprecated. Please use 'variable' instead.

# Plot the random effects
plot(bayesian_model, "Camera")

## Warning: Argument 'pars' is deprecated. Please use 'variable' instead.

plot(bayesian_model, "Authority")

## Warning: Argument 'pars' is deprecated. Please use 'variable' instead.

# Get unique levels of Authority and Camera from original data
authority_levels <- levels(as.factor(diversity_data$Authority))
camera_levels <- levels(as.factor(diversity_data$Camera))

# Create a new data frame for predictions
new_data <- expand.grid(
  Avg_Cogongrass_Cover = seq(
    min(diversity_data$Avg_Cogongrass_Cover),
    max(diversity_data$Avg_Cogongrass_Cover),
    length.out = 100
  ),
  Veg_shannon_index = quantile(
    diversity_data$Veg_shannon_index,
    probs = c(0.25, 0.5, 0.75)
  ),
  Cogon_Patch_Size = mean(diversity_data$Cogon_Patch_Size),
  CWD_presence = factor(c(0, 1)),
  Authority = NA,
  Camera = NA
)

predictions <- fitted(
  bayesian_model,
  newdata = new_data,
  re_formula = NA
)

new_data$fit   <- predictions[, "Estimate"]
new_data$lwr   <- predictions[, "Q2.5"]
new_data$upr   <- predictions[, "Q97.5"]



ggplot(
  new_data,
  aes(
    x = Avg_Cogongrass_Cover,
    y = fit,
    color = factor(Veg_shannon_index),
    fill = factor(Veg_shannon_index),
    linetype = CWD_presence
  )
) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.15, color = NA) +
  labs(
    x = "Average Cogongrass Cover (%)",
    y = "Predicted Wildlife Shannon Index",
    color = "Vegetation Shannon\n(quantiles)",
    fill = "Vegetation Shannon\n(quantiles)",
    linetype = "CWD Presence"
  ) +
  theme_minimal(base_size = 13)