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 = 3, 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: 3 chains, each with iter = 5000; warmup = 3000; thin = 1;
##          total post-warmup draws = 6000
## 
## 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.22     0.09     0.91 1.00     2635     3666
## 
## ~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.52 1.00     2098     3223
## 
## Regression Coefficients:
##                                         Estimate Est.Error l-95% CI u-95% CI
## Intercept                                   5.12      0.55     4.08     6.24
## StatusInvaded                              -4.68      0.51    -5.73    -3.72
## BehaviorForaging                            0.65      0.61    -0.55     1.88
## BehaviorTransit                             0.83      0.55    -0.22     1.91
## NameCanis_latrans                           1.23      0.80    -0.31     2.81
## NameDasypus_novemcinctus                    0.10      0.69    -1.26     1.45
## NameDidelphis_virginiana                    0.29      0.73    -1.14     1.72
## NameLynx_rufus                              1.16      0.81    -0.43     2.74
## NameMeleagris_gallopavo                     0.27      0.79    -1.30     1.79
## NameProcyon_lotor                           0.19      0.70    -1.14     1.57
## NameSciurus_carolinensis                    0.62      0.75    -0.85     2.11
## NameSylvilagus_floridanus                   0.44      0.74    -0.96     1.91
## StatusInvaded:BehaviorForaging             -0.52      0.61    -1.74     0.68
## StatusInvaded:BehaviorTransit              -1.34      0.56    -2.44    -0.28
## StatusInvaded:NameCanis_latrans             1.21      0.82    -0.36     2.85
## StatusInvaded:NameDasypus_novemcinctus     -0.03      0.69    -1.36     1.32
## StatusInvaded:NameDidelphis_virginiana      0.26      0.74    -1.18     1.72
## StatusInvaded:NameLynx_rufus                1.16      0.83    -0.46     2.84
## StatusInvaded:NameMeleagris_gallopavo       0.24      0.80    -1.31     1.83
## StatusInvaded:NameProcyon_lotor             0.05      0.70    -1.30     1.37
## StatusInvaded:NameSciurus_carolinensis      0.61      0.76    -0.89     2.12
## StatusInvaded:NameSylvilagus_floridanus     0.44      0.75    -1.01     1.85
##                                         Rhat Bulk_ESS Tail_ESS
## Intercept                               1.00     4228     4257
## StatusInvaded                           1.00     4667     4507
## BehaviorForaging                        1.00     3857     3925
## BehaviorTransit                         1.00     3650     3891
## NameCanis_latrans                       1.00     7205     4625
## NameDasypus_novemcinctus                1.00     4558     4035
## NameDidelphis_virginiana                1.00     5012     4435
## NameLynx_rufus                          1.00     7005     5059
## NameMeleagris_gallopavo                 1.00     6612     4741
## NameProcyon_lotor                       1.00     4511     4103
## NameSciurus_carolinensis                1.00     5186     4582
## NameSylvilagus_floridanus               1.00     4803     4634
## StatusInvaded:BehaviorForaging          1.00     3826     3673
## StatusInvaded:BehaviorTransit           1.00     3582     3879
## StatusInvaded:NameCanis_latrans         1.00     7425     4861
## StatusInvaded:NameDasypus_novemcinctus  1.00     4557     3978
## StatusInvaded:NameDidelphis_virginiana  1.00     5221     4238
## StatusInvaded:NameLynx_rufus            1.00     6789     4366
## StatusInvaded:NameMeleagris_gallopavo   1.00     6504     4557
## StatusInvaded:NameProcyon_lotor         1.00     4592     4207
## StatusInvaded:NameSciurus_carolinensis  1.00     5421     4351
## StatusInvaded:NameSylvilagus_floridanus 1.00     5085     4667
## 
## 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 = 3, cores = 4,
  control = list(adapt_delta = 0.99, max_treedepth = 15)
)

## Compiling Stan program...

## Start sampling

# Summarize model
summary(Loc_model)

##  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: 3 chains, each with iter = 5000; warmup = 3000; thin = 1;
##          total post-warmup draws = 6000
## 
## 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.20     0.01     0.74 1.00     1595     1748
## 
## ~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.47 1.00     2045     2502
## 
## Regression Coefficients:
##                                       Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                                 1.31      0.24     0.77     1.76 1.00
## BehLocPatch                              -2.04      0.24    -2.52    -1.59 1.00
## BehaviorForaging                          1.33      0.29     0.77     1.92 1.00
## BehaviorTransit                          -1.46      0.18    -1.81    -1.11 1.00
## NameCanis_latrans                         1.62      0.70     0.30     3.07 1.00
## NameDasypus_novemcinctus                  0.13      0.24    -0.34     0.60 1.00
## NameDidelphis_virginiana                  0.20      0.47    -0.72     1.12 1.00
## NameLynx_rufus                            1.78      0.65     0.59     3.13 1.00
## NameMeleagris_gallopavo                   0.25      0.71    -1.14     1.65 1.00
## NameProcyon_lotor                         0.59      0.25     0.10     1.09 1.00
## NameSciurus_carolinensis                  0.97      0.55    -0.08     2.08 1.00
## NameSylvilagus_floridanus                 0.73      0.54    -0.30     1.78 1.00
## BehLocPatch:BehaviorForaging             -2.58      0.41    -3.38    -1.79 1.00
## BehLocPatch:BehaviorTransit               2.41      0.25     1.90     2.91 1.00
## BehLocPatch:NameCanis_latrans             0.66      0.85    -0.97     2.38 1.00
## BehLocPatch:NameDasypus_novemcinctus     -0.19      0.33    -0.84     0.45 1.00
## BehLocPatch:NameDidelphis_virginiana      0.40      0.61    -0.78     1.63 1.00
## BehLocPatch:NameLynx_rufus                0.15      0.96    -1.65     2.04 1.00
## BehLocPatch:NameMeleagris_gallopavo       0.08      0.83    -1.55     1.73 1.00
## BehLocPatch:NameProcyon_lotor            -1.04      0.36    -1.75    -0.34 1.00
## BehLocPatch:NameSciurus_carolinensis     -0.72      0.79    -2.25     0.86 1.00
## BehLocPatch:NameSylvilagus_floridanus    -0.56      0.72    -1.99     0.84 1.00
##                                       Bulk_ESS Tail_ESS
## Intercept                                 2447     2401
## BehLocPatch                               2640     3378
## BehaviorForaging                          4152     4620
## BehaviorTransit                           3575     3327
## NameCanis_latrans                         5932     4357
## NameDasypus_novemcinctus                  4563     4623
## NameDidelphis_virginiana                  5427     4811
## NameLynx_rufus                            7620     4256
## NameMeleagris_gallopavo                   6271     4500
## NameProcyon_lotor                         4764     4633
## NameSciurus_carolinensis                  6190     4311
## NameSylvilagus_floridanus                 5864     4415
## BehLocPatch:BehaviorForaging              3652     3994
## BehLocPatch:BehaviorTransit               2696     3318
## BehLocPatch:NameCanis_latrans             6094     4707
## BehLocPatch:NameDasypus_novemcinctus      4835     4983
## BehLocPatch:NameDidelphis_virginiana      5468     4557
## BehLocPatch:NameLynx_rufus                6744     4586
## BehLocPatch:NameMeleagris_gallopavo       6543     4931
## BehLocPatch:NameProcyon_lotor             4856     3930
## BehLocPatch:NameSciurus_carolinensis      7099     4665
## BehLocPatch:NameSylvilagus_floridanus     5684     4621
## 
## 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