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)

# Tree PCQ Data
tree_data <- read_excel("C:/Users/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Tree_PCQ")

# Soil Data
fuel_data <- read_excel("C:/Users/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Fuel_Sampling")

# Veg Data
Veg_Cover <- read_excel("C:/Users/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Veg_Cover")

# Shrub Cover Data
shrub_data <- read_excel("C:/Users/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                         sheet = "Shrub_Cover")

# Site Data
CameraData <- read_excel("C:/Users/Kody Melissa/Downloads/CameraData.xlsx")

CameraLoc <- read_excel("C:/Users/Kody Melissa/Downloads/CameraLoc.xlsx",
                  sheet = "CameraLocations")

# First day of the study
study_start <- as.Date("2024-05-10")

# Effort Matrix
effort_matrix <- read_excel("C:/Users/Kody Melissa/Downloads/CameraLoc.xlsx",
                            sheet = "Daily_EM")

col_names <- names(effort_matrix)
date_cols <- col_names[-1]
converted_dates <- as.character(as.Date(as.numeric(date_cols), origin = "1899-12-30"))
names(effort_matrix) <- c("Plot", converted_dates)

# Pivot, filter, and add custom study week
effort_matrix <- effort_matrix %>%
  pivot_longer(
    cols = -Plot,
    names_to = "date",
    values_to = "effort"
  ) %>%
  filter(effort == 1) %>%
  mutate(
    date = as.Date(date),
    study_week = as.integer(floor((date - study_start) / 7)) + 1  # Week 1 starts on study_start
  ) %>%
  dplyr::select(Plot, date, study_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

Normal Turkey

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/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Tree_PCQ")

# Soil Data
fuel_data <- read_excel("C:/Users/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Fuel_Sampling")

# Veg Data
Veg_Cover <- read_excel("C:/Users/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Veg_Cover")

# Shrub Cover Data
shrub_data <- read_excel("C:/Users/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                         sheet = "Shrub_Cover")

# Site Data
CameraData <- read_excel("C:/Users/Kody Melissa/Downloads/CameraData.xlsx")

CameraLoc <- read_excel("C:/Users/Kody Melissa/Downloads/CameraLoc.xlsx",
                  sheet = "CameraLocations")

# First day of the study
study_start <- as.Date("2024-05-10")

# Effort Matrix
effort_matrix <- read_excel("C:/Users/Kody Melissa/Downloads/CameraLoc.xlsx",
                            sheet = "Daily_EM")

col_names <- names(effort_matrix)
date_cols <- col_names[-1]
converted_dates <- as.character(as.Date(as.numeric(date_cols), origin = "1899-12-30"))
names(effort_matrix) <- c("Plot", converted_dates)

# Pivot, filter, and add custom study week
effort_matrix <- effort_matrix %>%
  pivot_longer(
    cols = -Plot,
    names_to = "date",
    values_to = "effort"
  ) %>%
  filter(effort == 1) %>%
  mutate(
    date = as.Date(date),
    study_week = as.integer(floor((date - study_start) / 7)) + 1  # Week 1 starts on study_start
  ) %>%
  dplyr::select(Plot, date, study_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

Turkey- Minimum 10 observations

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/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Tree_PCQ")

# Soil Data
fuel_data <- read_excel("C:/Users/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Fuel_Sampling")

# Veg Data
Veg_Cover <- read_excel("C:/Users/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Veg_Cover")

# Shrub Cover Data
shrub_data <- read_excel("C:/Users/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                         sheet = "Shrub_Cover")

# Site Data
CameraData <- read_excel("C:/Users/Kody Melissa/Downloads/CameraData.xlsx")

CameraLoc <- read_excel("C:/Users/Kody Melissa/Downloads/CameraLoc.xlsx",
                  sheet = "CameraLocations")

# First day of the study
study_start <- as.Date("2024-05-10")

# Effort Matrix
effort_matrix <- read_excel("C:/Users/Kody Melissa/Downloads/CameraLoc.xlsx",
                            sheet = "Daily_EM")

col_names <- names(effort_matrix)
date_cols <- col_names[-1]
converted_dates <- as.character(as.Date(as.numeric(date_cols), origin = "1899-12-30"))
names(effort_matrix) <- c("Plot", converted_dates)

# Pivot, filter, and add custom study week
effort_matrix <- effort_matrix %>%
  pivot_longer(
    cols = -Plot,
    names_to = "date",
    values_to = "effort"
  ) %>%
  filter(effort == 1) %>%
  mutate(
    date = as.Date(date),
    study_week = as.integer(floor((date - study_start) / 7)) + 1  # Week 1 starts on study_start
  ) %>%
  dplyr::select(Plot, date, study_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

Turkey- Minimum 25 observations

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/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Tree_PCQ")

# Soil Data
fuel_data <- read_excel("C:/Users/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Fuel_Sampling")

# Veg Data
Veg_Cover <- read_excel("C:/Users/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Veg_Cover")

# Shrub Cover Data
shrub_data <- read_excel("C:/Users/Kody Melissa/Downloads/Field_Data_FL_AL_MS.xlsx",
                         sheet = "Shrub_Cover")

# Site Data
CameraData <- read_excel("C:/Users/Kody Melissa/Downloads/CameraData.xlsx")

CameraLoc <- read_excel("C:/Users/Kody Melissa/Downloads/CameraLoc.xlsx",
                  sheet = "CameraLocations")

# First day of the study
study_start <- as.Date("2024-05-10")

# Effort Matrix
effort_matrix <- read_excel("C:/Users/Kody Melissa/Downloads/CameraLoc.xlsx",
                            sheet = "Daily_EM")

col_names <- names(effort_matrix)
date_cols <- col_names[-1]
converted_dates <- as.character(as.Date(as.numeric(date_cols), origin = "1899-12-30"))
names(effort_matrix) <- c("Plot", converted_dates)

# Pivot, filter, and add custom study week
effort_matrix <- effort_matrix %>%
  pivot_longer(
    cols = -Plot,
    names_to = "date",
    values_to = "effort"
  ) %>%
  filter(effort == 1) %>%
  mutate(
    date = as.Date(date),
    study_week = as.integer(floor((date - study_start) / 7)) + 1  # Week 1 starts on study_start
  ) %>%
  dplyr::select(Plot, date, study_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