CWD_Behavior_Multinomial
Alan Ivory
2026-10-06
Set-up
knitr::opts_chunk$set(echo = TRUE)
pacman::p_load(tidyverse, readxl, lubridate, vegan,
brms, tidybayes, posterior, loo,
ggplot2, ggdist, patchwork, scales, forcats,
kableExtra, flextable, officer)
set.seed(97)
options(mc.cores = parallel::detectCores())
proj <- "C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/13_CWD_Wildlife"
# Camera Review Data
CameraData <- read_excel(file.path(proj, "02_Data", "CWD_Camera_Review.xlsx"),
sheet = "Alan")## Warning: Expecting numeric in H1065 / R1065C8: got 'none'## Warning: Expecting numeric in I1065 / R1065C9: got 'none'## Warning: Coercing numeric to date in D1066 / R1066C4## Warning: Expecting numeric in H1066 / R1066C8: got 'none'## Warning: Expecting numeric in I1066 / R1066C9: got 'none'## Warning: Expecting numeric in H1128 / R1128C8: got 'none'## Warning: Expecting numeric in I1128 / R1128C9: got 'none'## Warning: Expecting numeric in H1129 / R1129C8: got 'none'## Warning: Expecting numeric in I1129 / R1129C9: got 'none'## Warning: Expecting numeric in H1130 / R1130C8: got 'none'## Warning: Expecting numeric in I1130 / R1130C9: got 'none'## Warning: Expecting numeric in H1181 / R1181C8: got 'none'## Warning: Expecting numeric in I1181 / R1181C9: got 'none'## Warning: Expecting numeric in H1182 / R1182C8: got 'none'## Warning: Expecting numeric in I1182 / R1182C9: got 'none'## Warning: Expecting numeric in H1183 / R1183C8: got 'none'## Warning: Expecting numeric in I1183 / R1183C9: got 'none'## Warning: Expecting numeric in H2960 / R2960C8: got 'none'## Warning: Expecting numeric in I2960 / R2960C9: got 'none'## Warning: Expecting numeric in H3053 / R3053C8: got 'none'## Warning: Expecting numeric in I3053 / R3053C9: got 'none'## Warning: Expecting numeric in H3054 / R3054C8: got 'none'## Warning: Expecting numeric in I3054 / R3054C9: got 'none'## Warning: Expecting numeric in H5553 / R5553C8: got 'none'## Warning: Expecting numeric in I5553 / R5553C9: got 'none'## Warning: Expecting numeric in J5553 / R5553C10: got 'none'# Vegetation quadrat data (ground cover, soil moisture, canopy, decay class)
Master <- read_excel(file.path(proj, "02_Data", "CWD_Vegetation_Data.xlsx"),
sheet = "Master")
# Point-centered-quarter tree data
Point <- read_excel(file.path(proj, "02_Data", "CWD_Vegetation_Data.xlsx"),
sheet = "Point")
# Effort matrix: 1 = camera deployed that plot-day
effort_matrix <- read_excel(file.path(proj, "02_Data", "Cam_Dates_Tracker.xlsx"),
sheet = "Sheet1")## New names:
## • `` -> `...1`out_dir <- file.path(proj, "03_Output", "behavior")
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)Analysis settings
min_obs_species <- 20 # species need this many independent observations
min_grp_obs <- 30 # orders below this are lumped into "Other"
gap_mins <- 30 # independence filter
decay_mode <- "factor"
stopifnot(decay_mode %in% c("factor", "monotonic"))
decay_ordered <- decay_mode == "monotonic"
cat("Decay enters the model as:", decay_mode, "\n")## Decay enters the model as: factorData Preparation
CameraData <- CameraData %>%
mutate(
Plot = as.character(Plot),
Location = suppressWarnings(as.numeric(as.character(Location))),
Common_Name = str_to_lower(str_trim(as.character(Common_Name))),
Class = str_trim(as.character(Class)),
Order = str_trim(as.character(Order)),
Camera.Type = str_trim(as.character(Camera.Type)),
Date = as.Date(Date),
DateTime = as.POSIXct(paste(Date, format(Time, "%H:%M:%S")),
tz = "America/New_York"),
Year = year(Date),
Area = str_extract(Plot, "^[A-Z]+"),
Behavior = recode(str_to_lower(str_trim(as.character(Behavior))),
"search" = "Local_Search",
"foraging" = "Foraging",
"transit" = "Transit",
"none" = NA_character_)
)
effort_long <- effort_matrix %>%
rename(Date = 1) %>%
mutate(Date = as.Date(Date)) %>%
pivot_longer(-Date, names_to = "Plot", values_to = "deployed") %>%
mutate(deployed = suppressWarnings(as.numeric(as.character(deployed)))) %>%
filter(deployed == 1) %>%
dplyr::select(Plot, Date)
CamFilt <- CameraData %>%
filter(Location == 1, Year == 2026,
!is.na(Common_Name), !is.na(Behavior), !is.na(DateTime)) %>%
semi_join(effort_long, by = c("Plot", "Date"))
cat("Location codes in the raw review sheet, before filtering:\n")## Location codes in the raw review sheet, before filtering:print(table(CameraData$Location, useNA = "ifany"))##
## 0 1 2 <NA>
## 719 5127 1 33step_n <- c(
"Raw records" = nrow(CameraData),
"On the CWD (Location == 1)" = sum(CameraData$Location == 1, na.rm = TRUE),
" + 2026 only" = sum(CameraData$Location == 1 &
CameraData$Year == 2026, na.rm = TRUE),
" + species/behavior/time scored" = nrow(CameraData %>%
filter(Location == 1, Year == 2026, !is.na(Common_Name),
!is.na(Behavior), !is.na(DateTime))),
" + camera deployed that day" = nrow(CamFilt)
)
cat("\nAttrition:\n")##
## Attrition:print(data.frame(records = step_n))## records
## Raw records 5880
## On the CWD (Location == 1) 5127
## + 2026 only 4966
## + species/behavior/time scored 4966
## + camera deployed that day 4411Independence filter
One observation per species, per plot, per 30 minutes.
CamInd <- CamFilt %>%
arrange(Plot, Common_Name, DateTime) %>%
group_by(Plot, Common_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()
cat("Independent observations with a scored behavior:", nrow(CamInd), "\n")## Independent observations with a scored behavior: 3530print(table(CamInd$Behavior))##
## Foraging Local_Search Transit
## 267 1742 1521Subset to species with enough observations
species_counts <- CamInd %>%
group_by(Common_Name) %>%
summarize(count = n(), .groups = "drop")
species_subset <- species_counts %>%
filter(count >= min_obs_species) %>%
pull(Common_Name)
cat("Species retained:", length(species_subset), "covering",
sum(species_counts$count[species_counts$Common_Name %in% species_subset]),
"of", nrow(CamInd), "observations\n")## Species retained: 18 covering 3349 of 3530 observationsprint(as.data.frame(species_counts %>%
filter(Common_Name %in% species_subset) %>%
arrange(desc(count))), row.names = FALSE)## Common_Name count
## cotton_mouse 1357
## eastern_woodrat 413
## northern_cardinal 272
## nine_banded_armadillo 212
## eastern_chipmunk 187
## virginia_opossum 130
## carolina_wren 113
## raccoon 111
## hispid_cotton_rat 107
## eastern_gray_squirrel 104
## five_lined_skink 60
## eastern_cottontail 54
## white_tailed_deer 44
## eastern_towhee 43
## bobcat 41
## broad_headed_skink 36
## indigo_bunting 34
## chucks_wills_widow 31Functional guilds
guild_lookup <- tibble::tribble(
~Common_Name, ~Guild,
"cotton_mouse", "Small rodent",
"eastern_woodrat", "Small rodent",
"hispid_cotton_rat", "Small rodent",
"eastern_chipmunk", "Sciurid",
"eastern_gray_squirrel", "Sciurid",
"southern_flying_squirrel", "Sciurid",
"nine_banded_armadillo", "Ground omnivore",
"virginia_opossum", "Ground omnivore",
"raccoon", "Ground omnivore",
"striped_skunk", "Ground omnivore",
"bobcat", "Carnivore",
"coyote", "Carnivore",
"white_tailed_deer", "Herbivore",
"eastern_cottontail", "Herbivore",
"northern_cardinal", "Songbird",
"carolina_wren", "Songbird",
"eastern_towhee", "Songbird",
"indigo_bunting", "Songbird",
"brown_thrasher", "Songbird",
"gray_catbird", "Songbird",
"northern_house_wren", "Songbird",
"yellow_breasted_chat", "Songbird",
"blue_gray_gnatcatcher", "Songbird",
"chucks_wills_widow", "Nightjar",
"five_lined_skink", "Reptile",
"broad_headed_skink", "Reptile",
"six_lined_racerunner", "Reptile",
"southern_black_racer", "Reptile",
"northern_bobwhite", "Ground bird",
"wild_turkey", "Ground bird"
)Site covariates
Same derivation as the occupancy analysis: Daubenmire classes converted to midpoint percentages before averaging, tree metrics from the point-centered-quarter data.
daub_mid <- function(x) {
c(`0` = 0, `1` = 2.5, `2` = 15, `3` = 37.5,
`4` = 62.5, `5` = 85, `6` = 97.5)[as.character(x)]
}
cover_cols <- c("bg", "pl", "dll", "lg", "fw", "forbs", "pc",
"shrubs", "dh", "sdw", "vines")
site_struct <- Master %>%
mutate(across(all_of(cover_cols), ~ as.numeric(daub_mid(.x)),
.names = "{.col}_pct"),
sm_mean = rowMeans(across(c(sm1, sm2, sm3)), na.rm = TRUE)) %>%
group_by(cwd) %>%
summarise(decay = first(decay),
off_ground = mean(og > 0, na.rm = TRUE),
across(ends_with("_pct"), ~ mean(.x, na.rm = TRUE)),
soil_moist = mean(sm_mean, na.rm = TRUE),
canopy_cover = mean(cc, na.rm = TRUE),
.groups = "drop")
site_struct$gc_shannon <- vegan::diversity(
as.matrix(site_struct %>% dplyr::select(ends_with("_pct"))), index = "shannon")
tree_metrics <- Point %>%
group_by(cwd) %>%
summarise(tree_dens = 10000 / (mean(dist, na.rm = TRUE)^2),
basal_area = (10000 / (mean(dist, na.rm = TRUE)^2)) *
mean(pi * (dbh / 200)^2, na.rm = TRUE),
.groups = "drop")
site_env <- site_struct %>%
left_join(tree_metrics, by = "cwd") %>%
rename(Plot = cwd)Prepare covariates
cov_vars <- c("canopy_cover", "soil_moist", "gc_shannon",
"lg_pct", "shrubs_pct", "basal_area", "tree_dens", "off_ground")
site_scaled <- site_env %>% dplyr::select(Plot, decay, all_of(cov_vars))
sc <- scale(as.matrix(site_scaled[, cov_vars]))
cov_centers <- attr(sc, "scaled:center")
cov_scales <- attr(sc, "scaled:scale")
site_scaled[, cov_vars] <- as.data.frame(sc)
# Back-transform helpers for the figures
unscale <- function(z, v) z * cov_scales[[v]] + cov_centers[[v]]
rescale_to <- function(raw, v) (raw - cov_centers[[v]]) / cov_scales[[v]]Behavior Occurrence
dat <- CamInd %>%
filter(Common_Name %in% species_subset) %>%
left_join(site_scaled, by = "Plot") %>%
left_join(guild_lookup, by = "Common_Name") %>%
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,
# month factor to absorb temporal clustering
Month = factor(month(DateTime)),
# decay class as an ORDERED factor, 1 < 2 < 3 < 4 < 5.
decay_num = as.numeric(decay),
# ordered only in monotonic mode
decay = factor(as.integer(decay), levels = 1:5, ordered = decay_ordered),
# set factor baselines
Behavior = factor(Behavior, ordered = FALSE),
Behavior = relevel(Behavior, ref = "Local_Search"),
Plot = factor(Plot),
Area = factor(Area),
Common_Name = factor(Common_Name),
Class = factor(Class),
Guild = factor(Guild),
Camera.Type = factor(Camera.Type)
)
# Lump sparse orders so the Order grouping is estimable
ord_n <- dat %>% count(Order)
dat <- dat %>%
mutate(Order_grp = factor(ifelse(Order %in% ord_n$Order[ord_n$n >= min_grp_obs],
Order, "Other")))
stopifnot(!any(is.na(dat$Guild)), !any(is.na(dat$decay)), !any(is.na(dat$Behavior)),
nlevels(dat$decay) == 5,
is.ordered(dat$decay) == decay_ordered)
cat("Model data:", nrow(dat), "observations,", nlevels(dat$Common_Name),
"species,", nlevels(dat$Plot), "logs\n")## Model data: 3349 observations, 18 species, 50 logscat("\nBehavior x Guild:\n"); print(table(dat$Guild, dat$Behavior))##
## Behavior x Guild:##
## Local_Search Foraging Transit
## Carnivore 11 1 29
## Ground omnivore 95 84 274
## Herbivore 42 14 42
## Nightjar 27 0 4
## Reptile 51 0 45
## Sciurid 167 31 93
## Small rodent 862 60 955
## Songbird 361 63 38cat("\nBehavior x Order:\n"); print(table(dat$Order_grp, dat$Behavior))##
## Behavior x Order:##
## Local_Search Foraging Transit
## Artiodactyla 10 8 26
## Caprimulgiformes 27 0 4
## Carnivora 53 6 93
## Cingulata 29 78 105
## Didelphimorphia 24 1 105
## Lagomorpha 32 6 16
## Passeriformes 361 63 38
## Rodentia 1029 91 1048
## Squamata 51 0 45cat("\nBehavior x Class:\n"); print(table(dat$Class, dat$Behavior))##
## Behavior x Class:##
## Local_Search Foraging Transit
## Aves 388 63 42
## Mammalia 1177 190 1393
## Reptilia 51 0 45cat("\nObservations per decay class:\n")##
## Observations per decay class:print(table(dat$decay))##
## 1 2 3 4 5
## 228 432 869 974 846Objective: Do wildlife behaviors shift with CWD decay class?
behavior_formula <- if (decay_mode == "monotonic") {
bf(Behavior ~ mo(decay) * Guild + off_ground + canopy_cover + shrubs_pct +
(1 | Common_Name) + (1 | Plot) + (1 | Area) +
(1 | Camera.Type) + (1 | Month))
} else {
bf(Behavior ~ decay * Guild + off_ground + canopy_cover + shrubs_pct +
(1 | Common_Name) + (1 | Plot) + (1 | Area) +
(1 | Camera.Type) + (1 | Month))
}
print(behavior_formula)## Behavior ~ decay * Guild + off_ground + canopy_cover + shrubs_pct + (1 | Common_Name) + (1 | Plot) + (1 | Area) + (1 | Camera.Type) + (1 | Month)behavior_family <- categorical(link = "logit", refcat = "Local_Search")
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")
)
simo_slots <- as.data.frame(
get_prior(behavior_formula, family = behavior_family, data = dat)) %>%
filter(class == "simo", coef != "")
if (nrow(simo_slots) > 0) {
simo_priors <- do.call(c, unname(Map(
function(cf, dp) set_prior("dirichlet(1)", class = "simo",
coef = cf, dpar = dp),
simo_slots$coef, simo_slots$dpar)))
priors <- c(priors, simo_priors)
cat("Simplex parameters:", nrow(simo_slots), "\n")
print(as.data.frame(simo_slots[, c("dpar", "coef")]), row.names = FALSE)
} else {
cat("No simplex parameters (decay_mode = '", decay_mode, "').\n", sep = "")
}## No simplex parameters (decay_mode = 'factor').invisible(validate_prior(priors, behavior_formula,
family = behavior_family, data = dat))
cat("\nAll priors matched to model parameters.\n")##
## All priors matched to model parameters.Estimability check before fitting
cat("Observations per decay class:\n")## Observations per decay class:print(table(dat$decay))##
## 1 2 3 4 5
## 228 432 869 974 846cat("\nGuild x decay (all behaviors):\n")##
## Guild x decay (all behaviors):print(table(dat$Guild, dat$decay))##
## 1 2 3 4 5
## Carnivore 3 9 21 2 6
## Ground omnivore 27 45 110 177 94
## Herbivore 11 15 33 15 24
## Nightjar 0 1 0 14 16
## Reptile 9 6 17 53 11
## Sciurid 16 27 66 80 102
## Small rodent 125 240 491 504 517
## Songbird 37 89 131 129 76cat("\nGuild x decay (foraging only):\n")##
## Guild x decay (foraging only):print(table(droplevels(dat$Guild[dat$Behavior == "Foraging"]),
droplevels(dat$decay[dat$Behavior == "Foraging"])))##
## 1 2 3 4 5
## Carnivore 0 0 0 1 0
## Ground omnivore 2 3 13 49 17
## Herbivore 5 2 5 0 2
## Sciurid 0 1 6 6 18
## Small rodent 3 7 14 13 23
## Songbird 2 9 11 27 14# --- Empty cells: drop the offending guilds -------------------------------
cell_counts <- table(dat$Guild, dat$decay)
empty_guilds <- rownames(cell_counts)[apply(cell_counts == 0, 1, any)]
if (length(empty_guilds) > 0) {
dropped_n <- sum(dat$Guild %in% empty_guilds)
cat("\nGuilds missing at least one decay class, dropped because a monotonic",
"\nshape cannot be estimated for a guild absent from part of the gradient:\n")
for (gname in empty_guilds) {
missing_at <- colnames(cell_counts)[cell_counts[gname, ] == 0]
cat(sprintf(" %-18s no observations at decay class %s (%d obs, %s of data)\n",
gname, paste(missing_at, collapse = ", "),
sum(dat$Guild == gname),
scales::percent(mean(dat$Guild == gname), accuracy = 0.1)))
}
dat <- dat %>%
filter(!Guild %in% empty_guilds) %>%
droplevels()
cat(sprintf("\nRetained %d of %d observations (%d dropped), %d guilds, %d species.\n",
nrow(dat), nrow(dat) + dropped_n, dropped_n,
nlevels(dat$Guild), nlevels(dat$Common_Name)))
} else {
cat("\nNo empty Guild x decay cells.\n")
}##
## Guilds missing at least one decay class, dropped because a monotonic
## shape cannot be estimated for a guild absent from part of the gradient:
## Nightjar no observations at decay class 1, 3 (31 obs, 0.9% of data)
##
## Retained 3318 of 3349 observations (31 dropped), 7 guilds, 17 species.span <- table(dat$Guild, dat$decay)
classes_present <- rowSums(span > 0)
cat("\nDecay classes with data, per retained guild:\n")##
## Decay classes with data, per retained guild:print(classes_present)## Carnivore Ground omnivore Herbivore Reptile Sciurid
## 5 5 5 5 5
## Small rodent Songbird
## 5 5stopifnot(all(classes_present >= 3))
X_fac <- model.matrix(~ decay * Guild + canopy_cover + shrubs_pct,
data = dat %>% mutate(decay = factor(decay, ordered = FALSE)))
cat(sprintf("\nFixed-effect columns: %d as an unordered factor vs %d with mo().\n",
ncol(X_fac),
2 + nlevels(dat$Guild) + 1))##
## Fixed-effect columns: 37 as an unordered factor vs 10 with mo().fac_vars <- names(dat)[sapply(dat, is.factor)]
empty_lvls <- lapply(dat[fac_vars], function(x)
setdiff(levels(x), unique(as.character(x))))
empty_lvls <- empty_lvls[lengths(empty_lvls) > 0]
if (length(empty_lvls)) {
cat("\nFactors carrying unused levels after the drop:\n")
for (v in names(empty_lvls))
cat(sprintf(" %-14s %s\n", v, paste(empty_lvls[[v]], collapse = ", ")))
stop("Unused factor levels reached the model frame; droplevels() failed.")
} else {
cat("\nNo unused factor levels remain.\n")
}##
## No unused factor levels remain.sparse_forage <- dat %>%
group_by(Guild) %>%
summarise(n_forage = sum(Behavior == "Foraging"), n = n(), .groups = "drop") %>%
filter(n_forage < 5)
if (nrow(sparse_forage) > 0) {
cat("\nGuilds with fewer than five foraging records. These are kept, but",
"their\nforaging estimates are prior-driven and should not be",
"interpreted:\n")
print(as.data.frame(sparse_forage), row.names = FALSE)
}##
## Guilds with fewer than five foraging records. These are kept, but their
## foraging estimates are prior-driven and should not be interpreted:
## Guild n_forage n
## Carnivore 1 41
## Reptile 0 96Run the model
model_global <- brm(
behavior_formula, # Behavior ~ mo(decay) * Guild + off_ground + canopy_cover +
# shrubs_pct +
# (1|Common_Name) + (1|Plot) + (1|Area) +
# (1|Camera.Type) + (1|Month)
family = behavior_family,
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 5 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_global)## Warning: There were 5 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 ~ decay * Guild + off_ground + canopy_cover + shrubs_pct + (1 | Common_Name) + (1 | Plot) + (1 | Area) + (1 | Camera.Type) + (1 | Month)
## Data: dat (Number of observations: 3318)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Area (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 0.73 0.51 0.14 2.12 1.00 3904
## sd(muTransit_Intercept) 0.47 0.42 0.03 1.55 1.00 3430
## Tail_ESS
## sd(muForaging_Intercept) 2787
## sd(muTransit_Intercept) 3293
##
## ~Camera.Type (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 0.68 0.65 0.02 2.34 1.00 3307
## sd(muTransit_Intercept) 0.68 0.53 0.07 2.09 1.00 4314
## Tail_ESS
## sd(muForaging_Intercept) 4295
## sd(muTransit_Intercept) 3951
##
## ~Common_Name (Number of levels: 17)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 1.31 0.35 0.77 2.13 1.00 4343
## sd(muTransit_Intercept) 0.53 0.17 0.27 0.93 1.00 3416
## Tail_ESS
## sd(muForaging_Intercept) 4636
## sd(muTransit_Intercept) 4662
##
## ~Month (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 0.19 0.21 0.01 0.76 1.00 3340
## sd(muTransit_Intercept) 0.11 0.13 0.00 0.44 1.00 3661
## Tail_ESS
## sd(muForaging_Intercept) 4783
## sd(muTransit_Intercept) 5061
##
## ~Plot (Number of levels: 50)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(muForaging_Intercept) 0.49 0.18 0.13 0.86 1.00 2173
## sd(muTransit_Intercept) 0.57 0.09 0.40 0.77 1.00 3467
## Tail_ESS
## sd(muForaging_Intercept) 2067
## sd(muTransit_Intercept) 4888
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## muForaging_Intercept -1.94 0.99 -3.76 0.19 1.00
## muTransit_Intercept 0.09 0.75 -1.45 1.51 1.00
## muForaging_decay2 -0.43 0.57 -1.54 0.69 1.00
## muForaging_decay3 -0.35 0.52 -1.37 0.68 1.00
## muForaging_decay4 0.08 0.53 -0.96 1.12 1.00
## muForaging_decay5 -0.03 0.53 -1.06 1.02 1.00
## muForaging_GuildGroundomnivore 0.41 0.75 -1.07 1.85 1.00
## muForaging_GuildHerbivore 0.92 0.77 -0.62 2.40 1.00
## muForaging_GuildReptile -0.94 0.90 -2.71 0.85 1.00
## muForaging_GuildSciurid 0.02 0.78 -1.47 1.57 1.00
## muForaging_GuildSmallrodent -0.10 0.73 -1.51 1.35 1.00
## muForaging_GuildSongbird -0.12 0.70 -1.49 1.22 1.00
## muForaging_off_ground -0.23 0.15 -0.54 0.06 1.00
## muForaging_canopy_cover 0.20 0.14 -0.08 0.48 1.00
## muForaging_shrubs_pct 0.09 0.12 -0.16 0.33 1.00
## muForaging_decay2:GuildGroundomnivore -0.65 0.70 -2.05 0.73 1.00
## muForaging_decay3:GuildGroundomnivore 0.58 0.61 -0.64 1.78 1.00
## muForaging_decay4:GuildGroundomnivore 0.68 0.60 -0.49 1.87 1.00
## muForaging_decay5:GuildGroundomnivore 0.36 0.61 -0.86 1.54 1.00
## muForaging_decay2:GuildHerbivore 0.13 0.76 -1.41 1.58 1.00
## muForaging_decay3:GuildHerbivore 0.03 0.68 -1.30 1.35 1.00
## muForaging_decay4:GuildHerbivore -1.26 0.78 -2.83 0.24 1.00
## muForaging_decay5:GuildHerbivore -0.79 0.73 -2.21 0.65 1.00
## muForaging_decay2:GuildReptile -0.08 1.00 -2.03 1.87 1.00
## muForaging_decay3:GuildReptile -0.21 0.94 -2.10 1.62 1.00
## muForaging_decay4:GuildReptile -0.43 0.92 -2.25 1.36 1.00
## muForaging_decay5:GuildReptile -0.16 0.95 -2.01 1.68 1.00
## muForaging_decay2:GuildSciurid -0.18 0.79 -1.77 1.29 1.00
## muForaging_decay3:GuildSciurid 0.32 0.66 -0.98 1.61 1.00
## muForaging_decay4:GuildSciurid -0.18 0.66 -1.46 1.12 1.00
## muForaging_decay5:GuildSciurid 0.82 0.63 -0.41 2.08 1.00
## muForaging_decay2:GuildSmallrodent 0.42 0.63 -0.81 1.66 1.00
## muForaging_decay3:GuildSmallrodent -0.24 0.59 -1.40 0.91 1.00
## muForaging_decay4:GuildSmallrodent -0.33 0.58 -1.47 0.82 1.00
## muForaging_decay5:GuildSmallrodent -0.25 0.57 -1.37 0.89 1.00
## muForaging_decay2:GuildSongbird 0.09 0.62 -1.10 1.31 1.00
## muForaging_decay3:GuildSongbird -0.33 0.60 -1.49 0.85 1.00
## muForaging_decay4:GuildSongbird 0.66 0.57 -0.46 1.79 1.00
## muForaging_decay5:GuildSongbird 0.27 0.58 -0.88 1.41 1.00
## muTransit_decay2 -0.12 0.48 -1.05 0.82 1.00
## muTransit_decay3 0.05 0.44 -0.82 0.92 1.00
## muTransit_decay4 -0.51 0.48 -1.47 0.43 1.00
## muTransit_decay5 0.12 0.48 -0.81 1.07 1.00
## muTransit_GuildGroundomnivore 0.51 0.54 -0.53 1.57 1.00
## muTransit_GuildHerbivore -0.24 0.59 -1.37 0.96 1.00
## muTransit_GuildReptile 0.07 0.60 -1.09 1.27 1.00
## muTransit_GuildSciurid -0.46 0.57 -1.56 0.68 1.00
## muTransit_GuildSmallrodent 0.22 0.51 -0.74 1.22 1.00
## muTransit_GuildSongbird -2.24 0.55 -3.29 -1.10 1.00
## muTransit_off_ground -0.08 0.12 -0.31 0.16 1.00
## muTransit_canopy_cover 0.10 0.12 -0.14 0.33 1.00
## muTransit_shrubs_pct -0.16 0.10 -0.37 0.04 1.00
## muTransit_decay2:GuildGroundomnivore 0.32 0.55 -0.74 1.39 1.00
## muTransit_decay3:GuildGroundomnivore 0.24 0.48 -0.72 1.18 1.00
## muTransit_decay4:GuildGroundomnivore 0.55 0.52 -0.47 1.58 1.00
## muTransit_decay5:GuildGroundomnivore -0.09 0.52 -1.09 0.96 1.00
## muTransit_decay2:GuildHerbivore 0.88 0.67 -0.46 2.18 1.00
## muTransit_decay3:GuildHerbivore -0.16 0.58 -1.32 0.99 1.00
## muTransit_decay4:GuildHerbivore 0.59 0.66 -0.70 1.90 1.00
## muTransit_decay5:GuildHerbivore -0.18 0.61 -1.38 1.01 1.00
## muTransit_decay2:GuildReptile -0.42 0.72 -1.84 1.00 1.00
## muTransit_decay3:GuildReptile -0.39 0.62 -1.59 0.84 1.00
## muTransit_decay4:GuildReptile -0.54 0.59 -1.68 0.60 1.00
## muTransit_decay5:GuildReptile 0.30 0.69 -1.01 1.67 1.00
## muTransit_decay2:GuildSciurid 0.41 0.58 -0.73 1.52 1.00
## muTransit_decay3:GuildSciurid -0.12 0.52 -1.14 0.91 1.00
## muTransit_decay4:GuildSciurid -0.62 0.56 -1.71 0.47 1.00
## muTransit_decay5:GuildSciurid -1.08 0.54 -2.12 -0.02 1.00
## muTransit_decay2:GuildSmallrodent 0.19 0.45 -0.69 1.08 1.00
## muTransit_decay3:GuildSmallrodent -0.46 0.41 -1.26 0.36 1.00
## muTransit_decay4:GuildSmallrodent 0.25 0.46 -0.65 1.14 1.00
## muTransit_decay5:GuildSmallrodent -0.53 0.45 -1.39 0.35 1.00
## muTransit_decay2:GuildSongbird -0.56 0.62 -1.79 0.63 1.00
## muTransit_decay3:GuildSongbird 0.15 0.53 -0.90 1.17 1.00
## muTransit_decay4:GuildSongbird 0.10 0.57 -1.02 1.25 1.00
## muTransit_decay5:GuildSongbird -0.11 0.58 -1.23 1.03 1.00
## Bulk_ESS Tail_ESS
## muForaging_Intercept 5378 5717
## muTransit_Intercept 5691 5495
## muForaging_decay2 7746 6068
## muForaging_decay3 7215 6009
## muForaging_decay4 6632 5680
## muForaging_decay5 7463 5610
## muForaging_GuildGroundomnivore 9555 6172
## muForaging_GuildHerbivore 9112 5682
## muForaging_GuildReptile 12356 6178
## muForaging_GuildSciurid 8851 5137
## muForaging_GuildSmallrodent 7674 6004
## muForaging_GuildSongbird 6715 5873
## muForaging_off_ground 8192 5927
## muForaging_canopy_cover 8604 5950
## muForaging_shrubs_pct 8626 6033
## muForaging_decay2:GuildGroundomnivore 9170 6534
## muForaging_decay3:GuildGroundomnivore 9328 6409
## muForaging_decay4:GuildGroundomnivore 7748 5612
## muForaging_decay5:GuildGroundomnivore 8149 5882
## muForaging_decay2:GuildHerbivore 11393 6994
## muForaging_decay3:GuildHerbivore 9272 6738
## muForaging_decay4:GuildHerbivore 11882 5822
## muForaging_decay5:GuildHerbivore 11032 5790
## muForaging_decay2:GuildReptile 16456 5778
## muForaging_decay3:GuildReptile 16335 5598
## muForaging_decay4:GuildReptile 15654 5673
## muForaging_decay5:GuildReptile 16153 5931
## muForaging_decay2:GuildSciurid 11552 5970
## muForaging_decay3:GuildSciurid 9072 6551
## muForaging_decay4:GuildSciurid 8454 6553
## muForaging_decay5:GuildSciurid 8690 6152
## muForaging_decay2:GuildSmallrodent 7921 6067
## muForaging_decay3:GuildSmallrodent 7813 6123
## muForaging_decay4:GuildSmallrodent 7597 5893
## muForaging_decay5:GuildSmallrodent 7795 5996
## muForaging_decay2:GuildSongbird 8214 5858
## muForaging_decay3:GuildSongbird 7621 6253
## muForaging_decay4:GuildSongbird 7095 6168
## muForaging_decay5:GuildSongbird 7868 5619
## muTransit_decay2 6229 6075
## muTransit_decay3 5351 5407
## muTransit_decay4 5105 5815
## muTransit_decay5 5507 5583
## muTransit_GuildGroundomnivore 6108 5383
## muTransit_GuildHerbivore 8259 6462
## muTransit_GuildReptile 7736 5837
## muTransit_GuildSciurid 7025 5805
## muTransit_GuildSmallrodent 5988 6116
## muTransit_GuildSongbird 6314 6184
## muTransit_off_ground 5994 5197
## muTransit_canopy_cover 4999 5129
## muTransit_shrubs_pct 6179 6206
## muTransit_decay2:GuildGroundomnivore 8542 6695
## muTransit_decay3:GuildGroundomnivore 7229 6593
## muTransit_decay4:GuildGroundomnivore 6662 6441
## muTransit_decay5:GuildGroundomnivore 6681 5923
## muTransit_decay2:GuildHerbivore 8686 6338
## muTransit_decay3:GuildHerbivore 7612 6190
## muTransit_decay4:GuildHerbivore 7817 6522
## muTransit_decay5:GuildHerbivore 8588 6360
## muTransit_decay2:GuildReptile 11639 6453
## muTransit_decay3:GuildReptile 9512 6539
## muTransit_decay4:GuildReptile 7296 6582
## muTransit_decay5:GuildReptile 10611 6589
## muTransit_decay2:GuildSciurid 8137 6751
## muTransit_decay3:GuildSciurid 6468 5998
## muTransit_decay4:GuildSciurid 6375 6128
## muTransit_decay5:GuildSciurid 6981 6414
## muTransit_decay2:GuildSmallrodent 7132 6438
## muTransit_decay3:GuildSmallrodent 5827 5611
## muTransit_decay4:GuildSmallrodent 5738 5979
## muTransit_decay5:GuildSmallrodent 6121 5772
## muTransit_decay2:GuildSongbird 8142 7031
## muTransit_decay3:GuildSongbird 6469 5871
## muTransit_decay4:GuildSongbird 7583 6515
## muTransit_decay5:GuildSongbird 7546 6628
##
## 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).shared_formula <- bf(
Behavior ~ mo(decay, id = "d") * Guild + off_ground + canopy_cover + shrubs_pct +
(1 | Common_Name) + (1 | Plot) + (1 | Area) +
(1 | Camera.Type) + (1 | Month)
)
shared_simo <- as.data.frame(
get_prior(shared_formula, family = behavior_family, data = dat)) %>%
filter(class == "simo", coef != "")
shared_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"),
do.call(c, unname(Map(function(cf, dp) set_prior("dirichlet(1)", class = "simo",
coef = cf, dpar = dp),
shared_simo$coef, shared_simo$dpar))))
model_shared <- brm(
shared_formula,
family = behavior_family,
data = dat,
prior = shared_priors,
chains = 4, cores = 4, iter = 4000, warmup = 2000,
control = list(adapt_delta = 0.95, max_treedepth = 12)
)
summary(model_shared)Convergence
cat("Max Rhat: ", round(max(rhat(model_global), na.rm = TRUE), 4), "\n")## Max Rhat: 1.0038cat("Min bulk ESS:", round(min(posterior::summarise_draws(model_global)$ess_bulk,
na.rm = TRUE)), "\n")## Min bulk ESS: 1964cat("Divergent transitions:",
sum(subset(nuts_params(model_global), Parameter == "divergent__")$Value), "\n")## Divergent transitions: 5Helper Functions
# Probability of direction
pd_fn <- function(x) max(mean(x > 0), mean(x < 0))
# Display labels
common_names <- c(
"cotton_mouse" = "Cotton Mouse",
"eastern_woodrat" = "Eastern Woodrat",
"hispid_cotton_rat" = "Hispid Cotton Rat",
"eastern_chipmunk" = "Eastern Chipmunk",
"eastern_gray_squirrel" = "Eastern Gray Squirrel",
"nine_banded_armadillo" = "Nine-Banded Armadillo",
"virginia_opossum" = "Virginia Opossum",
"raccoon" = "Raccoon",
"bobcat" = "Bobcat",
"white_tailed_deer" = "White-Tailed Deer",
"eastern_cottontail" = "Eastern Cottontail",
"northern_cardinal" = "Northern Cardinal",
"carolina_wren" = "Carolina Wren",
"eastern_towhee" = "Eastern Towhee",
"indigo_bunting" = "Indigo Bunting",
"chucks_wills_widow" = "Chuck-Will's-Widow",
"five_lined_skink" = "Five-Lined Skink",
"broad_headed_skink" = "Broad-Headed Skink"
)
behavior_colors <- c("Local Search" = "#0072B2",
"Foraging" = "#009E73",
"Transit" = "#D55E00")
beh_factor <- function(x) {
factor(x,
levels = c("Local_Search", "Foraging", "Transit"),
labels = c("Local Search", "Foraging", "Transit"))
}
save_effect_docx <- function(df, file_name, caption_text) {
if (".category" %in% names(df)) {
df <- df %>%
mutate(Behavior = beh_factor(.category)) %>%
dplyr::select(-.category)
}
if ("Common_Name" %in% names(df)) {
df <- df %>%
mutate(Species = recode(as.character(Common_Name), !!!common_names)) %>%
dplyr::select(-Common_Name)
}
bold_flag <- df$pd >= 0.95
out <- df %>%
mutate(`Median (pp)` = round(median_pp, 2),
pd = round(pd, 2),
`Lower 95% CrI` = round(lower, 2),
`Upper 95% CrI` = round(upper, 2)) %>%
dplyr::select(any_of(c("Species", "Guild", "Order_grp", "Class", "Behavior")),
`Median (pp)`, pd, `Lower 95% CrI`, `Upper 95% CrI`)
ft <- flextable(out) %>%
set_caption(caption_text) %>%
bold(i = which(bold_flag), bold = TRUE) %>%
autofit()
save_as_docx(ft, path = file.path(out_dir, file_name))
list(table = out, flextable = ft)
}Parameter tables
make_pd_table <- function(model, file_name, caption_text) {
draws <- as_draws_df(model)
beta <- draws[, grep("^b_", colnames(draws)), drop = FALSE]
tab <- tibble::tibble(
parameter = gsub("^b_", "", colnames(beta)),
mean = colMeans(beta),
pd = apply(beta, 2, pd_fn),
lower = apply(beta, 2, quantile, probs = 0.025),
upper = apply(beta, 2, quantile, probs = 0.975)
) %>%
mutate(pd_95 = pd >= 0.95)
disp <- tab %>%
transmute(Parameter = parameter,
Mean = round(mean, 2),
pd = round(pd, 2),
`Lower 95% CrI` = round(lower, 2),
`Upper 95% CrI` = round(upper, 2))
ft <- flextable(disp) %>%
set_caption(caption_text) %>%
bold(i = which(tab$pd_95), bold = TRUE) %>%
autofit()
save_as_docx(ft, path = file.path(out_dir, file_name))
list(table = tab, flextable = ft)
}
pd_global <- make_pd_table(
model_global, "PD_Table_Behavior_Global.docx",
"Probability of Direction - decay (factor) x Guild + canopy + shrubs")## Warning: Dropping 'draws_df' class as required metadata was removed.Effect tables
newdata_profile <- dat %>%
distinct(Common_Name, Guild, Class, Order_grp) %>%
crossing(decay = factor(1:5, levels = 1:5, ordered = decay_ordered)) %>%
mutate(canopy_cover = 0, shrubs_pct = 0, off_ground = 0,
Plot = NA, Area = NA, Camera.Type = NA, Month = NA)
prof_draws <- add_epred_draws(newdata_profile, model_global,
re_formula = NA) %>%
ungroup()
summarise_prob <- function(df, ...) {
df %>%
group_by(.draw, .category, decay, ...) %>%
summarise(.epred = mean(.epred), .groups = "drop") %>%
group_by(.category, decay, ...) %>%
summarise(prob_pct = median(.epred) * 100,
lower = quantile(.epred, .025) * 100,
upper = quantile(.epred, .975) * 100,
.groups = "drop")
}
community_profile <- summarise_prob(prof_draws)
guild_profile <- summarise_prob(prof_draws, Guild)
print(community_profile)## # A tibble: 15 × 5
## .category decay prob_pct lower upper
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 Local_Search 1 49.0 25.2 73.3
## 2 Local_Search 2 49.9 28.0 72.5
## 3 Local_Search 3 52.2 28.7 74.1
## 4 Local_Search 4 54.0 28.2 74.9
## 5 Local_Search 5 50.4 26.1 72.9
## 6 Foraging 1 8.43 1.74 39.0
## 7 Foraging 2 5.86 1.21 29.8
## 8 Foraging 3 6.73 1.46 31.9
## 9 Foraging 4 11.0 2.75 43.0
## 10 Foraging 5 9.96 2.34 40.3
## 11 Transit 1 39.7 14.8 64.5
## 12 Transit 2 42.0 17.9 64.2
## 13 Transit 3 38.8 15.0 63.8
## 14 Transit 4 32.3 11.1 55.9
## 15 Transit 5 36.9 13.9 61.0save_profile_docx <- function(df, file_name, caption_text) {
out <- df %>%
mutate(Behavior = beh_factor(.category),
`Decay class` = as.integer(as.character(decay)),
`Probability (%)`= round(prob_pct, 1),
`Lower 95% CrI` = round(lower, 1),
`Upper 95% CrI` = round(upper, 1)) %>%
dplyr::select(any_of("Guild"), Behavior, `Decay class`,
`Probability (%)`, `Lower 95% CrI`, `Upper 95% CrI`) %>%
arrange(across(any_of("Guild")), Behavior, `Decay class`)
ft <- flextable(out) %>% set_caption(caption_text) %>% autofit()
save_as_docx(ft, path = file.path(out_dir, file_name))
out
}
save_profile_docx(community_profile, "Profile_Community_Decay.docx",
"Community-level predicted behavior probability by CWD decay class")## # A tibble: 15 × 5
## Behavior `Decay class` `Probability (%)` `Lower 95% CrI` `Upper 95% CrI`
## <fct> <int> <dbl> <dbl> <dbl>
## 1 Local Search 1 49 25.2 73.3
## 2 Local Search 2 49.9 28 72.5
## 3 Local Search 3 52.2 28.7 74.1
## 4 Local Search 4 54 28.2 74.9
## 5 Local Search 5 50.4 26.1 72.9
## 6 Foraging 1 8.4 1.7 39
## 7 Foraging 2 5.9 1.2 29.8
## 8 Foraging 3 6.7 1.5 31.9
## 9 Foraging 4 11 2.8 43
## 10 Foraging 5 10 2.3 40.3
## 11 Transit 1 39.7 14.8 64.5
## 12 Transit 2 42 17.9 64.2
## 13 Transit 3 38.8 15 63.8
## 14 Transit 4 32.3 11.1 55.9
## 15 Transit 5 36.9 13.9 61save_profile_docx(guild_profile, "Profile_Guild_Decay.docx",
"Guild-level predicted behavior probability by CWD decay class")## # A tibble: 105 × 6
## Guild Behavior `Decay class` `Probability (%)` `Lower 95% CrI`
## <fct> <fct> <int> <dbl> <dbl>
## 1 Carnivore Local Search 1 42.3 16.3
## 2 Carnivore Local Search 2 46.4 16.8
## 3 Carnivore Local Search 3 42.5 15.6
## 4 Carnivore Local Search 4 51.9 20.2
## 5 Carnivore Local Search 5 40 13.8
## 6 Carnivore Foraging 1 5.8 0.9
## 7 Carnivore Foraging 2 4.2 0.5
## 8 Carnivore Foraging 3 4.1 0.6
## 9 Carnivore Foraging 4 7.8 1.1
## 10 Carnivore Foraging 5 5.4 0.7
## # ℹ 95 more rows
## # ℹ 1 more variable: `Upper 95% CrI` <dbl>newdata_effect <- dat %>%
distinct(Common_Name, Guild, Class, Order_grp) %>%
crossing(decay_level = c("low", "high")) %>%
mutate(decay = factor(ifelse(decay_level == "low", 1L, 5L),
levels = 1:5, ordered = decay_ordered),
canopy_cover = 0, shrubs_pct = 0, off_ground = 0,
Plot = NA, Area = NA, Camera.Type = NA, Month = NA)
eff_draws <- add_epred_draws(newdata_effect, model_global, re_formula = NA) %>%
ungroup()
# Community-level (averaged over species)
community_effect <- eff_draws %>%
group_by(.draw, .category, decay_level) %>%
summarise(.epred = mean(.epred), .groups = "drop") %>%
pivot_wider(names_from = decay_level, values_from = .epred) %>%
mutate(delta = high - low) %>%
group_by(.category) %>%
summarise(median_pp = median(delta) * 100,
lower = quantile(delta, .025) * 100,
upper = quantile(delta, .975) * 100,
pd = pd_fn(delta), .groups = "drop")
print(community_effect)## # A tibble: 3 × 5
## .category median_pp lower upper pd
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Local_Search 1.11 -11.6 14.3 0.574
## 2 Foraging 1.24 -7.36 11.0 0.656
## 3 Transit -2.45 -16.8 10.8 0.648# Per guild
guild_effect <- eff_draws %>%
group_by(.draw, .category, Guild, decay_level) %>%
summarise(.epred = mean(.epred), .groups = "drop") %>%
pivot_wider(names_from = decay_level, values_from = .epred) %>%
mutate(delta = high - low) %>%
group_by(Guild, .category) %>%
summarise(median_pp = median(delta) * 100,
lower = quantile(delta, .025) * 100,
upper = quantile(delta, .975) * 100,
pd = pd_fn(delta), .groups = "drop")
order_effect <- eff_draws %>%
group_by(.draw, .category, Order_grp, decay_level) %>%
summarise(.epred = mean(.epred), .groups = "drop") %>%
pivot_wider(names_from = decay_level, values_from = .epred) %>%
mutate(delta = high - low) %>%
group_by(Order_grp, .category) %>%
summarise(median_pp = median(delta) * 100,
lower = quantile(delta, .025) * 100,
upper = quantile(delta, .975) * 100,
pd = pd_fn(delta), .groups = "drop")
class_effect <- eff_draws %>%
group_by(.draw, .category, Class, decay_level) %>%
summarise(.epred = mean(.epred), .groups = "drop") %>%
pivot_wider(names_from = decay_level, values_from = .epred) %>%
mutate(delta = high - low) %>%
group_by(Class, .category) %>%
summarise(median_pp = median(delta) * 100,
lower = quantile(delta, .025) * 100,
upper = quantile(delta, .975) * 100,
pd = pd_fn(delta), .groups = "drop")
# Per species, with the species random effects included
sp_draws <- add_epred_draws(newdata_effect, model_global,
re_formula = ~ (1 | Common_Name)) %>%
ungroup()
species_effect <- sp_draws %>%
pivot_wider(id_cols = c(.draw, .category, Common_Name),
names_from = decay_level, values_from = .epred) %>%
mutate(delta = high - low) %>%
group_by(Common_Name, .category) %>%
summarise(median_pp = median(delta) * 100,
lower = quantile(delta, .025) * 100,
upper = quantile(delta, .975) * 100,
pd = pd_fn(delta), .groups = "drop")
save_effect_docx(community_effect, "Effect_Community_Decay.docx",
"Community-level change in behavior probability from decay class 1 to 5")## $table
## # A tibble: 3 × 5
## Behavior `Median (pp)` pd `Lower 95% CrI` `Upper 95% CrI`
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Local Search 1.11 0.57 -11.6 14.3
## 2 Foraging 1.24 0.66 -7.36 11.0
## 3 Transit -2.45 0.65 -16.8 10.8
##
## $flextable
## a flextable object.
## col_keys: `Behavior`, `Median (pp)`, `pd`, `Lower 95% CrI`, `Upper 95% CrI`
## header has 1 row(s)
## body has 3 row(s)
## original dataset sample:
## 'data.frame': 3 obs. of 5 variables:
## $ Behavior : Factor w/ 3 levels "Local Search",..: 1 2 3
## $ Median (pp) : num 1.11 1.24 -2.45
## $ pd : num 0.57 0.66 0.65
## $ Lower 95% CrI: num -11.59 -7.36 -16.82
## $ Upper 95% CrI: num 14.3 11 10.8save_effect_docx(guild_effect, "Effect_Guild_Decay.docx",
"Guild-level change in behavior probability from decay class 1 to 5")## $table
## # A tibble: 21 × 6
## Guild Behavior `Median (pp)` pd `Lower 95% CrI` `Upper 95% CrI`
## <fct> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Carnivore Local Se… -2.07 0.59 -21.4 16.6
## 2 Carnivore Foraging -0.36 0.57 -11.6 10.4
## 3 Carnivore Transit 2.37 0.59 -18.8 24.1
## 4 Ground omnivore Local Se… -1.47 0.57 -21.2 17.9
## 5 Ground omnivore Foraging 1.45 0.67 -10.5 20.9
## 6 Ground omnivore Transit -1.12 0.54 -24.1 22.0
## 7 Herbivore Local Se… 5.7 0.67 -20.4 32.0
## 8 Herbivore Foraging -5.52 0.84 -34.4 9.83
## 9 Herbivore Transit 2.01 0.56 -26.2 29.7
## 10 Reptile Local Se… -7.9 0.71 -36.6 21.4
## # ℹ 11 more rows
##
## $flextable
## a flextable object.
## col_keys: `Guild`, `Behavior`, `Median (pp)`, `pd`, `Lower 95% CrI`, `Upper 95% CrI`
## header has 1 row(s)
## body has 21 row(s)
## original dataset sample:
## 'data.frame': 21 obs. of 6 variables:
## $ Guild : Factor w/ 7 levels "Carnivore","Ground omnivore",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ Behavior : Factor w/ 3 levels "Local Search",..: 1 2 3 1 2 3 1 2 3 1 ...
## $ Median (pp) : num -2.07 -0.36 2.37 -1.47 1.45 -1.12 5.7 -5.52 2.01 -7.9 ...
## $ pd : num 0.59 0.57 0.59 0.57 0.67 0.54 0.67 0.84 0.56 0.71 ...
## $ Lower 95% CrI: num -21.4 -11.6 -18.8 -21.2 -10.5 ...
## $ Upper 95% CrI: num 16.6 10.4 24.1 17.9 20.9 ...save_effect_docx(order_effect, "Effect_Order_Decay.docx",
paste("Order-level change in behavior probability from decay class 1 to 5.",
"Marginal summary of the guild model, not a separate fit; orders with",
"identical guild membership return identical values."))## $table
## # A tibble: 24 × 6
## Order_grp Behavior `Median (pp)` pd `Lower 95% CrI` `Upper 95% CrI`
## <fct> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Artiodactyla Local Se… 5.7 0.67 -20.4 32.0
## 2 Artiodactyla Foraging -5.52 0.84 -34.4 9.83
## 3 Artiodactyla Transit 2.01 0.56 -26.2 29.7
## 4 Carnivora Local Se… -1.76 0.59 -18.5 14.2
## 5 Carnivora Foraging 0.63 0.6 -9.02 13.3
## 6 Carnivora Transit 0.7 0.53 -18.3 19.9
## 7 Cingulata Local Se… -1.47 0.57 -21.2 17.9
## 8 Cingulata Foraging 1.45 0.67 -10.5 20.9
## 9 Cingulata Transit -1.12 0.54 -24.1 22.0
## 10 Didelphimorphia Local Se… -1.47 0.57 -21.2 17.9
## # ℹ 14 more rows
##
## $flextable
## a flextable object.
## col_keys: `Order_grp`, `Behavior`, `Median (pp)`, `pd`, `Lower 95% CrI`, `Upper 95% CrI`
## header has 1 row(s)
## body has 24 row(s)
## original dataset sample:
## 'data.frame': 24 obs. of 6 variables:
## $ Order_grp : Factor w/ 8 levels "Artiodactyla",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ Behavior : Factor w/ 3 levels "Local Search",..: 1 2 3 1 2 3 1 2 3 1 ...
## $ Median (pp) : num 5.7 -5.52 2.01 -1.76 0.63 0.7 -1.47 1.45 -1.12 -1.47 ...
## $ pd : num 0.67 0.84 0.56 0.59 0.6 0.53 0.57 0.67 0.54 0.57 ...
## $ Lower 95% CrI: num -20.39 -34.43 -26.2 -18.49 -9.02 ...
## $ Upper 95% CrI: num 32.01 9.83 29.7 14.15 13.28 ...save_effect_docx(class_effect, "Effect_Class_Decay.docx",
paste("Class-level change in behavior probability from decay class 1 to 5.",
"Marginal summary of the guild model, not a separate fit."))## $table
## # A tibble: 9 × 6
## Class Behavior `Median (pp)` pd `Lower 95% CrI` `Upper 95% CrI`
## <fct> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Aves Local Search -2.08 0.61 -19.3 14.5
## 2 Aves Foraging 1.74 0.67 -11.6 19.5
## 3 Aves Transit -0.17 0.52 -13.3 12.1
## 4 Mammalia Local Search 3.95 0.72 -10.0 18.3
## 5 Mammalia Foraging 1.21 0.64 -8.42 12.0
## 6 Mammalia Transit -5.47 0.75 -22.0 10.2
## 7 Reptilia Local Search -7.9 0.71 -36.6 21.4
## 8 Reptilia Foraging -0.41 0.65 -12.9 13.3
## 9 Reptilia Transit 8.39 0.7 -24.2 39.3
##
## $flextable
## a flextable object.
## col_keys: `Class`, `Behavior`, `Median (pp)`, `pd`, `Lower 95% CrI`, `Upper 95% CrI`
## header has 1 row(s)
## body has 9 row(s)
## original dataset sample:
## 'data.frame': 9 obs. of 6 variables:
## $ Class : Factor w/ 3 levels "Aves","Mammalia",..: 1 1 1 2 2 2 3 3 3
## $ Behavior : Factor w/ 3 levels "Local Search",..: 1 2 3 1 2 3 1 2 3
## $ Median (pp) : num -2.08 1.74 -0.17 3.95 1.21 -5.47 -7.9 -0.41 8.39
## $ pd : num 0.61 0.67 0.52 0.72 0.64 0.75 0.71 0.65 0.7
## $ Lower 95% CrI: num -19.33 -11.65 -13.31 -10.04 -8.42 ...
## $ Upper 95% CrI: num 14.5 19.5 12.1 18.3 12 ...save_effect_docx(species_effect, "Effect_Species_Decay.docx",
"Species-specific change in behavior probability from decay class 1 to 5")## $table
## # A tibble: 51 × 6
## Species Behavior `Median (pp)` pd `Lower 95% CrI` `Upper 95% CrI`
## <chr> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Bobcat Local S… -1.87 0.59 -20.4 16.1
## 2 Bobcat Foraging -0.2 0.57 -9.94 8.22
## 3 Bobcat Transit 2.25 0.59 -18.4 23.5
## 4 Broad-Headed Sk… Local S… -7.93 0.71 -37.6 21.8
## 5 Broad-Headed Sk… Foraging -0.18 0.65 -10.2 7.55
## 6 Broad-Headed Sk… Transit 8.53 0.71 -23.4 39.5
## 7 Carolina Wren Local S… -2.83 0.63 -21.4 16.1
## 8 Carolina Wren Foraging 3.15 0.67 -15.4 22.8
## 9 Carolina Wren Transit -0.25 0.53 -13.2 10.6
## 10 Cotton Mouse Local S… 9.18 0.87 -7.08 25.8
## # ℹ 41 more rows
##
## $flextable
## a flextable object.
## col_keys: `Species`, `Behavior`, `Median (pp)`, `pd`, `Lower 95% CrI`, `Upper 95% CrI`
## header has 1 row(s)
## body has 51 row(s)
## original dataset sample:
## 'data.frame': 51 obs. of 6 variables:
## $ Species : chr "Bobcat" "Bobcat" "Bobcat" "Broad-Headed Skink" ...
## $ Behavior : Factor w/ 3 levels "Local Search",..: 1 2 3 1 2 3 1 2 3 1 ...
## $ Median (pp) : num -1.87 -0.2 2.25 -7.93 -0.18 8.53 -2.83 3.15 -0.25 9.18 ...
## $ pd : num 0.59 0.57 0.59 0.71 0.65 0.71 0.63 0.67 0.53 0.87 ...
## $ Lower 95% CrI: num -20.43 -9.94 -18.35 -37.55 -10.17 ...
## $ Upper 95% CrI: num 16.12 8.22 23.47 21.78 7.55 ...Where along the gradient does the change happen?
if (decay_mode != "monotonic") {
cat("Skipped: simplexes exist only when decay_mode = 'monotonic'.\n")
} else {
guild_key <- setNames(levels(dat$Guild),
gsub("[^A-Za-z0-9]", "", levels(dat$Guild)))
ref_guild <- levels(dat$Guild)[1]
simo_draws <- as_draws_df(model_global) %>%
dplyr::select(starts_with("simo_")) %>%
pivot_longer(everything(), names_to = "param", values_to = "share") %>%
mutate(
dpar = sub("^simo_(mu[A-Za-z]+)_.*$", "\\1", param),
step = as.integer(gsub("^.*\\[|\\]$", "", param)),
term = sub("^simo_mu[A-Za-z]+_(.*[^0-9])[0-9]+\\[[0-9]+\\]$", "\\1", param),
guild_raw = ifelse(grepl("Guild", term), sub("^.*Guild", "", term), NA_character_),
Guild = ifelse(is.na(guild_raw),
paste0(ref_guild, " (reference)"),
unname(guild_key[guild_raw])),
Behavior = beh_factor(sub("^mu", "", dpar)),
Step = factor(step, levels = 1:4,
labels = c("1→2", "2→3", "3→4", "4→5")))
stopifnot(!any(is.na(simo_draws$Guild)), !any(is.na(simo_draws$step)),
all(simo_draws$step %in% 1:4))
simplex_summary <- simo_draws %>%
group_by(Behavior, Guild, Step) %>%
summarise(share = median(share),
lower = quantile(share, 0.025),
upper = quantile(share, 0.975),
.groups = "drop")
cat("Share of the total decay effect occurring at each step:\n")
print(as.data.frame(simplex_summary %>%
mutate(across(where(is.numeric), ~ round(.x, 3)))), row.names = FALSE)
bsp_summary <- as_draws_df(model_global) %>%
dplyr::select(starts_with("bsp_")) %>%
pivot_longer(everything(), names_to = "param", values_to = "value") %>%
mutate(dpar = sub("^bsp_(mu[A-Za-z]+)_.*$", "\\1", param),
term = sub("^bsp_mu[A-Za-z]+_", "", param),
Behavior = beh_factor(sub("^mu", "", dpar))) %>%
group_by(Behavior, term) %>%
summarise(mean = mean(value),
pd = pd_fn(value),
lower = quantile(value, 0.025),
upper = quantile(value, 0.975),
.groups = "drop")
cat("\nTotal monotonic decay effect (logit scale):\n")
print(as.data.frame(bsp_summary %>%
mutate(across(where(is.numeric), ~ round(.x, 2)))), row.names = FALSE)
ft_simplex <- flextable(simplex_summary %>%
mutate(Share = round(share, 3),
`Lower 95% CrI` = round(lower, 3),
`Upper 95% CrI` = round(upper, 3)) %>%
dplyr::select(Guild, Behavior, Step, Share,
`Lower 95% CrI`, `Upper 95% CrI`)) %>%
set_caption(paste("Share of each guild's total decay effect occurring",
"between adjacent decay classes.")) %>%
autofit()
save_as_docx(ft_simplex, path = file.path(out_dir, "Simplex_Guild_Decay.docx"))
ft_bsp <- flextable(bsp_summary %>%
mutate(Mean = round(mean, 2), pd = round(pd, 2),
`Lower 95% CrI` = round(lower, 2),
`Upper 95% CrI` = round(upper, 2)) %>%
dplyr::select(Term = term, Behavior, Mean, pd,
`Lower 95% CrI`, `Upper 95% CrI`)) %>%
set_caption("Total monotonic effect of decay class, logit scale.") %>%
bold(i = which(bsp_summary$pd >= 0.95), bold = TRUE) %>%
autofit()
save_as_docx(ft_bsp, path = file.path(out_dir, "Monotonic_Effect_Decay.docx"))}## Skipped: simplexes exist only when decay_mode = 'monotonic'.if (decay_mode != "monotonic") {
cat("Skipped: simplexes exist only when decay_mode = 'monotonic'.\n")
} else {
p_simplex <- ggplot(simplex_summary,
aes(Step, share, color = Behavior, group = Behavior)) +
geom_hline(yintercept = 0.25, linetype = 2, color = "grey60",
linewidth = 0.4) +
geom_line(position = position_dodge(width = 0.4), alpha = 0.5,
linewidth = 0.5) +
geom_linerange(aes(ymin = lower, ymax = upper),
position = position_dodge(width = 0.4), linewidth = 0.8) +
geom_point(position = position_dodge(width = 0.4), size = 2.1,
shape = 21, fill = "white", stroke = 0.7) +
facet_wrap(~ Guild) +
scale_color_manual(values = behavior_colors, name = "Behavior") +
scale_y_continuous(labels = label_percent(accuracy = 1)) +
labs(x = "Step between adjacent decay classes",
y = "Share of total decay effect",
caption = paste("Dashed line is 25%, the share expected if the effect",
"accrued evenly across the gradient.")) +
theme_classic(base_size = 10) +
theme(axis.title = element_text(face = "bold"),
plot.caption = element_text(hjust = 0, color = "grey35"))
print(p_simplex)
ggsave(file.path(out_dir, "Figure_Decay_Simplex.png"), p_simplex,
width = 9, height = 6, dpi = 600, bg = "white")}## Skipped: simplexes exist only when decay_mode = 'monotonic'.Pairwise contrasts between decay classes
# Posterior draws at every decay class
sp_prof_draws <- add_epred_draws(newdata_profile, model_global,
re_formula = ~ (1 | Common_Name)) %>%
ungroup()
# All pairwise differences among decay classes, hi - lo, in percentage points.
pairwise_decay <- function(dr, group_vars = character(0)) {
keys <- c(".draw", ".category", group_vars)
base <- dr %>%
group_by(across(all_of(c(keys, "decay")))) %>%
summarise(.epred = mean(.epred), .groups = "drop")
combos <- utils::combn(levels(droplevels(dr$decay)), 2, simplify = FALSE)
purrr::map_dfr(combos, function(p) {
lo <- p[1]; hi <- p[2]
a <- base %>% filter(decay == lo) %>%
dplyr::select(all_of(keys), .lo = .epred)
b <- base %>% filter(decay == hi) %>%
dplyr::select(all_of(keys), .hi = .epred)
dplyr::inner_join(a, b, by = keys) %>%
mutate(delta = .hi - .lo) %>%
group_by(across(all_of(c(".category", group_vars)))) %>%
summarise(class_hi = hi,
class_lo = lo,
median_pp = median(delta) * 100,
lower = quantile(delta, 0.025) * 100,
upper = quantile(delta, 0.975) * 100,
pd = pd_fn(delta),
.groups = "drop")
}) %>%
mutate(contrast = paste(class_hi, "vs", class_lo))
}
banded_decay <- function(dr, group_vars = character(0),
sound = c("1", "2"), advanced = c("4", "5")) {
keys <- c(".draw", ".category", group_vars)
dr %>%
filter(decay %in% c(sound, advanced)) %>%
mutate(band = ifelse(decay %in% sound, "sound", "advanced")) %>%
group_by(across(all_of(c(keys, "band")))) %>%
summarise(.epred = mean(.epred), .groups = "drop") %>%
pivot_wider(names_from = band, values_from = .epred) %>%
mutate(delta = advanced - sound) %>%
group_by(across(all_of(c(".category", group_vars)))) %>%
summarise(median_pp = median(delta) * 100,
lower = quantile(delta, 0.025) * 100,
upper = quantile(delta, 0.975) * 100,
pd = pd_fn(delta),
.groups = "drop")
}
decay_spread <- function(dr, group_vars = character(0)) {
keys <- c(".draw", ".category", group_vars)
dr %>%
group_by(across(all_of(c(keys, "decay")))) %>%
summarise(.epred = mean(.epred), .groups = "drop") %>%
group_by(across(all_of(keys))) %>%
summarise(rng = (max(.epred) - min(.epred)) * 100, .groups = "drop") %>%
group_by(across(all_of(c(".category", group_vars)))) %>%
summarise(median_range = median(rng),
lower = quantile(rng, 0.025),
upper = quantile(rng, 0.975),
.groups = "drop")
}How much data is behind each cell
decay_support <- dat %>%
group_by(Guild, decay) %>%
summarise(
n_obs = n(),
n_plots = n_distinct(Plot),
n_forage = sum(Behavior == "Foraging"),
n_plots_forage = n_distinct(Plot[Behavior == "Foraging"]),
top_plot_share = if (sum(Behavior == "Foraging") > 0) {
max(table(droplevels(factor(Plot[Behavior == "Foraging"])))) /
sum(Behavior == "Foraging")
} else NA_real_,
.groups = "drop")
cat("CWDs per decay class (the replication unit for decay):\n")## CWDs per decay class (the replication unit for decay):print(dat %>% distinct(Plot, decay) %>% count(decay, name = "n_logs"))## # A tibble: 5 × 2
## decay n_logs
## <fct> <int>
## 1 1 6
## 2 2 9
## 3 3 14
## 4 4 11
## 5 5 10cat("\nSupport per guild x decay cell:\n")##
## Support per guild x decay cell:print(as.data.frame(decay_support %>%
mutate(top_plot_share = round(top_plot_share, 2))), row.names = FALSE)## Guild decay n_obs n_plots n_forage n_plots_forage top_plot_share
## Carnivore 1 3 2 0 0 NA
## Carnivore 2 9 4 0 0 NA
## Carnivore 3 21 8 0 0 NA
## Carnivore 4 2 2 1 1 1.00
## Carnivore 5 6 3 0 0 NA
## Ground omnivore 1 27 6 2 2 0.50
## Ground omnivore 2 45 7 3 3 0.33
## Ground omnivore 3 110 13 13 8 0.23
## Ground omnivore 4 177 11 49 8 0.31
## Ground omnivore 5 94 9 17 6 0.41
## Herbivore 1 11 5 5 3 0.60
## Herbivore 2 15 7 2 2 0.50
## Herbivore 3 33 10 5 4 0.40
## Herbivore 4 15 8 0 0 NA
## Herbivore 5 24 9 2 2 0.50
## Reptile 1 9 1 0 0 NA
## Reptile 2 6 3 0 0 NA
## Reptile 3 17 4 0 0 NA
## Reptile 4 53 5 0 0 NA
## Reptile 5 11 4 0 0 NA
## Sciurid 1 16 4 0 0 NA
## Sciurid 2 27 6 1 1 1.00
## Sciurid 3 66 11 6 3 0.50
## Sciurid 4 80 7 6 4 0.33
## Sciurid 5 102 6 18 4 0.44
## Small rodent 1 125 6 3 3 0.33
## Small rodent 2 240 9 7 4 0.57
## Small rodent 3 491 13 14 6 0.29
## Small rodent 4 504 11 13 5 0.31
## Small rodent 5 517 10 23 4 0.61
## Songbird 1 37 6 2 2 0.50
## Songbird 2 89 9 9 5 0.33
## Songbird 3 131 14 11 6 0.36
## Songbird 4 129 11 27 9 0.30
## Songbird 5 76 10 14 7 0.21thin <- decay_support %>%
filter(n_forage == 0 | n_plots_forage < 3 |
(!is.na(top_plot_share) & top_plot_share > 0.5))
if (nrow(thin) > 0) {
cat("\nThinly supported foraging cells -- contrasts touching these should be",
"\nreported with the caveat, not as findings:\n")
print(as.data.frame(thin), row.names = FALSE)
}##
## Thinly supported foraging cells -- contrasts touching these should be
## reported with the caveat, not as findings:
## Guild decay n_obs n_plots n_forage n_plots_forage top_plot_share
## Carnivore 1 3 2 0 0 NA
## Carnivore 2 9 4 0 0 NA
## Carnivore 3 21 8 0 0 NA
## Carnivore 4 2 2 1 1 1.0000000
## Carnivore 5 6 3 0 0 NA
## Ground omnivore 1 27 6 2 2 0.5000000
## Herbivore 1 11 5 5 3 0.6000000
## Herbivore 2 15 7 2 2 0.5000000
## Herbivore 4 15 8 0 0 NA
## Herbivore 5 24 9 2 2 0.5000000
## Reptile 1 9 1 0 0 NA
## Reptile 2 6 3 0 0 NA
## Reptile 3 17 4 0 0 NA
## Reptile 4 53 5 0 0 NA
## Reptile 5 11 4 0 0 NA
## Sciurid 1 16 4 0 0 NA
## Sciurid 2 27 6 1 1 1.0000000
## Small rodent 2 240 9 7 4 0.5714286
## Small rodent 5 517 10 23 4 0.6086957
## Songbird 1 37 6 2 2 0.5000000ft_support <- flextable(decay_support %>%
mutate(top_plot_share = round(top_plot_share, 2)) %>%
rename(`Decay class` = decay, `Obs` = n_obs, `Logs` = n_plots,
`Foraging obs` = n_forage, `Logs w/ foraging` = n_plots_forage,
`Top log share` = top_plot_share)) %>%
set_caption(paste("Sampling support per guild x decay class.",
"CWDs, not detections, are the replication unit.")) %>%
autofit()
save_as_docx(ft_support, path = file.path(out_dir, "Support_Guild_Decay.docx"))All pairwise contrasts
pw_community <- pairwise_decay(prof_draws)
pw_guild <- pairwise_decay(prof_draws, "Guild")
pw_species <- pairwise_decay(sp_prof_draws, "Common_Name")
cat("Community-level pairwise contrasts:\n")## Community-level pairwise contrasts:print(as.data.frame(pw_community %>%
mutate(Behavior = beh_factor(.category)) %>%
dplyr::select(Behavior, contrast, median_pp, pd, lower, upper) %>%
mutate(across(where(is.numeric), ~ round(.x, 2)))), row.names = FALSE)## Behavior contrast median_pp pd lower upper
## Local Search 2 vs 1 1.03 0.56 -11.90 14.40
## Foraging 2 vs 1 -2.25 0.80 -15.37 3.98
## Transit 2 vs 1 2.13 0.62 -12.07 15.68
## Local Search 3 vs 1 3.10 0.69 -9.16 15.10
## Foraging 3 vs 1 -1.55 0.73 -12.79 4.65
## Transit 3 vs 1 -0.76 0.54 -13.61 12.00
## Local Search 4 vs 1 4.51 0.77 -8.00 16.96
## Foraging 4 vs 1 2.20 0.76 -5.86 13.12
## Transit 4 vs 1 -6.99 0.87 -20.36 5.10
## Local Search 5 vs 1 1.11 0.57 -11.59 14.29
## Foraging 5 vs 1 1.24 0.66 -7.36 11.04
## Transit 5 vs 1 -2.45 0.65 -16.82 10.82
## Local Search 3 vs 2 2.06 0.62 -10.93 14.85
## Foraging 3 vs 2 0.66 0.62 -6.66 9.56
## Transit 3 vs 2 -2.93 0.67 -16.35 10.60
## Local Search 4 vs 2 3.55 0.70 -10.41 16.32
## Foraging 4 vs 2 4.66 0.94 -1.56 19.36
## Transit 4 vs 2 -9.14 0.92 -22.33 3.45
## Local Search 5 vs 2 0.23 0.51 -14.04 14.06
## Foraging 5 vs 2 3.64 0.89 -2.83 17.47
## Transit 5 vs 2 -4.72 0.75 -19.12 9.54
## Local Search 4 vs 3 1.49 0.59 -11.09 13.58
## Foraging 4 vs 3 3.97 0.92 -1.90 16.34
## Transit 4 vs 3 -6.20 0.87 -17.89 4.68
## Local Search 5 vs 3 -1.85 0.62 -14.14 10.12
## Foraging 5 vs 3 2.89 0.87 -2.64 14.13
## Transit 5 vs 3 -1.89 0.62 -14.00 10.09
## Local Search 5 vs 4 -3.29 0.72 -14.72 8.31
## Foraging 5 vs 4 -0.93 0.64 -10.40 6.74
## Transit 5 vs 4 4.38 0.78 -6.61 16.01cat("\nGuild-level contrasts reaching pd >= 0.95:\n")##
## Guild-level contrasts reaching pd >= 0.95:print(as.data.frame(pw_guild %>%
filter(pd >= 0.95) %>%
mutate(Behavior = beh_factor(.category)) %>%
dplyr::select(Guild, Behavior, contrast, median_pp, pd, lower, upper) %>%
mutate(across(where(is.numeric), ~ round(.x, 2))) %>%
arrange(Guild, Behavior)), row.names = FALSE)## Guild Behavior contrast median_pp pd lower upper
## Ground omnivore Foraging 3 vs 2 4.16 0.95 -0.67 27.97
## Ground omnivore Foraging 4 vs 2 10.07 1.00 1.14 45.20
## Ground omnivore Foraging 5 vs 2 6.01 0.98 0.18 36.38
## Reptile Local Search 5 vs 4 -28.82 0.96 -57.35 3.16
## Reptile Transit 4 vs 1 -21.41 0.96 -46.52 2.72
## Reptile Transit 5 vs 4 30.37 0.97 -0.79 59.52
## Sciurid Foraging 5 vs 2 13.32 0.98 0.48 50.87
## Sciurid Foraging 5 vs 3 10.01 0.97 -0.82 37.97
## Sciurid Transit 4 vs 1 -18.63 0.97 -43.68 0.14
## Sciurid Transit 5 vs 1 -18.76 0.97 -44.32 0.25
## Sciurid Transit 4 vs 2 -26.34 0.99 -52.39 -4.31
## Sciurid Transit 5 vs 2 -26.57 0.99 -53.66 -3.67
## Sciurid Transit 4 vs 3 -17.24 0.98 -38.45 -1.37
## Sciurid Transit 5 vs 3 -17.33 0.98 -38.88 -1.62
## Songbird Foraging 4 vs 2 9.94 0.98 0.28 35.29
## Songbird Foraging 4 vs 3 13.00 1.00 2.02 41.24
## Songbird Foraging 5 vs 3 6.17 0.97 -0.26 30.09save_pairwise_docx <- function(df, file_name, caption_text) {
out <- df %>%
mutate(Behavior = beh_factor(.category), Contrast = contrast) %>%
{ if ("Common_Name" %in% names(.))
mutate(., Species = recode(as.character(Common_Name), !!!common_names)) else . } %>%
mutate(`Median (pp)` = round(median_pp, 2),
pd = round(pd, 2),
`Lower 95% CrI` = round(lower, 2),
`Upper 95% CrI` = round(upper, 2)) %>%
dplyr::select(any_of(c("Species", "Guild")), Behavior, Contrast,
`Median (pp)`, pd, `Lower 95% CrI`, `Upper 95% CrI`) %>%
arrange(across(any_of(c("Species", "Guild"))), Behavior, Contrast)
ft <- flextable(out) %>%
set_caption(caption_text) %>%
bold(i = which(out$pd >= 0.95), bold = TRUE) %>%
autofit()
save_as_docx(ft, path = file.path(out_dir, file_name))
out
}
save_pairwise_docx(pw_community, "Pairwise_Community_Decay.docx",
paste("Community-level pairwise contrasts between decay classes.",
"Ten comparisons per behavior; exploratory."))## # A tibble: 30 × 6
## Behavior Contrast `Median (pp)` pd `Lower 95% CrI` `Upper 95% CrI`
## <fct> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Local Search 2 vs 1 1.03 0.56 -11.9 14.4
## 2 Local Search 3 vs 1 3.1 0.69 -9.16 15.1
## 3 Local Search 3 vs 2 2.06 0.62 -10.9 14.8
## 4 Local Search 4 vs 1 4.51 0.77 -8 17.0
## 5 Local Search 4 vs 2 3.55 0.7 -10.4 16.3
## 6 Local Search 4 vs 3 1.49 0.59 -11.1 13.6
## 7 Local Search 5 vs 1 1.11 0.57 -11.6 14.3
## 8 Local Search 5 vs 2 0.23 0.51 -14.0 14.1
## 9 Local Search 5 vs 3 -1.85 0.62 -14.1 10.1
## 10 Local Search 5 vs 4 -3.29 0.72 -14.7 8.31
## # ℹ 20 more rowssave_pairwise_docx(pw_guild, "Pairwise_Guild_Decay.docx",
paste("Guild-level pairwise contrasts between decay classes.",
"Ten comparisons per behavior per guild; exploratory."))## # A tibble: 210 × 7
## Guild Behavior Contrast `Median (pp)` pd `Lower 95% CrI` `Upper 95% CrI`
## <fct> <fct> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Carniv… Local S… 2 vs 1 3.56 0.64 -15.9 23.0
## 2 Carniv… Local S… 3 vs 1 -0.01 0.5 -17.8 18.3
## 3 Carniv… Local S… 3 vs 2 -3.43 0.63 -25.3 19.2
## 4 Carniv… Local S… 4 vs 1 8.32 0.81 -10.5 28.1
## 5 Carniv… Local S… 4 vs 2 4.58 0.65 -18.6 29.0
## 6 Carniv… Local S… 4 vs 3 8.22 0.76 -14.3 31.1
## 7 Carniv… Local S… 5 vs 1 -2.07 0.59 -21.4 16.6
## 8 Carniv… Local S… 5 vs 2 -5.79 0.68 -29.5 18.2
## 9 Carniv… Local S… 5 vs 3 -2.2 0.58 -23.4 18.7
## 10 Carniv… Local S… 5 vs 4 -10.4 0.81 -34.2 12.5
## # ℹ 200 more rowssave_pairwise_docx(pw_species, "Pairwise_Species_Decay.docx",
"Species-level pairwise contrasts between decay classes.")## # A tibble: 510 × 7
## Species Behavior Contrast `Median (pp)` pd `Lower 95% CrI` `Upper 95% CrI`
## <chr> <fct> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Bobcat Local S… 2 vs 1 2.83 0.63 -15.2 22.4
## 2 Bobcat Local S… 3 vs 1 -0.42 0.52 -17.4 16.6
## 3 Bobcat Local S… 3 vs 2 -3.22 0.62 -25.1 17.9
## 4 Bobcat Local S… 4 vs 1 8.73 0.83 -9 29.2
## 5 Bobcat Local S… 4 vs 2 5.76 0.69 -16.8 30.2
## 6 Bobcat Local S… 4 vs 3 8.91 0.79 -11.8 32.8
## 7 Bobcat Local S… 5 vs 1 -1.87 0.59 -20.4 16.1
## 8 Bobcat Local S… 5 vs 2 -4.94 0.67 -28.2 18.3
## 9 Bobcat Local S… 5 vs 3 -1.61 0.57 -21.5 18.2
## 10 Bobcat Local S… 5 vs 4 -10.6 0.83 -34.7 11.2
## # ℹ 500 more rowsAdvanced versus sound logs, and overall spread
band_community <- banded_decay(prof_draws)
band_guild <- banded_decay(prof_draws, "Guild")
band_species <- banded_decay(sp_prof_draws, "Common_Name")
cat("Advanced (4-5) minus sound (1-2), community:\n")## Advanced (4-5) minus sound (1-2), community:print(as.data.frame(band_community %>%
mutate(Behavior = beh_factor(.category)) %>%
dplyr::select(Behavior, median_pp, pd, lower, upper) %>%
mutate(across(where(is.numeric), ~ round(.x, 2)))), row.names = FALSE)## Behavior median_pp pd lower upper
## Local Search 2.39 0.68 -7.68 12.16
## Foraging 2.96 0.90 -2.00 12.99
## Transit -5.94 0.88 -16.23 4.01cat("\nAdvanced (4-5) minus sound (1-2), by guild:\n")##
## Advanced (4-5) minus sound (1-2), by guild:print(as.data.frame(band_guild %>%
mutate(Behavior = beh_factor(.category)) %>%
dplyr::select(Guild, Behavior, median_pp, pd, lower, upper) %>%
mutate(across(where(is.numeric), ~ round(.x, 2))) %>%
arrange(Guild, Behavior)), row.names = FALSE)## Guild Behavior median_pp pd lower upper
## Carnivore Local Search 1.31 0.56 -14.15 16.85
## Carnivore Foraging 1.29 0.74 -4.86 13.97
## Carnivore Transit -3.40 0.66 -20.84 13.22
## Ground omnivore Local Search -1.91 0.60 -17.68 12.95
## Ground omnivore Foraging 5.74 0.97 -0.07 29.29
## Ground omnivore Transit -5.78 0.75 -24.51 10.89
## Herbivore Local Search 10.04 0.84 -9.89 30.49
## Herbivore Foraging -4.10 0.84 -26.63 6.30
## Herbivore Transit -4.04 0.64 -26.72 18.57
## Reptile Local Search 0.96 0.53 -21.19 24.10
## Reptile Foraging 0.07 0.53 -12.47 12.17
## Reptile Transit -1.07 0.53 -25.78 21.32
## Sciurid Local Search 14.20 0.88 -10.80 35.25
## Sciurid Foraging 7.27 0.96 -1.22 31.90
## Sciurid Transit -22.88 1.00 -42.78 -5.60
## Small rodent Local Search 7.93 0.92 -3.41 20.22
## Small rodent Foraging -0.23 0.57 -9.31 6.56
## Small rodent Transit -7.37 0.88 -20.83 5.48
## Songbird Local Search -6.67 0.89 -22.38 3.96
## Songbird Foraging 5.86 0.95 -1.22 23.68
## Songbird Transit 0.02 0.50 -8.12 8.06save_effect_docx(band_community, "Band_Community_Decay.docx",
"Community-level change in behavior probability, advanced (4-5) vs sound (1-2) logs")## $table
## # A tibble: 3 × 5
## Behavior `Median (pp)` pd `Lower 95% CrI` `Upper 95% CrI`
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Local Search 2.39 0.68 -7.68 12.2
## 2 Foraging 2.96 0.9 -2 13.0
## 3 Transit -5.94 0.88 -16.2 4.01
##
## $flextable
## a flextable object.
## col_keys: `Behavior`, `Median (pp)`, `pd`, `Lower 95% CrI`, `Upper 95% CrI`
## header has 1 row(s)
## body has 3 row(s)
## original dataset sample:
## 'data.frame': 3 obs. of 5 variables:
## $ Behavior : Factor w/ 3 levels "Local Search",..: 1 2 3
## $ Median (pp) : num 2.39 2.96 -5.94
## $ pd : num 0.68 0.9 0.88
## $ Lower 95% CrI: num -7.68 -2 -16.23
## $ Upper 95% CrI: num 12.16 12.99 4.01save_effect_docx(band_guild, "Band_Guild_Decay.docx",
"Guild-level change in behavior probability, advanced (4-5) vs sound (1-2) logs")## $table
## # A tibble: 21 × 6
## Guild Behavior `Median (pp)` pd `Lower 95% CrI` `Upper 95% CrI`
## <fct> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Carnivore Local Se… 1.31 0.56 -14.2 16.8
## 2 Ground omnivore Local Se… -1.91 0.6 -17.7 13.0
## 3 Herbivore Local Se… 10.0 0.84 -9.89 30.5
## 4 Reptile Local Se… 0.96 0.53 -21.2 24.1
## 5 Sciurid Local Se… 14.2 0.88 -10.8 35.2
## 6 Small rodent Local Se… 7.93 0.92 -3.41 20.2
## 7 Songbird Local Se… -6.67 0.89 -22.4 3.96
## 8 Carnivore Foraging 1.29 0.74 -4.86 14.0
## 9 Ground omnivore Foraging 5.74 0.97 -0.07 29.3
## 10 Herbivore Foraging -4.1 0.84 -26.6 6.3
## # ℹ 11 more rows
##
## $flextable
## a flextable object.
## col_keys: `Guild`, `Behavior`, `Median (pp)`, `pd`, `Lower 95% CrI`, `Upper 95% CrI`
## header has 1 row(s)
## body has 21 row(s)
## original dataset sample:
## 'data.frame': 21 obs. of 6 variables:
## $ Guild : Factor w/ 7 levels "Carnivore","Ground omnivore",..: 1 2 3 4 5 6 7 1 2 3 ...
## $ Behavior : Factor w/ 3 levels "Local Search",..: 1 1 1 1 1 1 1 2 2 2 ...
## $ Median (pp) : num 1.31 -1.91 10.04 0.96 14.2 ...
## $ pd : num 0.56 0.6 0.84 0.53 0.88 0.92 0.89 0.74 0.97 0.84 ...
## $ Lower 95% CrI: num -14.15 -17.68 -9.89 -21.19 -10.8 ...
## $ Upper 95% CrI: num 16.9 12.9 30.5 24.1 35.2 ...save_effect_docx(band_species, "Band_Species_Decay.docx",
"Species-level change in behavior probability, advanced (4-5) vs sound (1-2) logs")## $table
## # A tibble: 51 × 6
## Species Behavior `Median (pp)` pd `Lower 95% CrI` `Upper 95% CrI`
## <chr> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Bobcat Local S… 1.97 0.6 -12.9 17.6
## 2 Broad-Headed Sk… Local S… 1.12 0.54 -20.9 24.4
## 3 Carolina Wren Local S… -9.37 0.93 -24.7 3.22
## 4 Cotton Mouse Local S… 8.18 0.92 -3.51 20.4
## 5 Eastern Chipmunk Local S… 19.8 0.97 -1.58 38.5
## 6 Eastern Cottont… Local S… 11.8 0.87 -8.54 31.7
## 7 Eastern Gray Sq… Local S… 8.12 0.76 -15.6 29.0
## 8 Eastern Towhee Local S… -5.79 0.89 -21 4
## 9 Eastern Woodrat Local S… 8.39 0.92 -3.5 20.7
## 10 Five-Lined Skink Local S… 1.12 0.54 -20.8 24.8
## # ℹ 41 more rows
##
## $flextable
## a flextable object.
## col_keys: `Species`, `Behavior`, `Median (pp)`, `pd`, `Lower 95% CrI`, `Upper 95% CrI`
## header has 1 row(s)
## body has 51 row(s)
## original dataset sample:
## 'data.frame': 51 obs. of 6 variables:
## $ Species : chr "Bobcat" "Broad-Headed Skink" "Carolina Wren" "Cotton Mouse" ...
## $ Behavior : Factor w/ 3 levels "Local Search",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Median (pp) : num 1.97 1.12 -9.37 8.18 19.75 ...
## $ pd : num 0.6 0.54 0.93 0.92 0.97 0.87 0.76 0.89 0.92 0.54 ...
## $ Lower 95% CrI: num -12.89 -20.93 -24.72 -3.51 -1.58 ...
## $ Upper 95% CrI: num 17.55 24.42 3.22 20.43 38.54 ...spread_guild <- decay_spread(prof_draws, "Guild")
cat("\nPosterior spread across the five decay classes (pp).",
"\nRead the LOWER bound: it is the smallest spread the data are consistent with.\n")##
## Posterior spread across the five decay classes (pp).
## Read the LOWER bound: it is the smallest spread the data are consistent with.print(as.data.frame(spread_guild %>%
mutate(Behavior = beh_factor(.category)) %>%
dplyr::select(Guild, Behavior, median_range, lower, upper) %>%
mutate(across(where(is.numeric), ~ round(.x, 1))) %>%
arrange(Behavior, desc(median_range))), row.names = FALSE)## Guild Behavior median_range lower upper
## Reptile Local Search 36.3 13.9 59.9
## Sciurid Local Search 27.5 10.6 48.2
## Herbivore Local Search 26.9 9.0 50.0
## Carnivore Local Search 19.4 6.9 37.5
## Small rodent Local Search 16.8 6.3 30.2
## Songbird Local Search 16.2 4.9 38.7
## Ground omnivore Local Search 15.4 4.8 32.4
## Sciurid Foraging 15.4 2.4 52.1
## Herbivore Foraging 14.6 1.8 49.1
## Songbird Foraging 14.3 2.6 43.1
## Ground omnivore Foraging 11.3 1.5 46.4
## Carnivore Foraging 6.2 0.9 30.4
## Reptile Foraging 4.6 0.3 39.6
## Small rodent Foraging 4.0 0.5 22.5
## Reptile Transit 37.8 14.0 62.2
## Sciurid Transit 32.8 11.0 56.0
## Herbivore Transit 29.5 9.0 54.4
## Carnivore Transit 22.5 7.6 43.2
## Ground omnivore Transit 19.7 6.5 39.6
## Small rodent Transit 16.4 5.6 31.0
## Songbird Transit 9.1 1.9 27.2ft_spread <- flextable(spread_guild %>%
mutate(Behavior = beh_factor(.category),
`Median range (pp)` = round(median_range, 1),
`Lower 95% CrI` = round(lower, 1),
`Upper 95% CrI` = round(upper, 1)) %>%
dplyr::select(Guild, Behavior, `Median range (pp)`,
`Lower 95% CrI`, `Upper 95% CrI`)) %>%
set_caption(paste("Posterior spread (max - min) in behavior probability",
"across the five decay classes.")) %>%
autofit()
save_as_docx(ft_spread, path = file.path(out_dir, "Spread_Guild_Decay.docx"))Figures
Pairwise contrast matrices
pw_tile <- function(df, facet_var, fill_limit = NULL, ncol = NULL) {
d <- df %>%
mutate(Behavior = beh_factor(.category),
x = factor(class_lo, levels = levels(dat$decay)),
y = factor(class_hi, levels = rev(levels(dat$decay))),
credible = pd >= 0.95,
lab = sprintf("%+.0f\n%.2f", median_pp, pd))
lim <- if (is.null(fill_limit)) max(abs(d$median_pp), na.rm = TRUE) else fill_limit
ggplot(d, aes(x, y, fill = median_pp)) +
geom_tile(color = "white", linewidth = 0.6) +
geom_tile(data = dplyr::filter(d, credible),
color = "black", linewidth = 0.8, fill = NA) +
geom_text(aes(label = lab, fontface = ifelse(credible, "bold", "plain")),
size = 2.5, color = "grey10", lineheight = 0.9) +
scale_fill_gradient2(low = "#0072B2", mid = "#F0EFEC", high = "#D55E00",
midpoint = 0, limits = c(-lim, lim),
name = "Difference\n(pp)") +
facet_wrap(facet_var, ncol = ncol) +
labs(x = "Decay class (subtracted)", y = "Decay class") +
coord_equal() +
theme_minimal(base_size = 10) +
theme(panel.grid = element_blank(),
strip.text = element_text(face = "bold"))
}
p_pw_community <- pw_tile(pw_community, "Behavior")
print(p_pw_community)
ggsave(file.path(out_dir, "Figure_Pairwise_Community.png"), p_pw_community,
width = 9, height = 3.6, dpi = 600, bg = "white")
credible_guilds <- function(cat) {
pw_guild %>% filter(.category == cat, pd >= 0.95) %>%
pull(Guild) %>% unique() %>% sort()
}
all_guilds <- sort(unique(as.character(pw_guild$Guild)))
kf <- credible_guilds("Foraging"); if (!length(kf)) kf <- all_guilds
kt <- credible_guilds("Transit"); if (!length(kt)) kt <- all_guilds
cat("Foraging guilds retained:", paste(kf, collapse = ", "),
if (identical(kf, all_guilds)) " (no credible cell; showing all)" else "", "\n")## Foraging guilds retained: Ground omnivore, Sciurid, Songbirdcat("Transit guilds retained:", paste(kt, collapse = ", "),
if (identical(kt, all_guilds)) " (no credible cell; showing all)" else "", "\n")## Transit guilds retained: Reptile, Sciuridp_pw_forage <- pw_tile(pw_guild %>% filter(.category == "Foraging",
Guild %in% kf), "Guild",
ncol = max(1, length(kf)))
print(p_pw_forage)
ggsave(file.path(out_dir, "Figure_Pairwise_Guild_Foraging.png"), p_pw_forage,
width = max(4, 3 * length(kf)), height = 3.6, dpi = 600, bg = "white")
p_pw_transit <- pw_tile(pw_guild %>% filter(.category == "Transit",
Guild %in% kt), "Guild",
ncol = max(1, length(kt)))
print(p_pw_transit)
ggsave(file.path(out_dir, "Figure_Pairwise_Guild_Transit.png"), p_pw_transit,
width = max(4, 3.3 * length(kt)), height = 3.6, dpi = 600, bg = "white")
p_pw_forage_all <- pw_tile(pw_guild %>% filter(.category == "Foraging"), "Guild",
ncol = 4)
ggsave(file.path(out_dir, "FigureS_Pairwise_Guild_Foraging_all.png"),
p_pw_forage_all, width = 11, height = 6, dpi = 600, bg = "white")
p_pw_transit_all <- pw_tile(pw_guild %>% filter(.category == "Transit"), "Guild",
ncol = 4)
ggsave(file.path(out_dir, "FigureS_Pairwise_Guild_Transit_all.png"),
p_pw_transit_all, width = 11, height = 6, dpi = 600, bg = "white")Decay gradient, community
pred_grid <- dat %>%
distinct(Common_Name, Guild, Class, Order_grp) %>%
crossing(decay = factor(1:5, levels = 1:5, ordered = decay_ordered)) %>%
mutate(canopy_cover = 0, shrubs_pct = 0, off_ground = 0,
Plot = NA, Area = NA, Camera.Type = NA, Month = NA,
decay_x = as.integer(as.character(decay)))
behavior_summary <- add_epred_draws(pred_grid, model_global,
re_formula = NA, ndraws = 2000) %>%
ungroup() %>%
group_by(decay_x, .draw, .category) %>%
summarise(.epred = mean(.epred), .groups = "drop") %>%
group_by(decay_x, .category) %>%
median_qi(.epred, .width = 0.95) %>%
mutate(Behavior = beh_factor(.category))
plot_proportions <- dat %>%
mutate(decay_x = as.integer(as.character(decay))) %>%
group_by(Plot, decay_x, Behavior) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(Plot, decay_x) %>%
mutate(site_prob = n / sum(n)) %>%
ungroup() %>%
mutate(Behavior = beh_factor(Behavior))
pd_x <- position_dodge(width = 0.34)
p_community <- ggplot(behavior_summary,
aes(decay_x, .epred, color = Behavior, fill = Behavior)) +
geom_point(data = plot_proportions, inherit.aes = FALSE,
aes(decay_x, site_prob, fill = Behavior),
shape = 21, color = "black", alpha = 0.25, size = 2.2, stroke = 0.5,
position = position_jitter(width = 0.14, height = 0)) +
geom_line(linewidth = 0.5, alpha = 0.45, position = pd_x) +
geom_linerange(aes(ymin = .lower, ymax = .upper), linewidth = 0.9,
position = pd_x) +
geom_point(size = 2.6, shape = 21, color = "black", stroke = 0.6,
position = pd_x) +
scale_color_manual(values = behavior_colors, name = "Behavior") +
scale_fill_manual(values = behavior_colors, name = "Behavior") +
scale_y_continuous(limits = c(0, 1), labels = label_percent(accuracy = 1)) +
scale_x_continuous(breaks = 1:5) +
labs(x = "CWD Decay Class", y = "Predicted Probability") +
theme_classic(base_size = 12) +
theme(axis.title = element_text(face = "bold"))
print(p_community)
ggsave(file.path(out_dir, "Figure_Behavior_by_Decay.png"), p_community,
width = 7, height = 5, dpi = 600, bg = "white")Decay gradient, by guild
behavior_summary_guild <- add_epred_draws(pred_grid, model_global,
re_formula = NA, ndraws = 2000) %>%
ungroup() %>%
group_by(Guild, decay_x, .draw, .category) %>%
summarise(.epred = mean(.epred), .groups = "drop") %>%
group_by(Guild, decay_x, .category) %>%
median_qi(.epred, .width = 0.95) %>%
mutate(Behavior = beh_factor(.category))
p_guild <- ggplot(behavior_summary_guild,
aes(decay_x, .epred, color = Behavior, fill = Behavior)) +
geom_line(linewidth = 0.5, alpha = 0.45, position = pd_x) +
geom_linerange(aes(ymin = .lower, ymax = .upper), linewidth = 0.8,
position = pd_x) +
geom_point(size = 2.1, shape = 21, color = "black", stroke = 0.5,
position = pd_x) +
facet_wrap(~ Guild) +
scale_color_manual(values = behavior_colors, name = "Behavior") +
scale_fill_manual(values = behavior_colors, name = "Behavior") +
scale_y_continuous(limits = c(0, 1), labels = label_percent()) +
scale_x_continuous(breaks = 1:5) +
labs(x = "CWD Decay Class", y = "Behavior Probability") +
theme_classic(base_size = 11)
print(p_guild)
ggsave(file.path(out_dir, "Figure_Behavior_by_Guild.png"), p_guild,
width = 9, height = 6, dpi = 600, bg = "white")Decay gradient, by species
behavior_summary_sp <- add_epred_draws(pred_grid, model_global,
re_formula = ~ (1 | Common_Name),
ndraws = 2000) %>%
ungroup() %>%
group_by(Common_Name, decay_x, .category) %>%
median_qi(.epred, .width = 0.95) %>%
mutate(Behavior = beh_factor(.category),
Species = recode(as.character(Common_Name), !!!common_names))
p_species <- ggplot(behavior_summary_sp,
aes(decay_x, .epred, color = Behavior, fill = Behavior)) +
geom_line(linewidth = 0.45, alpha = 0.45, position = pd_x) +
geom_linerange(aes(ymin = .lower, ymax = .upper), linewidth = 0.6,
position = pd_x) +
geom_point(size = 1.6, shape = 21, color = "black", stroke = 0.4,
position = pd_x) +
facet_wrap(~ Species, ncol = 4) +
scale_color_manual(values = behavior_colors, name = "Behavior") +
scale_fill_manual(values = behavior_colors, name = "Behavior") +
scale_y_continuous(limits = c(0, 1), labels = label_percent()) +
scale_x_continuous(breaks = 1:5) +
labs(x = "CWD Decay Class", y = "Behavior Probability") +
theme_classic(base_size = 10)
print(p_species)
ggsave(file.path(out_dir, "Figure_Behavior_by_Species.png"), p_species,
width = 11, height = 9, dpi = 600, bg = "white")Decay class 5 minus decay class 1
contrast_draws <- eff_draws %>%
group_by(.draw, .category, decay_level) %>%
summarise(.epred = mean(.epred), .groups = "drop") %>%
pivot_wider(names_from = decay_level, values_from = .epred) %>%
mutate(delta_pp = (high - low) * 100,
Behavior = beh_factor(.category))
p_contrast <- ggplot(contrast_draws,
aes(x = delta_pp, y = fct_rev(Behavior), fill = Behavior)) +
geom_vline(xintercept = 0, linetype = 2, color = "grey40") +
stat_halfeye(.width = 0.95, point_interval = "median_qi",
slab_alpha = 0.75, normalize = "xy") +
geom_text(data = community_effect %>% mutate(Behavior = beh_factor(.category)),
aes(x = upper, y = fct_rev(Behavior),
label = sprintf("pd = %.2f", pd)),
inherit.aes = FALSE, hjust = -0.15, size = 3.3) +
scale_fill_manual(values = behavior_colors, guide = "none") +
scale_x_continuous(expand = expansion(mult = c(0.05, 0.18))) +
labs(x = "Decay class 5 \u2212 decay class 1 (percentage points)", y = NULL) +
theme_classic(base_size = 12) +
theme(axis.title = element_text(face = "bold"))
print(p_contrast)
ggsave(file.path(out_dir, "Figure_Decay_Contrast.png"), p_contrast,
width = 7, height = 4.5, dpi = 600, bg = "white")Species-level contrasts
p_sp_contrast <- species_effect %>%
mutate(Behavior = beh_factor(.category),
Species = recode(as.character(Common_Name), !!!common_names),
credible = factor(ifelse(pd >= 0.95, "yes", "no"),
levels = c("no", "yes"))) %>%
ggplot(aes(median_pp, fct_reorder(Species, median_pp),
color = Behavior, alpha = credible)) +
geom_vline(xintercept = 0, linetype = 2, color = "grey40") +
geom_pointrange(aes(xmin = lower, xmax = upper), size = 0.35) +
facet_wrap(~ Behavior, scales = "free_x") +
scale_color_manual(values = behavior_colors, guide = "none") +
scale_alpha_manual(values = c(no = 0.35, yes = 1), drop = FALSE,
name = "pd >= 0.95") +
labs(x = "Decay class 5 \u2212 decay class 1 (percentage points)", y = NULL) +
theme_classic(base_size = 10)
print(p_sp_contrast)
ggsave(file.path(out_dir, "Figure_Species_Decay_Contrast.png"), p_sp_contrast,
width = 11, height = 6, dpi = 600, bg = "white")Combined figure
combined_plot <- (p_community / p_guild) +
plot_layout(heights = c(1, 1.6)) +
plot_annotation(tag_levels = "A")
print(combined_plot)
ggsave(file.path(out_dir, "Figure_Behavior_Combined.png"), combined_plot,
width = 9, height = 11, dpi = 600, bg = "white")Sensitivity: air temperature
dat_temp <- dat %>%
filter(!is.na(Air.TemperatureF)) %>%
mutate(temp_z = as.numeric(scale(Air.TemperatureF)))
cat("Temperature subset:", nrow(dat_temp), "of", nrow(dat), "observations\n")
print(table(dat_temp$Camera.Type))
model_temp <- brm(
bf(Behavior ~ mo(decay) * Guild + canopy_cover + shrubs_pct + temp_z +
(1 | Common_Name) +(1 | Plot) + (1 | Area) +
(1 | Camera.Type) + (1 | Month)),
family = categorical(link = "logit", refcat = "Local_Search"),
data = dat_temp,
prior = priors,
chains = 4, cores = 4, iter = 4000, warmup = 2000,
control = list(adapt_delta = 0.95, max_treedepth = 12)
)
summary(model_temp)Save
saveRDS(
list(
model_name = paste0("decay (", decay_mode, ") x Guild + off_ground + ",
"canopy_cover + shrubs_pct"),
model_global = model_global,
dat = dat,
community_profile = community_profile,
guild_profile = guild_profile,
pw_community = pw_community,
pw_guild = pw_guild,
pw_species = pw_species,
band_community = band_community,
band_guild = band_guild,
band_species = band_species,
spread_guild = spread_guild,
decay_support = decay_support,
community_effect = community_effect,
guild_effect = guild_effect,
order_effect = order_effect,
class_effect = class_effect,
species_effect = species_effect,
cov_centers = cov_centers,
cov_scales = cov_scales,
guild_lookup = guild_lookup
),
file.path(out_dir, "behavior_results.rds")
)