C4_O2_Behavior
Alan Ivory
2026-03-24
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, knitr, kableExtra, officer, flextable, purrr)
# 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 rowsVegetation 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.2742307Objective 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")), # 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 9 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 9 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.34 0.02 1.30 1.00 2600
## sd(muTransit_Intercept) 0.21 0.21 0.01 0.77 1.00 3551
## Tail_ESS
## sd(muForaging_Intercept) 3660
## sd(muTransit_Intercept) 4692
##
## ~Month (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 0.66 0.22 0.34 1.21 1.00 3410
## sd(muTransit_Intercept) 0.17 0.11 0.01 0.43 1.00 2431
## Tail_ESS
## sd(muForaging_Intercept) 5074
## sd(muTransit_Intercept) 3422
##
## ~Site (Number of levels: 5)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 0.16 0.17 0.00 0.58 1.00 3974
## sd(muTransit_Intercept) 0.48 0.25 0.20 1.14 1.00 4322
## Tail_ESS
## sd(muForaging_Intercept) 5026
## sd(muTransit_Intercept) 4970
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI
## muForaging_Intercept -0.50 0.52 -1.51
## muTransit_Intercept 1.51 0.42 0.66
## muForaging_StatusInvaded -0.21 0.44 -1.07
## muForaging_NameDasypus_novemcinctus 1.84 0.53 0.80
## muForaging_NameDidelphis_virginiana 0.90 0.72 -0.50
## muForaging_NameLynx_rufus -0.44 0.89 -2.24
## muForaging_NameMeleagris_gallopavo 0.22 0.68 -1.14
## muForaging_NameOdocoileus_virginianus 1.12 0.37 0.40
## muForaging_NameProcyon_lotor -0.71 0.55 -1.82
## muForaging_NameSciurus_carolinensis 0.08 0.69 -1.31
## muForaging_NameSylvilagus_floridanus 1.54 0.62 0.35
## muForaging_StatusInvaded:NameDasypus_novemcinctus 0.53 0.61 -0.65
## muForaging_StatusInvaded:NameDidelphis_virginiana 0.06 0.79 -1.50
## muForaging_StatusInvaded:NameLynx_rufus -0.19 0.96 -2.10
## muForaging_StatusInvaded:NameMeleagris_gallopavo 0.13 0.85 -1.56
## muForaging_StatusInvaded:NameOdocoileus_virginianus -0.14 0.45 -1.03
## muForaging_StatusInvaded:NameProcyon_lotor -0.17 0.67 -1.50
## muForaging_StatusInvaded:NameSciurus_carolinensis 0.16 0.77 -1.34
## muForaging_StatusInvaded:NameSylvilagus_floridanus 0.43 0.75 -1.01
## muTransit_StatusInvaded 0.19 0.38 -0.55
## muTransit_NameDasypus_novemcinctus 0.69 0.49 -0.23
## muTransit_NameDidelphis_virginiana 0.22 0.66 -1.03
## muTransit_NameLynx_rufus 1.45 0.74 0.07
## muTransit_NameMeleagris_gallopavo 0.03 0.57 -1.07
## muTransit_NameOdocoileus_virginianus -0.17 0.29 -0.75
## muTransit_NameProcyon_lotor -0.50 0.41 -1.29
## muTransit_NameSciurus_carolinensis -0.87 0.60 -2.02
## muTransit_NameSylvilagus_floridanus -0.11 0.58 -1.22
## muTransit_StatusInvaded:NameDasypus_novemcinctus 0.03 0.58 -1.12
## muTransit_StatusInvaded:NameDidelphis_virginiana 0.35 0.74 -1.07
## muTransit_StatusInvaded:NameLynx_rufus 0.63 0.87 -1.04
## muTransit_StatusInvaded:NameMeleagris_gallopavo -0.22 0.80 -1.77
## muTransit_StatusInvaded:NameOdocoileus_virginianus -0.25 0.39 -1.01
## muTransit_StatusInvaded:NameProcyon_lotor -0.06 0.50 -1.03
## muTransit_StatusInvaded:NameSciurus_carolinensis -0.49 0.69 -1.86
## muTransit_StatusInvaded:NameSylvilagus_floridanus 0.13 0.71 -1.23
## u-95% CI Rhat Bulk_ESS
## muForaging_Intercept 0.51 1.00 5088
## muTransit_Intercept 2.32 1.00 4381
## muForaging_StatusInvaded 0.65 1.00 6109
## muForaging_NameDasypus_novemcinctus 2.88 1.00 7577
## muForaging_NameDidelphis_virginiana 2.32 1.00 11031
## muForaging_NameLynx_rufus 1.28 1.00 14124
## muForaging_NameMeleagris_gallopavo 1.54 1.00 11467
## muForaging_NameOdocoileus_virginianus 1.85 1.00 7063
## muForaging_NameProcyon_lotor 0.38 1.00 8418
## muForaging_NameSciurus_carolinensis 1.40 1.00 9721
## muForaging_NameSylvilagus_floridanus 2.78 1.00 9033
## muForaging_StatusInvaded:NameDasypus_novemcinctus 1.72 1.00 7514
## muForaging_StatusInvaded:NameDidelphis_virginiana 1.59 1.00 12435
## muForaging_StatusInvaded:NameLynx_rufus 1.67 1.00 13650
## muForaging_StatusInvaded:NameMeleagris_gallopavo 1.81 1.00 13873
## muForaging_StatusInvaded:NameOdocoileus_virginianus 0.74 1.00 6067
## muForaging_StatusInvaded:NameProcyon_lotor 1.13 1.00 9111
## muForaging_StatusInvaded:NameSciurus_carolinensis 1.65 1.00 10217
## muForaging_StatusInvaded:NameSylvilagus_floridanus 1.90 1.00 10561
## muTransit_StatusInvaded 0.94 1.00 5097
## muTransit_NameDasypus_novemcinctus 1.68 1.00 7850
## muTransit_NameDidelphis_virginiana 1.56 1.00 10351
## muTransit_NameLynx_rufus 2.97 1.00 11449
## muTransit_NameMeleagris_gallopavo 1.16 1.00 10260
## muTransit_NameOdocoileus_virginianus 0.40 1.00 6702
## muTransit_NameProcyon_lotor 0.31 1.00 7424
## muTransit_NameSciurus_carolinensis 0.33 1.00 9088
## muTransit_NameSylvilagus_floridanus 1.06 1.00 10286
## muTransit_StatusInvaded:NameDasypus_novemcinctus 1.17 1.00 7673
## muTransit_StatusInvaded:NameDidelphis_virginiana 1.80 1.00 10624
## muTransit_StatusInvaded:NameLynx_rufus 2.36 1.00 12496
## muTransit_StatusInvaded:NameMeleagris_gallopavo 1.34 1.00 11833
## muTransit_StatusInvaded:NameOdocoileus_virginianus 0.52 1.00 5243
## muTransit_StatusInvaded:NameProcyon_lotor 0.90 1.00 6921
## muTransit_StatusInvaded:NameSciurus_carolinensis 0.85 1.00 9260
## muTransit_StatusInvaded:NameSylvilagus_floridanus 1.54 1.00 9303
## Tail_ESS
## muForaging_Intercept 5540
## muTransit_Intercept 5056
## muForaging_StatusInvaded 5252
## muForaging_NameDasypus_novemcinctus 5989
## muForaging_NameDidelphis_virginiana 6645
## muForaging_NameLynx_rufus 5488
## muForaging_NameMeleagris_gallopavo 6405
## muForaging_NameOdocoileus_virginianus 5964
## muForaging_NameProcyon_lotor 6516
## muForaging_NameSciurus_carolinensis 6523
## muForaging_NameSylvilagus_floridanus 6429
## muForaging_StatusInvaded:NameDasypus_novemcinctus 6640
## muForaging_StatusInvaded:NameDidelphis_virginiana 6811
## muForaging_StatusInvaded:NameLynx_rufus 5128
## muForaging_StatusInvaded:NameMeleagris_gallopavo 6050
## muForaging_StatusInvaded:NameOdocoileus_virginianus 5681
## muForaging_StatusInvaded:NameProcyon_lotor 5845
## muForaging_StatusInvaded:NameSciurus_carolinensis 6014
## muForaging_StatusInvaded:NameSylvilagus_floridanus 6569
## muTransit_StatusInvaded 5483
## muTransit_NameDasypus_novemcinctus 6297
## muTransit_NameDidelphis_virginiana 6245
## muTransit_NameLynx_rufus 5884
## muTransit_NameMeleagris_gallopavo 6385
## muTransit_NameOdocoileus_virginianus 6073
## muTransit_NameProcyon_lotor 6354
## muTransit_NameSciurus_carolinensis 6373
## muTransit_NameSylvilagus_floridanus 6347
## muTransit_StatusInvaded:NameDasypus_novemcinctus 6221
## muTransit_StatusInvaded:NameDidelphis_virginiana 6629
## muTransit_StatusInvaded:NameLynx_rufus 6018
## muTransit_StatusInvaded:NameMeleagris_gallopavo 6124
## muTransit_StatusInvaded:NameOdocoileus_virginianus 5431
## muTransit_StatusInvaded:NameProcyon_lotor 5899
## muTransit_StatusInvaded:NameSciurus_carolinensis 6328
## muTransit_StatusInvaded:NameSylvilagus_floridanus 6243
##
## 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.1 31.4
## p_waic 43.5 1.9
## waic 5078.1 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.2 31.4
## p_loo 43.7 2.0
## looic 5078.5 62.8
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 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 samplingsummary(model_loc)## 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.63 0.45 0.04 1.78 1.00 2652
## sd(muTransit_Intercept) 0.56 0.37 0.13 1.57 1.00 4248
## Tail_ESS
## sd(muForaging_Intercept) 2514
## sd(muTransit_Intercept) 4428
##
## ~Month (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 0.44 0.21 0.12 0.94 1.00 3118
## sd(muTransit_Intercept) 0.24 0.16 0.01 0.62 1.00 2386
## Tail_ESS
## sd(muForaging_Intercept) 4168
## sd(muTransit_Intercept) 2806
##
## ~Site (Number of levels: 5)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 0.69 0.40 0.11 1.70 1.00 2956
## sd(muTransit_Intercept) 0.61 0.33 0.21 1.44 1.00 4656
## Tail_ESS
## sd(muForaging_Intercept) 2418
## sd(muTransit_Intercept) 5868
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI
## muForaging_Intercept -0.06 0.72 -1.48
## muTransit_Intercept 0.92 0.62 -0.32
## muForaging_BehLocPatch -1.30 0.60 -2.47
## muForaging_NameDasypus_novemcinctus 1.61 0.55 0.55
## muForaging_NameDidelphis_virginiana 0.31 0.75 -1.18
## muForaging_NameLynx_rufus -0.60 0.88 -2.36
## muForaging_NameMeleagris_gallopavo 0.13 0.85 -1.56
## muForaging_NameOdocoileus_virginianus 0.29 0.44 -0.58
## muForaging_NameProcyon_lotor -1.21 0.57 -2.33
## muForaging_NameSciurus_carolinensis -0.13 0.66 -1.44
## muForaging_NameSylvilagus_floridanus 1.06 0.69 -0.29
## muForaging_BehLocPatch:NameDasypus_novemcinctus -0.32 0.81 -1.94
## muForaging_BehLocPatch:NameDidelphis_virginiana -0.16 0.96 -2.07
## muForaging_BehLocPatch:NameLynx_rufus -0.01 0.99 -1.95
## muForaging_BehLocPatch:NameMeleagris_gallopavo -0.07 0.98 -2.01
## muForaging_BehLocPatch:NameOdocoileus_virginianus -0.27 0.62 -1.48
## muForaging_BehLocPatch:NameProcyon_lotor -0.14 0.94 -1.99
## muForaging_BehLocPatch:NameSciurus_carolinensis -0.04 0.97 -1.95
## muForaging_BehLocPatch:NameSylvilagus_floridanus -0.13 0.97 -2.08
## muTransit_BehLocPatch 1.85 0.50 0.89
## muTransit_NameDasypus_novemcinctus 0.49 0.51 -0.50
## muTransit_NameDidelphis_virginiana 0.14 0.67 -1.15
## muTransit_NameLynx_rufus 1.47 0.74 0.05
## muTransit_NameMeleagris_gallopavo -0.34 0.79 -1.90
## muTransit_NameOdocoileus_virginianus -0.63 0.40 -1.44
## muTransit_NameProcyon_lotor -0.41 0.46 -1.30
## muTransit_NameSciurus_carolinensis -0.79 0.59 -1.94
## muTransit_NameSylvilagus_floridanus -0.13 0.65 -1.41
## muTransit_BehLocPatch:NameDasypus_novemcinctus -0.38 0.71 -1.72
## muTransit_BehLocPatch:NameDidelphis_virginiana 0.51 0.88 -1.18
## muTransit_BehLocPatch:NameLynx_rufus 0.04 0.98 -1.82
## muTransit_BehLocPatch:NameMeleagris_gallopavo 0.18 0.95 -1.64
## muTransit_BehLocPatch:NameOdocoileus_virginianus 0.08 0.51 -0.93
## muTransit_BehLocPatch:NameProcyon_lotor 0.39 0.74 -1.03
## muTransit_BehLocPatch:NameSciurus_carolinensis 0.13 0.94 -1.72
## muTransit_BehLocPatch:NameSylvilagus_floridanus 0.24 0.91 -1.51
## u-95% CI Rhat Bulk_ESS
## muForaging_Intercept 1.37 1.00 5899
## muTransit_Intercept 2.10 1.00 5143
## muForaging_BehLocPatch -0.16 1.00 8497
## muForaging_NameDasypus_novemcinctus 2.73 1.00 8576
## muForaging_NameDidelphis_virginiana 1.75 1.00 10630
## muForaging_NameLynx_rufus 1.09 1.00 11963
## muForaging_NameMeleagris_gallopavo 1.78 1.00 13302
## muForaging_NameOdocoileus_virginianus 1.16 1.00 7465
## muForaging_NameProcyon_lotor -0.11 1.00 9192
## muForaging_NameSciurus_carolinensis 1.14 1.00 9507
## muForaging_NameSylvilagus_floridanus 2.43 1.00 10905
## muForaging_BehLocPatch:NameDasypus_novemcinctus 1.26 1.00 12697
## muForaging_BehLocPatch:NameDidelphis_virginiana 1.70 1.00 14081
## muForaging_BehLocPatch:NameLynx_rufus 1.94 1.00 15700
## muForaging_BehLocPatch:NameMeleagris_gallopavo 1.84 1.00 15097
## muForaging_BehLocPatch:NameOdocoileus_virginianus 0.94 1.00 9143
## muForaging_BehLocPatch:NameProcyon_lotor 1.65 1.00 13385
## muForaging_BehLocPatch:NameSciurus_carolinensis 1.83 1.00 15079
## muForaging_BehLocPatch:NameSylvilagus_floridanus 1.76 1.00 14281
## muTransit_BehLocPatch 2.83 1.00 7287
## muTransit_NameDasypus_novemcinctus 1.51 1.00 8026
## muTransit_NameDidelphis_virginiana 1.51 1.00 10259
## muTransit_NameLynx_rufus 2.95 1.00 11589
## muTransit_NameMeleagris_gallopavo 1.22 1.00 12656
## muTransit_NameOdocoileus_virginianus 0.16 1.00 6160
## muTransit_NameProcyon_lotor 0.47 1.00 6586
## muTransit_NameSciurus_carolinensis 0.36 1.00 8991
## muTransit_NameSylvilagus_floridanus 1.19 1.00 11070
## muTransit_BehLocPatch:NameDasypus_novemcinctus 1.03 1.00 10107
## muTransit_BehLocPatch:NameDidelphis_virginiana 2.27 1.00 13336
## muTransit_BehLocPatch:NameLynx_rufus 1.94 1.00 15199
## muTransit_BehLocPatch:NameMeleagris_gallopavo 2.07 1.00 14630
## muTransit_BehLocPatch:NameOdocoileus_virginianus 1.08 1.00 7492
## muTransit_BehLocPatch:NameProcyon_lotor 1.89 1.00 10421
## muTransit_BehLocPatch:NameSciurus_carolinensis 2.00 1.00 14115
## muTransit_BehLocPatch:NameSylvilagus_floridanus 2.07 1.00 12622
## Tail_ESS
## muForaging_Intercept 5231
## muTransit_Intercept 4687
## muForaging_BehLocPatch 5933
## muForaging_NameDasypus_novemcinctus 6390
## muForaging_NameDidelphis_virginiana 6542
## muForaging_NameLynx_rufus 5486
## muForaging_NameMeleagris_gallopavo 6014
## muForaging_NameOdocoileus_virginianus 6439
## muForaging_NameProcyon_lotor 6626
## muForaging_NameSciurus_carolinensis 6065
## muForaging_NameSylvilagus_floridanus 6560
## muForaging_BehLocPatch:NameDasypus_novemcinctus 6507
## muForaging_BehLocPatch:NameDidelphis_virginiana 6129
## muForaging_BehLocPatch:NameLynx_rufus 6005
## muForaging_BehLocPatch:NameMeleagris_gallopavo 5772
## muForaging_BehLocPatch:NameOdocoileus_virginianus 6356
## muForaging_BehLocPatch:NameProcyon_lotor 5915
## muForaging_BehLocPatch:NameSciurus_carolinensis 5755
## muForaging_BehLocPatch:NameSylvilagus_floridanus 5346
## muTransit_BehLocPatch 6335
## muTransit_NameDasypus_novemcinctus 6122
## muTransit_NameDidelphis_virginiana 6268
## muTransit_NameLynx_rufus 6539
## muTransit_NameMeleagris_gallopavo 5530
## muTransit_NameOdocoileus_virginianus 5406
## muTransit_NameProcyon_lotor 5145
## muTransit_NameSciurus_carolinensis 5778
## muTransit_NameSylvilagus_floridanus 6642
## muTransit_BehLocPatch:NameDasypus_novemcinctus 6490
## muTransit_BehLocPatch:NameDidelphis_virginiana 6222
## muTransit_BehLocPatch:NameLynx_rufus 5983
## muTransit_BehLocPatch:NameMeleagris_gallopavo 5611
## muTransit_BehLocPatch:NameOdocoileus_virginianus 5729
## muTransit_BehLocPatch:NameProcyon_lotor 6493
## muTransit_BehLocPatch:NameSciurus_carolinensis 5942
## muTransit_BehLocPatch:NameSylvilagus_floridanus 5884
##
## 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:
## 5 (0.4%) 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.7 24.7
## p_waic 36.9 1.9
## waic 2017.4 49.4
##
## 5 (0.4%) 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 -1008.8 24.7
## p_loo 37.0 1.9
## looic 2017.6 49.4
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.7]).
##
## 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) +
# 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.3,
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.15, 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.3, "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) +
# 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.3,
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.3, "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.05840118 0.538125 FALSE
## b_muTransit_Intercept 0.91968444 0.932500 FALSE
## b_muForaging_BehLocPatch -1.29734100 0.984625 TRUE
## b_muForaging_NameDasypus_novemcinctus 1.61040189 0.998625 TRUE
## b_muForaging_NameDidelphis_virginiana 0.31277574 0.664625 FALSE
## b_muForaging_NameLynx_rufus -0.59856004 0.755125 FALSEpd_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.5010418 0.845000 FALSE
## b_muTransit_Intercept 1.5109298 0.998125 TRUE
## b_muForaging_StatusInvaded -0.2106001 0.682250 FALSE
## b_muForaging_NameDasypus_novemcinctus 1.8363116 0.999625 TRUE
## b_muForaging_NameDidelphis_virginiana 0.9031730 0.893375 FALSE
## b_muForaging_NameLynx_rufus -0.4446561 0.688625 FALSEresults_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.058401179 0.538125 FALSE
## b_muTransit_Intercept 0.919684441 0.932500 FALSE
## b_muForaging_BehLocPatch -1.297341004 0.984625 TRUE
## b_muForaging_NameDasypus_novemcinctus 1.610401888 0.998625 TRUE
## b_muForaging_NameDidelphis_virginiana 0.312775738 0.664625 FALSE
## b_muForaging_NameLynx_rufus -0.598560039 0.755125 FALSE
## b_muForaging_NameMeleagris_gallopavo 0.125684888 0.564250 FALSE
## b_muForaging_NameOdocoileus_virginianus 0.291804710 0.752125 FALSE
## b_muForaging_NameProcyon_lotor -1.207954766 0.986250 TRUE
## b_muForaging_NameSciurus_carolinensis -0.131622216 0.580000 FALSE
## b_muForaging_NameSylvilagus_floridanus 1.062730967 0.937625 FALSE
## b_muForaging_BehLocPatch:NameDasypus_novemcinctus -0.315627348 0.649125 FALSE
## b_muForaging_BehLocPatch:NameDidelphis_virginiana -0.156714709 0.559000 FALSE
## b_muForaging_BehLocPatch:NameLynx_rufus -0.009619574 0.502750 FALSE
## b_muForaging_BehLocPatch:NameMeleagris_gallopavo -0.068273025 0.522375 FALSE
## b_muForaging_BehLocPatch:NameOdocoileus_virginianus -0.265221085 0.657875 FALSE
## b_muForaging_BehLocPatch:NameProcyon_lotor -0.144141232 0.555875 FALSE
## b_muForaging_BehLocPatch:NameSciurus_carolinensis -0.042576445 0.510500 FALSE
## b_muForaging_BehLocPatch:NameSylvilagus_floridanus -0.126320321 0.549125 FALSE
## b_muTransit_BehLocPatch 1.845941976 1.000000 TRUE
## b_muTransit_NameDasypus_novemcinctus 0.486853304 0.832750 FALSE
## b_muTransit_NameDidelphis_virginiana 0.144899299 0.582625 FALSE
## b_muTransit_NameLynx_rufus 1.467439145 0.979625 TRUE
## b_muTransit_NameMeleagris_gallopavo -0.339232006 0.667125 FALSE
## b_muTransit_NameOdocoileus_virginianus -0.626341411 0.940375 FALSE
## b_muTransit_NameProcyon_lotor -0.413787003 0.818125 FALSE
## b_muTransit_NameSciurus_carolinensis -0.793357088 0.911625 FALSE
## b_muTransit_NameSylvilagus_floridanus -0.128928881 0.578750 FALSE
## b_muTransit_BehLocPatch:NameDasypus_novemcinctus -0.377947655 0.709000 FALSE
## b_muTransit_BehLocPatch:NameDidelphis_virginiana 0.505121014 0.713625 FALSE
## b_muTransit_BehLocPatch:NameLynx_rufus 0.040860257 0.518000 FALSE
## b_muTransit_BehLocPatch:NameMeleagris_gallopavo 0.179041659 0.574000 FALSE
## b_muTransit_BehLocPatch:NameOdocoileus_virginianus 0.076774962 0.559250 FALSE
## b_muTransit_BehLocPatch:NameProcyon_lotor 0.387262671 0.699625 FALSE
## b_muTransit_BehLocPatch:NameSciurus_carolinensis 0.132139381 0.551375 FALSE
## b_muTransit_BehLocPatch:NameSylvilagus_floridanus 0.243169241 0.601375 FALSEprint(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.50104175 0.845000
## b_muTransit_Intercept 1.51092983 0.998125
## b_muForaging_StatusInvaded -0.21060007 0.682250
## b_muForaging_NameDasypus_novemcinctus 1.83631165 0.999625
## b_muForaging_NameDidelphis_virginiana 0.90317304 0.893375
## b_muForaging_NameLynx_rufus -0.44465613 0.688625
## b_muForaging_NameMeleagris_gallopavo 0.21924714 0.627875
## b_muForaging_NameOdocoileus_virginianus 1.11654790 0.999500
## b_muForaging_NameProcyon_lotor -0.71485764 0.909500
## b_muForaging_NameSciurus_carolinensis 0.07684438 0.550875
## b_muForaging_NameSylvilagus_floridanus 1.53702054 0.995375
## b_muForaging_StatusInvaded:NameDasypus_novemcinctus 0.53092208 0.806750
## b_muForaging_StatusInvaded:NameDidelphis_virginiana 0.05997647 0.535500
## b_muForaging_StatusInvaded:NameLynx_rufus -0.19331483 0.576750
## b_muForaging_StatusInvaded:NameMeleagris_gallopavo 0.12636956 0.559625
## b_muForaging_StatusInvaded:NameOdocoileus_virginianus -0.13987203 0.619750
## b_muForaging_StatusInvaded:NameProcyon_lotor -0.16679537 0.601500
## b_muForaging_StatusInvaded:NameSciurus_carolinensis 0.16273762 0.585375
## b_muForaging_StatusInvaded:NameSylvilagus_floridanus 0.42539745 0.713000
## b_muTransit_StatusInvaded 0.18555792 0.686125
## b_muTransit_NameDasypus_novemcinctus 0.69339642 0.924000
## b_muTransit_NameDidelphis_virginiana 0.21785632 0.626625
## b_muTransit_NameLynx_rufus 1.45450841 0.980625
## b_muTransit_NameMeleagris_gallopavo 0.02887784 0.510625
## b_muTransit_NameOdocoileus_virginianus -0.17225504 0.722625
## b_muTransit_NameProcyon_lotor -0.50352009 0.891000
## b_muTransit_NameSciurus_carolinensis -0.86596152 0.922625
## b_muTransit_NameSylvilagus_floridanus -0.10944797 0.575500
## b_muTransit_StatusInvaded:NameDasypus_novemcinctus 0.02653855 0.513750
## b_muTransit_StatusInvaded:NameDidelphis_virginiana 0.35428071 0.678125
## b_muTransit_StatusInvaded:NameLynx_rufus 0.62657790 0.768000
## b_muTransit_StatusInvaded:NameMeleagris_gallopavo -0.22200445 0.614000
## b_muTransit_StatusInvaded:NameOdocoileus_virginianus -0.24626439 0.733000
## b_muTransit_StatusInvaded:NameProcyon_lotor -0.06069070 0.545000
## b_muTransit_StatusInvaded:NameSciurus_carolinensis -0.48782989 0.759750
## b_muTransit_StatusInvaded:NameSylvilagus_floridanus 0.13218805 0.568625
## 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 FALSEPublication Tables
make_pub_table_with_ci <- function(model, pd_table, file_name, caption_text) {
sum_mod <- summary(model)$fixed # regression coefficients summary
ci_df <- as.data.frame(sum_mod) %>%
tibble::rownames_to_column(var = "parameter") %>%
dplyr::select(parameter, Estimate, `l-95% CI`, `u-95% CI`) %>%
rename(
mean = Estimate,
l_95 = `l-95% CI`,
u_95 = `u-95% CI`
)
pd_table <- pd_table %>%
dplyr::mutate(pd = as.numeric(pd))
# Merge with pd_table
df_pub <- pd_table %>%
dplyr::select(parameter, mean_pd = mean, pd) %>%
left_join(ci_df, by = "parameter") %>%
dplyr::mutate(
mean = round(mean, 3),
pd = round(pd, 3),
l_95 = round(l_95, 3),
u_95 = round(u_95, 3),
bold_flag = pd >= 0.95
) %>%
dplyr::select(parameter, mean, pd, l_95, u_95, bold_flag)
bold_rows <- which(df_pub$bold_flag == TRUE)
ft <- flextable(df_pub %>% dplyr::select(-bold_flag)) %>%
set_caption(caption_text) %>%
bold(i = bold_rows, j = 1:5, bold = TRUE) %>%
autofit()
save_as_docx(ft, path = file_name)
return(ft)
}
# Status model
ft_status <- make_pub_table_with_ci(model_status, results_status,
file_name = "Status_Model_Parameters_CI.docx",
caption_text = "Status Model Parameter Estimates with 95% Credible Intervals")## Warning: There were 9 divergent transitions after warmup. Increasing
## adapt_delta above 0.95 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup# Location model
ft_loc <- make_pub_table_with_ci(model_loc, results_loc,
file_name = "Location_Model_Parameters_CI.docx",
caption_text = "Location Model Parameter Estimates with 95% Credible Intervals")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("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("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("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")