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/alanivory34428/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/alanivory34428/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/alanivory34428/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/alanivory34428/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/alanivory34428/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/04_Wildlife/02_Data/CameraData.xlsx")

CameraLoc <- read_excel("C:/Users/alanivory34428/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/alanivory34428/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: 174 × 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
## # ℹ 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

#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 = 4000, warmup = 1000, 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: 327) 
##   Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
##          total post-warmup draws = 12000
## 
## Multilevel Hyperparameters:
## ~Plot (Number of levels: 32) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.24      0.11     0.03     0.47 1.00     2988     2979
## 
## Regression Coefficients:
##                                          Estimate Est.Error l-95% CI u-95% CI
## Intercept                                   -0.02      0.43    -0.86     0.85
## StatusInvaded                                0.08      0.53    -0.95     1.13
## BehaviorForaging                             0.58      0.30    -0.00     1.17
## BehaviorTransit                              1.36      0.27     0.84     1.88
## NameCanis_latrans                            0.18      0.47    -0.76     1.11
## NameDidelphis_virginiana                     0.69      0.70    -0.61     2.12
## NameLynx_rufus                               1.16      0.63    -0.02     2.43
## NameOdocoileus_virginianus                  -1.05      0.40    -1.84    -0.26
## NameProcyon_lotor                            0.38      0.46    -0.52     1.28
## NameSciurus_carolinensis                    -0.13      0.65    -1.36     1.16
## NameSciurus_niger                            1.54      0.77     0.11     3.16
## NameSus_scrofa                               0.11      0.62    -1.08     1.35
## NameSylvilagus_floridanus                    1.16      0.53     0.13     2.22
## NameVulpes_vulpes                            1.32      0.93    -0.32     3.33
## StatusInvaded:BehaviorForaging              -0.63      0.38    -1.38     0.11
## StatusInvaded:BehaviorTransit               -0.04      0.33    -0.68     0.59
## StatusInvaded:NameCanis_latrans              1.05      0.63    -0.18     2.28
## StatusInvaded:NameDidelphis_virginiana       0.05      0.82    -1.61     1.64
## StatusInvaded:NameLynx_rufus                 0.10      0.79    -1.49     1.63
## StatusInvaded:NameOdocoileus_virginianus     0.36      0.48    -0.59     1.31
## StatusInvaded:NameProcyon_lotor              0.04      0.56    -1.06     1.13
## StatusInvaded:NameSciurus_carolinensis       0.09      0.77    -1.44     1.57
## StatusInvaded:NameSciurus_niger             -1.05      1.01    -3.08     0.90
## StatusInvaded:NameSus_scrofa                 1.75      1.48    -0.81     5.10
## StatusInvaded:NameSylvilagus_floridanus     -0.33      0.67    -1.62     0.97
## StatusInvaded:NameVulpes_vulpes              0.32      1.66    -2.64     3.92
##                                          Rhat Bulk_ESS Tail_ESS
## Intercept                                1.00     2313     4413
## StatusInvaded                            1.00     2266     4654
## BehaviorForaging                         1.00     6561     8245
## BehaviorTransit                          1.00     6315     7963
## NameCanis_latrans                        1.00     2839     5598
## NameDidelphis_virginiana                 1.00     4317     6644
## NameLynx_rufus                           1.00     4141     6623
## NameOdocoileus_virginianus               1.00     2419     4446
## NameProcyon_lotor                        1.00     2700     5061
## NameSciurus_carolinensis                 1.00     4496     7333
## NameSciurus_niger                        1.00     4894     7473
## NameSus_scrofa                           1.00     4396     7461
## NameSylvilagus_floridanus                1.00     3466     6224
## NameVulpes_vulpes                        1.00     6439     7381
## StatusInvaded:BehaviorForaging           1.00     6593     8253
## StatusInvaded:BehaviorTransit            1.00     6340     8595
## StatusInvaded:NameCanis_latrans          1.00     3368     5790
## StatusInvaded:NameDidelphis_virginiana   1.00     4342     6962
## StatusInvaded:NameLynx_rufus             1.00     4273     6094
## StatusInvaded:NameOdocoileus_virginianus 1.00     2320     4337
## StatusInvaded:NameProcyon_lotor          1.00     2765     4765
## StatusInvaded:NameSciurus_carolinensis   1.00     4501     6566
## StatusInvaded:NameSciurus_niger          1.00     5656     7872
## StatusInvaded:NameSus_scrofa             1.00     8634     7304
## StatusInvaded:NameSylvilagus_floridanus  1.00     3727     7104
## StatusInvaded:NameVulpes_vulpes          1.00     7799     7690
## 
## Further Distributional Parameters:
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi     2.00      0.16     1.69     2.33 1.00     9654     8086
## 
## 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 = 4000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  control = list(adapt_delta = 0.99, max_treedepth = 15)
)
## Compiling Stan program...
## Start sampling
## Warning: There were 8000 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 3.06, 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.
##  Family: beta 
##   Links: mu = logit; phi = identity 
## Formula: Proportion ~ BehLoc * Behavior + BehLoc * Name + (1 | Plot) + (1 | Camera.Type) 
##    Data: invaded_data (Number of observations: 197) 
##   Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## 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.72      0.35     0.15     1.52 1.50        7       12
## 
## ~Plot (Number of levels: 17) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.49      0.21     0.19     0.93 1.66        6       19
## 
## Regression Coefficients:
##                                           Estimate   Est.Error     l-95% CI
## Intercept                                     0.24        0.33        -0.34
## BehLocPatch                                  -0.45        0.43        -1.26
## BehaviorForaging                              0.19        0.24        -0.44
## BehaviorTransit                               0.69        0.20         0.29
## NameCanis_latrans                             1.36        0.45         0.54
## NameDidelphis_virginiana                      0.66        0.29         0.08
## NameLynx_rufus                                1.88        0.47         0.93
## NameOdocoileus_virginianus                   -0.91        0.26        -1.48
## NameProcyon_lotor                             0.32        0.24        -0.25
## NameSciurus_carolinensis                      0.01        0.33        -0.58
## NameSciurus_niger                             0.03        0.55        -0.96
## NameSus_scrofa                           -80555.26    88295.16   -206047.01
## NameSylvilagus_floridanus                     0.99        0.37         0.26
## NameVulpes_vulpes                             2.59        1.18         0.42
## BehLocPatch:BehaviorForaging                 -0.31        0.40        -1.25
## BehLocPatch:BehaviorTransit                   1.79        0.27         1.33
## BehLocPatch:NameCanis_latrans                -0.84        0.86        -2.59
## BehLocPatch:NameDidelphis_virginiana         -0.17        0.59        -1.36
## BehLocPatch:NameLynx_rufus                   -0.98        1.06        -2.94
## BehLocPatch:NameOdocoileus_virginianus        0.03        0.52        -0.67
## BehLocPatch:NameProcyon_lotor                -0.15        0.67        -1.00
## BehLocPatch:NameSciurus_carolinensis          1.19        0.99        -0.59
## BehLocPatch:NameSciurus_niger                 1.36        1.56        -0.80
## BehLocPatch:NameSus_scrofa                80556.89    88295.29   -112798.56
## BehLocPatch:NameSylvilagus_floridanus        -1.22        0.83        -2.54
## BehLocPatch:NameVulpes_vulpes          -3801242.66 10462543.96 -33077548.10
##                                          u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept                                    0.85 1.66        6       22
## BehLocPatch                                  0.35 1.93        6       21
## BehaviorForaging                             0.66 1.67        7       12
## BehaviorTransit                              1.06 1.46        8       45
## NameCanis_latrans                            2.27 1.71        6       13
## NameDidelphis_virginiana                     1.24 1.32       11       39
## NameLynx_rufus                               2.75 1.53        7       22
## NameOdocoileus_virginianus                  -0.54 1.48        8       20
## NameProcyon_lotor                            0.72 1.38       10       32
## NameSciurus_carolinensis                     0.65 1.14       23       58
## NameSciurus_niger                            1.08 1.74        6       20
## NameSus_scrofa                          112798.85 2.36        5       17
## NameSylvilagus_floridanus                    1.79 1.52        7       18
## NameVulpes_vulpes                            4.98 1.39        9       35
## BehLocPatch:BehaviorForaging                 0.36 1.44        8       12
## BehLocPatch:BehaviorTransit                  2.46 1.31       10       25
## BehLocPatch:NameCanis_latrans                0.95 1.96        6       25
## BehLocPatch:NameDidelphis_virginiana         0.99 1.37        9       39
## BehLocPatch:NameLynx_rufus                   1.16 1.36        9       34
## BehLocPatch:NameOdocoileus_virginianus       1.11 1.87        6       17
## BehLocPatch:NameProcyon_lotor                1.27 1.63        7       23
## BehLocPatch:NameSciurus_carolinensis         3.11 1.37       10       23
## BehLocPatch:NameSciurus_niger                4.44 1.86        6       27
## BehLocPatch:NameSus_scrofa              206048.59 2.36        5       17
## BehLocPatch:NameSylvilagus_floridanus        0.37 1.48        8       43
## BehLocPatch:NameVulpes_vulpes          6422678.83 3.06        7       12
## 
## Further Distributional Parameters:
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi     2.91      0.33     2.44     3.58 1.42        8       29
## 
## 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/alanivory34428/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/04_Wildlife/02_Data/CameraLoc.xlsx",
                  sheet = "CameraLocations")

CameraData <- read_excel("C:/Users/alanivory34428/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.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
## # ℹ 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.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       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
)
## 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 + (1 | Authority) + (1 | Camera)
##    Data: diversity_data
## 
## REML criterion at convergence: 38.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7385 -0.6306 -0.1460  0.4177  2.7655 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Authority (Intercept) 0.0000   0.0000  
##  Camera    (Intercept) 0.0260   0.1613  
##  Residual              0.1279   0.3577  
## Number of obs: 32, groups:  Authority, 4; Camera, 4
## 
## Fixed effects:
##                                         Estimate Std. Error        df t value
## (Intercept)                             0.450256   0.121460  4.613881   3.707
## Avg_Cogongrass_Cover                   -0.166255   0.155059 26.891653  -1.072
## Veg_shannon_index                       0.010820   0.108964 25.628316   0.099
## Cogon_Patch_Size                       -0.005934   0.078788 26.390634  -0.075
## Avg_Cogongrass_Cover:Veg_shannon_index -0.222671   0.080042 26.921473  -2.782
##                                        Pr(>|t|)   
## (Intercept)                             0.01606 * 
## Avg_Cogongrass_Cover                    0.29316   
## Veg_shannon_index                       0.92168   
## Cogon_Patch_Size                        0.94053   
## Avg_Cogongrass_Cover:Veg_shannon_index  0.00975 **
## ---
## 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.269                     
## Vg_shnnn_nd -0.021  0.611              
## Cgn_Ptch_Sz -0.148 -0.477 -0.058       
## Avg_C_C:V__  0.434  0.685  0.074 -0.403
## 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")
## [[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 = 4000,          
  warmup = 2000,       
  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
summary(bayesian_model)
##  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 = 4000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Multilevel Hyperparameters:
## ~Authority (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.19      0.20     0.01     0.70 1.00     3288     3581
## 
## ~Camera (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.26      0.23     0.01     0.85 1.00     2323     2857
## 
## Regression Coefficients:
##                                                  Estimate Est.Error l-95% CI
## Intercept                                            0.41      0.23    -0.09
## scaleAvg_Cogongrass_Cover                           -0.13      0.18    -0.46
## scaleVeg_shannon_index                               0.02      0.13    -0.22
## scaleCogon_Patch_Size                               -0.03      0.09    -0.21
## scaleAvg_Cogongrass_Cover:scaleVeg_shannon_index    -0.20      0.09    -0.38
##                                                  u-95% CI Rhat Bulk_ESS
## Intercept                                            0.83 1.00     3847
## scaleAvg_Cogongrass_Cover                            0.23 1.00     3653
## scaleVeg_shannon_index                               0.27 1.00     5746
## scaleCogon_Patch_Size                                0.15 1.00     4898
## scaleAvg_Cogongrass_Cover:scaleVeg_shannon_index    -0.01 1.00     4398
##                                                  Tail_ESS
## Intercept                                            3061
## scaleAvg_Cogongrass_Cover                            4344
## scaleVeg_shannon_index                               5148
## scaleCogon_Patch_Size                                5132
## scaleAvg_Cogongrass_Cover:scaleVeg_shannon_index     4934
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.38      0.06     0.29     0.52 1.00     5356     4650
## 
## 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")