C4_O2_Behavior
Alan Ivory
2025-12-18
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)
# Install dependencies
#install.packages(c("posterior", "RcppParallel", "jsonlite"))
# Then install cmdstanr from the Stan developers’ repository
#install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
# Tree PCQ Data
tree_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
sheet = "Tree_PCQ")
# Soil Data
fuel_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
sheet = "Fuel_Sampling")
# Veg Data
Veg_Cover <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
sheet = "Veg_Cover")
# Shrub Cover Data
shrub_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
sheet = "Shrub_Cover")
# Site Data
CameraData <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/04_Wildlife/02_Data/CameraData.xlsx")
CameraLoc <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/04_Wildlife/02_Data/CameraLoc.xlsx",
sheet = "CameraLocations")
# Add effort data
effort_matrix <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/04_Wildlife/02_Data/CameraLoc.xlsx",
sheet = "Effort_Matrix_Full") %>%
pivot_longer(cols = matches("^202[4-5]-"), names_to = "week", values_to = "days") %>%
filter(days == "7") %>%
dplyr::select(Plot, week)Number of quadrats sampled per plot
I moved this from a later section because the filtering process removed quadrats that did not capture any species. Rows labeled as “None” were removed, suggesting that the number of quadrats sampled per plot is not consistent across all plots.
# Count the total number of quadrats per plot
quadrat_count <- Veg_Cover %>%
group_by(Plot) %>%
summarize(total_quadrats = n_distinct(Quadrat), .groups = "drop")Filter All data to only include specified species (Per PLANTS database)
#Filter tree data to only include trees with "tree" in the growth column
tree_data <- dplyr::filter(tree_data, Growth == "Tree")
#Filter Veg Cover to exclude Shrubs and Trees
Veg_Cover <- dplyr::filter(Veg_Cover, Growth != "Shrub" & Growth != "Tree")
#Filter Shrub Cover to only include Shrubs and Trees
shrub_data <- dplyr::filter(shrub_data, Growth == "Shrub" | Growth == "Tree")Filter all data to only include species found at 3% of all sites
# Calculate the total number of sites
total_sites <- nrow(CameraLoc)
# Function to filter data by frequency
filter_by_frequency <- function(df) {
# Group data by species and calculate the frequency
freq <- df %>%
group_by(Species) %>%
summarise(Frequency = n_distinct(Plot) / nrow(CameraLoc) * 100) %>%
filter(Frequency >= 3)
# Filter the original data to include only species with frequency >= 3%
filtered_df <- df %>%
filter(Species %in% freq$Species)
return(filtered_df)
}
# Filter tree data by frequency
tree_data <- filter_by_frequency(tree_data)
# Filter Veg Cover data by frequency
Veg_Cover <- filter_by_frequency(Veg_Cover)
# Filter Shrub Cover data by frequency
shrub_data <- filter_by_frequency(shrub_data)Shrub Cover Conversion
# Total length of Shrub cover at a site
shrub_cover <- shrub_data %>%
mutate(Cover = Line_End - Line_Start) %>%
group_by(Species_Name, Plot) %>%
summarise(Shrub_Total_Cover = sum(Cover, na.rm = TRUE), .groups = "drop") %>%
mutate(Shrub_Percent_Cover = Shrub_Total_Cover / 3000 * 100)
# Summed length of shrub over at a site
shrub_cover_summed <- shrub_cover %>%
group_by(Plot) %>%
summarize(total_shrub_cover = sum(Shrub_Total_Cover, na.rm = TRUE), .groups = "drop")Herbacous Cover Conversion
# Combine Plot and Quadrat columns
Veg_Cover <- Veg_Cover %>%
mutate(Plot_Quadrat = paste(Plot, Quadrat, sep = '_'))
# Join with CogonSites to get site information
Veg_Cover <- Veg_Cover %>%
left_join(CameraLoc, by = "Plot")
# Sum species cover across quadrats for each species at each plot
veg_cover_summed <- Veg_Cover %>%
group_by(Plot, Species_Name) %>%
summarize(total_cover = sum(Cover_Per, na.rm = TRUE), .groups = "drop")
# Calculate average herbaceous species cover
avg_species_cover <- veg_cover_summed %>%
left_join(quadrat_count, by = "Plot") %>%
mutate(avg_cover = total_cover / total_quadrats)Merging Herb cover with Shrub
This species matrix includes herbaceous and shrub species
# Merge shrub cover with herbaceous average cover
combined_cover <- avg_species_cover %>%
full_join(
shrub_cover %>%
dplyr::select(Plot, Species_Name, Shrub_Percent_Cover),
by = c("Plot", "Species_Name")
) %>%
mutate(
overlap_flag = ifelse(!is.na(avg_cover) & !is.na(Shrub_Percent_Cover), TRUE, FALSE), # Flag overlaps
final_cover = case_when(
!is.na(avg_cover) & is.na(Shrub_Percent_Cover) ~ avg_cover, # Use herbaceous cover if no shrub data
is.na(avg_cover) & !is.na(Shrub_Percent_Cover) ~ Shrub_Percent_Cover, # Use shrub cover if no herbaceous data
TRUE ~ NA_real_ # Leave as NA where overlaps exist
)
)
# Species Matrix
species_matrix <- combined_cover %>%
dplyr::select(Plot, Species_Name, final_cover) %>%
pivot_wider(
names_from = Species_Name,
values_from = final_cover,
values_fill = 0
)Summarize Cogongrass Cover
avg_cogongrass_cover <- species_matrix %>%
group_by(Plot) %>%
summarize(Avg_Cogongrass_Cover = sum(Imperata_cylindrica, na.rm = TRUE) / n(), .groups = "drop")Herbacous Shannon Diversity Index
# Summarize species cover by site
site_species_cover <- Veg_Cover %>%
group_by(Plot, Species_Name) %>%
summarize(total_cover = sum(Cover_Per, na.rm = TRUE)) %>%
ungroup()## `summarise()` has grouped output by 'Plot'. You can override using the
## `.groups` argument.# Calculate Shannon diversity per site
Veg_shannon_diversity <- site_species_cover %>%
group_by(Plot) %>%
mutate(proportion = total_cover / sum(total_cover)) %>%
summarize(Veg_shannon_index = -sum(proportion * log(proportion), na.rm = TRUE))
print(Veg_shannon_diversity)## # A tibble: 206 × 2
## Plot Veg_shannon_index
## <chr> <dbl>
## 1 BI200 2.28
## 2 BI201 2.20
## 3 BI202 1.50
## 4 BI97 1.82
## 5 BI99 3.06
## 6 BN210 2.97
## 7 BN211 2.43
## 8 BN212 2.22
## 9 BN96 3.05
## 10 BN98 2.79
## # ℹ 196 more 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 (optional)
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(student_t(3, 0, 2.5), class = "Intercept", dpar = "muForaging"),
prior(student_t(3, 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 | Plot) + (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...## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install the appropriate version of Rtools for 4.5.1 from
## https://cran.r-project.org/bin/windows/Rtools/.## Trying to compile a simple C file## Running "C:/PROGRA~1/R/R-45~1.1/bin/x64/Rcmd.exe" SHLIB foo.c
## using C compiler: 'gcc.exe (GCC) 14.2.0'
## gcc -I"C:/PROGRA~1/R/R-45~1.1/include" -DNDEBUG -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/Rcpp/include/" -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/RcppEigen/include/" -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/RcppEigen/include/unsupported" -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/BH/include" -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/StanHeaders/include/src/" -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/StanHeaders/include/" -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/RcppParallel/include/" -DRCPP_PARALLEL_USE_TBB=1 -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include "C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp" -std=c++1y -I"c:/rtools45/x86_64-w64-mingw32.static.posix/include" -O2 -Wall -std=gnu2x -mfpmath=sse -msse2 -mstackrealign -c foo.c -o foo.o
## cc1.exe: warning: command-line option '-std=c++14' is valid for C++/ObjC++ but not for C
## In file included from C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Core:19,
## from C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Dense:1,
## from C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [C:/PROGRA~1/R/R-45~1.1/etc/x64/Makeconf:289: foo.o] Error 1## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install the appropriate version of Rtools for 4.5.1 from
## https://cran.r-project.org/bin/windows/Rtools/.## Start sampling## Warning: There were 24 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 fit diagnostics / model comparison
summary(model_status)## Warning: There were 24 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 | Plot) + (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.30 0.31 0.01 1.13 1.00 3171
## sd(muTransit_Intercept) 0.17 0.18 0.00 0.66 1.00 3669
## Tail_ESS
## sd(muForaging_Intercept) 4475
## sd(muTransit_Intercept) 3744
##
## ~Month (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 0.69 0.24 0.36 1.26 1.00 3595
## sd(muTransit_Intercept) 0.16 0.11 0.01 0.41 1.00 3025
## Tail_ESS
## sd(muForaging_Intercept) 3969
## sd(muTransit_Intercept) 3523
##
## ~Plot (Number of levels: 32)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 0.59 0.13 0.37 0.89 1.00 3486
## sd(muTransit_Intercept) 0.43 0.10 0.26 0.65 1.00 3379
## Tail_ESS
## sd(muForaging_Intercept) 5327
## sd(muTransit_Intercept) 4649
##
## ~Site (Number of levels: 5)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 0.22 0.21 0.01 0.77 1.00 4037
## sd(muTransit_Intercept) 0.47 0.29 0.08 1.20 1.00 2974
## Tail_ESS
## sd(muForaging_Intercept) 5511
## sd(muTransit_Intercept) 2338
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI
## muForaging_Intercept -0.53 0.53 -1.61
## muTransit_Intercept 1.50 0.42 0.65
## muForaging_StatusInvaded -0.23 0.48 -1.17
## muForaging_NameDasypus_novemcinctus 1.75 0.53 0.72
## muForaging_NameDidelphis_virginiana 0.72 0.73 -0.72
## muForaging_NameLynx_rufus -0.44 0.89 -2.21
## muForaging_NameMeleagris_gallopavo 0.25 0.69 -1.12
## muForaging_NameOdocoileus_virginianus 1.21 0.38 0.46
## muForaging_NameProcyon_lotor -0.85 0.56 -1.98
## muForaging_NameSciurus_carolinensis -0.06 0.68 -1.39
## muForaging_NameSylvilagus_floridanus 1.68 0.60 0.51
## muForaging_StatusInvaded:NameDasypus_novemcinctus 0.67 0.63 -0.56
## muForaging_StatusInvaded:NameDidelphis_virginiana 0.08 0.81 -1.48
## muForaging_StatusInvaded:NameLynx_rufus -0.19 0.96 -2.08
## muForaging_StatusInvaded:NameMeleagris_gallopavo 0.11 0.86 -1.61
## muForaging_StatusInvaded:NameOdocoileus_virginianus -0.30 0.46 -1.19
## muForaging_StatusInvaded:NameProcyon_lotor -0.14 0.67 -1.44
## muForaging_StatusInvaded:NameSciurus_carolinensis 0.27 0.74 -1.18
## muForaging_StatusInvaded:NameSylvilagus_floridanus 0.40 0.75 -1.03
## muTransit_StatusInvaded 0.27 0.41 -0.54
## muTransit_NameDasypus_novemcinctus 0.79 0.48 -0.12
## muTransit_NameDidelphis_virginiana 0.30 0.68 -1.03
## muTransit_NameLynx_rufus 1.48 0.74 0.09
## muTransit_NameMeleagris_gallopavo 0.15 0.55 -0.91
## muTransit_NameOdocoileus_virginianus -0.14 0.30 -0.73
## muTransit_NameProcyon_lotor -0.55 0.41 -1.35
## muTransit_NameSciurus_carolinensis -0.78 0.61 -1.97
## muTransit_NameSylvilagus_floridanus -0.12 0.58 -1.22
## muTransit_StatusInvaded:NameDasypus_novemcinctus 0.00 0.58 -1.14
## muTransit_StatusInvaded:NameDidelphis_virginiana 0.37 0.74 -1.09
## muTransit_StatusInvaded:NameLynx_rufus 0.59 0.86 -1.06
## muTransit_StatusInvaded:NameMeleagris_gallopavo -0.19 0.78 -1.72
## muTransit_StatusInvaded:NameOdocoileus_virginianus -0.31 0.40 -1.12
## muTransit_StatusInvaded:NameProcyon_lotor 0.04 0.50 -0.95
## muTransit_StatusInvaded:NameSciurus_carolinensis -0.57 0.68 -1.92
## muTransit_StatusInvaded:NameSylvilagus_floridanus 0.10 0.72 -1.28
## u-95% CI Rhat Bulk_ESS
## muForaging_Intercept 0.51 1.00 4586
## muTransit_Intercept 2.33 1.00 5058
## muForaging_StatusInvaded 0.70 1.00 5082
## muForaging_NameDasypus_novemcinctus 2.81 1.00 7091
## muForaging_NameDidelphis_virginiana 2.16 1.00 9580
## muForaging_NameLynx_rufus 1.27 1.00 13301
## muForaging_NameMeleagris_gallopavo 1.61 1.00 10211
## muForaging_NameOdocoileus_virginianus 1.96 1.00 6673
## muForaging_NameProcyon_lotor 0.26 1.00 7924
## muForaging_NameSciurus_carolinensis 1.24 1.00 10366
## muForaging_NameSylvilagus_floridanus 2.87 1.00 8745
## muForaging_StatusInvaded:NameDasypus_novemcinctus 1.89 1.00 8079
## muForaging_StatusInvaded:NameDidelphis_virginiana 1.67 1.00 11079
## muForaging_StatusInvaded:NameLynx_rufus 1.62 1.00 13415
## muForaging_StatusInvaded:NameMeleagris_gallopavo 1.78 1.00 12787
## muForaging_StatusInvaded:NameOdocoileus_virginianus 0.60 1.00 5671
## muForaging_StatusInvaded:NameProcyon_lotor 1.18 1.00 7993
## muForaging_StatusInvaded:NameSciurus_carolinensis 1.71 1.00 10316
## muForaging_StatusInvaded:NameSylvilagus_floridanus 1.87 1.00 11021
## muTransit_StatusInvaded 1.10 1.00 4695
## muTransit_NameDasypus_novemcinctus 1.76 1.00 7266
## muTransit_NameDidelphis_virginiana 1.65 1.00 10030
## muTransit_NameLynx_rufus 2.97 1.00 10972
## muTransit_NameMeleagris_gallopavo 1.26 1.00 10595
## muTransit_NameOdocoileus_virginianus 0.43 1.00 6571
## muTransit_NameProcyon_lotor 0.28 1.00 6749
## muTransit_NameSciurus_carolinensis 0.43 1.00 9867
## muTransit_NameSylvilagus_floridanus 1.02 1.00 9827
## muTransit_StatusInvaded:NameDasypus_novemcinctus 1.17 1.00 7442
## muTransit_StatusInvaded:NameDidelphis_virginiana 1.78 1.00 10226
## muTransit_StatusInvaded:NameLynx_rufus 2.32 1.00 11077
## muTransit_StatusInvaded:NameMeleagris_gallopavo 1.34 1.00 11249
## muTransit_StatusInvaded:NameOdocoileus_virginianus 0.46 1.00 4740
## muTransit_StatusInvaded:NameProcyon_lotor 1.00 1.00 5711
## muTransit_StatusInvaded:NameSciurus_carolinensis 0.78 1.00 8822
## muTransit_StatusInvaded:NameSylvilagus_floridanus 1.54 1.00 10621
## Tail_ESS
## muForaging_Intercept 5713
## muTransit_Intercept 4715
## muForaging_StatusInvaded 4977
## muForaging_NameDasypus_novemcinctus 6332
## muForaging_NameDidelphis_virginiana 5888
## muForaging_NameLynx_rufus 6098
## muForaging_NameMeleagris_gallopavo 6726
## muForaging_NameOdocoileus_virginianus 6193
## muForaging_NameProcyon_lotor 5554
## muForaging_NameSciurus_carolinensis 6516
## muForaging_NameSylvilagus_floridanus 6388
## muForaging_StatusInvaded:NameDasypus_novemcinctus 6869
## muForaging_StatusInvaded:NameDidelphis_virginiana 5301
## muForaging_StatusInvaded:NameLynx_rufus 6365
## muForaging_StatusInvaded:NameMeleagris_gallopavo 5836
## muForaging_StatusInvaded:NameOdocoileus_virginianus 4989
## muForaging_StatusInvaded:NameProcyon_lotor 5863
## muForaging_StatusInvaded:NameSciurus_carolinensis 5385
## muForaging_StatusInvaded:NameSylvilagus_floridanus 6712
## muTransit_StatusInvaded 5190
## muTransit_NameDasypus_novemcinctus 5929
## muTransit_NameDidelphis_virginiana 6089
## muTransit_NameLynx_rufus 5920
## muTransit_NameMeleagris_gallopavo 6405
## muTransit_NameOdocoileus_virginianus 6079
## muTransit_NameProcyon_lotor 4682
## muTransit_NameSciurus_carolinensis 6483
## muTransit_NameSylvilagus_floridanus 6125
## muTransit_StatusInvaded:NameDasypus_novemcinctus 6074
## muTransit_StatusInvaded:NameDidelphis_virginiana 6144
## muTransit_StatusInvaded:NameLynx_rufus 6258
## muTransit_StatusInvaded:NameMeleagris_gallopavo 5305
## muTransit_StatusInvaded:NameOdocoileus_virginianus 5613
## muTransit_StatusInvaded:NameProcyon_lotor 5949
## muTransit_StatusInvaded:NameSciurus_carolinensis 5628
## muTransit_StatusInvaded:NameSylvilagus_floridanus 6692
##
## 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:
## 7 (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 -2488.1 32.0
## p_waic 75.0 2.2
## waic 4976.2 64.0
##
## 7 (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 -2488.3 32.0
## p_loo 75.2 2.2
## looic 4976.6 64.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.6]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.# Posterior predictions contrasting Invaded vs Non_Invaded, marginalizing random effects
# Need to work on this later
# Visualize the interaction
#ce_status_name <- conditional_effects(model_status,
# effects = "Status:Name",
# categorical = TRUE)
#plot(ce_status_name)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(student_t(3, 0, 2.5), class = "Intercept", dpar = "muForaging"),
prior(student_t(3, 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 | Plot) + (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...## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install the appropriate version of Rtools for 4.5.1 from
## https://cran.r-project.org/bin/windows/Rtools/.## Trying to compile a simple C file## Running "C:/PROGRA~1/R/R-45~1.1/bin/x64/Rcmd.exe" SHLIB foo.c
## using C compiler: 'gcc.exe (GCC) 14.2.0'
## gcc -I"C:/PROGRA~1/R/R-45~1.1/include" -DNDEBUG -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/Rcpp/include/" -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/RcppEigen/include/" -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/RcppEigen/include/unsupported" -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/BH/include" -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/StanHeaders/include/src/" -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/StanHeaders/include/" -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/RcppParallel/include/" -DRCPP_PARALLEL_USE_TBB=1 -I"C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include "C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp" -std=c++1y -I"c:/rtools45/x86_64-w64-mingw32.static.posix/include" -O2 -Wall -std=gnu2x -mfpmath=sse -msse2 -mstackrealign -c foo.c -o foo.o
## cc1.exe: warning: command-line option '-std=c++14' is valid for C++/ObjC++ but not for C
## In file included from C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Core:19,
## from C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Dense:1,
## from C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## C:/Users/DrewIvory/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [C:/PROGRA~1/R/R-45~1.1/etc/x64/Makeconf:289: foo.o] Error 1## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install the appropriate version of Rtools for 4.5.1 from
## https://cran.r-project.org/bin/windows/Rtools/.## Start sampling## Warning: There were 11 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 problemssummary(model_loc)## Warning: There were 11 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 | Plot) + (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.39 0.38 0.01 1.44 1.00 4227
## sd(muTransit_Intercept) 0.38 0.36 0.01 1.35 1.00 3697
## Tail_ESS
## sd(muForaging_Intercept) 4516
## sd(muTransit_Intercept) 4283
##
## ~Month (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 0.45 0.21 0.11 0.95 1.00 2553
## sd(muTransit_Intercept) 0.32 0.18 0.03 0.72 1.00 2655
## Tail_ESS
## sd(muForaging_Intercept) 2147
## sd(muTransit_Intercept) 2783
##
## ~Plot (Number of levels: 17)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 0.90 0.30 0.42 1.58 1.00 2677
## sd(muTransit_Intercept) 0.60 0.19 0.28 1.01 1.00 3014
## Tail_ESS
## sd(muForaging_Intercept) 3547
## sd(muTransit_Intercept) 4056
##
## ~Site (Number of levels: 5)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 0.45 0.40 0.01 1.48 1.00 2895
## sd(muTransit_Intercept) 0.43 0.33 0.02 1.27 1.00 3009
## Tail_ESS
## sd(muForaging_Intercept) 3516
## sd(muTransit_Intercept) 3907
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI
## muForaging_Intercept -0.16 0.66 -1.48
## muTransit_Intercept 1.20 0.58 0.01
## muForaging_BehLocPatch -1.31 0.61 -2.50
## muForaging_NameDasypus_novemcinctus 1.52 0.55 0.47
## muForaging_NameDidelphis_virginiana 0.22 0.76 -1.30
## muForaging_NameLynx_rufus -0.56 0.85 -2.26
## muForaging_NameMeleagris_gallopavo 0.14 0.86 -1.58
## muForaging_NameOdocoileus_virginianus 0.29 0.45 -0.61
## muForaging_NameProcyon_lotor -1.25 0.58 -2.40
## muForaging_NameSciurus_carolinensis -0.11 0.67 -1.44
## muForaging_NameSylvilagus_floridanus 0.99 0.69 -0.39
## muForaging_BehLocPatch:NameDasypus_novemcinctus -0.32 0.83 -1.98
## muForaging_BehLocPatch:NameDidelphis_virginiana -0.18 0.96 -2.08
## muForaging_BehLocPatch:NameLynx_rufus -0.03 1.00 -1.97
## muForaging_BehLocPatch:NameMeleagris_gallopavo -0.06 0.97 -1.96
## muForaging_BehLocPatch:NameOdocoileus_virginianus -0.25 0.63 -1.50
## muForaging_BehLocPatch:NameProcyon_lotor -0.15 0.93 -2.00
## muForaging_BehLocPatch:NameSciurus_carolinensis -0.05 0.98 -1.93
## muForaging_BehLocPatch:NameSylvilagus_floridanus -0.12 0.98 -2.08
## muTransit_BehLocPatch 1.92 0.51 0.93
## muTransit_NameDasypus_novemcinctus 0.59 0.51 -0.41
## muTransit_NameDidelphis_virginiana 0.27 0.67 -1.01
## muTransit_NameLynx_rufus 1.38 0.74 -0.01
## muTransit_NameMeleagris_gallopavo -0.29 0.80 -1.82
## muTransit_NameOdocoileus_virginianus -0.77 0.40 -1.57
## muTransit_NameProcyon_lotor -0.36 0.45 -1.22
## muTransit_NameSciurus_carolinensis -0.78 0.59 -1.90
## muTransit_NameSylvilagus_floridanus -0.10 0.65 -1.37
## muTransit_BehLocPatch:NameDasypus_novemcinctus -0.33 0.72 -1.72
## muTransit_BehLocPatch:NameDidelphis_virginiana 0.52 0.88 -1.16
## muTransit_BehLocPatch:NameLynx_rufus 0.03 0.97 -1.83
## muTransit_BehLocPatch:NameMeleagris_gallopavo 0.21 0.95 -1.60
## muTransit_BehLocPatch:NameOdocoileus_virginianus 0.07 0.52 -0.98
## muTransit_BehLocPatch:NameProcyon_lotor 0.45 0.75 -0.99
## muTransit_BehLocPatch:NameSciurus_carolinensis 0.13 0.96 -1.70
## muTransit_BehLocPatch:NameSylvilagus_floridanus 0.27 0.92 -1.49
## u-95% CI Rhat Bulk_ESS
## muForaging_Intercept 1.14 1.00 4613
## muTransit_Intercept 2.31 1.00 4391
## muForaging_BehLocPatch -0.11 1.00 6621
## muForaging_NameDasypus_novemcinctus 2.63 1.00 6691
## muForaging_NameDidelphis_virginiana 1.67 1.00 9390
## muForaging_NameLynx_rufus 1.07 1.00 10812
## muForaging_NameMeleagris_gallopavo 1.80 1.00 11403
## muForaging_NameOdocoileus_virginianus 1.19 1.00 5687
## muForaging_NameProcyon_lotor -0.12 1.00 6814
## muForaging_NameSciurus_carolinensis 1.17 1.00 8937
## muForaging_NameSylvilagus_floridanus 2.34 1.00 9080
## muForaging_BehLocPatch:NameDasypus_novemcinctus 1.28 1.00 9262
## muForaging_BehLocPatch:NameDidelphis_virginiana 1.66 1.00 11814
## muForaging_BehLocPatch:NameLynx_rufus 1.93 1.00 10989
## muForaging_BehLocPatch:NameMeleagris_gallopavo 1.86 1.00 11146
## muForaging_BehLocPatch:NameOdocoileus_virginianus 0.98 1.00 6303
## muForaging_BehLocPatch:NameProcyon_lotor 1.72 1.00 12429
## muForaging_BehLocPatch:NameSciurus_carolinensis 1.85 1.00 12241
## muForaging_BehLocPatch:NameSylvilagus_floridanus 1.77 1.00 11654
## muTransit_BehLocPatch 2.93 1.00 4788
## muTransit_NameDasypus_novemcinctus 1.59 1.00 6555
## muTransit_NameDidelphis_virginiana 1.63 1.00 8727
## muTransit_NameLynx_rufus 2.82 1.00 9156
## muTransit_NameMeleagris_gallopavo 1.30 1.00 10909
## muTransit_NameOdocoileus_virginianus -0.00 1.00 5384
## muTransit_NameProcyon_lotor 0.53 1.00 5816
## muTransit_NameSciurus_carolinensis 0.41 1.00 8561
## muTransit_NameSylvilagus_floridanus 1.19 1.00 8437
## muTransit_BehLocPatch:NameDasypus_novemcinctus 1.11 1.00 6981
## muTransit_BehLocPatch:NameDidelphis_virginiana 2.28 1.00 9362
## muTransit_BehLocPatch:NameLynx_rufus 1.94 1.00 11554
## muTransit_BehLocPatch:NameMeleagris_gallopavo 2.10 1.00 11190
## muTransit_BehLocPatch:NameOdocoileus_virginianus 1.08 1.00 4832
## muTransit_BehLocPatch:NameProcyon_lotor 1.94 1.00 9620
## muTransit_BehLocPatch:NameSciurus_carolinensis 2.01 1.00 10704
## muTransit_BehLocPatch:NameSylvilagus_floridanus 2.10 1.00 11201
## Tail_ESS
## muForaging_Intercept 4927
## muTransit_Intercept 4106
## muForaging_BehLocPatch 5697
## muForaging_NameDasypus_novemcinctus 5879
## muForaging_NameDidelphis_virginiana 5915
## muForaging_NameLynx_rufus 6156
## muForaging_NameMeleagris_gallopavo 5387
## muForaging_NameOdocoileus_virginianus 5934
## muForaging_NameProcyon_lotor 6506
## muForaging_NameSciurus_carolinensis 6516
## muForaging_NameSylvilagus_floridanus 6171
## muForaging_BehLocPatch:NameDasypus_novemcinctus 5897
## muForaging_BehLocPatch:NameDidelphis_virginiana 6143
## muForaging_BehLocPatch:NameLynx_rufus 4873
## muForaging_BehLocPatch:NameMeleagris_gallopavo 5662
## muForaging_BehLocPatch:NameOdocoileus_virginianus 5805
## muForaging_BehLocPatch:NameProcyon_lotor 6052
## muForaging_BehLocPatch:NameSciurus_carolinensis 5995
## muForaging_BehLocPatch:NameSylvilagus_floridanus 6218
## muTransit_BehLocPatch 4715
## muTransit_NameDasypus_novemcinctus 6105
## muTransit_NameDidelphis_virginiana 6788
## muTransit_NameLynx_rufus 6271
## muTransit_NameMeleagris_gallopavo 6035
## muTransit_NameOdocoileus_virginianus 6197
## muTransit_NameProcyon_lotor 5993
## muTransit_NameSciurus_carolinensis 6007
## muTransit_NameSylvilagus_floridanus 6251
## muTransit_BehLocPatch:NameDasypus_novemcinctus 6014
## muTransit_BehLocPatch:NameDidelphis_virginiana 6126
## muTransit_BehLocPatch:NameLynx_rufus 5915
## muTransit_BehLocPatch:NameMeleagris_gallopavo 6048
## muTransit_BehLocPatch:NameOdocoileus_virginianus 4540
## muTransit_BehLocPatch:NameProcyon_lotor 5831
## muTransit_BehLocPatch:NameSciurus_carolinensis 6070
## muTransit_BehLocPatch:NameSylvilagus_floridanus 5741
##
## 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:
## 8 (0.6%) p_waic estimates greater than 0.4. We recommend trying loo instead.##
## Computed from 8000 by 1300 log-likelihood matrix.
##
## Estimate SE
## elpd_waic -981.9 25.3
## p_waic 49.2 2.2
## waic 1963.9 50.7
##
## 8 (0.6%) 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 -982.1 25.3
## p_loo 49.4 2.2
## looic 1964.3 50.7
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.7]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.# Posterior predictions contrasting Patch vs Outside, marginalizing random effects
# Need new code hereWithin Patch Behavior Plots
# Basic conditional effects for all predictors (must include categorical = TRUE)
conditional_effects(model_loc, categorical = TRUE)## Warning: Interactions cannot be plotted directly if 'categorical' is TRUE.
## Please use argument 'conditions' instead.
# For specific predictors of interest:
# 1) BehLoc effect (inside vs outside patches) - your main question
ce_behloc <- conditional_effects(model_loc, effects = "BehLoc", categorical = TRUE)
plot(ce_behloc)
# Customize a specific plot:
plot(ce_behloc, plot = FALSE)[[1]] +
labs(title = "Behavior by Patch Location (Invaded Sites Only)",
x = "Location",
y = "Predicted Probability") +
theme_minimal()
# Multiple effects at once:
ce_multiple <- conditional_effects(model_loc,
effects = c("BehLoc"),
categorical = TRUE)
plot(ce_multiple)
# Conditional effects showing interaction
#ce_int <- conditional_effects(model_loc,
# effects = "BehLoc:Name",
# categorical = TRUE)
#plot(ce_int)