Normal

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)

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

# 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

This is not needed for non-ordination analysis. Moving the threshold down to 0% to keep the option, but to ensure it has no effect for now.

# 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 >= 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.
## 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

Plot Species level relationship with cogongrass cover

species_index <- which(dimnames(data_list$y)$species == "Procyon_lotor")
# Create a set of values across the range of observed cogongrass values
cogon.pred.vals <- seq(min(data_list$occ.covs$Avg_Cogongrass_Cover),
                       max(data_list$occ.covs$Avg_Cogongrass_Cover),
                       length.out = 100)

# Scale predicted values by mean and standard deviation used to fit the model
cogon.pred.vals.scale <- (cogon.pred.vals - mean(data_list$occ.covs$Avg_Cogongrass_Cover)) /
  sd(data_list$occ.covs$Avg_Cogongrass_Cover)

# Predict occupancy across cogongrass cover values at mean values of all other variables
pred.df <- as.matrix(data.frame(intercept = 1, Avg_Cogongrass_Cover =
                                cogon.pred.vals.scale, 'I(Avg_Cogongrass_Cover^2)' = 0,
                                Cogon_Patch_Size = 0, Veg_shannon_index = 0,
                                total_shrub_cover = 0, Tree_Density = 0,
                                Avg_Canopy_Cover = 0, avg_veg_height = 0, Auth = 0))

out.pred <- predict(ms_fullQ_fullQ, pred.df, species = species_index)
## Warning in predict.msPGOcc(ms_fullQ_fullQ, pred.df, species = species_index):
## 'species' is not an argument
str(out.pred)
## List of 3
##  $ psi.0.samples: num [1:4000, 1:8, 1:100] 0.4047 0.2256 0.297 0.0707 0.2208 ...
##  $ z.0.samples  : int [1:4000, 1:8, 1:100] 1 0 0 0 0 1 0 1 1 1 ...
##  $ call         : language predict.msPGOcc(object = ms_fullQ_fullQ, X.0 = pred.df, species = species_index)
##  - attr(*, "class")= chr "predict.msPGOcc"
str(out.pred$psi.0.samples)
##  num [1:4000, 1:8, 1:100] 0.4047 0.2256 0.297 0.0707 0.2208 ...
psi.0.quants <- apply(out.pred$psi.0.samples, c(3), function(x) quantile(x, prob = c(0.025, 0.5, 0.975)))

str(psi.0.quants)
##  num [1:3, 1:100] 0.0017 0.12 0.8617 0.0017 0.1196 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:3] "2.5%" "50%" "97.5%"
##   ..$ : NULL
psi.plot.dat <- data.frame(
  psi.med = psi.0.quants[2, ],
  psi.low = psi.0.quants[1, ],
  psi.high = psi.0.quants[3, ],
  Avg_Cogongrass_Cover = cogon.pred.vals
)

str(psi.plot.dat)
## 'data.frame':    100 obs. of  4 variables:
##  $ psi.med             : num  0.12 0.12 0.123 0.12 0.121 ...
##  $ psi.low             : num  0.0017 0.0017 0.00195 0.00194 0.00212 ...
##  $ psi.high            : num  0.862 0.848 0.843 0.839 0.839 ...
##  $ Avg_Cogongrass_Cover: num  -0.708 -0.675 -0.641 -0.608 -0.575 ...
ggplot(psi.plot.dat, aes(x = Avg_Cogongrass_Cover, y = psi.med)) +
  geom_ribbon(aes(ymin = psi.low, ymax = psi.high), fill = "grey70") +
  geom_line() +
  theme_bw() +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = "Average Cogongrass Cover", y = "Occupancy Probability") +
  ggtitle("Procyon_lotor Occupancy vs Cogongrass Cover")