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.33      0.24     0.09     0.98 1.00     3291     4301
## 
## ~Site (Number of levels: 5) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.16      0.16     0.00     0.55 1.00     2797     3887
## 
## Regression Coefficients:
##                                         Estimate Est.Error l-95% CI u-95% CI
## Intercept                                   5.11      0.55     4.05     6.24
## StatusInvaded                              -4.67      0.51    -5.73    -3.72
## BehaviorForaging                            0.67      0.61    -0.49     1.89
## BehaviorTransit                             0.84      0.56    -0.29     1.97
## NameCanis_latrans                           1.22      0.80    -0.31     2.85
## NameDasypus_novemcinctus                    0.08      0.69    -1.25     1.41
## NameDidelphis_virginiana                    0.28      0.72    -1.15     1.72
## NameLynx_rufus                              1.15      0.82    -0.43     2.75
## NameMeleagris_gallopavo                     0.28      0.78    -1.23     1.85
## NameProcyon_lotor                           0.18      0.68    -1.12     1.55
## NameSciurus_carolinensis                    0.62      0.76    -0.90     2.12
## NameSylvilagus_floridanus                   0.46      0.74    -0.98     1.89
## StatusInvaded:BehaviorForaging             -0.54      0.61    -1.77     0.64
## StatusInvaded:BehaviorTransit              -1.34      0.57    -2.49    -0.22
## StatusInvaded:NameCanis_latrans             1.22      0.82    -0.39     2.83
## StatusInvaded:NameDasypus_novemcinctus     -0.02      0.69    -1.36     1.32
## StatusInvaded:NameDidelphis_virginiana      0.28      0.72    -1.14     1.70
## StatusInvaded:NameLynx_rufus                1.16      0.82    -0.45     2.79
## StatusInvaded:NameMeleagris_gallopavo       0.24      0.79    -1.33     1.80
## StatusInvaded:NameProcyon_lotor             0.07      0.68    -1.28     1.37
## StatusInvaded:NameSciurus_carolinensis      0.61      0.76    -0.88     2.11
## StatusInvaded:NameSylvilagus_floridanus     0.40      0.74    -1.05     1.85
##                                         Rhat Bulk_ESS Tail_ESS
## Intercept                               1.00     6106     5718
## StatusInvaded                           1.00     7441     5571
## BehaviorForaging                        1.00     6204     5167
## BehaviorTransit                         1.00     6075     5456
## NameCanis_latrans                       1.00     9804     6383
## NameDasypus_novemcinctus                1.00     7007     6181
## NameDidelphis_virginiana                1.00     7329     5596
## NameLynx_rufus                          1.00     9114     6226
## NameMeleagris_gallopavo                 1.00     9287     6316
## NameProcyon_lotor                       1.00     6697     5893
## NameSciurus_carolinensis                1.00     7846     6271
## NameSylvilagus_floridanus               1.00     7791     6295
## StatusInvaded:BehaviorForaging          1.00     6275     5464
## StatusInvaded:BehaviorTransit           1.00     6114     5423
## StatusInvaded:NameCanis_latrans         1.00     9669     6109
## StatusInvaded:NameDasypus_novemcinctus  1.00     6740     5746
## StatusInvaded:NameDidelphis_virginiana  1.00     7200     6241
## StatusInvaded:NameLynx_rufus            1.00     9840     6334
## StatusInvaded:NameMeleagris_gallopavo   1.00     9552     6296
## StatusInvaded:NameProcyon_lotor         1.00     6773     5520
## StatusInvaded:NameSciurus_carolinensis  1.00     7248     6120
## StatusInvaded:NameSylvilagus_floridanus 1.00     7993     6152
## 
## 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 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: Examine the pairs() plot to diagnose sampling problems

# Summarize model
summary(Loc_model)

## 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: 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.22      0.22     0.01     0.75 1.00     2317     2640
## 
## ~Site (Number of levels: 5) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.14      0.13     0.00     0.48 1.00     3101     3950
## 
## Regression Coefficients:
##                                       Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                                 1.32      0.26     0.81     1.76 1.00
## BehLocPatch                              -2.04      0.23    -2.50    -1.59 1.00
## BehaviorForaging                          1.32      0.29     0.78     1.89 1.00
## BehaviorTransit                          -1.46      0.18    -1.80    -1.13 1.00
## NameCanis_latrans                         1.60      0.69     0.33     3.01 1.00
## NameDasypus_novemcinctus                  0.14      0.25    -0.35     0.63 1.00
## NameDidelphis_virginiana                  0.18      0.47    -0.72     1.15 1.00
## NameLynx_rufus                            1.79      0.66     0.54     3.14 1.00
## NameMeleagris_gallopavo                   0.25      0.71    -1.12     1.65 1.00
## NameProcyon_lotor                         0.59      0.25     0.11     1.08 1.00
## NameSciurus_carolinensis                  0.98      0.56    -0.08     2.12 1.00
## NameSylvilagus_floridanus                 0.72      0.53    -0.30     1.79 1.00
## BehLocPatch:BehaviorForaging             -2.57      0.41    -3.39    -1.80 1.00
## BehLocPatch:BehaviorTransit               2.40      0.25     1.91     2.89 1.00
## BehLocPatch:NameCanis_latrans             0.69      0.85    -0.92     2.40 1.00
## BehLocPatch:NameDasypus_novemcinctus     -0.18      0.34    -0.87     0.48 1.00
## BehLocPatch:NameDidelphis_virginiana      0.41      0.61    -0.77     1.63 1.00
## BehLocPatch:NameLynx_rufus                0.16      0.97    -1.69     2.06 1.00
## BehLocPatch:NameMeleagris_gallopavo       0.09      0.83    -1.51     1.75 1.00
## BehLocPatch:NameProcyon_lotor            -1.05      0.36    -1.74    -0.34 1.00
## BehLocPatch:NameSciurus_carolinensis     -0.74      0.82    -2.34     0.85 1.00
## BehLocPatch:NameSylvilagus_floridanus    -0.55      0.74    -2.01     0.90 1.00
##                                       Bulk_ESS Tail_ESS
## Intercept                                 3835     2554
## BehLocPatch                               4953     5415
## BehaviorForaging                          6439     5998
## BehaviorTransit                           5609     5745
## NameCanis_latrans                        10610     5626
## NameDasypus_novemcinctus                  8792     6741
## NameDidelphis_virginiana                  8351     5673
## NameLynx_rufus                           12299     5115
## NameMeleagris_gallopavo                  10849     6072
## NameProcyon_lotor                         7715     6484
## NameSciurus_carolinensis                 11266     6046
## NameSylvilagus_floridanus                11007     6259
## BehLocPatch:BehaviorForaging              6190     6138
## BehLocPatch:BehaviorTransit               5063     5139
## BehLocPatch:NameCanis_latrans            10619     5583
## BehLocPatch:NameDasypus_novemcinctus      8931     6306
## BehLocPatch:NameDidelphis_virginiana      7847     5839
## BehLocPatch:NameLynx_rufus               12768     5468
## BehLocPatch:NameMeleagris_gallopavo      11479     6512
## BehLocPatch:NameProcyon_lotor             8192     5929
## BehLocPatch:NameSciurus_carolinensis     11130     5341
## BehLocPatch:NameSylvilagus_floridanus    11391     6059
## 
## 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`?

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