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

# 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: 174 × 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
## # ℹ 164 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.2836106

Objective 2: Determine whether wildlife behaviors differ in and around cogongrass patches, specifically as it relates to whether specific taxa avoid cogongrass invaded areas.

Behavior Occurence

CameraData <- CameraData%>%
  dplyr::select(-Status)

O2_data <- CameraData %>%
  left_join(CameraLoc_O2, by = "Plot")

# Creating proportions of observations for each behavior type
behavior_proportions <- O2_data %>%
  group_by(Plot, Name, Behavior, Status, Camera.Type, BehLoc) %>%
  summarise(ObservationCount = n()) %>%
  ungroup() %>%
  group_by(Plot, Name, Status, Camera.Type, BehLoc) %>%
  mutate(Proportion = ObservationCount / sum(ObservationCount)) %>%
  ungroup()
## `summarise()` has grouped output by 'Plot', 'Name', 'Behavior', 'Status',
## 'Camera.Type'. You can override using the `.groups` argument.

Behavior Model

This model is out of data, compared to the one below, because I wanted local search to be the reference variable. That is because local search was a catch-all for non-foraging or transit behaviors.

#formula <- bf(Proportion ~ Status * Behavior + Status * Name + (1 | Plot),
#              family = Beta(link = "logit"))

# Adjust Proportion for Beta distribution
#n <- nrow(behavior_proportions)
#behavior_proportions$Proportion <- (behavior_proportions$Proportion * (n - 1) + 0.5) / n

# Fit the model
#brms_model <- brm(formula = formula,
#                  data = behavior_proportions,
#                  iter = 4000, warmup = 1000, chains = 4, cores = 4,
#                  control = list(adapt_delta = 0.95))

# Diagnostic Plots
#plot(brms_model)

# Summary Statistics
#summary(brms_model)

Behavior Model- Switched to local search and dasypus novemcinctus being the reference variables

# Ensure Behavior is an unordered factor
if(!is.factor(behavior_proportions$Behavior) || is.ordered(behavior_proportions$Behavior)) {
  behavior_proportions$Behavior <- factor(behavior_proportions$Behavior, ordered = FALSE)
}

# Make Local_Search the reference category
behavior_proportions$Behavior <- relevel(behavior_proportions$Behavior, ref = "Local_Search")

# Ensure Name is an unordered factor
if(!is.factor(behavior_proportions$Name) || is.ordered(behavior_proportions$Name)) {
  behavior_proportions$Name <- factor(behavior_proportions$Name, ordered = FALSE)
}

# Make Dasypus_novemcinctus the reference category
behavior_proportions$Name <- relevel(behavior_proportions$Name, ref = "Dasypus_novemcinctus")

# Ensure Status is an unordered factor
if(!is.factor(behavior_proportions$Status) || is.ordered(behavior_proportions$Status)) {
  behavior_proportions$Status <- factor(behavior_proportions$Status, ordered = FALSE)
}

# Make "Non_Invaded" the reference category
behavior_proportions$Status <- relevel(behavior_proportions$Status, ref = "Non_Invaded")

# Adjust Proportion for Beta distribution
n <- nrow(behavior_proportions)
behavior_proportions$Proportion <- (behavior_proportions$Proportion * (n - 1) + 0.5) / n

# Fit the Bayesian GLMM with new reference levels
brms_model_releveled <- brm(
  formula = bf(Proportion ~ Status * Behavior + Status * Name + (1 | Plot), 
               family = Beta(link = "logit")),
  data = behavior_proportions,
  iter = 10000, warmup = 3000, chains = 4, cores = 4,
  control = list(adapt_delta = 0.95)
)
## Compiling Stan program...
## Start sampling
plot(brms_model_releveled)

# Summary
summary(brms_model_releveled)
##  Family: beta 
##   Links: mu = logit; phi = identity 
## Formula: Proportion ~ Status * Behavior + Status * Name + (1 | Plot) 
##    Data: behavior_proportions (Number of observations: 323) 
##   Draws: 4 chains, each with iter = 10000; warmup = 3000; thin = 1;
##          total post-warmup draws = 28000
## 
## Multilevel Hyperparameters:
## ~Plot (Number of levels: 32) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.19      0.11     0.01     0.42 1.00     6538     8122
## 
## Regression Coefficients:
##                                          Estimate Est.Error l-95% CI u-95% CI
## Intercept                                   -0.02      0.43    -0.84     0.84
## StatusInvaded                                0.02      0.52    -1.00     1.04
## BehaviorForaging                             0.58      0.30     0.00     1.16
## BehaviorTransit                              1.32      0.27     0.80     1.84
## NameCanis_latrans                            0.17      0.47    -0.77     1.08
## NameDidelphis_virginiana                     0.67      0.70    -0.66     2.10
## NameLynx_rufus                               1.14      0.63    -0.06     2.39
## NameOdocoileus_virginianus                  -1.02      0.40    -1.84    -0.26
## NameProcyon_lotor                            0.37      0.47    -0.56     1.27
## NameSciurus_carolinensis                    -0.15      0.65    -1.43     1.13
## NameSciurus_niger                            1.50      0.78     0.05     3.11
## NameSus_scrofa                               0.12      0.63    -1.09     1.38
## NameSylvilagus_floridanus                    1.11      0.53     0.08     2.16
## NameVulpes_vulpes                            1.28      0.92    -0.37     3.27
## StatusInvaded:BehaviorForaging              -0.70      0.38    -1.44     0.05
## StatusInvaded:BehaviorTransit                0.03      0.33    -0.61     0.68
## StatusInvaded:NameCanis_latrans              1.03      0.63    -0.19     2.27
## StatusInvaded:NameDidelphis_virginiana       0.02      0.82    -1.62     1.60
## StatusInvaded:NameLynx_rufus                 0.11      0.79    -1.46     1.66
## StatusInvaded:NameOdocoileus_virginianus     0.38      0.48    -0.55     1.33
## StatusInvaded:NameProcyon_lotor              0.13      0.56    -0.98     1.25
## StatusInvaded:NameSciurus_carolinensis       0.11      0.78    -1.44     1.65
## StatusInvaded:NameSciurus_niger             -0.99      1.01    -3.00     0.95
## StatusInvaded:NameSus_scrofa                 1.63      1.54    -0.96     5.15
## StatusInvaded:NameSylvilagus_floridanus      0.14      0.69    -1.21     1.48
## StatusInvaded:NameVulpes_vulpes              0.27      1.65    -2.70     3.84
##                                          Rhat Bulk_ESS Tail_ESS
## Intercept                                1.00     4600     7845
## StatusInvaded                            1.00     4603     7952
## BehaviorForaging                         1.00    13454    19291
## BehaviorTransit                          1.00    14608    19504
## NameCanis_latrans                        1.00     5536    10161
## NameDidelphis_virginiana                 1.00     8408    14518
## NameLynx_rufus                           1.00     7918    13089
## NameOdocoileus_virginianus               1.00     4594     8458
## NameProcyon_lotor                        1.00     5453     8953
## NameSciurus_carolinensis                 1.00     7855    13517
## NameSciurus_niger                        1.00     9604    16223
## NameSus_scrofa                           1.00     8459    14427
## NameSylvilagus_floridanus                1.00     7046    11367
## NameVulpes_vulpes                        1.00    13508    17550
## StatusInvaded:BehaviorForaging           1.00    13588    18960
## StatusInvaded:BehaviorTransit            1.00    14301    19628
## StatusInvaded:NameCanis_latrans          1.00     6451    11972
## StatusInvaded:NameDidelphis_virginiana   1.00     7842    12898
## StatusInvaded:NameLynx_rufus             1.00     8387    14058
## StatusInvaded:NameOdocoileus_virginianus 1.00     4706     8169
## StatusInvaded:NameProcyon_lotor          1.00     5543     9369
## StatusInvaded:NameSciurus_carolinensis   1.00     8359    14687
## StatusInvaded:NameSciurus_niger          1.00    10874    16076
## StatusInvaded:NameSus_scrofa             1.00    16620    15825
## StatusInvaded:NameSylvilagus_floridanus  1.00     7699    12175
## StatusInvaded:NameVulpes_vulpes          1.00    17888    17130
## 
## Further Distributional Parameters:
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi     1.89      0.15     1.61     2.20 1.00    23309    21410
## 
## 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).
# Conditional Effects for interactions
ce_behavior <- conditional_effects(brms_model_releveled, effects = "Status:Behavior")
plot(ce_behavior)

ce_name <- conditional_effects(brms_model_releveled, effects = "Status:Name")
plot(ce_name)

Behavior Model Within cogongrass patches

# Filter data for invaded status
invaded_data <- behavior_proportions %>% 
  filter(Status == "Invaded")

# Adjust Proportion for Beta distribution
n <- nrow(invaded_data)
invaded_data$Proportion <- (invaded_data$Proportion * (n - 1) + 0.5) / n

Loc_model <- brm(
  Proportion ~ BehLoc * Behavior + BehLoc * Name + (1 | Plot) + (1 | Camera.Type), 
  data = invaded_data, 
  family = Beta(link = "logit"), 
  iter = 10000,
  warmup = 3000,
  chains = 4,
  cores = 4,
  control = list(adapt_delta = 0.99, max_treedepth = 15)
)
## Compiling Stan program...
## Start sampling
## Warning: There were 1 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 27999 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 15. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: The largest R-hat is 2.99, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
# Summarize model
summary(Loc_model)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 1 divergent transitions after warmup. Increasing
## adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: beta 
##   Links: mu = logit; phi = identity 
## Formula: Proportion ~ BehLoc * Behavior + BehLoc * Name + (1 | Plot) + (1 | Camera.Type) 
##    Data: invaded_data (Number of observations: 193) 
##   Draws: 4 chains, each with iter = 10000; warmup = 3000; thin = 1;
##          total post-warmup draws = 28000
## 
## Multilevel Hyperparameters:
## ~Camera.Type (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.61      0.33     0.15     1.36 1.18       18       75
## 
## ~Plot (Number of levels: 17) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.43      0.20     0.13     0.89 1.27       13       40
## 
## Regression Coefficients:
##                                            Estimate    Est.Error     l-95% CI
## Intercept                                     -0.02         0.41        -0.75
## BehLocPatch                                   -0.35         0.53        -1.27
## BehaviorForaging                               0.23         0.28        -0.31
## BehaviorTransit                                0.66         0.27         0.03
## NameCanis_latrans                              1.48         0.48         0.52
## NameDidelphis_virginiana                       0.82         0.46        -0.08
## NameLynx_rufus                                 1.89         0.47         1.18
## NameOdocoileus_virginianus                    -0.78         0.28        -1.28
## NameProcyon_lotor                              0.54         0.34        -0.08
## NameSciurus_carolinensis                       0.05         0.37        -0.69
## NameSciurus_niger                              0.45         0.68        -0.95
## NameSus_scrofa                           -163517.21    411291.23   -972432.68
## NameSylvilagus_floridanus                      1.34         0.42         0.47
## NameVulpes_vulpes                              2.84         1.61         0.38
## BehLocPatch:BehaviorForaging                  -0.98         0.54        -2.11
## BehLocPatch:BehaviorTransit                    2.08         0.48         1.26
## BehLocPatch:NameCanis_latrans                 -1.11         0.71        -2.54
## BehLocPatch:NameDidelphis_virginiana          -0.55         0.84        -2.54
## BehLocPatch:NameLynx_rufus                    -1.23         1.33        -3.50
## BehLocPatch:NameOdocoileus_virginianus        -0.01         0.50        -0.97
## BehLocPatch:NameProcyon_lotor                 -0.52         0.67        -1.89
## BehLocPatch:NameSciurus_carolinensis           1.26         1.43        -1.25
## BehLocPatch:NameSciurus_niger                  0.75         1.99        -2.64
## BehLocPatch:NameSus_scrofa                163518.52    411291.16   -521882.80
## BehLocPatch:NameSylvilagus_floridanus         -0.79         0.80        -2.34
## BehLocPatch:NameVulpes_vulpes          221498769.57 416566857.87 -49104752.71
##                                             u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept                                       0.82 1.10       29      111
## BehLocPatch                                     0.66 1.22       14       64
## BehaviorForaging                                0.76 1.16       19       46
## BehaviorTransit                                 1.12 1.31       11       15
## NameCanis_latrans                               2.41 1.09       43      153
## NameDidelphis_virginiana                        1.78 1.07       48       44
## NameLynx_rufus                                  2.95 1.05       78       81
## NameOdocoileus_virginianus                     -0.22 1.04       47      160
## NameProcyon_lotor                               1.19 1.12       24      106
## NameSciurus_carolinensis                        0.75 1.10       65      138
## NameSciurus_niger                               1.73 1.08       45       99
## NameSus_scrofa                             521882.22 2.88        5       11
## NameSylvilagus_floridanus                       2.17 1.09       44       93
## NameVulpes_vulpes                               6.83 1.33       10       16
## BehLocPatch:BehaviorForaging                   -0.04 1.14       23       82
## BehLocPatch:BehaviorTransit                     3.08 1.24       14       52
## BehLocPatch:NameCanis_latrans                   0.22 1.05       67      136
## BehLocPatch:NameDidelphis_virginiana            1.17 1.13       22       19
## BehLocPatch:NameLynx_rufus                      1.81 1.05       40       98
## BehLocPatch:NameOdocoileus_virginianus          1.01 1.16       19       50
## BehLocPatch:NameProcyon_lotor                   0.82 1.12       24       50
## BehLocPatch:NameSciurus_carolinensis            4.64 1.11       33       90
## BehLocPatch:NameSciurus_niger                   5.28 1.20       16       57
## BehLocPatch:NameSus_scrofa                 972435.20 2.88        5       11
## BehLocPatch:NameSylvilagus_floridanus           0.95 1.06       56       61
## BehLocPatch:NameVulpes_vulpes          1177831277.90 2.99        5       14
## 
## Further Distributional Parameters:
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi     2.85      0.32     2.29     3.53 1.16       18       91
## 
## 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).
# Diagnostics and posterior checks
plot(Loc_model)

pp_check(Loc_model, type = 'dens_overlay')
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

pp_check(Loc_model, type = 'hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(Loc_model, type = 'boxplot')
## Using 10 posterior draws for ppc type 'boxplot' by default.

pp_check(Loc_model, type = 'intervals')
## Using all posterior draws for ppc type 'intervals' by default.

pp_check(Loc_model, type = 'scatter_avg')
## Using all posterior draws for ppc type 'scatter_avg' by default.

Within Patch Behavior Plots

# Compute marginal effects
marginal_effects_data <- conditional_effects(Loc_model)

# Plot marginal effects
plot(marginal_effects_data)

# Obtain predictions
predicted_values <- posterior_epred(Loc_model, newdata = invaded_data)

# Plot predicted vs. observed
plot(invaded_data$Proportion, colMeans(predicted_values))
abline(a = 0, b = 1, col = "red")

# Example interaction plot for BehLoc and Behavior
interaction_plot_data <- conditional_effects(Loc_model, effects = "BehLoc:Behavior")

# Convert to data frame for ggplot
df <- as.data.frame(interaction_plot_data$`BehLoc:Behavior`)

ggplot(df, aes(x = BehLoc, y = estimate__, color = Behavior, group = Behavior)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower__, ymax = upper__), alpha = 0.2) +
  theme_minimal() +
  labs(title = "Interaction Plot of BehLoc and Behavior",
       y = "Proportion",
       x = "BehLoc")

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

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 × 14
##    Plot  Wild_shannon_index Veg_shannon_index   Lat  Long Status     
##    <chr>              <dbl>             <dbl> <dbl> <dbl> <chr>      
##  1 BI201              0.305              2.70  30.8 -86.9 Invaded    
##  2 BN211              0.303              2.43  30.8 -86.9 Non_Invaded
##  3 EI100              0.292              2.72  31.0 -87.1 Invaded    
##  4 EI102              0.412              2.15  31.0 -87.1 Invaded    
##  5 EI104              0.310              2.90  31.0 -87.1 Invaded    
##  6 EI106              0.446              1.63  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              2.94  30.8 -87.1 Invaded    
## # ℹ 22 more rows
## # ℹ 8 more variables: Start_Date <dttm>, Camera <chr>, Cogon_Patch_Size <dbl>,
## #   VegetationDiversity <dbl>, PostTreatmentDensities <dbl>, Authority <chr>,
## #   Auth <dbl>, Avg_Cogongrass_Cover <dbl>
summary(diversity_data)
##      Plot           Wild_shannon_index Veg_shannon_index      Lat       
##  Length:32          Min.   :0.04539    Min.   :1.183     Min.   :28.68  
##  Class :character   1st Qu.:0.29022    1st Qu.:2.230     1st Qu.:30.76  
##  Mode  :character   Median :0.44221    Median :2.446     Median :30.77  
##                     Mean   :0.56021    Mean   :2.411     Mean   :30.45  
##                     3rd Qu.:0.75653    3rd Qu.:2.713     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       Avg_Cogongrass_Cover
##  Min.   :0.0000         Length:32          Min.   :1.000   Min.   : 0.000      
##  1st Qu.:0.0000         Class :character   1st Qu.:2.000   1st Qu.: 0.000      
##  Median :0.0000         Mode  :character   Median :3.000   Median : 2.143      
##  Mean   :0.8978                            Mean   :2.844   Mean   :13.683      
##  3rd Qu.:1.4650                            3rd Qu.:3.000   3rd Qu.:20.089      
##  Max.   :3.7100                            Max.   :4.000   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 Model

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

# 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 + (1 | Authority) + (1 | Camera)
##    Data: diversity_data
## 
## REML criterion at convergence: 44.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4432 -0.5454 -0.2552  0.3598  2.7706 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Authority (Intercept) 0.006203 0.07876 
##  Camera    (Intercept) 0.010002 0.10001 
##  Residual              0.151297 0.38897 
## Number of obs: 32, groups:  Authority, 5; Camera, 4
## 
## Fixed effects:
##                                        Estimate Std. Error       df t value
## (Intercept)                             0.54514    0.10252  2.49580   5.317
## Avg_Cogongrass_Cover                    0.04079    0.09692 25.48601   0.421
## Veg_shannon_index                       0.04254    0.08810 19.36180   0.483
## Cogon_Patch_Size                       -0.05781    0.09147 23.50451  -0.632
## Avg_Cogongrass_Cover:Veg_shannon_index -0.10232    0.06567 20.69791  -1.558
##                                        Pr(>|t|)  
## (Intercept)                              0.0204 *
## Avg_Cogongrass_Cover                     0.6774  
## Veg_shannon_index                        0.6346  
## Cogon_Patch_Size                         0.5335  
## Avg_Cogongrass_Cover:Veg_shannon_index   0.1344  
## ---
## 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
## Avg_Cgngr_C  0.066                     
## Vg_shnnn_nd -0.099  0.221              
## Cgn_Ptch_Sz -0.005 -0.594 -0.215       
## Avg_C_C:V__  0.204  0.414 -0.375 -0.321

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")
## [[1]]

## 
## [[2]]

Bayesian Diversity Model

bayesian_model <- brm(
  Wild_shannon_index ~ scale(Avg_Cogongrass_Cover) * scale(Veg_shannon_index) + scale(Cogon_Patch_Size) +
    (1 | Authority) + (1 | Camera),
  data = diversity_data,
  family = gaussian(),
  chains = 4,          
  iter = 10000,          
  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(cauchy(0, 1), class = "Intercept"),  # Weakly informative cauchy prior for intercept
    prior(cauchy(0, 1), class = "sd")       # For random effect standard deviations
  )
)
## Compiling Stan program...
## Start sampling
## Warning: There were 9 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
summary(bayesian_model)
## Warning: There were 9 divergent transitions after warmup. Increasing
## adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: Wild_shannon_index ~ scale(Avg_Cogongrass_Cover) * scale(Veg_shannon_index) + scale(Cogon_Patch_Size) + (1 | Authority) + (1 | Camera) 
##    Data: diversity_data (Number of observations: 32) 
##   Draws: 4 chains, each with iter = 10000; warmup = 3000; thin = 1;
##          total post-warmup draws = 28000
## 
## Multilevel Hyperparameters:
## ~Authority (Number of levels: 5) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.20      0.17     0.01     0.64 1.00    10132    12670
## 
## ~Camera (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.22      0.21     0.01     0.79 1.00     9455    12044
## 
## Regression Coefficients:
##                                                  Estimate Est.Error l-95% CI
## Intercept                                            0.51      0.21     0.06
## scaleAvg_Cogongrass_Cover                            0.05      0.11    -0.15
## scaleVeg_shannon_index                               0.04      0.10    -0.15
## scaleCogon_Patch_Size                               -0.07      0.10    -0.27
## scaleAvg_Cogongrass_Cover:scaleVeg_shannon_index    -0.10      0.07    -0.24
##                                                  u-95% CI Rhat Bulk_ESS
## Intercept                                            0.91 1.00    11666
## scaleAvg_Cogongrass_Cover                            0.27 1.00    17174
## scaleVeg_shannon_index                               0.24 1.00    22196
## scaleCogon_Patch_Size                                0.13 1.00    17396
## scaleAvg_Cogongrass_Cover:scaleVeg_shannon_index     0.05 1.00    19306
##                                                  Tail_ESS
## Intercept                                            8653
## scaleAvg_Cogongrass_Cover                           17875
## scaleVeg_shannon_index                              19408
## scaleCogon_Patch_Size                               18892
## scaleAvg_Cogongrass_Cover:scaleVeg_shannon_index    18484
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.41      0.06     0.31     0.55 1.00    21872    18820
## 
## 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),  # Use mean value for this predictor
  Authority = authority_levels[1],  # Assign a representative level for Authority
  Camera = camera_levels[1]  # Assign a representative level for Camera
)

# Repeat new_data to match for each combination of Camera levels
new_data <- as.data.frame(new_data)
new_data <- new_data[rep(seq_len(nrow(new_data)), each = length(authority_levels) * length(camera_levels)), ]
new_data$Authority <- rep(authority_levels, each = 100 * length(camera_levels))
new_data$Camera <- rep(camera_levels, times = 100 * length(authority_levels))

# Generate predictions including random effects
predictions <- fitted(bayesian_model, newdata = new_data, re_formula = NULL)
new_data$Wild_shannon_index <- rowMeans(predictions)

ggplot(new_data, aes(x = Avg_Cogongrass_Cover, y = Wild_shannon_index, color = factor(Veg_shannon_index))) +
  geom_line() +
  labs(x = "Average Cogongrass Cover", y = "Wild Shannon Index", color = "Vegetation Shannon Index") +
  theme_minimal() +
  ggtitle("Effect of Cogongrass Cover and Veg. Shannon Index on Wild Shannon Index")