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, forcats, tidybayes, ggdist)

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

# Install cmdstanr
#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

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

Place holder in case we need to filter vegetation by frequency

# Total number of sites
total_sites <- nrow(CameraLoc)

# Function to filter data by frequency
filter_by_frequency <- function(df) {
  # Group data by species and calculate the frequency
  freq <- df %>%
    group_by(Species) %>%
    summarise(Frequency = n_distinct(Plot) / nrow(CameraLoc) * 100) %>%
    filter(Frequency >= 0)
  
  # Filter the original data to include only species with frequency >= XX%
  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

# 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 and turkeys with greater than 10 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 10
species_subset <- species_counts %>%
  filter(count > 10) %>%
  pull(Name)

# Filter CameraData to only include species with count > 10
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 later
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)

# 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()

dat <- O2_data %>%
  filter(!is.na(Behavior), !is.na(Status), !is.na(BehLoc)) %>%
  mutate(
    # time-of-day as proportion of a 24-hour day, then z-score and quadratic
    time_prop = (hour(DateTime)*3600 + minute(DateTime)*60 + second(DateTime)) / 86400,
    time_z = as.numeric(scale(time_prop)),
    time_z2 = time_z^2,
    # temperature standardized
    temp_z = as.numeric(scale(Air.TemperatureC)),
    # month factor to absorb temporal clustering
    Month = factor(month(DateTime)),
    # set factor baselines
    Behavior = factor(Behavior, ordered = FALSE),
    Behavior = relevel(Behavior, ref = "Local_Search"),
    Status   = factor(Status, levels = c("Non_Invaded","Invaded")),
    BehLoc   = factor(BehLoc, levels = c("Non_Patch","Patch")),  # ensure Non_Patch vs Patch
    Plot     = factor(Plot),
    Name     = factor(Name),
    Site     = factor(Site),
    Camera.Type = factor(Camera.Type)
  )

Behavior Model

# Priors (weakly informative, regularizing)
priors <- c(
  prior(normal(0, 1), class = "b", dpar = "muForaging"),
  prior(normal(0, 1), class = "b", dpar = "muTransit"),
  prior(normal(0, 2.5), class = "Intercept", dpar = "muForaging"),
  prior(normal(0, 2.5), class = "Intercept", dpar = "muTransit"),
  prior(exponential(1), class = "sd", dpar = "muForaging"),
  prior(exponential(1), class = "sd", dpar = "muTransit")
)

# 1) BETWEEN SITES: Do behaviors differ at invaded vs non-invaded sites?
# Multilevel multinomial (categorical) model with population-level effects and random intercepts
model_status <- brm(
  bf(Behavior ~ Status * Name +
       (1 | Site) + (1 | Camera.Type) + (1 | Month)),
  family = categorical(link = "logit", refcat = "Local_Search"),
  data = dat,
  prior = priors,
  chains = 4, cores = 4, iter = 4000, warmup = 2000,
  control = list(adapt_delta = 0.95, max_treedepth = 12)
)

## Compiling Stan program...

## Start sampling

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

# Quick diagnostics
summary(model_status)

## Warning: There were 12 divergent transitions after warmup. Increasing
## adapt_delta above 0.95 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

##  Family: categorical 
##   Links: muForaging = logit; muTransit = logit 
## Formula: Behavior ~ Status * Name + (1 | Site) + (1 | Camera.Type) + (1 | Month) 
##    Data: dat (Number of observations: 2815) 
##   Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Multilevel Hyperparameters:
## ~Camera.Type (Number of levels: 4) 
##                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept)     0.40      0.35     0.02     1.29 1.00     2803
## sd(muTransit_Intercept)      0.21      0.21     0.01     0.76 1.00     3526
##                          Tail_ESS
## sd(muForaging_Intercept)     3875
## sd(muTransit_Intercept)      4034
## 
## ~Month (Number of levels: 9) 
##                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept)     0.67      0.23     0.35     1.22 1.00     3114
## sd(muTransit_Intercept)      0.17      0.11     0.01     0.44 1.00     2445
##                          Tail_ESS
## sd(muForaging_Intercept)     4265
## sd(muTransit_Intercept)      3680
## 
## ~Site (Number of levels: 5) 
##                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept)     0.17      0.17     0.00     0.61 1.00     4102
## sd(muTransit_Intercept)      0.47      0.23     0.20     1.07 1.00     4653
##                          Tail_ESS
## sd(muForaging_Intercept)     4329
## sd(muTransit_Intercept)      5508
## 
## Regression Coefficients:
##                                                     Estimate Est.Error l-95% CI
## muForaging_Intercept                                   -0.51      0.53    -1.55
## muTransit_Intercept                                     1.51      0.42     0.69
## muForaging_StatusInvaded                               -0.21      0.44    -1.07
## muForaging_NameDasypus_novemcinctus                     1.84      0.54     0.81
## muForaging_NameDidelphis_virginiana                     0.89      0.73    -0.56
## muForaging_NameLynx_rufus                              -0.44      0.89    -2.26
## muForaging_NameMeleagris_gallopavo                      0.23      0.67    -1.09
## muForaging_NameOdocoileus_virginianus                   1.12      0.38     0.39
## muForaging_NameProcyon_lotor                           -0.71      0.55    -1.81
## muForaging_NameSciurus_carolinensis                     0.07      0.69    -1.28
## muForaging_NameSylvilagus_floridanus                    1.54      0.61     0.36
## muForaging_StatusInvaded:NameDasypus_novemcinctus       0.53      0.62    -0.69
## muForaging_StatusInvaded:NameDidelphis_virginiana       0.08      0.80    -1.51
## muForaging_StatusInvaded:NameLynx_rufus                -0.19      0.94    -2.04
## muForaging_StatusInvaded:NameMeleagris_gallopavo        0.12      0.87    -1.65
## muForaging_StatusInvaded:NameOdocoileus_virginianus    -0.14      0.45    -1.02
## muForaging_StatusInvaded:NameProcyon_lotor             -0.17      0.67    -1.47
## muForaging_StatusInvaded:NameSciurus_carolinensis       0.17      0.75    -1.32
## muForaging_StatusInvaded:NameSylvilagus_floridanus      0.42      0.74    -1.02
## muTransit_StatusInvaded                                 0.20      0.39    -0.56
## muTransit_NameDasypus_novemcinctus                      0.70      0.49    -0.26
## muTransit_NameDidelphis_virginiana                      0.21      0.66    -1.05
## muTransit_NameLynx_rufus                                1.47      0.73     0.08
## muTransit_NameMeleagris_gallopavo                       0.05      0.56    -1.02
## muTransit_NameOdocoileus_virginianus                   -0.17      0.30    -0.76
## muTransit_NameProcyon_lotor                            -0.49      0.41    -1.29
## muTransit_NameSciurus_carolinensis                     -0.86      0.62    -2.05
## muTransit_NameSylvilagus_floridanus                    -0.11      0.58    -1.23
## muTransit_StatusInvaded:NameDasypus_novemcinctus        0.01      0.59    -1.14
## muTransit_StatusInvaded:NameDidelphis_virginiana        0.36      0.74    -1.08
## muTransit_StatusInvaded:NameLynx_rufus                  0.62      0.87    -1.03
## muTransit_StatusInvaded:NameMeleagris_gallopavo        -0.24      0.79    -1.76
## muTransit_StatusInvaded:NameOdocoileus_virginianus     -0.26      0.40    -1.06
## muTransit_StatusInvaded:NameProcyon_lotor              -0.08      0.50    -1.08
## muTransit_StatusInvaded:NameSciurus_carolinensis       -0.49      0.68    -1.81
## muTransit_StatusInvaded:NameSylvilagus_floridanus       0.12      0.70    -1.24
##                                                     u-95% CI Rhat Bulk_ESS
## muForaging_Intercept                                    0.52 1.00     4346
## muTransit_Intercept                                     2.31 1.00     5564
## muForaging_StatusInvaded                                0.64 1.00     5923
## muForaging_NameDasypus_novemcinctus                     2.90 1.00     8038
## muForaging_NameDidelphis_virginiana                     2.27 1.00     9883
## muForaging_NameLynx_rufus                               1.27 1.00    12736
## muForaging_NameMeleagris_gallopavo                      1.55 1.00    11025
## muForaging_NameOdocoileus_virginianus                   1.86 1.00     6383
## muForaging_NameProcyon_lotor                            0.34 1.00     8589
## muForaging_NameSciurus_carolinensis                     1.43 1.00     9417
## muForaging_NameSylvilagus_floridanus                    2.77 1.00     8546
## muForaging_StatusInvaded:NameDasypus_novemcinctus       1.73 1.00     7915
## muForaging_StatusInvaded:NameDidelphis_virginiana       1.66 1.00    12312
## muForaging_StatusInvaded:NameLynx_rufus                 1.61 1.00    13458
## muForaging_StatusInvaded:NameMeleagris_gallopavo        1.77 1.00    13431
## muForaging_StatusInvaded:NameOdocoileus_virginianus     0.74 1.00     5934
## muForaging_StatusInvaded:NameProcyon_lotor              1.14 1.00     9001
## muForaging_StatusInvaded:NameSciurus_carolinensis       1.65 1.00    10506
## muForaging_StatusInvaded:NameSylvilagus_floridanus      1.88 1.00    10009
## muTransit_StatusInvaded                                 0.97 1.00     5088
## muTransit_NameDasypus_novemcinctus                      1.67 1.00     8280
## muTransit_NameDidelphis_virginiana                      1.49 1.00    10185
## muTransit_NameLynx_rufus                                2.93 1.00    10870
## muTransit_NameMeleagris_gallopavo                       1.18 1.00    10680
## muTransit_NameOdocoileus_virginianus                    0.40 1.00     6503
## muTransit_NameProcyon_lotor                             0.30 1.00     7606
## muTransit_NameSciurus_carolinensis                      0.36 1.00     8873
## muTransit_NameSylvilagus_floridanus                     1.07 1.00     9221
## muTransit_StatusInvaded:NameDasypus_novemcinctus        1.17 1.00     7186
## muTransit_StatusInvaded:NameDidelphis_virginiana        1.81 1.00     9982
## muTransit_StatusInvaded:NameLynx_rufus                  2.35 1.00    12094
## muTransit_StatusInvaded:NameMeleagris_gallopavo         1.34 1.00    11778
## muTransit_StatusInvaded:NameOdocoileus_virginianus      0.53 1.00     5041
## muTransit_StatusInvaded:NameProcyon_lotor               0.90 1.00     6334
## muTransit_StatusInvaded:NameSciurus_carolinensis        0.84 1.00     8259
## muTransit_StatusInvaded:NameSylvilagus_floridanus       1.53 1.00    10095
##                                                     Tail_ESS
## muForaging_Intercept                                    4482
## muTransit_Intercept                                     5452
## muForaging_StatusInvaded                                5919
## muForaging_NameDasypus_novemcinctus                     6254
## muForaging_NameDidelphis_virginiana                     5956
## muForaging_NameLynx_rufus                               6115
## muForaging_NameMeleagris_gallopavo                      6623
## muForaging_NameOdocoileus_virginianus                   5826
## muForaging_NameProcyon_lotor                            6222
## muForaging_NameSciurus_carolinensis                     6373
## muForaging_NameSylvilagus_floridanus                    6355
## muForaging_StatusInvaded:NameDasypus_novemcinctus       5804
## muForaging_StatusInvaded:NameDidelphis_virginiana       5954
## muForaging_StatusInvaded:NameLynx_rufus                 5907
## muForaging_StatusInvaded:NameMeleagris_gallopavo        6430
## muForaging_StatusInvaded:NameOdocoileus_virginianus     5510
## muForaging_StatusInvaded:NameProcyon_lotor              6610
## muForaging_StatusInvaded:NameSciurus_carolinensis       6028
## muForaging_StatusInvaded:NameSylvilagus_floridanus      6753
## muTransit_StatusInvaded                                 5106
## muTransit_NameDasypus_novemcinctus                      6625
## muTransit_NameDidelphis_virginiana                      5804
## muTransit_NameLynx_rufus                                5196
## muTransit_NameMeleagris_gallopavo                       6396
## muTransit_NameOdocoileus_virginianus                    6266
## muTransit_NameProcyon_lotor                             5925
## muTransit_NameSciurus_carolinensis                      5526
## muTransit_NameSylvilagus_floridanus                     6663
## muTransit_StatusInvaded:NameDasypus_novemcinctus        6404
## muTransit_StatusInvaded:NameDidelphis_virginiana        5879
## muTransit_StatusInvaded:NameLynx_rufus                  6394
## muTransit_StatusInvaded:NameMeleagris_gallopavo         5760
## muTransit_StatusInvaded:NameOdocoileus_virginianus      4954
## muTransit_StatusInvaded:NameProcyon_lotor               6174
## muTransit_StatusInvaded:NameSciurus_carolinensis        5541
## muTransit_StatusInvaded:NameSylvilagus_floridanus       6119
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

plot(model_status)

waic(model_status)

## Warning: 
## 5 (0.2%) p_waic estimates greater than 0.4. We recommend trying loo instead.

## 
## Computed from 8000 by 2815 log-likelihood matrix.
## 
##           Estimate   SE
## elpd_waic  -2539.2 31.4
## p_waic        43.6  1.9
## waic        5078.3 62.8
## 
## 5 (0.2%) p_waic estimates greater than 0.4. We recommend trying loo instead.

loo(model_status)

## 
## Computed from 8000 by 2815 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -2539.3 31.4
## p_loo        43.7  2.0
## looic      5078.6 62.9
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.6]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

Behavior Model Within cogongrass patches

# 2) WITHIN INVADED: Do behaviors differ inside vs outside cogongrass patches?
dat_inv <- dat %>% filter(Status == "Invaded")


# Priors for the within-invaded model (same structure, applies to both non-reference categories)
priors_inv <- c(
  prior(normal(0, 1), class = "b", dpar = "muForaging"),
  prior(normal(0, 1), class = "b", dpar = "muTransit"),
  prior(normal(0, 2.5), class = "Intercept", dpar = "muForaging"),
  prior(normal(0, 2.5), class = "Intercept", dpar = "muTransit"),
  prior(exponential(1), class = "sd", dpar = "muForaging"),
  prior(exponential(1), class = "sd", dpar = "muTransit")
)

model_loc <- brm(
  bf(Behavior ~ BehLoc * Name +
       (1 | Site) + (1 | Camera.Type) + (1 | Month)),
  family = categorical(link = "logit", refcat = "Local_Search"),
  data = dat_inv,
  prior = priors_inv,
  chains = 4, cores = 4, iter = 4000, warmup = 2000,
  control = list(adapt_delta = 0.95, max_treedepth = 12)
)

## 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

summary(model_loc)

## Warning: There were 1 divergent transitions after warmup. Increasing
## adapt_delta above 0.95 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

##  Family: categorical 
##   Links: muForaging = logit; muTransit = logit 
## Formula: Behavior ~ BehLoc * Name + (1 | Site) + (1 | Camera.Type) + (1 | Month) 
##    Data: dat_inv (Number of observations: 1300) 
##   Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Multilevel Hyperparameters:
## ~Camera.Type (Number of levels: 4) 
##                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept)     0.62      0.45     0.05     1.77 1.00     2711
## sd(muTransit_Intercept)      0.55      0.35     0.11     1.47 1.00     2916
##                          Tail_ESS
## sd(muForaging_Intercept)     2576
## sd(muTransit_Intercept)      2384
## 
## ~Month (Number of levels: 9) 
##                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept)     0.43      0.20     0.11     0.91 1.00     2145
## sd(muTransit_Intercept)      0.25      0.16     0.01     0.63 1.00     2448
##                          Tail_ESS
## sd(muForaging_Intercept)     1856
## sd(muTransit_Intercept)      3102
## 
## ~Site (Number of levels: 5) 
##                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept)     0.70      0.40     0.14     1.69 1.00     2963
## sd(muTransit_Intercept)      0.61      0.34     0.20     1.46 1.00     3087
##                          Tail_ESS
## sd(muForaging_Intercept)     2649
## sd(muTransit_Intercept)      4535
## 
## Regression Coefficients:
##                                                   Estimate Est.Error l-95% CI
## muForaging_Intercept                                 -0.03      0.70    -1.42
## muTransit_Intercept                                   0.93      0.60    -0.32
## muForaging_BehLocPatch                               -1.29      0.59    -2.46
## muForaging_NameDasypus_novemcinctus                   1.60      0.55     0.56
## muForaging_NameDidelphis_virginiana                   0.34      0.75    -1.14
## muForaging_NameLynx_rufus                            -0.59      0.87    -2.34
## muForaging_NameMeleagris_gallopavo                    0.11      0.85    -1.58
## muForaging_NameOdocoileus_virginianus                 0.27      0.44    -0.59
## muForaging_NameProcyon_lotor                         -1.23      0.57    -2.35
## muForaging_NameSciurus_carolinensis                  -0.14      0.67    -1.48
## muForaging_NameSylvilagus_floridanus                  1.06      0.68    -0.26
## muForaging_BehLocPatch:NameDasypus_novemcinctus      -0.32      0.80    -1.92
## muForaging_BehLocPatch:NameDidelphis_virginiana      -0.17      0.93    -2.04
## muForaging_BehLocPatch:NameLynx_rufus                -0.02      1.00    -1.99
## muForaging_BehLocPatch:NameMeleagris_gallopavo       -0.05      0.99    -2.02
## muForaging_BehLocPatch:NameOdocoileus_virginianus    -0.28      0.62    -1.47
## muForaging_BehLocPatch:NameProcyon_lotor             -0.14      0.92    -1.96
## muForaging_BehLocPatch:NameSciurus_carolinensis      -0.05      0.98    -1.97
## muForaging_BehLocPatch:NameSylvilagus_floridanus     -0.12      0.96    -2.07
## muTransit_BehLocPatch                                 1.84      0.50     0.89
## muTransit_NameDasypus_novemcinctus                    0.49      0.51    -0.49
## muTransit_NameDidelphis_virginiana                    0.16      0.67    -1.12
## muTransit_NameLynx_rufus                              1.47      0.74     0.08
## muTransit_NameMeleagris_gallopavo                    -0.32      0.80    -1.91
## muTransit_NameOdocoileus_virginianus                 -0.63      0.40    -1.42
## muTransit_NameProcyon_lotor                          -0.41      0.45    -1.29
## muTransit_NameSciurus_carolinensis                   -0.79      0.59    -1.94
## muTransit_NameSylvilagus_floridanus                  -0.12      0.64    -1.36
## muTransit_BehLocPatch:NameDasypus_novemcinctus       -0.36      0.71    -1.72
## muTransit_BehLocPatch:NameDidelphis_virginiana        0.53      0.89    -1.16
## muTransit_BehLocPatch:NameLynx_rufus                  0.05      0.99    -1.85
## muTransit_BehLocPatch:NameMeleagris_gallopavo         0.17      0.93    -1.65
## muTransit_BehLocPatch:NameOdocoileus_virginianus      0.09      0.52    -0.96
## muTransit_BehLocPatch:NameProcyon_lotor               0.41      0.75    -1.02
## muTransit_BehLocPatch:NameSciurus_carolinensis        0.13      0.94    -1.67
## muTransit_BehLocPatch:NameSylvilagus_floridanus       0.25      0.92    -1.51
##                                                   u-95% CI Rhat Bulk_ESS
## muForaging_Intercept                                  1.35 1.00     3261
## muTransit_Intercept                                   2.04 1.00     3259
## muForaging_BehLocPatch                               -0.14 1.00     5243
## muForaging_NameDasypus_novemcinctus                   2.68 1.00     5996
## muForaging_NameDidelphis_virginiana                   1.80 1.00     9393
## muForaging_NameLynx_rufus                             1.05 1.00     9712
## muForaging_NameMeleagris_gallopavo                    1.74 1.00     9396
## muForaging_NameOdocoileus_virginianus                 1.12 1.00     4860
## muForaging_NameProcyon_lotor                         -0.14 1.00     6259
## muForaging_NameSciurus_carolinensis                   1.14 1.00     8496
## muForaging_NameSylvilagus_floridanus                  2.41 1.00     8607
## muForaging_BehLocPatch:NameDasypus_novemcinctus       1.26 1.00     9559
## muForaging_BehLocPatch:NameDidelphis_virginiana       1.66 1.00    12085
## muForaging_BehLocPatch:NameLynx_rufus                 1.94 1.00    13532
## muForaging_BehLocPatch:NameMeleagris_gallopavo        1.88 1.00    12813
## muForaging_BehLocPatch:NameOdocoileus_virginianus     0.97 1.00     5296
## muForaging_BehLocPatch:NameProcyon_lotor              1.66 1.00    10540
## muForaging_BehLocPatch:NameSciurus_carolinensis       1.87 1.00    11686
## muForaging_BehLocPatch:NameSylvilagus_floridanus      1.73 1.00    12069
## muTransit_BehLocPatch                                 2.86 1.00     4861
## muTransit_NameDasypus_novemcinctus                    1.52 1.00     5387
## muTransit_NameDidelphis_virginiana                    1.55 1.00     7795
## muTransit_NameLynx_rufus                              2.99 1.00     9806
## muTransit_NameMeleagris_gallopavo                     1.26 1.00    10597
## muTransit_NameOdocoileus_virginianus                  0.14 1.00     4215
## muTransit_NameProcyon_lotor                           0.49 1.00     4772
## muTransit_NameSciurus_carolinensis                    0.38 1.00     7009
## muTransit_NameSylvilagus_floridanus                   1.18 1.00     7883
## muTransit_BehLocPatch:NameDasypus_novemcinctus        1.05 1.00     7602
## muTransit_BehLocPatch:NameDidelphis_virginiana        2.31 1.00    10611
## muTransit_BehLocPatch:NameLynx_rufus                  2.00 1.00    11829
## muTransit_BehLocPatch:NameMeleagris_gallopavo         2.00 1.00    12704
## muTransit_BehLocPatch:NameOdocoileus_virginianus      1.05 1.00     4916
## muTransit_BehLocPatch:NameProcyon_lotor               1.94 1.00     8505
## muTransit_BehLocPatch:NameSciurus_carolinensis        1.96 1.00    13409
## muTransit_BehLocPatch:NameSylvilagus_floridanus       2.11 1.00    11597
##                                                   Tail_ESS
## muForaging_Intercept                                  4632
## muTransit_Intercept                                   4297
## muForaging_BehLocPatch                                4892
## muForaging_NameDasypus_novemcinctus                   5961
## muForaging_NameDidelphis_virginiana                   6417
## muForaging_NameLynx_rufus                             6155
## muForaging_NameMeleagris_gallopavo                    6196
## muForaging_NameOdocoileus_virginianus                 5848
## muForaging_NameProcyon_lotor                          6095
## muForaging_NameSciurus_carolinensis                   6327
## muForaging_NameSylvilagus_floridanus                  6205
## muForaging_BehLocPatch:NameDasypus_novemcinctus       6029
## muForaging_BehLocPatch:NameDidelphis_virginiana       5809
## muForaging_BehLocPatch:NameLynx_rufus                 5983
## muForaging_BehLocPatch:NameMeleagris_gallopavo        5914
## muForaging_BehLocPatch:NameOdocoileus_virginianus     5649
## muForaging_BehLocPatch:NameProcyon_lotor              5740
## muForaging_BehLocPatch:NameSciurus_carolinensis       5283
## muForaging_BehLocPatch:NameSylvilagus_floridanus      6207
## muTransit_BehLocPatch                                 5160
## muTransit_NameDasypus_novemcinctus                    5846
## muTransit_NameDidelphis_virginiana                    5884
## muTransit_NameLynx_rufus                              5937
## muTransit_NameMeleagris_gallopavo                     5640
## muTransit_NameOdocoileus_virginianus                  5209
## muTransit_NameProcyon_lotor                           5735
## muTransit_NameSciurus_carolinensis                    6009
## muTransit_NameSylvilagus_floridanus                   6339
## muTransit_BehLocPatch:NameDasypus_novemcinctus        5953
## muTransit_BehLocPatch:NameDidelphis_virginiana        6187
## muTransit_BehLocPatch:NameLynx_rufus                  5861
## muTransit_BehLocPatch:NameMeleagris_gallopavo         6165
## muTransit_BehLocPatch:NameOdocoileus_virginianus      5566
## muTransit_BehLocPatch:NameProcyon_lotor               5732
## muTransit_BehLocPatch:NameSciurus_carolinensis        5716
## muTransit_BehLocPatch:NameSylvilagus_floridanus       5727
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

plot(model_loc)

waic(model_loc)

## Warning: 
## 7 (0.5%) p_waic estimates greater than 0.4. We recommend trying loo instead.

## 
## Computed from 8000 by 1300 log-likelihood matrix.
## 
##           Estimate   SE
## elpd_waic  -1008.8 24.7
## p_waic        37.1  1.9
## waic        2017.7 49.5
## 
## 7 (0.5%) p_waic estimates greater than 0.4. We recommend trying loo instead.

loo(model_loc)

## 
## Computed from 8000 by 1300 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1009.0 24.7
## p_loo        37.3  1.9
## looic      2018.0 49.5
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.9]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

Between Site Behavior Plots

# Prediction grid for invasion status
newdata <- expand_grid(
  Status = c("Non_Invaded", "Invaded"),
  Name = unique(dat$Name)  # Include all species to get marginal effect
) %>%
  mutate(
    Site = NA,
    Camera.Type = NA,
    Month = NA
  )

# Get posterior predictions
pred_draws <- add_epred_draws(
  newdata,
  model_status,
  re_formula = NA,  # Marginalizing over random effects
  ndraws = 2000
) %>%
  ungroup()

# Marginal predictions across species (average effect)
marginal_predictions <- pred_draws %>%
  group_by(Status, .draw, .category) %>%
  summarise(.epred = mean(.epred), .groups = "drop")

# Summary statistics
behavior_summary <- marginal_predictions %>%
  group_by(Status, .category) %>%
  median_qi(.epred, .width = 0.95) %>%
  mutate(
    Behavior = .category,
    Probability = .epred,
    Lower = .lower,
    Upper = .upper
  ) %>%
  dplyr::select(Status, Behavior, Probability, Lower, Upper)

# Reorder behaviors
behavior_summary <- behavior_summary %>%
  mutate(
    Behavior = factor(Behavior, 
                      levels = c("Local_Search", "Foraging", "Transit"),
                      labels = c("Local Search", "Foraging", "Transit")),
    Status = factor(Status, 
                    levels = c("Non_Invaded", "Invaded"),
                    labels = c("Non-Invaded", "Invaded"))
  )

# Color palette
status_colors <- c("Non-Invaded" = "#0072B2", "Invaded" = "#D55E00")

# Plot-level proportions
plot_proportions <- dat %>%
  group_by(Plot, Status, Behavior) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(Plot, Status) %>%
  mutate(
    total = sum(n),
    site_prob = n / total
  ) %>%
  ungroup() %>%
  mutate(
    Status = factor(Status, 
                    levels = c("Non_Invaded", "Invaded"),
                    labels = c("Non-Invaded", "Invaded")),
    Behavior = factor(Behavior, 
                      levels = c("Local_Search", "Foraging", "Transit"),
                      labels = c("Local Search", "Foraging", "Transit"))
  )

print(plot_proportions)

## # A tibble: 94 × 6
##    Plot  Status      Behavior         n total site_prob
##    <fct> <fct>       <fct>        <int> <int>     <dbl>
##  1 BI201 Invaded     Local Search     6    58    0.103 
##  2 BI201 Invaded     Foraging        13    58    0.224 
##  3 BI201 Invaded     Transit         39    58    0.672 
##  4 BN211 Non-Invaded Local Search     8    83    0.0964
##  5 BN211 Non-Invaded Foraging        21    83    0.253 
##  6 BN211 Non-Invaded Transit         54    83    0.651 
##  7 EI100 Invaded     Local Search     9    69    0.130 
##  8 EI100 Invaded     Foraging        17    69    0.246 
##  9 EI100 Invaded     Transit         43    69    0.623 
## 10 EI102 Invaded     Local Search    16    59    0.271 
## # ℹ 84 more rows

# Plot
p <- ggplot(behavior_summary, aes(x = Behavior, y = Probability, fill = Status)) +
  geom_col(position = position_dodge(width = 0.8), 
           width = 0.7,
           color = "black",
           linewidth = 0.3) +
  geom_errorbar(aes(ymin = Lower, ymax = Upper),
                position = position_dodge(width = 0.8),
                width = 0.25,
                linewidth = 0.5) +
  # Add plot-level data points
  geom_point(data = plot_proportions,
             aes(x = Behavior, y = site_prob, fill = Status),
             position = position_dodge(width = 0.8),
             size = 2.5,
             shape = 21,
             color = "black",
             alpha = 0.5,
             stroke = 0.8) +
  scale_fill_manual(values = status_colors,
                    name = "Site Status") +
  scale_y_continuous(
    limits = c(0, 1),
    breaks = seq(0, 1, 0.1),
    expand = expansion(mult = c(0, 0.02)),
    labels = scales::label_percent(accuracy = 1)
  ) +
  labs(
    x = "Behavior Category",
    y = "Predicted Probability",
    title = NULL
  ) +
  theme_classic(base_size = 12) +
  theme(
    axis.text.x = element_text(size = 11, color = "black"),
    axis.text.y = element_text(size = 10, color = "black"),
    axis.title = element_text(size = 12, face = "bold"),
    axis.line = element_line(linewidth = 0.5),
    axis.ticks = element_line(linewidth = 0.5),
    legend.position = c(0.10, 0.85),
    legend.title = element_text(size = 11, face = "bold"),
    legend.text = element_text(size = 10),
    legend.background = element_rect(fill = "white", color = "black", linewidth = 0.3),
    legend.key.size = unit(1.5, "cm"),
    panel.grid.major.y = element_line(color = "gray90", linewidth = 0.3),
    plot.margin = margin(10, 10, 10, 10)
  )

print(p)

ggsave("Figure_Behavior_by_Invasion_Status.png", 
       plot = p,
       width = 7, 
       height = 5, 
       dpi = 600,
       bg = "white")

Predicted Probabilities

# Create prediction grid for BehLoc
newdata_loc <- expand_grid(
  BehLoc = c("Non_Patch", "Patch"),
  Name = unique(dat_inv$Name)  # Include all species to get marginal effect
) %>%
  mutate(
    Site = NA,
    Camera.Type = NA,
    Month = NA
  )

# Posterior predictions
pred_draws_loc <- add_epred_draws(
  newdata_loc,
  model_loc,
  re_formula = NA,  # Marginalizing over random effects
  ndraws = 2000
) %>%
  ungroup()

marginal_predictions_loc <- pred_draws_loc %>%
  group_by(BehLoc, .draw, .category) %>%
  summarise(.epred = mean(.epred), .groups = "drop")

# Summary statistics
behavior_summary_loc <- marginal_predictions_loc %>%
  group_by(BehLoc, .category) %>%
  median_qi(.epred, .width = 0.95) %>%
  mutate(
    Behavior = .category,
    Probability = .epred,
    Lower = .lower,
    Upper = .upper
  ) %>%
  dplyr::select(BehLoc, Behavior, Probability, Lower, Upper)

# Reorder factors
behavior_summary_loc <- behavior_summary_loc %>%
  mutate(
    Behavior = factor(Behavior, 
                      levels = c("Local_Search", "Foraging", "Transit"),
                      labels = c("Local Search", "Foraging", "Transit")),
    BehLoc = factor(BehLoc, 
                    levels = c("Non_Patch", "Patch"),
                    labels = c("Non-Patch", "Patch"))
  )


# Plot-level proportions for overlay points

plot_proportions_loc <- dat_inv %>%
  group_by(Plot, BehLoc, Behavior) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(Plot, BehLoc) %>%
  mutate(
    total = sum(n),
    site_prob = n / total
  ) %>%
  ungroup() %>%
  mutate(
    BehLoc = factor(BehLoc, 
                    levels = c("Non_Patch", "Patch"),
                    labels = c("Non-Patch", "Patch")),
    Behavior = factor(Behavior, 
                      levels = c("Local_Search", "Foraging", "Transit"),
                      labels = c("Local Search", "Foraging", "Transit"))
  )


# Color palette
location_colors <- c("Non-Patch" = "#0072B2", "Patch" = "#D55E00")

p_loc <- ggplot(behavior_summary_loc, aes(x = Behavior, y = Probability, fill = BehLoc)) +
  geom_col(position = position_dodge(width = 0.8), 
           width = 0.7,
           color = "black",
           linewidth = 0.3) +
  geom_errorbar(aes(ymin = Lower, ymax = Upper),
                position = position_dodge(width = 0.8),
                width = 0.25,
                linewidth = 0.5) +
  # Add site-level data points
  geom_point(data = plot_proportions_loc,
             aes(x = Behavior, y = site_prob, fill = BehLoc),
             position = position_dodge(width = 0.8),
             size = 2.5,
             shape = 21,
             color = "black",
             alpha = 0.5,
             stroke = 0.8) +
  scale_fill_manual(values = location_colors,
                    name = "Location") +
  scale_y_continuous(
    limits = c(0, 1),
    breaks = seq(0, 1, 0.1),
    expand = expansion(mult = c(0, 0.02)),
    labels = scales::label_percent(accuracy = 1)
  ) +
  labs(
    x = "Behavior Category",
    y = "Predicted Probability",
    title = NULL
  ) +
  theme_classic(base_size = 12) +
  theme(
    axis.text.x = element_text(size = 11, color = "black"),
    axis.text.y = element_text(size = 10, color = "black"),
    axis.title = element_text(size = 12, face = "bold"),
    axis.line = element_line(linewidth = 0.5),
    axis.ticks = element_line(linewidth = 0.5),
    legend.position = c(0.14, 0.83),
    legend.title = element_text(size = 11, face = "bold"),
    legend.text = element_text(size = 10),
    legend.background = element_rect(fill = "white", color = "black", linewidth = 0.3),
    legend.key.size = unit(1.7, "cm"),
    panel.grid.major.y = element_line(color = "gray90", linewidth = 0.3),
    plot.margin = margin(10, 10, 10, 10)
  )

print(p_loc)

ggsave("Figure_Behavior_by_Location.png", 
       plot = p_loc,
       width = 7, 
       height = 5, 
       dpi = 600,
       bg = "white")

Probability of Direction

pd <- function(x) {
  max(mean(x > 0), mean(x < 0))
}


draws_loc <- as_draws_df(model_loc)
draws_status <- as_draws_df(model_status)

beta_loc <- draws_loc[, grep("^b_", colnames(draws_loc))]

## Warning: Dropping 'draws_df' class as required metadata was removed.

beta_status <- draws_status[, grep("^b_", colnames(draws_status))]

## Warning: Dropping 'draws_df' class as required metadata was removed.

pd_vals_loc <- apply(beta_loc, 2, pd)
beta_means_loc <- colMeans(beta_loc)

results_loc <- data.frame(
  parameter = colnames(beta_loc),
  mean = beta_means_loc,
  pd = pd_vals_loc,
  pd_95 = pd_vals_loc >= 0.95
)

head(results_loc)

##                                                                   parameter
## b_muForaging_Intercept                               b_muForaging_Intercept
## b_muTransit_Intercept                                 b_muTransit_Intercept
## b_muForaging_BehLocPatch                           b_muForaging_BehLocPatch
## b_muForaging_NameDasypus_novemcinctus b_muForaging_NameDasypus_novemcinctus
## b_muForaging_NameDidelphis_virginiana b_muForaging_NameDidelphis_virginiana
## b_muForaging_NameLynx_rufus                     b_muForaging_NameLynx_rufus
##                                              mean       pd pd_95
## b_muForaging_Intercept                -0.03080426 0.513000 FALSE
## b_muTransit_Intercept                  0.93393204 0.939625 FALSE
## b_muForaging_BehLocPatch              -1.28570639 0.986375  TRUE
## b_muForaging_NameDasypus_novemcinctus  1.60448011 0.999250  TRUE
## b_muForaging_NameDidelphis_virginiana  0.33537520 0.674750 FALSE
## b_muForaging_NameLynx_rufus           -0.59398428 0.742500 FALSE

pd_vals_status <- apply(beta_status, 2, pd)
beta_means_status <- colMeans(beta_status)

results_status <- data.frame(
  parameter = colnames(beta_status),
  mean = beta_means_status,
  pd = pd_vals_status,
  pd_95 = pd_vals_status >= 0.95
)

head(results_status)

##                                                                   parameter
## b_muForaging_Intercept                               b_muForaging_Intercept
## b_muTransit_Intercept                                 b_muTransit_Intercept
## b_muForaging_StatusInvaded                       b_muForaging_StatusInvaded
## b_muForaging_NameDasypus_novemcinctus b_muForaging_NameDasypus_novemcinctus
## b_muForaging_NameDidelphis_virginiana b_muForaging_NameDidelphis_virginiana
## b_muForaging_NameLynx_rufus                     b_muForaging_NameLynx_rufus
##                                             mean       pd pd_95
## b_muForaging_Intercept                -0.5083134 0.841625 FALSE
## b_muTransit_Intercept                  1.5135275 0.998625  TRUE
## b_muForaging_StatusInvaded            -0.2122971 0.686500 FALSE
## b_muForaging_NameDasypus_novemcinctus  1.8384012 1.000000  TRUE
## b_muForaging_NameDidelphis_virginiana  0.8868583 0.885875 FALSE
## b_muForaging_NameLynx_rufus           -0.4405670 0.683875 FALSE

results_loc$parameter <- gsub("^b_", "", results_loc$parameter)
results_status$parameter <- gsub("^b_", "", results_status$parameter)

print(results_loc)

##                                                                                             parameter
## b_muForaging_Intercept                                                           muForaging_Intercept
## b_muTransit_Intercept                                                             muTransit_Intercept
## b_muForaging_BehLocPatch                                                       muForaging_BehLocPatch
## b_muForaging_NameDasypus_novemcinctus                             muForaging_NameDasypus_novemcinctus
## b_muForaging_NameDidelphis_virginiana                             muForaging_NameDidelphis_virginiana
## b_muForaging_NameLynx_rufus                                                 muForaging_NameLynx_rufus
## b_muForaging_NameMeleagris_gallopavo                               muForaging_NameMeleagris_gallopavo
## b_muForaging_NameOdocoileus_virginianus                         muForaging_NameOdocoileus_virginianus
## b_muForaging_NameProcyon_lotor                                           muForaging_NameProcyon_lotor
## b_muForaging_NameSciurus_carolinensis                             muForaging_NameSciurus_carolinensis
## b_muForaging_NameSylvilagus_floridanus                           muForaging_NameSylvilagus_floridanus
## b_muForaging_BehLocPatch:NameDasypus_novemcinctus     muForaging_BehLocPatch:NameDasypus_novemcinctus
## b_muForaging_BehLocPatch:NameDidelphis_virginiana     muForaging_BehLocPatch:NameDidelphis_virginiana
## b_muForaging_BehLocPatch:NameLynx_rufus                         muForaging_BehLocPatch:NameLynx_rufus
## b_muForaging_BehLocPatch:NameMeleagris_gallopavo       muForaging_BehLocPatch:NameMeleagris_gallopavo
## b_muForaging_BehLocPatch:NameOdocoileus_virginianus muForaging_BehLocPatch:NameOdocoileus_virginianus
## b_muForaging_BehLocPatch:NameProcyon_lotor                   muForaging_BehLocPatch:NameProcyon_lotor
## b_muForaging_BehLocPatch:NameSciurus_carolinensis     muForaging_BehLocPatch:NameSciurus_carolinensis
## b_muForaging_BehLocPatch:NameSylvilagus_floridanus   muForaging_BehLocPatch:NameSylvilagus_floridanus
## b_muTransit_BehLocPatch                                                         muTransit_BehLocPatch
## b_muTransit_NameDasypus_novemcinctus                               muTransit_NameDasypus_novemcinctus
## b_muTransit_NameDidelphis_virginiana                               muTransit_NameDidelphis_virginiana
## b_muTransit_NameLynx_rufus                                                   muTransit_NameLynx_rufus
## b_muTransit_NameMeleagris_gallopavo                                 muTransit_NameMeleagris_gallopavo
## b_muTransit_NameOdocoileus_virginianus                           muTransit_NameOdocoileus_virginianus
## b_muTransit_NameProcyon_lotor                                             muTransit_NameProcyon_lotor
## b_muTransit_NameSciurus_carolinensis                               muTransit_NameSciurus_carolinensis
## b_muTransit_NameSylvilagus_floridanus                             muTransit_NameSylvilagus_floridanus
## b_muTransit_BehLocPatch:NameDasypus_novemcinctus       muTransit_BehLocPatch:NameDasypus_novemcinctus
## b_muTransit_BehLocPatch:NameDidelphis_virginiana       muTransit_BehLocPatch:NameDidelphis_virginiana
## b_muTransit_BehLocPatch:NameLynx_rufus                           muTransit_BehLocPatch:NameLynx_rufus
## b_muTransit_BehLocPatch:NameMeleagris_gallopavo         muTransit_BehLocPatch:NameMeleagris_gallopavo
## b_muTransit_BehLocPatch:NameOdocoileus_virginianus   muTransit_BehLocPatch:NameOdocoileus_virginianus
## b_muTransit_BehLocPatch:NameProcyon_lotor                     muTransit_BehLocPatch:NameProcyon_lotor
## b_muTransit_BehLocPatch:NameSciurus_carolinensis       muTransit_BehLocPatch:NameSciurus_carolinensis
## b_muTransit_BehLocPatch:NameSylvilagus_floridanus     muTransit_BehLocPatch:NameSylvilagus_floridanus
##                                                            mean       pd pd_95
## b_muForaging_Intercept                              -0.03080426 0.513000 FALSE
## b_muTransit_Intercept                                0.93393204 0.939625 FALSE
## b_muForaging_BehLocPatch                            -1.28570639 0.986375  TRUE
## b_muForaging_NameDasypus_novemcinctus                1.60448011 0.999250  TRUE
## b_muForaging_NameDidelphis_virginiana                0.33537520 0.674750 FALSE
## b_muForaging_NameLynx_rufus                         -0.59398428 0.742500 FALSE
## b_muForaging_NameMeleagris_gallopavo                 0.11261003 0.555375 FALSE
## b_muForaging_NameOdocoileus_virginianus              0.27248601 0.730125 FALSE
## b_muForaging_NameProcyon_lotor                      -1.22687431 0.987250  TRUE
## b_muForaging_NameSciurus_carolinensis               -0.14489733 0.577125 FALSE
## b_muForaging_NameSylvilagus_floridanus               1.06174356 0.940250 FALSE
## b_muForaging_BehLocPatch:NameDasypus_novemcinctus   -0.32451369 0.657375 FALSE
## b_muForaging_BehLocPatch:NameDidelphis_virginiana   -0.16719722 0.570375 FALSE
## b_muForaging_BehLocPatch:NameLynx_rufus             -0.01702145 0.500125 FALSE
## b_muForaging_BehLocPatch:NameMeleagris_gallopavo    -0.04922888 0.522375 FALSE
## b_muForaging_BehLocPatch:NameOdocoileus_virginianus -0.27906570 0.673875 FALSE
## b_muForaging_BehLocPatch:NameProcyon_lotor          -0.13568735 0.554250 FALSE
## b_muForaging_BehLocPatch:NameSciurus_carolinensis   -0.04544090 0.521125 FALSE
## b_muForaging_BehLocPatch:NameSylvilagus_floridanus  -0.12145032 0.547750 FALSE
## b_muTransit_BehLocPatch                              1.83700150 1.000000  TRUE
## b_muTransit_NameDasypus_novemcinctus                 0.48553058 0.827000 FALSE
## b_muTransit_NameDidelphis_virginiana                 0.15885013 0.587500 FALSE
## b_muTransit_NameLynx_rufus                           1.47436355 0.981000  TRUE
## b_muTransit_NameMeleagris_gallopavo                 -0.31625107 0.653750 FALSE
## b_muTransit_NameOdocoileus_virginianus              -0.62566676 0.946125 FALSE
## b_muTransit_NameProcyon_lotor                       -0.41329255 0.826250 FALSE
## b_muTransit_NameSciurus_carolinensis                -0.79313235 0.910250 FALSE
## b_muTransit_NameSylvilagus_floridanus               -0.12072143 0.574000 FALSE
## b_muTransit_BehLocPatch:NameDasypus_novemcinctus    -0.35919885 0.699375 FALSE
## b_muTransit_BehLocPatch:NameDidelphis_virginiana     0.52957471 0.714750 FALSE
## b_muTransit_BehLocPatch:NameLynx_rufus               0.04701014 0.517750 FALSE
## b_muTransit_BehLocPatch:NameMeleagris_gallopavo      0.17458212 0.575625 FALSE
## b_muTransit_BehLocPatch:NameOdocoileus_virginianus   0.08613186 0.570500 FALSE
## b_muTransit_BehLocPatch:NameProcyon_lotor            0.41040036 0.701250 FALSE
## b_muTransit_BehLocPatch:NameSciurus_carolinensis     0.12796945 0.554250 FALSE
## b_muTransit_BehLocPatch:NameSylvilagus_floridanus    0.24843130 0.603500 FALSE

print(results_status)

##                                                                                                 parameter
## b_muForaging_Intercept                                                               muForaging_Intercept
## b_muTransit_Intercept                                                                 muTransit_Intercept
## b_muForaging_StatusInvaded                                                       muForaging_StatusInvaded
## b_muForaging_NameDasypus_novemcinctus                                 muForaging_NameDasypus_novemcinctus
## b_muForaging_NameDidelphis_virginiana                                 muForaging_NameDidelphis_virginiana
## b_muForaging_NameLynx_rufus                                                     muForaging_NameLynx_rufus
## b_muForaging_NameMeleagris_gallopavo                                   muForaging_NameMeleagris_gallopavo
## b_muForaging_NameOdocoileus_virginianus                             muForaging_NameOdocoileus_virginianus
## b_muForaging_NameProcyon_lotor                                               muForaging_NameProcyon_lotor
## b_muForaging_NameSciurus_carolinensis                                 muForaging_NameSciurus_carolinensis
## b_muForaging_NameSylvilagus_floridanus                               muForaging_NameSylvilagus_floridanus
## b_muForaging_StatusInvaded:NameDasypus_novemcinctus     muForaging_StatusInvaded:NameDasypus_novemcinctus
## b_muForaging_StatusInvaded:NameDidelphis_virginiana     muForaging_StatusInvaded:NameDidelphis_virginiana
## b_muForaging_StatusInvaded:NameLynx_rufus                         muForaging_StatusInvaded:NameLynx_rufus
## b_muForaging_StatusInvaded:NameMeleagris_gallopavo       muForaging_StatusInvaded:NameMeleagris_gallopavo
## b_muForaging_StatusInvaded:NameOdocoileus_virginianus muForaging_StatusInvaded:NameOdocoileus_virginianus
## b_muForaging_StatusInvaded:NameProcyon_lotor                   muForaging_StatusInvaded:NameProcyon_lotor
## b_muForaging_StatusInvaded:NameSciurus_carolinensis     muForaging_StatusInvaded:NameSciurus_carolinensis
## b_muForaging_StatusInvaded:NameSylvilagus_floridanus   muForaging_StatusInvaded:NameSylvilagus_floridanus
## b_muTransit_StatusInvaded                                                         muTransit_StatusInvaded
## b_muTransit_NameDasypus_novemcinctus                                   muTransit_NameDasypus_novemcinctus
## b_muTransit_NameDidelphis_virginiana                                   muTransit_NameDidelphis_virginiana
## b_muTransit_NameLynx_rufus                                                       muTransit_NameLynx_rufus
## b_muTransit_NameMeleagris_gallopavo                                     muTransit_NameMeleagris_gallopavo
## b_muTransit_NameOdocoileus_virginianus                               muTransit_NameOdocoileus_virginianus
## b_muTransit_NameProcyon_lotor                                                 muTransit_NameProcyon_lotor
## b_muTransit_NameSciurus_carolinensis                                   muTransit_NameSciurus_carolinensis
## b_muTransit_NameSylvilagus_floridanus                                 muTransit_NameSylvilagus_floridanus
## b_muTransit_StatusInvaded:NameDasypus_novemcinctus       muTransit_StatusInvaded:NameDasypus_novemcinctus
## b_muTransit_StatusInvaded:NameDidelphis_virginiana       muTransit_StatusInvaded:NameDidelphis_virginiana
## b_muTransit_StatusInvaded:NameLynx_rufus                           muTransit_StatusInvaded:NameLynx_rufus
## b_muTransit_StatusInvaded:NameMeleagris_gallopavo         muTransit_StatusInvaded:NameMeleagris_gallopavo
## b_muTransit_StatusInvaded:NameOdocoileus_virginianus   muTransit_StatusInvaded:NameOdocoileus_virginianus
## b_muTransit_StatusInvaded:NameProcyon_lotor                     muTransit_StatusInvaded:NameProcyon_lotor
## b_muTransit_StatusInvaded:NameSciurus_carolinensis       muTransit_StatusInvaded:NameSciurus_carolinensis
## b_muTransit_StatusInvaded:NameSylvilagus_floridanus     muTransit_StatusInvaded:NameSylvilagus_floridanus
##                                                              mean       pd
## b_muForaging_Intercept                                -0.50831343 0.841625
## b_muTransit_Intercept                                  1.51352753 0.998625
## b_muForaging_StatusInvaded                            -0.21229710 0.686500
## b_muForaging_NameDasypus_novemcinctus                  1.83840119 1.000000
## b_muForaging_NameDidelphis_virginiana                  0.88685830 0.885875
## b_muForaging_NameLynx_rufus                           -0.44056697 0.683875
## b_muForaging_NameMeleagris_gallopavo                   0.23434054 0.637625
## b_muForaging_NameOdocoileus_virginianus                1.11766650 0.998750
## b_muForaging_NameProcyon_lotor                        -0.70880579 0.910375
## b_muForaging_NameSciurus_carolinensis                  0.07280573 0.541500
## b_muForaging_NameSylvilagus_floridanus                 1.53769281 0.993875
## b_muForaging_StatusInvaded:NameDasypus_novemcinctus    0.52815691 0.806875
## b_muForaging_StatusInvaded:NameDidelphis_virginiana    0.07727428 0.542250
## b_muForaging_StatusInvaded:NameLynx_rufus             -0.19125876 0.572750
## b_muForaging_StatusInvaded:NameMeleagris_gallopavo     0.12088302 0.559250
## b_muForaging_StatusInvaded:NameOdocoileus_virginianus -0.13828183 0.618875
## b_muForaging_StatusInvaded:NameProcyon_lotor          -0.16999033 0.597000
## b_muForaging_StatusInvaded:NameSciurus_carolinensis    0.16662018 0.590500
## b_muForaging_StatusInvaded:NameSylvilagus_floridanus   0.42200810 0.718500
## b_muTransit_StatusInvaded                              0.19764191 0.693250
## b_muTransit_NameDasypus_novemcinctus                   0.69934281 0.924000
## b_muTransit_NameDidelphis_virginiana                   0.20634420 0.616875
## b_muTransit_NameLynx_rufus                             1.46537849 0.982625
## b_muTransit_NameMeleagris_gallopavo                    0.04601522 0.522625
## b_muTransit_NameOdocoileus_virginianus                -0.16757613 0.713250
## b_muTransit_NameProcyon_lotor                         -0.49080071 0.885000
## b_muTransit_NameSciurus_carolinensis                  -0.86039418 0.914250
## b_muTransit_NameSylvilagus_floridanus                 -0.10988265 0.582250
## b_muTransit_StatusInvaded:NameDasypus_novemcinctus     0.01118142 0.501500
## b_muTransit_StatusInvaded:NameDidelphis_virginiana     0.36218959 0.689750
## b_muTransit_StatusInvaded:NameLynx_rufus               0.61847830 0.758375
## b_muTransit_StatusInvaded:NameMeleagris_gallopavo     -0.24241798 0.626000
## b_muTransit_StatusInvaded:NameOdocoileus_virginianus  -0.25833187 0.738750
## b_muTransit_StatusInvaded:NameProcyon_lotor           -0.07883214 0.557250
## b_muTransit_StatusInvaded:NameSciurus_carolinensis    -0.48509744 0.758125
## b_muTransit_StatusInvaded:NameSylvilagus_floridanus    0.11725317 0.561500
##                                                       pd_95
## b_muForaging_Intercept                                FALSE
## b_muTransit_Intercept                                  TRUE
## b_muForaging_StatusInvaded                            FALSE
## b_muForaging_NameDasypus_novemcinctus                  TRUE
## b_muForaging_NameDidelphis_virginiana                 FALSE
## b_muForaging_NameLynx_rufus                           FALSE
## b_muForaging_NameMeleagris_gallopavo                  FALSE
## b_muForaging_NameOdocoileus_virginianus                TRUE
## b_muForaging_NameProcyon_lotor                        FALSE
## b_muForaging_NameSciurus_carolinensis                 FALSE
## b_muForaging_NameSylvilagus_floridanus                 TRUE
## b_muForaging_StatusInvaded:NameDasypus_novemcinctus   FALSE
## b_muForaging_StatusInvaded:NameDidelphis_virginiana   FALSE
## b_muForaging_StatusInvaded:NameLynx_rufus             FALSE
## b_muForaging_StatusInvaded:NameMeleagris_gallopavo    FALSE
## b_muForaging_StatusInvaded:NameOdocoileus_virginianus FALSE
## b_muForaging_StatusInvaded:NameProcyon_lotor          FALSE
## b_muForaging_StatusInvaded:NameSciurus_carolinensis   FALSE
## b_muForaging_StatusInvaded:NameSylvilagus_floridanus  FALSE
## b_muTransit_StatusInvaded                             FALSE
## b_muTransit_NameDasypus_novemcinctus                  FALSE
## b_muTransit_NameDidelphis_virginiana                  FALSE
## b_muTransit_NameLynx_rufus                             TRUE
## b_muTransit_NameMeleagris_gallopavo                   FALSE
## b_muTransit_NameOdocoileus_virginianus                FALSE
## b_muTransit_NameProcyon_lotor                         FALSE
## b_muTransit_NameSciurus_carolinensis                  FALSE
## b_muTransit_NameSylvilagus_floridanus                 FALSE
## b_muTransit_StatusInvaded:NameDasypus_novemcinctus    FALSE
## b_muTransit_StatusInvaded:NameDidelphis_virginiana    FALSE
## b_muTransit_StatusInvaded:NameLynx_rufus              FALSE
## b_muTransit_StatusInvaded:NameMeleagris_gallopavo     FALSE
## b_muTransit_StatusInvaded:NameOdocoileus_virginianus  FALSE
## b_muTransit_StatusInvaded:NameProcyon_lotor           FALSE
## b_muTransit_StatusInvaded:NameSciurus_carolinensis    FALSE
## b_muTransit_StatusInvaded:NameSylvilagus_floridanus   FALSE

Species Specific Site Effects

common_names <- c(
  "Canis_latrans" = "Coyote",
  "Odocoileus_virginianus" = "White-Tailed Deer",
  "Dasypus_novemcinctus" = "Nine-Banded Armadillo",
  "Procyon_lotor" = "Raccoon",
  "Lynx_rufus" = "Bobcat",
  "Didelphis_virginiana" = "Virginia Opossum",
  "Sylvilagus_floridanus" = "Eastern Cottontail Rabbit",
  "Meleagris_gallopavo" = "Wild Turkey",
  "Sciurus_carolinensis" = "Gray Squirrel"
)

newdata_sp <- expand_grid(
  Status = c("Non_Invaded", "Invaded"),
  Name = unique(dat$Name)
) %>%
  mutate(
    Site = NA,
    Camera.Type = NA,
    Month = NA
  )

pred_draws_sp <- add_epred_draws(
  newdata_sp,
  model_status,
  re_formula = ~(1 | Name),   # ← ONLY species RE
  ndraws = 2000
)


behavior_summary_sp <- pred_draws_sp %>%
  group_by(Name, Status, .category) %>%
  median_qi(.epred, .width = 0.95) %>%
  mutate(
    Behavior = factor(.category,
                      levels = c("Local_Search", "Foraging", "Transit"),
                      labels = c("Local Search", "Foraging", "Transit")),
    Status = factor(Status,
                    levels = c("Non_Invaded", "Invaded"),
                    labels = c("Non-Invaded", "Invaded"))
  )

behavior_summary_sp <- behavior_summary_sp %>%
  mutate(
    Name = recode(Name, !!!common_names),
    Name = factor(Name, levels = common_names)
  )


p_species <- ggplot(behavior_summary_sp,
                    aes(x = Behavior, y = .epred, fill = Status)) +
  geom_col(position = position_dodge(0.8),
           width = 0.7,
           color = "black") +
  geom_errorbar(aes(ymin = .lower, ymax = .upper),
                position = position_dodge(0.8),
                width = 0.25) +
  facet_wrap(~ Name) +
  scale_y_continuous(limits = c(0,1),
                     labels = scales::label_percent()) +
  labs(
  x = "Behavior Category",
  y = "Behavior Probability"
  ) +

  theme_classic()

print(p_species)

ggsave("Figure_Behavior_by_Species.png", 
       plot = p_species,
       width = 7, 
       height = 5, 
       dpi = 600,
       bg = "white")

Species Specific Within Effects

newdata_loc_sp <- expand_grid(
  BehLoc = c("Non_Patch", "Patch"),
  Name = unique(dat_inv$Name)
) %>%
  mutate(
    Site = NA,
    Camera.Type = NA,
    Month = NA
  )

pred_draws_loc_sp <- add_epred_draws(
  newdata_loc_sp,
  model_loc,
  re_formula = ~(1 | Name),   # species RE
  ndraws = 2000
)


behavior_summary_loc_sp <- pred_draws_loc_sp %>%
  group_by(Name, BehLoc, .category) %>%
  median_qi(.epred, .width = 0.95) %>%
  mutate(
    Behavior = factor(.category,
                      levels = c("Local_Search", "Foraging", "Transit"),
                      labels = c("Local Search", "Foraging", "Transit")),
    BehLoc = factor(BehLoc,
                    levels = c("Non_Patch", "Patch"),
                    labels = c("Non-Patch", "Patch"))
  )

behavior_summary_loc_sp <- behavior_summary_loc_sp %>%
  mutate(
    Name = recode(Name, !!!common_names),
    Name = factor(Name, levels = common_names)
  )


plot_proportions_loc_sp <- dat_inv %>%
  group_by(Name, Plot, BehLoc, Behavior) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(Name, Plot, BehLoc) %>%
  mutate(site_prob = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    BehLoc = factor(BehLoc,
                    levels = c("Non_Patch", "Patch"),
                    labels = c("Non-Patch", "Patch")),
    Behavior = factor(Behavior,
                      levels = c("Local_Search", "Foraging", "Transit"),
                      labels = c("Local Search", "Foraging", "Transit"))
  )

location_colors <- c("Non-Patch" = "#0072B2",
                     "Patch" = "#D55E00")

p_loc_species <- ggplot(behavior_summary_loc_sp,
                        aes(x = Behavior, y = .epred, fill = BehLoc)) +
  geom_col(position = position_dodge(0.8),
           width = 0.7,
           color = "black") +
  geom_errorbar(aes(ymin = .lower, ymax = .upper),
                position = position_dodge(0.8),
                width = 0.25) +
# geom_point(data = plot_proportions_loc_sp,
#            aes(x = Behavior, y = site_prob, fill = BehLoc),
#             position = position_dodge(0.8),
#             size = 2.5,
#             shape = 21,
#             color = "black",
#             alpha = 0.5,
#             stroke = 0.8) +
  facet_wrap(~ Name) +
  scale_fill_manual(values = location_colors,
                    name = "Location") +
  scale_y_continuous(limits = c(0,1),
                     labels = scales::label_percent()) +
  labs(
  x = "Behavior Category",
  y = "Behavior Probability"
  ) +
  theme_classic(base_size = 12)

print(p_loc_species)

ggsave("Figure_Behavior_by_Species_Within.png", 
       plot = p_loc_species,
       width = 7, 
       height = 5, 
       dpi = 600,
       bg = "white")

C4_O2_Behavior

Alan Ivory

2026-02-19