C4_O2_Behavior

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, cmdstanr, lubridate)

# Install dependencies
#install.packages(c("posterior", "RcppParallel", "jsonlite"))

# Then install cmdstanr from the Stan developers’ repository
#install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))


# 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

# Calculate the total number of sites
total_sites <- nrow(CameraLoc)

# Function to filter data by frequency
filter_by_frequency <- function(df) {
  # Group data by species and calculate the frequency
  freq <- df %>%
    group_by(Species) %>%
    summarise(Frequency = n_distinct(Plot) / nrow(CameraLoc) * 100) %>%
    filter(Frequency >= 3)
  
  # Filter the original data to include only species with frequency >= 3%
  filtered_df <- df %>%
    filter(Species %in% freq$Species)
  
  return(filtered_df)
}

# Filter tree data by frequency
tree_data <- filter_by_frequency(tree_data)

# Filter Veg Cover data by frequency
Veg_Cover <- filter_by_frequency(Veg_Cover)

# Filter Shrub Cover data by frequency
shrub_data <- filter_by_frequency(shrub_data)

Shrub Cover Conversion

# Total length of Shrub cover at a site
shrub_cover <- shrub_data %>%
  mutate(Cover = Line_End - Line_Start) %>%
  group_by(Species_Name, Plot) %>%
  summarise(Shrub_Total_Cover = sum(Cover, na.rm = TRUE), .groups = "drop") %>%
  mutate(Shrub_Percent_Cover = Shrub_Total_Cover / 3000 * 100)

# Summed length of shrub over at a site
shrub_cover_summed <- shrub_cover %>%
  group_by(Plot) %>%
  summarize(total_shrub_cover = sum(Shrub_Total_Cover, na.rm = TRUE), .groups = "drop")

Herbacous Cover Conversion

# Combine Plot and Quadrat columns
Veg_Cover <- Veg_Cover %>%
  mutate(Plot_Quadrat = paste(Plot, Quadrat, sep = '_'))

# Join with CogonSites to get site information
Veg_Cover <- Veg_Cover %>%
  left_join(CameraLoc, by = "Plot")

# Sum species cover across quadrats for each species at each plot
veg_cover_summed <- Veg_Cover %>%
  group_by(Plot, Species_Name) %>%
  summarize(total_cover = sum(Cover_Per, na.rm = TRUE), .groups = "drop")

# Calculate average herbaceous species cover
avg_species_cover <- veg_cover_summed %>%
  left_join(quadrat_count, by = "Plot") %>%
  mutate(avg_cover = total_cover / total_quadrats)

Merging Herb cover with Shrub

This species matrix includes herbaceous and shrub species

# Merge shrub cover with herbaceous average cover
combined_cover <- avg_species_cover %>%
  full_join(
    shrub_cover %>%
      dplyr::select(Plot, Species_Name, Shrub_Percent_Cover),
    by = c("Plot", "Species_Name")
  ) %>%
  mutate(
    overlap_flag = ifelse(!is.na(avg_cover) & !is.na(Shrub_Percent_Cover), TRUE, FALSE), # Flag overlaps
    final_cover = case_when(
      !is.na(avg_cover) & is.na(Shrub_Percent_Cover) ~ avg_cover,  # Use herbaceous cover if no shrub data
      is.na(avg_cover) & !is.na(Shrub_Percent_Cover) ~ Shrub_Percent_Cover, # Use shrub cover if no herbaceous data
      TRUE ~ NA_real_ # Leave as NA where overlaps exist
    )
  )

# Species Matrix
species_matrix <- combined_cover %>%
  dplyr::select(Plot, Species_Name, final_cover) %>%
  pivot_wider(
    names_from = Species_Name,
    values_from = final_cover,
    values_fill = 0
  )

Summarize Cogongrass Cover

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

Herbacous Shannon Diversity Index

# Summarize species cover by site
site_species_cover <- Veg_Cover %>%
  group_by(Plot, Species_Name) %>%
  summarize(total_cover = sum(Cover_Per, na.rm = TRUE)) %>%
  ungroup()

## `summarise()` has grouped output by 'Plot'. You can override using the
## `.groups` argument.

# Calculate Shannon diversity per site
Veg_shannon_diversity <- site_species_cover %>%
  group_by(Plot) %>%
  mutate(proportion = total_cover / sum(total_cover)) %>%
  summarize(Veg_shannon_index = -sum(proportion * log(proportion), na.rm = TRUE))

print(Veg_shannon_diversity)

## # A tibble: 206 × 2
##    Plot  Veg_shannon_index
##    <chr>             <dbl>
##  1 BI200              2.28
##  2 BI201              2.20
##  3 BI202              1.50
##  4 BI97               1.82
##  5 BI99               3.06
##  6 BN210              2.97
##  7 BN211              2.43
##  8 BN212              2.22
##  9 BN96               3.05
## 10 BN98               2.79
## # ℹ 196 more rows

Vegetation Height

if (!is.numeric(fuel_data$Height)) {
  fuel_data$Height <- as.numeric(as.character(fuel_data$Height))
}

## Warning: NAs introduced by coercion

# Calculate average vegetation height per plot
veg_height <- fuel_data %>%
  group_by(Plot) %>%
  summarize(avg_veg_height = mean(Height, na.rm = TRUE), .groups = "drop")

Tree Metrics

# Tree density from point-centered quarter data
if (!is.numeric(tree_data$Distance)) {
  tree_data$Distance <- as.numeric(as.character(tree_data$Distance))
}

tree_density_data <- tree_data %>%
  group_by(Plot) %>%
  summarize(Average_Distance = mean(Distance) / 100,  # Convert to meters
            Tree_Density = 10000 / (Average_Distance^2))  # Convert to trees per hectare

# Average canopy cover from vegetation quadrats
tree_canopy_data <- Veg_Cover %>%
  distinct(Plot, Quadrat, .keep_all = TRUE) %>%  # Ensure each quadrat counts once per plot
  group_by(Plot) %>%
  summarize(Avg_Canopy_Cover = mean(Canopy_Cover, na.rm = TRUE), .groups = "drop") # Calculate the average canopy cover per plot

cor(tree_density_data$Tree_Density, tree_canopy_data$Avg_Canopy_Cover)

## [1] 0.2742307

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

Add Covariates

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

Subset to mammals with >2 observations

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

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

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

Create Detection Matrix

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

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

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

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

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

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

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

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

## [1] 32 36  9

Prepare Site Covariates

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

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

# Standardizing covariates
covariates_matrix <- scale(covariates_matrix)

Combine data into a list

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

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

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

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

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

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

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

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

str(data_list)

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

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

O2_data <- O2_data %>%
  mutate(
    DateTime = update(Date,
                      hour   = hour(Time),
                      minute = minute(Time),
                      second = second(Time))
  )

gap_mins <- 30

O2_data <- O2_data %>%
  filter(!is.na(DateTime)) %>%
  arrange(Plot, Name, DateTime) %>%
  group_by(Plot, Name) %>%
  group_modify(~{
    df <- .x
    keep <- logical(nrow(df))
    last_kept <- as.POSIXct(NA, tz = tz(df$DateTime[1]))
    for (i in seq_len(nrow(df))) {
      if (is.na(last_kept) || difftime(df$DateTime[i], last_kept, units = "mins") > gap_mins) {
        keep[i] <- TRUE
        last_kept <- df$DateTime[i]
      }
    }
    df[keep, , drop = FALSE]
  }) %>%
  ungroup()

# Creating proportions of observations for each behavior type
behavior_counts <- O2_data %>%
  group_by(Plot, Name, Behavior, Status, Site, Camera.Type, BehLoc) %>%
  summarise(ObservationCount = n(), .groups = "drop_last") %>%
  mutate(TotalObs = sum(ObservationCount)) %>%
  ungroup()

# collapse rare species (threshold = 10 rows; tune if needed)
keep_names <- behavior_counts %>%
  count(Name) %>% filter(n >= 10) %>% pull(Name)

behavior_counts <- behavior_counts %>%
  mutate(Name_group = if_else(Name %in% keep_names, Name, "Other"))

# quick check
behavior_counts %>% count(Name_group) %>% arrange(n)

## # A tibble: 9 × 2
##   Name_group                 n
##   <chr>                  <int>
## 1 Lynx_rufus                12
## 2 Didelphis_virginiana      13
## 3 Meleagris_gallopavo       14
## 4 Sciurus_carolinensis      15
## 5 Sylvilagus_floridanus     19
## 6 Canis_latrans             27
## 7 Dasypus_novemcinctus      31
## 8 Procyon_lotor             47
## 9 Odocoileus_virginianus   137

# 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

# Ensure Behavior is an unordered factor
if(!is.factor(behavior_counts$Behavior) || is.ordered(behavior_counts$Behavior)) {
  behavior_counts$Behavior <- factor(behavior_counts$Behavior, ordered = FALSE)
}
# Make Foraging the reference category
behavior_counts$Behavior <- relevel(behavior_counts$Behavior, ref = "Local_Search")

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

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

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

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

## Ensure Site is a factor
if(!is.factor(behavior_counts$Plot) || is.ordered(behavior_counts$Plot)) {
  behavior_counts$Plot <- factor(behavior_counts$Plot, ordered = FALSE)
}

priors <- c(
  prior(normal(0, 1), class = "b"),              # regularizing slopes
  prior(student_t(3, 0, 2.5), class = "Intercept"),
  prior(exponential(1), class = "sd")            # for random effects
)


# Fit the model
brms_model <- brm(
  ObservationCount | trials(TotalObs) ~ Status * Behavior + Status * Name + (1 | Site) + (1 | Camera.Type),
  data = behavior_counts,
  family = binomial(),
  prior = priors,
  iter = 5000, warmup = 3000, chains = 4, cores = 4,
  control = list(adapt_delta = 0.99, max_treedepth = 15)
)

## Compiling Stan program...

## Start sampling

# Diagnostic Plots
plot(brms_model)

# Summary Statistics
summary(brms_model)

##  Family: binomial 
##   Links: mu = logit 
## Formula: ObservationCount | trials(TotalObs) ~ Status * Behavior + Status * Name + (1 | Site) + (1 | Camera.Type) 
##    Data: behavior_counts (Number of observations: 315) 
##   Draws: 4 chains, each with iter = 5000; warmup = 3000; 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.32      0.23     0.09     0.93 1.00     3700     4409
## 
## ~Site (Number of levels: 5) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.15      0.15     0.00     0.54 1.00     2493     3293
## 
## Regression Coefficients:
##                                         Estimate Est.Error l-95% CI u-95% CI
## Intercept                                   5.12      0.54     4.10     6.23
## StatusInvaded                              -4.67      0.50    -5.73    -3.74
## BehaviorForaging                            0.66      0.59    -0.47     1.86
## BehaviorTransit                             0.83      0.56    -0.27     1.94
## NameCanis_latrans                           1.23      0.82    -0.36     2.83
## NameDasypus_novemcinctus                    0.08      0.69    -1.24     1.44
## NameDidelphis_virginiana                    0.28      0.74    -1.15     1.74
## NameLynx_rufus                              1.15      0.82    -0.42     2.76
## NameMeleagris_gallopavo                     0.28      0.80    -1.28     1.82
## NameProcyon_lotor                           0.18      0.69    -1.12     1.57
## NameSciurus_carolinensis                    0.62      0.76    -0.83     2.10
## NameSylvilagus_floridanus                   0.44      0.74    -1.02     1.87
## StatusInvaded:BehaviorForaging             -0.54      0.60    -1.73     0.58
## StatusInvaded:BehaviorTransit              -1.34      0.56    -2.45    -0.24
## StatusInvaded:NameCanis_latrans             1.20      0.82    -0.38     2.81
## StatusInvaded:NameDasypus_novemcinctus     -0.02      0.69    -1.38     1.33
## StatusInvaded:NameDidelphis_virginiana      0.28      0.74    -1.21     1.70
## StatusInvaded:NameLynx_rufus                1.15      0.83    -0.47     2.80
## StatusInvaded:NameMeleagris_gallopavo       0.24      0.81    -1.35     1.83
## StatusInvaded:NameProcyon_lotor             0.07      0.69    -1.32     1.37
## StatusInvaded:NameSciurus_carolinensis      0.63      0.76    -0.85     2.09
## StatusInvaded:NameSylvilagus_floridanus     0.43      0.73    -1.00     1.86
##                                         Rhat Bulk_ESS Tail_ESS
## Intercept                               1.00     5840     5773
## StatusInvaded                           1.00     7196     6052
## BehaviorForaging                        1.00     6259     5252
## BehaviorTransit                         1.00     6016     5204
## NameCanis_latrans                       1.00    10056     6420
## NameDasypus_novemcinctus                1.00     6833     5272
## NameDidelphis_virginiana                1.00     6777     6201
## NameLynx_rufus                          1.00    10298     6216
## NameMeleagris_gallopavo                 1.00     8536     6442
## NameProcyon_lotor                       1.00     7166     5849
## NameSciurus_carolinensis                1.00     7244     6328
## NameSylvilagus_floridanus               1.00     7521     5866
## StatusInvaded:BehaviorForaging          1.00     6234     5518
## StatusInvaded:BehaviorTransit           1.00     5942     5003
## StatusInvaded:NameCanis_latrans         1.00     9686     6386
## StatusInvaded:NameDasypus_novemcinctus  1.00     7060     5835
## StatusInvaded:NameDidelphis_virginiana  1.00     6798     5933
## StatusInvaded:NameLynx_rufus            1.00    10534     6299
## StatusInvaded:NameMeleagris_gallopavo   1.00     7911     6157
## StatusInvaded:NameProcyon_lotor         1.00     7098     5731
## StatusInvaded:NameSciurus_carolinensis  1.00     7315     5997
## StatusInvaded:NameSylvilagus_floridanus 1.00     7770     6366
## 
## 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).

Behavior Model Within cogongrass patches

# Create counts only for invaded sites
behavior_counts_invaded <- O2_data %>%
  filter(Status == "Invaded") %>% 
  group_by(Plot, Name, Behavior, Status, Site, Camera.Type, BehLoc) %>%
  summarise(ObservationCount = n(), .groups = "drop_last") %>%
  mutate(TotalObs = sum(ObservationCount)) %>%
  ungroup()

# Ensure Behavior is an unordered factor
if(!is.factor(behavior_counts_invaded$Behavior) || is.ordered(behavior_counts_invaded$Behavior)) {
  behavior_counts_invaded$Behavior <- factor(behavior_counts_invaded$Behavior, ordered = FALSE)
}
# Make Foraging the reference category
behavior_counts_invaded$Behavior <- relevel(behavior_counts_invaded$Behavior, ref = "Local_Search")

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

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

priors <- c(
  prior(normal(0, 1), class = "b"),              # regularizing slopes
  prior(student_t(3, 0, 2.5), class = "Intercept"),
  prior(exponential(1), class = "sd")            # for random effects
)



# Fit the model (binomial with trials)
Loc_model <- brm(
  ObservationCount | trials(TotalObs) ~ BehLoc * Behavior + BehLoc * Name +
    (1 | Site) + (1 | Camera.Type),
  data = behavior_counts_invaded,
  family = binomial(),
  prior = priors,
  iter = 5000, warmup = 3000, chains = 4, cores = 4,
  control = list(adapt_delta = 0.99, max_treedepth = 15)
)

## Compiling Stan program...

## Start sampling

## Warning: There were 5 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

# Summarize model
summary(Loc_model)

## Warning: There were 5 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: binomial 
##   Links: mu = logit 
## Formula: ObservationCount | trials(TotalObs) ~ BehLoc * Behavior + BehLoc * Name + (1 | Site) + (1 | Camera.Type) 
##    Data: behavior_counts_invaded (Number of observations: 186) 
##   Draws: 4 chains, each with iter = 5000; warmup = 3000; 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.21      0.21     0.01     0.72 1.00     2245     2840
## 
## ~Site (Number of levels: 5) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.14      0.14     0.00     0.49 1.00     3064     4314
## 
## Regression Coefficients:
##                                       Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                                 1.31      0.25     0.80     1.77 1.00
## BehLocPatch                              -2.04      0.23    -2.50    -1.59 1.00
## BehaviorForaging                          1.33      0.29     0.76     1.92 1.00
## BehaviorTransit                          -1.46      0.17    -1.80    -1.12 1.00
## NameCanis_latrans                         1.62      0.70     0.32     3.03 1.00
## NameDasypus_novemcinctus                  0.13      0.24    -0.36     0.61 1.00
## NameDidelphis_virginiana                  0.19      0.47    -0.73     1.12 1.00
## NameLynx_rufus                            1.79      0.66     0.55     3.16 1.00
## NameMeleagris_gallopavo                   0.24      0.71    -1.14     1.66 1.00
## NameProcyon_lotor                         0.59      0.25     0.11     1.08 1.00
## NameSciurus_carolinensis                  0.98      0.55    -0.06     2.13 1.00
## NameSylvilagus_floridanus                 0.73      0.54    -0.33     1.80 1.00
## BehLocPatch:BehaviorForaging             -2.59      0.41    -3.40    -1.79 1.00
## BehLocPatch:BehaviorTransit               2.40      0.25     1.92     2.91 1.00
## BehLocPatch:NameCanis_latrans             0.67      0.84    -0.94     2.32 1.00
## BehLocPatch:NameDasypus_novemcinctus     -0.19      0.33    -0.83     0.47 1.00
## BehLocPatch:NameDidelphis_virginiana      0.41      0.62    -0.75     1.63 1.00
## BehLocPatch:NameLynx_rufus                0.12      0.93    -1.68     1.99 1.00
## BehLocPatch:NameMeleagris_gallopavo       0.09      0.83    -1.52     1.72 1.00
## BehLocPatch:NameProcyon_lotor            -1.05      0.35    -1.72    -0.36 1.00
## BehLocPatch:NameSciurus_carolinensis     -0.75      0.82    -2.29     0.88 1.00
## BehLocPatch:NameSylvilagus_floridanus    -0.55      0.71    -1.94     0.84 1.00
##                                       Bulk_ESS Tail_ESS
## Intercept                                 3890     2631
## BehLocPatch                               4852     5575
## BehaviorForaging                          6476     5966
## BehaviorTransit                           5812     5888
## NameCanis_latrans                        11125     5750
## NameDasypus_novemcinctus                  7967     6293
## NameDidelphis_virginiana                  7951     5649
## NameLynx_rufus                           13330     5416
## NameMeleagris_gallopavo                  11399     6419
## NameProcyon_lotor                         8366     6414
## NameSciurus_carolinensis                 11826     6001
## NameSylvilagus_floridanus                10467     6027
## BehLocPatch:BehaviorForaging              5698     5422
## BehLocPatch:BehaviorTransit               4722     5337
## BehLocPatch:NameCanis_latrans            11130     6144
## BehLocPatch:NameDasypus_novemcinctus      8137     6092
## BehLocPatch:NameDidelphis_virginiana      8921     6311
## BehLocPatch:NameLynx_rufus               13274     5815
## BehLocPatch:NameMeleagris_gallopavo      11896     5898
## BehLocPatch:NameProcyon_lotor             8233     6488
## BehLocPatch:NameSciurus_carolinensis     12057     6400
## BehLocPatch:NameSylvilagus_floridanus     9981     5748
## 
## 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 `binwidth`.

pp_check(Loc_model, type = 'boxplot')

## Using 10 posterior draws for ppc type 'boxplot' by default.
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?Notch went outside hinges
## ℹ Do you want `notch = FALSE`?Notch went outside hinges
## ℹ Do you want `notch = FALSE`?Notch went outside hinges
## ℹ Do you want `notch = FALSE`?Notch went outside hinges
## ℹ Do you want `notch = FALSE`?Notch went outside hinges
## ℹ Do you want `notch = FALSE`?Notch went outside hinges
## ℹ Do you want `notch = FALSE`?Notch went outside hinges
## ℹ Do you want `notch = FALSE`?Notch went outside hinges
## ℹ Do you want `notch = FALSE`?Notch went outside hinges
## ℹ Do you want `notch = FALSE`?

pp_check(Loc_model, type = 'intervals')

## Using all posterior draws for ppc type 'intervals' by default.

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

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)

## Setting all 'trials' variables to 1 by default if not specified otherwise.

# Plot marginal effects
plot(marginal_effects_data)

## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"

## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"

## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"

## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"

# Add proportions for plotting predicted vs observed:
invaded_data <- behavior_counts_invaded %>%
  group_by(Plot, Name, Status, Camera.Type, BehLoc) %>%
  mutate(Proportion = ObservationCount / sum(ObservationCount)) %>%
  ungroup()

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

# Plot predicted (posterior mean) vs observed
plot(invaded_data$Proportion,
     colMeans(predicted_values),
     xlab = "Observed Proportion",
     ylab = "Predicted Proportion",
     pch = 19)
abline(a = 0, b = 1, col = "red", lwd = 2)

# Extract interaction marginal effect
interaction_plot_data <- conditional_effects(Loc_model, effects = "BehLoc:Behavior")

## Setting all 'trials' variables to 1 by default if not specified otherwise.

df <- as.data.frame(interaction_plot_data$`BehLoc:Behavior`)

ggplot(df, aes(x = BehLoc, y = estimate__, color = Behavior, group = Behavior)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = lower__, ymax = upper__, fill = Behavior),
              alpha = 0.2, color = NA) +
  theme_minimal() +
  labs(title = "Interaction Plot: BehLoc × Behavior",
       y = "Estimated Probability",
       x = "Behavior Location (Patch vs Non-Patch)")

AI Version