Wildlife_CWD_Multi_Species_Occupancy
Alan Ivory
2026-10-06
knitr::opts_chunk$set(echo = TRUE)
pacman::p_load(tidyverse, readxl, lubridate, vegan,
spOccupancy, coda, MCMCvis, statmod,
patchwork, ggplot2, ggrepel,
kableExtra, flextable, officer)
set.seed(97)
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", "occupancy")
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)Analysis settings
occasion_days <- 7 # length of a survey occasion, in days
require_full <- TRUE # keep a plot-occasion only if the camera ran all 7 days
min_det <- 10 # species need this many independent detections
min_plots <- 5 # species need to be present at this many plots
max_naive_occ <- 0.90 # species must not be found at more than this fraction
mammals_only <- FALSE # TRUE restricts the community to Mammalia
gap_mins <- 30 # independence filter
# MCMC settings
n.samples <- 10000
n.burn <- 3000
n.thin <- 2
n.chains <- 4
n.report <- 1000Data Preparation
Clean the camera observations
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)),
Camera.Type = str_trim(as.character(Camera.Type)),
Date = as.Date(Date),
Time_chr = format(Time, "%H:%M:%S"),
DateTime = as.POSIXct(paste(Date, Time_chr), tz = "America/New_York"),
Year = year(Date),
Area = str_extract(Plot, "^[A-Z]+")
)
# Effort matrix to long format: one row per deployed plot-day
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, ## only use observations from 2026 (exclude observations from 2020-2023)
!is.na(Common_Name), !is.na(DateTime)) %>%
semi_join(effort_long, by = c("Plot", "Date"))Independence 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 detections:", nrow(CamInd), "\n")## Independent detections: 3530Survey occasions
Weekly windows sampling “occurrences”.
start_date <- min(effort_long$Date)
end_date <- max(effort_long$Date)
n_occ <- floor(as.numeric(end_date - start_date + 1) / occasion_days)
assign_occ <- function(d) floor(as.numeric(d - start_date) / occasion_days) + 1
effort_occ <- effort_long %>%
mutate(occ = assign_occ(Date)) %>%
filter(occ <= n_occ) %>%
count(Plot, occ, name = "days_active")
cat("Occasions:", n_occ, "| window:", format(start_date), "to", format(end_date), "\n")## Occasions: 13 | window: 2026-04-06 to 2026-07-07print(table(effort_occ$days_active))##
## 1 2 3 4 5 6 7
## 1 8 18 11 28 17 544if (require_full) {
effort_occ <- effort_occ %>% filter(days_active == occasion_days)
}
cat("Plot-occasions retained:", nrow(effort_occ), "of", 50 * n_occ, "\n")## Plot-occasions retained: 544 of 650Species selection
# Naive occupancy on the retained (full-week) data
naive_tab <- CamInd %>%
mutate(occ = assign_occ(Date)) %>%
semi_join(effort_occ, by = c("Plot", "occ")) %>%
distinct(Common_Name, Plot) %>%
count(Common_Name, name = "n_plots_full") %>%
mutate(naive_occ = n_plots_full / length(unique(effort_occ$Plot)))
sp_summary <- CamInd %>%
{ if (mammals_only) filter(., Class == "Mammalia") else . } %>%
group_by(Common_Name) %>%
summarise(n_det = n(), n_plots = n_distinct(Plot),
Class = first(na.omit(Class)), .groups = "drop") %>%
left_join(naive_tab, by = "Common_Name") %>%
mutate(naive_occ = replace_na(naive_occ, 0)) %>%
arrange(desc(naive_occ))
species_subset <- sp_summary %>%
filter(n_det >= min_det, n_plots >= min_plots, naive_occ <= max_naive_occ) %>%
pull(Common_Name)
# Report what was dropped for being near-ubiquitous
sp_summary %>%
filter(n_det >= min_det, n_plots >= min_plots, naive_occ > max_naive_occ) %>%
dplyr::select(Common_Name, n_det, naive_occ) %>%
print()## # A tibble: 2 × 3
## Common_Name n_det naive_occ
## <chr> <int> <dbl>
## 1 northern_cardinal 272 0.94
## 2 cotton_mouse 1357 0.92cat("\nSpecies retained:", length(species_subset), "of", nrow(sp_summary), "\n")##
## Species retained: 21 of 54print(as.data.frame(sp_summary %>% filter(Common_Name %in% species_subset) %>%
dplyr::select(Common_Name, Class, n_det, naive_occ)),
row.names = FALSE)## Common_Name Class n_det naive_occ
## nine_banded_armadillo Mammalia 212 0.80
## carolina_wren Aves 113 0.66
## eastern_woodrat Mammalia 413 0.60
## virginia_opossum Mammalia 130 0.60
## hispid_cotton_rat Mammalia 107 0.54
## raccoon Mammalia 111 0.50
## eastern_chipmunk Mammalia 187 0.48
## eastern_gray_squirrel Mammalia 104 0.46
## white_tailed_deer Mammalia 44 0.46
## eastern_cottontail Mammalia 54 0.44
## bobcat Mammalia 41 0.38
## eastern_towhee Aves 43 0.30
## broad_headed_skink Reptilia 36 0.26
## five_lined_skink Reptilia 60 0.20
## indigo_bunting Aves 34 0.20
## northern_house_wren Aves 14 0.18
## southern_black_racer Reptilia 12 0.18
## brown_thrasher Aves 13 0.14
## southern_flying_squirrel Mammalia 10 0.14
## gray_catbird Aves 14 0.12
## six_lined_racerunner Reptilia 16 0.10Create Detection Matrix
site_names <- sort(unique(effort_long$Plot))
J <- length(site_names)
K <- n_occ
N <- length(species_subset)
# Grid of every plot x occasion, flagged retained / not
grid <- expand_grid(Plot = site_names, occ = 1:K) %>%
left_join(effort_occ, by = c("Plot", "occ")) %>%
mutate(days_active = replace_na(days_active, 0),
sampled = days_active > 0)
samp_mat <- matrix(grid$sampled[order(match(grid$Plot, site_names), grid$occ)],
nrow = J, ncol = K, byrow = TRUE,
dimnames = list(site_names, paste0("wk", 1:K)))
det_long <- CamInd %>%
filter(Common_Name %in% species_subset) %>%
mutate(occ = assign_occ(Date)) %>%
semi_join(effort_occ, by = c("Plot", "occ")) %>%
distinct(Common_Name, Plot, occ)
detection_array <- array(NA_integer_, dim = c(N, J, K),
dimnames = list(species = species_subset,
site = site_names,
week = paste0("wk", 1:K)))
for (s in seq_along(species_subset)) {
m <- matrix(0L, nrow = J, ncol = K,
dimnames = list(site_names, paste0("wk", 1:K)))
d <- det_long %>% filter(Common_Name == species_subset[s])
if (nrow(d)) m[cbind(match(d$Plot, site_names), d$occ)] <- 1L
m[!samp_mat] <- NA_integer_ # not deployed -> missing, not absent
detection_array[s, , ] <- m
}
dim(detection_array) # Should be num_species x num_sites x num_weeks## [1] 21 50 13cat("\nNaive occupancy:\n")##
## Naive occupancy:print(round(sort(apply(detection_array, 1,
function(x) mean(rowSums(x, na.rm = TRUE) > 0)),
decreasing = TRUE), 3))## nine_banded_armadillo carolina_wren eastern_woodrat
## 0.80 0.66 0.60
## virginia_opossum hispid_cotton_rat raccoon
## 0.60 0.54 0.50
## eastern_chipmunk eastern_gray_squirrel white_tailed_deer
## 0.48 0.46 0.46
## eastern_cottontail bobcat eastern_towhee
## 0.44 0.38 0.30
## broad_headed_skink five_lined_skink indigo_bunting
## 0.26 0.20 0.20
## northern_house_wren southern_black_racer brown_thrasher
## 0.18 0.18 0.14
## southern_flying_squirrel gray_catbird six_lined_racerunner
## 0.14 0.12 0.10cat("\nProportion of y that is NA (unsampled):",
round(mean(is.na(detection_array)), 3), "\n")##
## Proportion of y that is NA (unsampled): 0.163Site covariates
Daubenmire cover classes are converted to their midpoint percentages before averaging, and tree density and basal area come from the point-centered-quarter data by the Cottam and Curtis estimator.
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")
# Ground-cover heterogeneity
site_struct$gc_shannon <- vegan::diversity(
as.matrix(site_struct %>% dplyr::select(ends_with("_pct"))), index = "shannon")
# Tree metrics from point-centered quarter data
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) %>%
mutate(Area = str_extract(Plot, "^[A-Z]+")) %>%
filter(Plot %in% site_names) %>%
arrange(match(Plot, site_names))
stopifnot(identical(site_env$Plot, site_names))Distribution and correlation of covariates
occ_vars <- c("decay", "canopy_cover", "soil_moist", "gc_shannon",
"lg_pct", "shrubs_pct", "sdw_pct", "basal_area", "tree_dens",
"off_ground")
cat("CWDs by proportion of quadrats off the ground:\n")## CWDs by proportion of quadrats off the ground:print(table(round(site_env$off_ground, 3)))##
## 0 0.125 0.25 0.333 0.5 0.75 1
## 34 2 3 1 4 2 4cat(sprintf("CWDs with any off-ground quadrat: %d of %d\n",
sum(site_env$off_ground > 0), nrow(site_env)))## CWDs with any off-ground quadrat: 16 of 50og_vals <- sort(unique(Master$og))
if (any(!og_vals %in% c(0, 1))) {
cat("\nNOTE: 'og' contains undocumented value(s):", setdiff(og_vals, c(0, 1)),
"\n Logs affected:",
paste(unique(Master$cwd[!Master$og %in% c(0, 1)]), collapse = ", "),
"\n These are being treated as off-ground. Confirm that is correct.\n")
}##
## NOTE: 'og' contains undocumented value(s): 2
## Logs affected: CE10
## These are being treated as off-ground. Confirm that is correct.cor_tab <- print(round(cor(site_env[, occ_vars], use = "pairwise.complete.obs"), 2))## decay canopy_cover soil_moist gc_shannon lg_pct shrubs_pct sdw_pct
## decay 1.00 0.27 0.00 -0.20 0.00 0.10 0.12
## canopy_cover 0.27 1.00 -0.14 -0.62 -0.54 -0.05 0.14
## soil_moist 0.00 -0.14 1.00 0.08 0.13 -0.16 -0.05
## gc_shannon -0.20 -0.62 0.08 1.00 0.49 0.20 0.15
## lg_pct 0.00 -0.54 0.13 0.49 1.00 -0.19 -0.18
## shrubs_pct 0.10 -0.05 -0.16 0.20 -0.19 1.00 0.24
## sdw_pct 0.12 0.14 -0.05 0.15 -0.18 0.24 1.00
## basal_area -0.07 0.41 0.00 -0.14 -0.12 -0.33 -0.06
## tree_dens 0.10 0.36 -0.07 -0.43 -0.17 -0.21 -0.05
## off_ground -0.37 -0.12 0.15 0.06 -0.10 0.04 -0.22
## basal_area tree_dens off_ground
## decay -0.07 0.10 -0.37
## canopy_cover 0.41 0.36 -0.12
## soil_moist 0.00 -0.07 0.15
## gc_shannon -0.14 -0.43 0.06
## lg_pct -0.12 -0.17 -0.10
## shrubs_pct -0.33 -0.21 0.04
## sdw_pct -0.06 -0.05 -0.22
## basal_area 1.00 0.41 0.20
## tree_dens 0.41 1.00 -0.01
## off_ground 0.20 -0.01 1.00par(mfrow = c(3, 3))
hist(site_env$decay, main = "CWD Decay Class", xlab = "Decay class (1-5)")
hist(site_env$canopy_cover, main = "Canopy Cover", xlab = "Canopy cover (%)")
hist(site_env$soil_moist, main = "Soil Moisture", xlab = "Soil moisture")
hist(site_env$gc_shannon, main = "Ground Cover H'", xlab = "Shannon index")
hist(site_env$lg_pct, main = "Live Grass Cover", xlab = "Cover (%)")
hist(site_env$shrubs_pct, main = "Shrub Cover", xlab = "Cover (%)")
hist(site_env$basal_area, main = "Basal Area", xlab = "m2 / ha")
hist(site_env$tree_dens, main = "Tree Density", xlab = "Trees per ha")
hist(site_env$off_ground, main = "CWD Off Ground", xlab = "Prop. quadrats")
par(mfrow = c(1, 1))Prepare Site Covariates
Covariates are scaled and centered so the model treats them on a comparable scale. The centers and scales are kept so that the figures can put predictions back on real units (decay class 1-5, canopy cover 0-100%).
Random-effect grouping variables must be numeric in spOccupancy. The label lookups are kept so figures can show CE / CW / D
covariates_matrix <- scale(as.matrix(site_env[, occ_vars]))
occ_centers <- attr(covariates_matrix, "scaled:center")
occ_scales <- attr(covariates_matrix, "scaled:scale")
occ.covs <- as.data.frame(covariates_matrix)
occ.covs$decay_sq <- occ.covs$decay^2 # quadratic term
area_levels <- sort(unique(site_env$Area))
occ.covs$Area <- as.numeric(factor(site_env$Area, levels = area_levels))
rownames(occ.covs) <- site_names
str(occ.covs)## 'data.frame': 50 obs. of 12 variables:
## $ decay : num 1.391 0.618 1.391 -0.155 0.618 ...
## $ canopy_cover: num -0.104 0.707 1.519 1.438 1.722 ...
## $ soil_moist : num 2.392 0.186 -0.639 -1.456 0.312 ...
## $ gc_shannon : num 0.312 -2.424 -1.36 -2.692 -2.854 ...
## $ lg_pct : num -0.835 -0.438 -1.074 -1.153 -1.114 ...
## $ shrubs_pct : num 0.607 -1.086 1.125 -0.522 0.419 ...
## $ sdw_pct : num 3.389 -0.833 0.452 0.269 -0.236 ...
## $ basal_area : num -0.5186 0.3688 -0.2768 0.6845 -0.0356 ...
## $ tree_dens : num -0.356 -0.423 -0.737 4.387 1.12 ...
## $ off_ground : num -0.559 -0.559 -0.559 -0.559 0.232 ...
## $ decay_sq : num 1.9361 0.3824 1.9361 0.0239 0.3824 ...
## $ Area : num 1 1 1 1 1 1 1 1 1 1 ...Prepare Detection Covariates
Week enters as a linear and a quadratic term to allow a seasonal peak.
site_id is a site-level random intercept on detection.
It absorbs CWD-to-CWD variation in how intensively a resident animal
uses its CWD
week_mat <- matrix(rep(1:K, each = J), nrow = J, ncol = K)
week_z <- (week_mat - mean(1:K)) / sd(1:K)
# Dominant camera type per plot (CW03 and D11 each swapped mid-season)
cam_type <- CameraData %>%
filter(Plot %in% site_names, Year == 2026) %>%
count(Plot, Camera.Type) %>%
group_by(Plot) %>%
slice_max(n, n = 1, with_ties = FALSE) %>%
ungroup() %>%
arrange(match(Plot, site_names))
cam_levels <- sort(unique(cam_type$Camera.Type))
det.covs <- list(
week = week_z,
week_sq = week_z^2,
off_ground = as.numeric(occ.covs$off_ground),
cam = as.numeric(factor(cam_type$Camera.Type, levels = cam_levels)),
site_id = seq_len(J),
decay_det = as.numeric(occ.covs$decay),
canopy_det = as.numeric(occ.covs$canopy_cover)
)
str(det.covs)## List of 7
## $ week : num [1:50, 1:13] -1.54 -1.54 -1.54 -1.54 -1.54 ...
## $ week_sq : num [1:50, 1:13] 2.37 2.37 2.37 2.37 2.37 ...
## $ off_ground: num [1:50] -0.559 -0.559 -0.559 -0.559 0.232 ...
## $ cam : num [1:50] 2 1 2 2 2 1 1 1 2 2 ...
## $ site_id : int [1:50] 1 2 3 4 5 6 7 8 9 10 ...
## $ decay_det : num [1:50] 1.391 0.618 1.391 -0.155 0.618 ...
## $ canopy_det: num [1:50] -0.104 0.707 1.519 1.438 1.722 ...Combine data into a list
data_list <- list(
y = detection_array,
occ.covs = occ.covs,
det.covs = det.covs
)
stopifnot(dim(data_list$y)[2] == nrow(data_list$occ.covs))
stopifnot(all(sapply(det.covs[c("week", "week_sq")],
function(m) all(dim(m) == c(J, K)))))
stopifnot(all(sapply(det.covs[c("off_ground", "cam", "site_id",
"decay_det", "canopy_det")], length) == J))
stopifnot(all(apply(data_list$y, 1, function(x) sum(x, na.rm = TRUE)) > 0))
stopifnot(is.numeric(occ.covs$Area), is.numeric(det.covs$cam))
stopifnot(all(sapply(seq_len(N),
function(i) all(is.na(data_list$y[i, , ]) == !samp_mat))))
cat("All structural checks passed.\n")## All structural checks passed.str(data_list, max.level = 2)## List of 3
## $ y : int [1:21, 1:50, 1:13] NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "dimnames")=List of 3
## $ occ.covs:'data.frame': 50 obs. of 12 variables:
## ..$ decay : num [1:50] 1.391 0.618 1.391 -0.155 0.618 ...
## ..$ canopy_cover: num [1:50] -0.104 0.707 1.519 1.438 1.722 ...
## ..$ soil_moist : num [1:50] 2.392 0.186 -0.639 -1.456 0.312 ...
## ..$ gc_shannon : num [1:50] 0.312 -2.424 -1.36 -2.692 -2.854 ...
## ..$ lg_pct : num [1:50] -0.835 -0.438 -1.074 -1.153 -1.114 ...
## ..$ shrubs_pct : num [1:50] 0.607 -1.086 1.125 -0.522 0.419 ...
## ..$ sdw_pct : num [1:50] 3.389 -0.833 0.452 0.269 -0.236 ...
## ..$ basal_area : num [1:50] -0.5186 0.3688 -0.2768 0.6845 -0.0356 ...
## ..$ tree_dens : num [1:50] -0.356 -0.423 -0.737 4.387 1.12 ...
## ..$ off_ground : num [1:50] -0.559 -0.559 -0.559 -0.559 0.232 ...
## ..$ decay_sq : num [1:50] 1.9361 0.3824 1.9361 0.0239 0.3824 ...
## ..$ Area : num [1:50] 1 1 1 1 1 1 1 1 1 1 ...
## $ det.covs:List of 7
## ..$ week : num [1:50, 1:13] -1.54 -1.54 -1.54 -1.54 -1.54 ...
## ..$ week_sq : num [1:50, 1:13] 2.37 2.37 2.37 2.37 2.37 ...
## ..$ off_ground: num [1:50] -0.559 -0.559 -0.559 -0.559 0.232 ...
## ..$ cam : num [1:50] 2 1 2 2 2 1 1 1 2 2 ...
## ..$ site_id : int [1:50] 1 2 3 4 5 6 7 8 9 10 ...
## ..$ decay_det : num [1:50] 1.391 0.618 1.391 -0.155 0.618 ...
## ..$ canopy_det: num [1:50] -0.104 0.707 1.519 1.438 1.722 ...Objective: Do CWD decay class and nearby vegetation determine which species use a CWD?
Site Covariates
# Detection formulas
det.null <- ~ (1 | site_id)
det.week <- ~ week + week_sq + (1 | cam) + (1 | site_id)
det.ground <- ~ off_ground + (1 | cam) + (1 | site_id)
det.wk.ground <- ~ week + off_ground + (1 | cam) + (1 | site_id)
det.week.only <- ~ week + week_sq + (1 | site_id)
det.full <- ~ week + week_sq + off_ground + (1 | cam) + (1 | site_id)
# Occupancy formulas
occ.null <- ~ 1
occ.decay <- ~ decay + (1 | Area)
occ.decay.quad <- ~ decay + decay_sq + (1 | Area)
occ.canopy <- ~ canopy_cover + (1 | Area)
occ.log.stand <- ~ decay + canopy_cover + (1 | Area)
occ.conceal <- ~ shrubs_pct + lg_pct + gc_shannon + (1 | Area)
occ.structure <- ~ canopy_cover + basal_area + tree_dens + (1 | Area)
occ.global <- ~ decay + canopy_cover + soil_moist + gc_shannon +
shrubs_pct + basal_area + off_ground + (1 | Area)Priors
priors <- list(
# Community-level effects
beta.comm.normal = list(mean = 0, var = 2.73),
alpha.comm.normal = list(mean = 0, var = 2.73),
# Species-level variance (hierarchical shrinkage)
tau.sq.beta.ig = list(2, 1),
tau.sq.alpha.ig = list(2, 1),
# Random effect variances
sigma.sq.psi.ig = list(2, 1), # occupancy REs (Area)
sigma.sq.p.ig = list(2, 1) # detection REs (cam, site_id)
)Run the occupancy model
H0 – Null. Species use CWDs independently of any collected variable, and the detection process acts independently of collected variables.
H1 – Decay. Decay is the driver. As a CWD softens it gains cavities and invertebrate prey, so more species use it as decay advances.
H2 – Intermediate decay. Use peaks at intermediate decay: sound CWDs offer no interior refuge or food, while fully collapsed CWDs have lost the structure that makes a perch or a den site.
H3 – Stand openness. The CWD is incidental; what matters is the stand it sits in, along the open-pine to closed canopy gradient that canopy cover indexes.
H4 – CWD + stand. Decay and stand openness act additively: the log provides structure, the stand determines which species are present to use it.
H5 – Concealment. Understory cover governs use: species travel and forage on CWDs where shrubs and ground vegetation provide escape cover from predators.
H6 – Forest structure. Mature forest structure drives the community: dense, high-basal-area stands support the arboreal and closed-canopy species regardless of the log itself.
H7 – Decay as use intensity. Decay does not change which species occupy a CWD, it changes how often they use it, so decay belongs on detection rather than occupancy.
H8 – Global. CWD condition, stand openness, moisture and understory all contribute.
Model Comparison
waicOcc(ms_null_null, by.sp = FALSE)## elpd pD WAIC
## -2313.3233 248.9036 5124.4538waicOcc(ms_week_decay, by.sp = FALSE)## elpd pD WAIC
## -2204.6446 309.5596 5028.4083waicOcc(ms_ground_decayQ, by.sp = FALSE)## elpd pD WAIC
## -2241.1014 288.3186 5058.8401waicOcc(ms_wkground_canopy, by.sp = FALSE)## elpd pD WAIC
## -2232.5065 290.6617 5046.3364waicOcc(ms_weekonly_logstand, by.sp = FALSE)## elpd pD WAIC
## -2222.0075 304.3283 5052.6716waicOcc(ms_wkground_conceal, by.sp = FALSE)## elpd pD WAIC
## -2212.8065 305.0465 5035.7060waicOcc(ms_full_structure, by.sp = FALSE)## elpd pD WAIC
## -2182.0212 318.1414 5000.3252waicOcc(ms_full_canopy, by.sp = FALSE)## elpd pD WAIC
## -2209.3063 306.8852 5032.3831waicOcc(ms_full_global, by.sp = FALSE)## elpd pD WAIC
## -2150.889 333.352 4968.482WAIC table
model_objects <- c(
"H0: Null" = "ms_null_null",
"H1: Decay" = "ms_week_decay",
"H2: Intermediate decay" = "ms_ground_decayQ",
"H3: Stand openness" = "ms_wkground_canopy",
"H4: Log + stand" = "ms_weekonly_logstand",
"H5: Concealment" = "ms_wkground_conceal",
"H6: Forest structure" = "ms_full_structure",
"H7: Decay as use intensity" = "ms_full_canopy",
"H8: Global" = "ms_full_global"
)
have <- sapply(model_objects, exists)
mod_list <- mget(model_objects[have])
names(mod_list) <- names(model_objects)[have]
cat("Models available for comparison:", length(mod_list), "of",
length(model_objects), "\n")## Models available for comparison: 9 of 9waic_tab <- bind_rows(lapply(mod_list, function(m)
as.data.frame(t(waicOcc(m, by.sp = FALSE)))), .id = "Model") %>%
mutate(delta = WAIC - min(WAIC),
weight = exp(-0.5 * delta) / sum(exp(-0.5 * delta))) %>%
arrange(WAIC) %>%
mutate(across(c(elpd, pD, WAIC, delta), ~ round(.x, 1)),
weight = round(weight, 3))
print(as.data.frame(waic_tab), row.names = FALSE)## Model elpd pD WAIC delta weight
## H8: Global -2150.9 333.4 4968.5 0.0 1
## H6: Forest structure -2182.0 318.1 5000.3 31.8 0
## H1: Decay -2204.6 309.6 5028.4 59.9 0
## H7: Decay as use intensity -2209.3 306.9 5032.4 63.9 0
## H5: Concealment -2212.8 305.0 5035.7 67.2 0
## H3: Stand openness -2232.5 290.7 5046.3 77.9 0
## H4: Log + stand -2222.0 304.3 5052.7 84.2 0
## H2: Intermediate decay -2241.1 288.3 5058.8 90.4 0
## H0: Null -2313.3 248.9 5124.5 156.0 0write.csv(waic_tab, file.path(out_dir, "WAIC_hypotheses.csv"), row.names = FALSE)Per-species WAIC
waic_sp <- bind_rows(lapply(names(mod_list), function(nm) {
w <- waicOcc(mod_list[[nm]], by.sp = TRUE)
tibble(Model = nm, species = dimnames(data_list$y)[[1]], WAIC = w[, "WAIC"])
})) %>%
group_by(species) %>%
mutate(delta = WAIC - min(WAIC)) %>%
ungroup()
best_by_sp <- waic_sp %>%
group_by(species) %>%
slice_min(WAIC, n = 1) %>%
ungroup() %>%
arrange(species)
print(as.data.frame(best_by_sp %>% dplyr::select(species, Model)), row.names = FALSE)## species Model
## bobcat H8: Global
## broad_headed_skink H8: Global
## brown_thrasher H8: Global
## carolina_wren H8: Global
## eastern_chipmunk H6: Forest structure
## eastern_cottontail H2: Intermediate decay
## eastern_gray_squirrel H8: Global
## eastern_towhee H8: Global
## eastern_woodrat H8: Global
## five_lined_skink H8: Global
## gray_catbird H3: Stand openness
## hispid_cotton_rat H5: Concealment
## indigo_bunting H2: Intermediate decay
## nine_banded_armadillo H2: Intermediate decay
## northern_house_wren H3: Stand openness
## raccoon H8: Global
## six_lined_racerunner H3: Stand openness
## southern_black_racer H8: Global
## southern_flying_squirrel H4: Log + stand
## virginia_opossum H8: Global
## white_tailed_deer H2: Intermediate decaywrite.csv(waic_sp, file.path(out_dir, "WAIC_by_species.csv"), row.names = FALSE)Model Validation
A Bayesian p-value near 0.5 indicates adequate fit; values near 0 or 1 indicate the replicated data do not look like the observed data.
ppc_ft <- ppcOcc(
ms_full_global,
fit.stat = "freeman-tukey",
group = 1
)## Currently on species 1 out of 21## Currently on species 2 out of 21## Currently on species 3 out of 21## Currently on species 4 out of 21## Currently on species 5 out of 21## Currently on species 6 out of 21## Currently on species 7 out of 21## Currently on species 8 out of 21## Currently on species 9 out of 21## Currently on species 10 out of 21## Currently on species 11 out of 21## Currently on species 12 out of 21## Currently on species 13 out of 21## Currently on species 14 out of 21## Currently on species 15 out of 21## Currently on species 16 out of 21## Currently on species 17 out of 21## Currently on species 18 out of 21## Currently on species 19 out of 21## Currently on species 20 out of 21## Currently on species 21 out of 21ppc_chisq <- ppcOcc(
ms_full_global,
fit.stat = "chi-square",
group = 1
)## Currently on species 1 out of 21## Currently on species 2 out of 21## Currently on species 3 out of 21## Currently on species 4 out of 21## Currently on species 5 out of 21## Currently on species 6 out of 21## Currently on species 7 out of 21## Currently on species 8 out of 21## Currently on species 9 out of 21## Currently on species 10 out of 21## Currently on species 11 out of 21## Currently on species 12 out of 21## Currently on species 13 out of 21## Currently on species 14 out of 21## Currently on species 15 out of 21## Currently on species 16 out of 21## Currently on species 17 out of 21## Currently on species 18 out of 21## Currently on species 19 out of 21## Currently on species 20 out of 21## Currently on species 21 out of 21summary(ppc_ft)##
## Call:
## ppcOcc(object = ms_full_global, fit.stat = "freeman-tukey", group = 1)
##
## Samples per Chain: 10000
## Burn-in: 3000
## Thinning Rate: 2
## Number of Chains: 4
## Total Posterior Samples: 14000
##
## ----------------------------------------
## Community Level
## ----------------------------------------
## Bayesian p-value: 0.5116
##
## ----------------------------------------
## Species Level
## ----------------------------------------
## nine_banded_armadillo Bayesian p-value: 0.7621
## carolina_wren Bayesian p-value: 0.5974
## eastern_woodrat Bayesian p-value: 0.3321
## virginia_opossum Bayesian p-value: 0.7029
## hispid_cotton_rat Bayesian p-value: 0.5326
## raccoon Bayesian p-value: 0.4631
## eastern_chipmunk Bayesian p-value: 0.4331
## eastern_gray_squirrel Bayesian p-value: 0.6054
## white_tailed_deer Bayesian p-value: 0.6575
## eastern_cottontail Bayesian p-value: 0.5937
## bobcat Bayesian p-value: 0.4749
## eastern_towhee Bayesian p-value: 0.3616
## broad_headed_skink Bayesian p-value: 0.5882
## five_lined_skink Bayesian p-value: 0.4273
## indigo_bunting Bayesian p-value: 0.4039
## northern_house_wren Bayesian p-value: 0.4873
## southern_black_racer Bayesian p-value: 0.5404
## brown_thrasher Bayesian p-value: 0.5156
## southern_flying_squirrel Bayesian p-value: 0.4357
## gray_catbird Bayesian p-value: 0.4624
## six_lined_racerunner Bayesian p-value: 0.3666
## Fit statistic: freeman-tukeysummary(ppc_chisq)##
## Call:
## ppcOcc(object = ms_full_global, fit.stat = "chi-square", group = 1)
##
## Samples per Chain: 10000
## Burn-in: 3000
## Thinning Rate: 2
## Number of Chains: 4
## Total Posterior Samples: 14000
##
## ----------------------------------------
## Community Level
## ----------------------------------------
## Bayesian p-value: 0.4893
##
## ----------------------------------------
## Species Level
## ----------------------------------------
## nine_banded_armadillo Bayesian p-value: 0.4949
## carolina_wren Bayesian p-value: 0.4158
## eastern_woodrat Bayesian p-value: 0.5034
## virginia_opossum Bayesian p-value: 0.4899
## hispid_cotton_rat Bayesian p-value: 0.4622
## raccoon Bayesian p-value: 0.4265
## eastern_chipmunk Bayesian p-value: 0.4603
## eastern_gray_squirrel Bayesian p-value: 0.4629
## white_tailed_deer Bayesian p-value: 0.4256
## eastern_cottontail Bayesian p-value: 0.453
## bobcat Bayesian p-value: 0.4535
## eastern_towhee Bayesian p-value: 0.4411
## broad_headed_skink Bayesian p-value: 0.5269
## five_lined_skink Bayesian p-value: 0.6089
## indigo_bunting Bayesian p-value: 0.561
## northern_house_wren Bayesian p-value: 0.4984
## southern_black_racer Bayesian p-value: 0.5151
## brown_thrasher Bayesian p-value: 0.4985
## southern_flying_squirrel Bayesian p-value: 0.4879
## gray_catbird Bayesian p-value: 0.4886
## six_lined_racerunner Bayesian p-value: 0.6003
## Fit statistic: chi-squareRandom effect on detection
fit_compare <- bind_rows(
as.data.frame(t(waicOcc(ms_full_global_nosite, by.sp = FALSE))) %>%
mutate(Model = "H8: Global (no site RE)"),
as.data.frame(t(waicOcc(ms_full_global, by.sp = FALSE))) %>%
mutate(Model = "H8: Global + (1|site_id) on p")
) %>%
mutate(delta = WAIC - min(WAIC)) %>%
dplyr::select(Model, elpd, pD, WAIC, delta)
print(fit_compare)## Model elpd pD WAIC delta
## 1 H8: Global (no site RE) -2515.502 223.7604 5478.526 510.0435
## 2 H8: Global + (1|site_id) on p -2150.889 333.3520 4968.482 0.0000ppc_nosite <- ppcOcc(ms_full_global_nosite, fit.stat = "freeman-tukey", group = 1)## Currently on species 1 out of 21## Currently on species 2 out of 21## Currently on species 3 out of 21## Currently on species 4 out of 21## Currently on species 5 out of 21## Currently on species 6 out of 21## Currently on species 7 out of 21## Currently on species 8 out of 21## Currently on species 9 out of 21## Currently on species 10 out of 21## Currently on species 11 out of 21## Currently on species 12 out of 21## Currently on species 13 out of 21## Currently on species 14 out of 21## Currently on species 15 out of 21## Currently on species 16 out of 21## Currently on species 17 out of 21## Currently on species 18 out of 21## Currently on species 19 out of 21## Currently on species 20 out of 21## Currently on species 21 out of 21bp_compare <- tibble(
species = dimnames(data_list$y)[[1]],
bp_before = sapply(seq_len(N), function(i)
mean(ppc_nosite$fit.y.rep[, i] > ppc_nosite$fit.y[, i])),
bp_after = sapply(seq_len(N), function(i)
mean(ppc_ft$fit.y.rep[, i] > ppc_ft$fit.y[, i]))
) %>%
mutate(across(where(is.numeric), ~ round(.x, 3))) %>%
arrange(bp_before)
print(as.data.frame(bp_compare), row.names = FALSE)## species bp_before bp_after
## eastern_woodrat 0.000 0.332
## eastern_chipmunk 0.020 0.433
## nine_banded_armadillo 0.059 0.762
## five_lined_skink 0.060 0.427
## virginia_opossum 0.067 0.703
## bobcat 0.138 0.475
## eastern_gray_squirrel 0.187 0.605
## carolina_wren 0.230 0.597
## six_lined_racerunner 0.260 0.367
## southern_flying_squirrel 0.295 0.436
## hispid_cotton_rat 0.306 0.533
## northern_house_wren 0.310 0.487
## eastern_cottontail 0.313 0.594
## raccoon 0.332 0.463
## indigo_bunting 0.338 0.404
## broad_headed_skink 0.401 0.588
## eastern_towhee 0.423 0.362
## southern_black_racer 0.455 0.540
## white_tailed_deer 0.474 0.657
## gray_catbird 0.504 0.462
## brown_thrasher 0.571 0.516write.csv(bp_compare, file.path(out_dir, "bayesian_p_before_after.csv"),
row.names = FALSE)Occupancy Identifiability
sd_site <- sqrt(mean(ms_full_global$sigma.sq.p.samples[, "site_id"]))
gh <- gauss.quad.prob(60, dist = "normal", mu = 0, sigma = sd_site)
sp.names <- dimnames(data_list$y)[[1]]
beta0 <- colMeans(ms_full_global$beta.samples)[paste0("(Intercept)-", sp.names)]
alpha0 <- colMeans(ms_full_global$alpha.samples)[paste0("(Intercept)-", sp.names)]
naive <- apply(data_list$y, 1, function(x) mean(rowSums(x, na.rm = TRUE) > 0))
ident <- tibble(
species = sp.names,
psi = plogis(beta0),
p_week = plogis(alpha0),
p_det = sapply(alpha0, function(a)
sum(gh$weights * (1 - (1 - plogis(a + gh$nodes))^K))),
naive_occ = naive
) %>%
mutate(expected_naive = psi * p_det,
gap = expected_naive - naive_occ,
unobserved = psi - naive_occ) %>%
arrange(desc(gap))
print(as.data.frame(ident %>% mutate(across(where(is.numeric), ~ round(.x, 3)))),
row.names = FALSE)## species psi p_week p_det naive_occ expected_naive gap
## five_lined_skink 0.943 0.041 0.469 0.20 0.442 0.242
## indigo_bunting 0.970 0.036 0.440 0.20 0.426 0.226
## eastern_towhee 0.946 0.052 0.523 0.30 0.495 0.195
## southern_flying_squirrel 0.940 0.024 0.354 0.14 0.333 0.193
## six_lined_racerunner 0.924 0.020 0.313 0.10 0.290 0.190
## broad_headed_skink 0.968 0.040 0.464 0.26 0.449 0.189
## eastern_gray_squirrel 0.918 0.113 0.706 0.46 0.648 0.188
## eastern_chipmunk 0.928 0.114 0.708 0.48 0.657 0.177
## southern_black_racer 0.968 0.025 0.364 0.18 0.352 0.172
## brown_thrasher 0.907 0.023 0.344 0.14 0.312 0.172
## raccoon 0.964 0.102 0.683 0.50 0.659 0.159
## northern_house_wren 0.957 0.024 0.353 0.18 0.337 0.157
## bobcat 0.975 0.054 0.534 0.38 0.521 0.141
## gray_catbird 0.931 0.016 0.279 0.12 0.260 0.140
## eastern_cottontail 0.969 0.060 0.557 0.44 0.540 0.100
## eastern_woodrat 0.970 0.120 0.719 0.60 0.698 0.098
## hispid_cotton_rat 0.973 0.072 0.599 0.54 0.583 0.043
## virginia_opossum 0.965 0.090 0.652 0.60 0.629 0.029
## white_tailed_deer 0.984 0.046 0.496 0.46 0.488 0.028
## carolina_wren 0.985 0.095 0.666 0.66 0.656 -0.004
## nine_banded_armadillo 0.988 0.117 0.713 0.80 0.704 -0.096
## unobserved
## 0.743
## 0.770
## 0.646
## 0.800
## 0.824
## 0.708
## 0.458
## 0.448
## 0.788
## 0.767
## 0.464
## 0.777
## 0.595
## 0.811
## 0.529
## 0.370
## 0.433
## 0.365
## 0.524
## 0.325
## 0.188write.csv(ident, file.path(out_dir, "psi_identifiability.csv"), row.names = FALSE)
cat("\nmean psi across species:", round(mean(ident$psi), 3), "\n")##
## mean psi across species: 0.956cat("mean gap (expected - observed naive occupancy):", round(mean(ident$gap), 3), "\n")## mean gap (expected - observed naive occupancy): 0.13cat("species over-predicted by > 0.05:", sum(ident$gap > 0.05), "of", N, "\n")## species over-predicted by > 0.05: 16 of 21if (mean(ident$psi) > 0.9 || sum(ident$gap > 0.05) > N / 2) {
cat("\n*** WARNING: occupancy is saturating and the model over-predicts what",
"\n was seen. The detection random effect is likely absorbing the",
"\n occupancy signal. Do not report psi or the occupancy coefficients",
"\n from this fit without addressing it.\n")
}##
## *** WARNING: occupancy is saturating and the model over-predicts what
## was seen. The detection random effect is likely absorbing the
## occupancy signal. Do not report psi or the occupancy coefficients
## from this fit without addressing it.Convergence
Rhat should be below 1.01 and effective sample sizes should be 400+.
cat("Max Rhat (occupancy betas): ", round(max(ms_full_global$rhat$beta), 4), "\n")## Max Rhat (occupancy betas): 1.0606cat("Max Rhat (detection alphas):", round(max(ms_full_global$rhat$alpha), 4), "\n")## Max Rhat (detection alphas): 1.0493cat("Min ESS (occupancy betas): ", round(min(ms_full_global$ESS$beta)), "\n")## Min ESS (occupancy betas): 222cat("Min ESS (detection alphas): ", round(min(ms_full_global$ESS$alpha)), "\n")## Min ESS (detection alphas): 582Probability of Direction - Occupancy
pd <- function(x) {
max(mean(x > 0), mean(x < 0))
}
get_quantile <- function(samples, prob) {
apply(samples, 2, quantile, probs = prob)
}
beta_samps <- ms_full_global$beta.samples
results_species <- data.frame(
parameter = colnames(beta_samps),
mean = colMeans(beta_samps),
pd = apply(beta_samps, 2, pd),
lower_95_cri = get_quantile(beta_samps, 0.025),
upper_95_cri = get_quantile(beta_samps, 0.975)
)
results_species$pd_95 <- results_species$pd >= 0.95
head(results_species)## parameter mean
## (Intercept)-nine_banded_armadillo (Intercept)-nine_banded_armadillo 4.443899
## (Intercept)-carolina_wren (Intercept)-carolina_wren 4.154803
## (Intercept)-eastern_woodrat (Intercept)-eastern_woodrat 3.477198
## (Intercept)-virginia_opossum (Intercept)-virginia_opossum 3.309150
## (Intercept)-hispid_cotton_rat (Intercept)-hispid_cotton_rat 3.575092
## (Intercept)-raccoon (Intercept)-raccoon 3.296434
## pd lower_95_cri upper_95_cri pd_95
## (Intercept)-nine_banded_armadillo 1.0000000 2.212897 7.731388 TRUE
## (Intercept)-carolina_wren 0.9999286 2.017093 7.245362 TRUE
## (Intercept)-eastern_woodrat 0.9996429 1.424397 6.157945 TRUE
## (Intercept)-virginia_opossum 0.9992857 1.337921 5.577245 TRUE
## (Intercept)-hispid_cotton_rat 0.9995000 1.335175 6.435944 TRUE
## (Intercept)-raccoon 0.9987143 1.233616 5.725193 TRUEbeta_comm <- ms_full_global$beta.comm.samples
results_comm <- data.frame(
parameter = colnames(beta_comm),
mean = colMeans(beta_comm),
pd = apply(beta_comm, 2, pd),
lower_95_cri = get_quantile(beta_comm, 0.025),
upper_95_cri = get_quantile(beta_comm, 0.975)
)
results_comm$pd_95 <- results_comm$pd >= 0.95
results_comm## parameter mean pd lower_95_cri upper_95_cri pd_95
## (Intercept) (Intercept) 3.1358723 1.0000000 1.6609935 4.6577846 TRUE
## decay decay -0.1826473 0.6435714 -1.1699618 0.9153587 FALSE
## canopy_cover canopy_cover 1.9272648 0.9990714 0.6690618 3.2691442 TRUE
## soil_moist soil_moist -0.8772601 0.9177143 -2.1780060 0.3074653 FALSE
## gc_shannon gc_shannon 0.7538890 0.9271429 -0.3143801 1.8356924 FALSE
## shrubs_pct shrubs_pct 0.2351838 0.6778571 -0.7900393 1.2492708 FALSE
## basal_area basal_area -1.1497069 0.9914286 -2.2606721 -0.1744629 TRUE
## off_ground off_ground 0.1486143 0.6044286 -0.7268230 1.1899594 FALSEdoc <- read_docx()
ft_species <- flextable(results_species) %>%
colformat_double(digits = 4) %>%
set_header_labels(
parameter = "Parameter",
mean = "Mean",
pd = "pd",
lower_95_cri = "Lower 95% CrI",
upper_95_cri = "Upper 95% CrI"
) %>%
theme_booktabs() %>%
autofit() %>%
bold(i = ~ pd >= 0.95, bold = TRUE, part = "body")
doc <- doc %>%
body_add_par("Species-Level Results", style = "heading 1") %>%
body_add_par("", style = "Normal") %>%
body_add_flextable(ft_species)
ft_comm <- flextable(results_comm) %>%
colformat_double(digits = 4) %>%
set_header_labels(
parameter = "Parameter",
mean = "Mean",
pd = "pd",
lower_95_cri = "Lower 95% CrI",
upper_95_cri = "Upper 95% CrI"
) %>%
theme_booktabs() %>%
autofit() %>%
bold(i = ~ pd >= 0.95, bold = TRUE, part = "body")
doc <- doc %>%
body_add_break() %>%
body_add_par("Community-Level Results", style = "heading 1") %>%
body_add_par("", style = "Normal") %>%
body_add_flextable(ft_comm)
print(doc, target = file.path(out_dir, "occ_model_results.docx"))
write.csv(results_species, file.path(out_dir, "species_occupancy_params.csv"),
row.names = FALSE)
write.csv(results_comm, file.path(out_dir, "community_occupancy_params.csv"),
row.names = FALSE)H2: do rodents and birds respond differently?
h2_model <- ms_full_global
guild_lookup_occ <- tibble::tribble(
~Name, ~Guild,
"eastern_woodrat", "Rodent",
"hispid_cotton_rat", "Rodent",
"eastern_chipmunk", "Rodent",
"eastern_gray_squirrel", "Rodent",
"southern_flying_squirrel", "Rodent",
"carolina_wren", "Bird",
"eastern_towhee", "Bird",
"indigo_bunting", "Bird",
"northern_house_wren", "Bird",
"brown_thrasher", "Bird",
"gray_catbird", "Bird",
"nine_banded_armadillo", "Ground omnivore",
"virginia_opossum", "Ground omnivore",
"raccoon", "Ground omnivore",
"white_tailed_deer", "Herbivore",
"eastern_cottontail", "Herbivore",
"bobcat", "Carnivore",
"broad_headed_skink", "Reptile",
"five_lined_skink", "Reptile",
"southern_black_racer", "Reptile",
"six_lined_racerunner", "Reptile"
)
beta_h2 <- h2_model$beta.samples
# spOccupancy names these "<term>-<species>"
sp_names <- sub("^[^-]*-", "", colnames(beta_h2))
term_name <- sub("-.*$", "", colnames(beta_h2))
missing_sp <- setdiff(unique(sp_names), guild_lookup_occ$Name)
if (length(missing_sp)) {
stop("Species in the model with no guild assigned: ",
paste(missing_sp, collapse = ", "))
}
guild_of <- setNames(guild_lookup_occ$Guild, guild_lookup_occ$Name)[sp_names]
cat("Species per guild in the occupancy model:\n")## Species per guild in the occupancy model:print(table(guild_of[!duplicated(sp_names)]))##
## Bird Carnivore Ground omnivore Herbivore Reptile
## 6 1 3 2 4
## Rodent
## 5# Posterior of the guild mean coefficient, per term x guild
guild_posterior <- function(term, guild) {
cols <- which(term_name == term & guild_of == guild)
if (!length(cols)) return(NULL)
rowMeans(beta_h2[, cols, drop = FALSE])
}
rest_posterior <- function(term, guild) {
cols <- which(term_name == term & guild_of != guild)
if (!length(cols)) return(NULL)
rowMeans(beta_h2[, cols, drop = FALSE])
}
pd_fn2 <- function(x) max(mean(x > 0), mean(x < 0))
terms_available <- unique(term_name)
h2_terms <- intersect(c("shrubs_pct", "sdw_pct"), terms_available)
h2_guild <- bind_rows(lapply(h2_terms, function(tm) {
bind_rows(lapply(sort(unique(guild_of)), function(g) {
v <- guild_posterior(tm, g)
if (is.null(v)) return(NULL)
tibble(Term = tm, Guild = g,
n_species = sum(term_name == tm & guild_of == g),
mean = mean(v), pd = pd_fn2(v),
lower = quantile(v, 0.025), upper = quantile(v, 0.975))
}))
}))
cat("\nGuild-level occupancy coefficients (H2 terms):\n")##
## Guild-level occupancy coefficients (H2 terms):print(as.data.frame(h2_guild %>%
mutate(across(where(is.numeric), ~ round(.x, 3)))), row.names = FALSE)## Term Guild n_species mean pd lower upper
## shrubs_pct Bird 6 0.217 0.655 -0.966 1.339
## shrubs_pct Carnivore 1 0.471 0.685 -1.536 2.586
## shrubs_pct Ground omnivore 3 0.018 0.512 -1.205 1.212
## shrubs_pct Herbivore 2 0.078 0.557 -1.453 1.431
## shrubs_pct Reptile 4 0.374 0.725 -0.866 1.709
## shrubs_pct Rodent 5 0.303 0.720 -0.723 1.375# Contrast: this guild versus every other species in the model
h2_contrast <- bind_rows(lapply(h2_terms, function(tm) {
bind_rows(lapply(sort(unique(guild_of)), function(g) {
a <- guild_posterior(tm, g); b <- rest_posterior(tm, g)
if (is.null(a) || is.null(b)) return(NULL)
d <- a - b
tibble(Term = tm, Guild = g, contrast = "guild - all others",
mean = mean(d), pd = pd_fn2(d),
lower = quantile(d, 0.025), upper = quantile(d, 0.975))
}))
}))
cat("\nGuild versus all other species:\n")##
## Guild versus all other species:print(as.data.frame(h2_contrast %>%
mutate(across(where(is.numeric), ~ round(.x, 3)))), row.names = FALSE)## Term Guild contrast mean pd lower upper
## shrubs_pct Bird guild - all others -0.030 0.522 -0.828 0.699
## shrubs_pct Carnivore guild - all others 0.245 0.610 -1.461 2.094
## shrubs_pct Ground omnivore guild - all others -0.256 0.710 -1.282 0.584
## shrubs_pct Herbivore guild - all others -0.176 0.610 -1.541 0.900
## shrubs_pct Reptile guild - all others 0.168 0.633 -0.671 1.218
## shrubs_pct Rodent guild - all others 0.086 0.581 -0.639 0.926# The two headline predictions, stated as pass/fail
verdict <- function(term, guild, label) {
if (!term %in% terms_available) {
return(tibble(Prediction = label, Result = "not testable in this model",
mean = NA_real_, pd = NA_real_))
}
v <- guild_posterior(term, guild)
tibble(Prediction = label,
Result = ifelse(pd_fn2(v) >= 0.95 & mean(v) > 0, "supported",
ifelse(pd_fn2(v) >= 0.95 & mean(v) < 0, "opposite direction",
"not supported")),
mean = round(mean(v), 3), pd = round(pd_fn2(v), 3))
}
h2_verdict <- bind_rows(
verdict("shrubs_pct", "Rodent", "H2a: rodent occupancy rises with shrub cover"),
verdict("sdw_pct", "Rodent", "H2a (alt): rodent occupancy rises with standing dead woody cover"),
verdict("off_ground", "Bird", "H2b: bird occupancy rises with log height off ground"))
cat("\nH2 verdicts:\n"); print(as.data.frame(h2_verdict), row.names = FALSE)##
## H2 verdicts:## Prediction
## H2a: rodent occupancy rises with shrub cover
## H2a (alt): rodent occupancy rises with standing dead woody cover
## H2b: bird occupancy rises with log height off ground
## Result mean pd
## not supported 0.303 0.720
## not testable in this model NA NA
## not supported 0.179 0.621cat("\nReminder: shrub and standing-dead cover are proxies for vegetation",
"\nHEIGHT, which was not measured. A null here is a null against the proxy.\n")##
## Reminder: shrub and standing-dead cover are proxies for vegetation
## HEIGHT, which was not measured. A null here is a null against the proxy.write.csv(h2_guild, file.path(out_dir, "H2_guild_occupancy_effects.csv"), row.names = FALSE)
write.csv(h2_contrast, file.path(out_dir, "H2_guild_contrasts.csv"), row.names = FALSE)
write.csv(h2_verdict, file.path(out_dir, "H2_verdicts.csv"), row.names = FALSE)h2_plot <- h2_guild %>%
mutate(Term = recode(Term,
off_ground = "Log off ground",
shrubs_pct = "Shrub cover",
sdw_pct = "Standing dead woody cover"),
credible = ifelse(pd >= 0.95, "yes", "no"),
Guild = forcats::fct_reorder(Guild, mean))
p_h2 <- ggplot(h2_plot, aes(mean, Guild, colour = credible)) +
geom_vline(xintercept = 0, linetype = 2, colour = "grey40") +
geom_pointrange(aes(xmin = lower, xmax = upper), linewidth = 0.7, size = 0.4) +
facet_wrap(~ Term, scales = "free_x") +
scale_colour_manual(values = c(no = "grey65", yes = "#D55E00"),
name = "pd >= 0.95") +
labs(x = "Occupancy coefficient (logit scale, per SD)", y = NULL) +
theme_classic(base_size = 11) +
theme(strip.text = element_text(face = "bold"),
plot.caption = element_text(hjust = 0, colour = "grey35"))
print(p_h2)
ggsave(file.path(out_dir, "Figure_H2_guild_effects.png"), p_h2,
width = 9, height = 4.5, dpi = 300, bg = "white")Probability of Direction - Detection
alpha_samps <- ms_full_global$alpha.samples
results_species_det <- data.frame(
parameter = colnames(alpha_samps),
mean = colMeans(alpha_samps),
pd = apply(alpha_samps, 2, pd),
lower_95_cri = get_quantile(alpha_samps, 0.025),
upper_95_cri = get_quantile(alpha_samps, 0.975)
)
results_species_det$pd_95 <- results_species_det$pd >= 0.95
head(results_species_det)## parameter mean
## (Intercept)-nine_banded_armadillo (Intercept)-nine_banded_armadillo -2.023847
## (Intercept)-carolina_wren (Intercept)-carolina_wren -2.252078
## (Intercept)-eastern_woodrat (Intercept)-eastern_woodrat -1.990582
## (Intercept)-virginia_opossum (Intercept)-virginia_opossum -2.317678
## (Intercept)-hispid_cotton_rat (Intercept)-hispid_cotton_rat -2.562905
## (Intercept)-raccoon (Intercept)-raccoon -2.169783
## pd lower_95_cri upper_95_cri pd_95
## (Intercept)-nine_banded_armadillo 1 -3.022815 -1.0357981 TRUE
## (Intercept)-carolina_wren 1 -3.316993 -1.1964551 TRUE
## (Intercept)-eastern_woodrat 1 -3.044049 -0.9707985 TRUE
## (Intercept)-virginia_opossum 1 -3.311330 -1.3094656 TRUE
## (Intercept)-hispid_cotton_rat 1 -3.602157 -1.5305513 TRUE
## (Intercept)-raccoon 1 -3.226496 -1.1269673 TRUEalpha_comm <- ms_full_global$alpha.comm.samples
results_comm_det <- data.frame(
parameter = colnames(alpha_comm),
mean = colMeans(alpha_comm),
pd = apply(alpha_comm, 2, pd),
lower_95_cri = get_quantile(alpha_comm, 0.025),
upper_95_cri = get_quantile(alpha_comm, 0.975)
)
results_comm_det$pd_95 <- results_comm_det$pd >= 0.95
results_comm_det## parameter mean pd lower_95_cri upper_95_cri pd_95
## (Intercept) (Intercept) -2.88385569 1.0000000 -3.4120719 -2.3697790 TRUE
## week week -0.09420850 0.7800000 -0.3390157 0.1446147 FALSE
## week_sq week_sq -0.28073492 0.9898571 -0.5283216 -0.0448024 TRUE
## off_ground off_ground -0.02567156 0.5792857 -0.2891873 0.2420160 FALSEwrite.csv(results_comm_det, file.path(out_dir, "community_detection_params.csv"),
row.names = FALSE)Figures: occupancy across the decay gradient
Prediction helper
# Raw decay gradient, converted to the scaled units the model was fit on
decay.pred.vals <- seq(1, 5, length.out = 100)
decay.mean <- occ_centers[["decay"]]
decay.sd <- occ_scales[["decay"]]
decay.pred.scale <- (decay.pred.vals - decay.mean) / decay.sd
n_pred <- length(decay.pred.vals)
# Predicted occupancy across the decay gradient from one set of posterior draws.
# suffix = "" for the community, or "-<species>" for a single species.
predict_psi_decay <- function(samples, suffix = "", canopy_scaled = 0) {
b_int <- samples[, paste0("(Intercept)", suffix)]
b_decay <- samples[, paste0("decay", suffix)]
b_canopy <- samples[, paste0("canopy_cover", suffix)]
n_samples <- length(b_int)
logit_psi <- matrix(NA, nrow = n_samples, ncol = n_pred)
for (i in 1:n_samples) {
logit_psi[i, ] <- b_int[i] +
b_decay[i] * decay.pred.scale +
b_canopy[i] * canopy_scaled
# soil_moist, gc_shannon, shrubs_pct and basal_area are all held at their
# means, which are 0 on the scaled axis, so they drop out of the sum
}
psi_samples <- plogis(logit_psi)
psi.quants <- apply(psi_samples, 2, quantile, probs = c(0.025, 0.5, 0.975))
data.frame(
decay = decay.pred.vals,
psi.low = psi.quants[1, ],
psi.med = psi.quants[2, ],
psi.high = psi.quants[3, ]
)
}
# Observed reference: for each log, the proportion of the modelled species
# actually detected there, plotted against that log's decay class
observed_by_log <- tibble(
Plot = dimnames(ms_full_global$y)[[2]],
obs = apply(ms_full_global$y, 2, function(m) mean(rowSums(m, na.rm = TRUE) > 0))
) %>%
left_join(site_env %>% dplyr::select(Plot, decay), by = "Plot")Community occupancy across decay class
psi.comm.dat <- predict_psi_decay(ms_full_global$beta.comm.samples)
community_plot <- ggplot(psi.comm.dat, aes(x = decay, y = psi.med)) +
geom_ribbon(aes(ymin = psi.low, ymax = psi.high),
fill = "grey70", alpha = 0.5) +
geom_line(linewidth = 1.3, color = "#0072B2") +
geom_point(data = observed_by_log, inherit.aes = FALSE,
aes(x = decay, y = obs),
shape = 21, fill = "grey40", color = "black", alpha = 0.5,
size = 2.4, position = position_jitter(width = 0.12, height = 0)) +
scale_y_continuous(limits = c(0, 1)) +
scale_x_continuous(breaks = 1:5) +
labs(x = "CWD Decay Class",
y = "Community Occupancy Probability") +
theme_bw(base_size = 14) +
theme(axis.title = element_text(face = "bold"))
community_plot
ggsave(file.path(out_dir, "Figure_Community_Decay.png"),
community_plot, width = 8, height = 6, dpi = 600, bg = "white")Species Lookup
SpeciesLookup <- tibble::tribble(
~Name, ~Scientific_Name, ~Common_Name,
"nine_banded_armadillo", "Dasypus novemcinctus", "Nine-Banded Armadillo",
"carolina_wren", "Thryothorus ludovicianus", "Carolina Wren",
"eastern_woodrat", "Neotoma floridana", "Eastern Woodrat",
"virginia_opossum", "Didelphis virginiana", "Virginia Opossum",
"hispid_cotton_rat", "Sigmodon hispidus", "Hispid Cotton Rat",
"raccoon", "Procyon lotor", "Raccoon",
"eastern_chipmunk", "Tamias striatus", "Eastern Chipmunk",
"eastern_gray_squirrel", "Sciurus carolinensis", "Eastern Gray Squirrel",
"white_tailed_deer", "Odocoileus virginianus", "White-Tailed Deer",
"eastern_cottontail", "Sylvilagus floridanus", "Eastern Cottontail",
"bobcat", "Lynx rufus", "Bobcat",
"eastern_towhee", "Pipilo erythrophthalmus", "Eastern Towhee",
"broad_headed_skink", "Plestiodon laticeps", "Broad-Headed Skink",
"five_lined_skink", "Plestiodon fasciatus", "Five-Lined Skink",
"indigo_bunting", "Passerina cyanea", "Indigo Bunting",
"northern_house_wren", "Troglodytes aedon", "Northern House Wren",
"southern_black_racer", "Coluber constrictor", "Southern Black Racer",
"brown_thrasher", "Toxostoma rufum", "Brown Thrasher",
"southern_flying_squirrel", "Glaucomys volans", "Southern Flying Squirrel",
"gray_catbird", "Dumetella carolinensis", "Gray Catbird",
"six_lined_racerunner", "Aspidoscelis sexlineata", "Six-Lined Racerunner"
)Species-level occupancy across decay class
sp.names <- ms_full_global$sp.names
species_predictions <- list()
for (sp in sp.names) {
species_predictions[[sp]] <- predict_psi_decay(
ms_full_global$beta.samples, suffix = paste0("-", sp)
) %>%
mutate(species = sp)
}
all_predictions <- do.call(rbind, species_predictions) %>%
left_join(SpeciesLookup %>% dplyr::select(Name, Common_Name),
by = c("species" = "Name")) %>%
mutate(Common_Name = ifelse(is.na(Common_Name),
gsub("_", " ", species), Common_Name))
# Flag species whose decay effect reaches pd >= 0.95, for shading
decay_pd <- results_species %>%
filter(grepl("^decay-", parameter)) %>%
transmute(species = sub("^decay-", "", parameter),
decay_pd = pd,
credible = factor(ifelse(pd >= 0.95, "yes", "no"),
levels = c("no", "yes")))
all_predictions <- all_predictions %>% left_join(decay_pd, by = "species")
faceted_plot <- ggplot(all_predictions,
aes(x = decay, y = psi.med)) +
geom_ribbon(aes(ymin = psi.low, ymax = psi.high, fill = credible),
alpha = 0.4) +
geom_line(aes(color = credible), linewidth = 1) +
facet_wrap(~ Common_Name, ncol = 4) +
scale_color_manual(values = c(no = "grey35", yes = "#D55E00"),
drop = FALSE, name = "decay pd >= 0.95") +
scale_fill_manual(values = c(no = "grey75", yes = "#E9B7A0"),
drop = FALSE, name = "decay pd >= 0.95") +
scale_y_continuous(limits = c(0, 1)) +
scale_x_continuous(breaks = 1:5) +
labs(x = "CWD Decay Class", y = "Occupancy Probability") +
theme_bw(base_size = 11) +
theme(strip.text = element_text(face = "bold", size = 9),
legend.position = "bottom")
faceted_plot
ggsave(file.path(out_dir, "Figure_Species_Decay.png"),
faceted_plot, width = 12, height = 9, dpi = 600, bg = "white")Decay across levels of canopy cover
canopy_scenarios <- data.frame(
scenario = factor(c("Low", "Average", "High"),
levels = c("Low", "Average", "High")),
canopy_scaled = c(-1, 0, 1),
color = c("#D55E00", "#009E73", "#0072B2")
)
community_predictions <- list()
for (j in 1:nrow(canopy_scenarios)) {
community_predictions[[j]] <- predict_psi_decay(
ms_full_global$beta.comm.samples,
canopy_scaled = canopy_scenarios$canopy_scaled[j]
) %>%
mutate(scenario = canopy_scenarios$scenario[j])
}
psi.canopy.dat <- do.call(rbind, community_predictions) %>%
mutate(scenario = factor(scenario, levels = c("Low", "Average", "High")))
canopy_plot <- ggplot(psi.canopy.dat,
aes(x = decay, y = psi.med,
color = scenario, fill = scenario)) +
geom_ribbon(aes(ymin = psi.low, ymax = psi.high), alpha = 0.2, color = NA) +
geom_line(linewidth = 1.3) +
scale_color_manual(values = canopy_scenarios$color, name = "Canopy Cover") +
scale_fill_manual(values = canopy_scenarios$color, name = "Canopy Cover") +
scale_y_continuous(limits = c(0, 1)) +
scale_x_continuous(breaks = 1:5) +
labs(x = "CWD Decay Class", y = "Community Occupancy Probability") +
theme_bw(base_size = 14) +
theme(legend.position = "bottom",
axis.title = element_text(face = "bold"))
canopy_plot
ggsave(file.path(out_dir, "Figure_Decay_by_Canopy.png"),
canopy_plot, width = 8, height = 6, dpi = 600, bg = "white")
write.csv(psi.canopy.dat,
file.path(out_dir, "community_predictions_by_decay.csv"),
row.names = FALSE)Combined figure
combined_plot <- (community_plot / faceted_plot) +
plot_layout(heights = c(1, 2)) +
plot_annotation(tag_levels = "A")
combined_plot
ggsave(file.path(out_dir, "Figure_Decay_Combined.png"),
combined_plot, width = 12, height = 14, dpi = 600, bg = "white")Batched Parameter Estimates
Caterpillar plots of the posterior, in chunks so the labels stay readable.
MCMCplot(ms_full_global$beta.samples, ref_ovl = TRUE, ci = c(50, 95)) # Occupancy
MCMCplot(ms_full_global$alpha.samples, ref_ovl = TRUE, ci = c(50, 95)) # Detection
## Occupancy
n_params <- ncol(ms_full_global$beta.samples)
chunk_size <- 10
for (i in seq(1, n_params, by = chunk_size)) {
end <- min(i + chunk_size - 1, n_params)
param_names <- colnames(ms_full_global$beta.samples)[i:end]
file_name <- file.path(out_dir,
paste0("MCMCplot_Occupancy_Params_", i, "_to_", end, ".png"))
png(filename = file_name, width = 1200, height = 800, res = 150)
MCMCplot(ms_full_global$beta.samples[, param_names, drop = FALSE],
ref_ovl = TRUE,
ci = c(50, 95),
main = paste0("Occupancy Parameters: ", i, " to ", end))
dev.off()
}
## Detection
n_params <- ncol(ms_full_global$alpha.samples)
for (i in seq(1, n_params, by = chunk_size)) {
end <- min(i + chunk_size - 1, n_params)
param_names <- colnames(ms_full_global$alpha.samples)[i:end]
file_name <- file.path(out_dir,
paste0("MCMCplot_Detection_Params_", i, "_to_", end, ".png"))
png(filename = file_name, width = 1200, height = 800, res = 150)
MCMCplot(ms_full_global$alpha.samples[, param_names, drop = FALSE],
ref_ovl = TRUE,
ci = c(50, 95),
main = paste0("Detection Parameters: ", i, " to ", end))
dev.off()
}Facet Plots
plot_mcmc_faceted_groups <- function(samples_matrix,
chunk_size = 8,
title_text = "Parameters",
pd_cutoff = 0.95) {
param_names <- colnames(samples_matrix)
df <- as.data.frame(samples_matrix)
df$Iteration <- 1:nrow(df)
df_long <- df %>%
pivot_longer(-Iteration,
names_to = "Parameter",
values_to = "Value") %>%
mutate(
ParamIndex = match(Parameter, param_names),
Group = ceiling(ParamIndex / chunk_size)
)
# Posterior summaries + probability of direction
summary_df <- df_long %>%
group_by(Parameter, Group) %>%
summarise(
mean = mean(Value),
l95 = quantile(Value, 0.025),
u95 = quantile(Value, 0.975),
l50 = quantile(Value, 0.25),
u50 = quantile(Value, 0.75),
pd = max(mean(Value > 0), mean(Value < 0)),
.groups = "drop"
) %>%
mutate(
Credible = pd >= pd_cutoff
)
summary_df <- summary_df %>%
group_by(Group) %>%
mutate(Parameter = factor(Parameter,
levels = rev(unique(Parameter)))) %>%
ungroup()
ggplot(summary_df,
aes(x = mean, y = Parameter)) +
geom_vline(xintercept = 0, linetype = "dashed") +
# 95% CI
geom_errorbarh(aes(xmin = l95,
xmax = u95,
size = Credible),
height = 0) +
# 50% CI
geom_errorbarh(aes(xmin = l50,
xmax = u50,
size = Credible),
height = 0) +
# Posterior mean
geom_point(aes(fill = Credible),
shape = 21,
size = 3,
stroke = 1) +
facet_wrap(~Group, scales = "free_y", ncol = 2) +
scale_size_manual(values = c("TRUE" = 1.2,
"FALSE" = 0.5),
guide = "none") +
scale_fill_manual(values = c("TRUE" = "black",
"FALSE" = "white"),
guide = "none") +
labs(x = "Parameter Estimate",
y = "",
title = title_text) +
theme_bw() +
theme(
strip.background = element_blank(),
strip.text = element_blank(),
axis.text.y = element_text(size = 10)
)
}
plot_mcmc_faceted_groups(
ms_full_global$beta.samples,
chunk_size = 8,
title_text = "",
pd_cutoff = 0.95
)## Warning: `geom_errobarh()` was deprecated in ggplot2 4.0.0.
## ℹ Please use the `orientation` argument of `geom_errorbar()` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggplot2 package.
## Please report the issue at <https://github.com/tidyverse/ggplot2/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `height` was translated to `width`.
## `height` was translated to `width`.
ggsave(file.path(out_dir, "mcmc_faceted_groups_occupancy.png"),
width = 10, height = 8, dpi = 300)## `height` was translated to `width`.
## `height` was translated to `width`.plot_mcmc_faceted_groups(
ms_full_global$alpha.samples,
chunk_size = 8,
title_text = "",
pd_cutoff = 0.95
)## `height` was translated to `width`.
## `height` was translated to `width`.
ggsave(file.path(out_dir, "mcmc_faceted_groups_detection.png"),
width = 10, height = 8, dpi = 300)## `height` was translated to `width`.
## `height` was translated to `width`.Save
saveRDS(
list(
best_name = "H8: Global",
ms_full_global = ms_full_global,
mod_list = mod_list,
waic_tab = waic_tab,
waic_sp = waic_sp,
fit_compare = fit_compare,
bp_compare = bp_compare,
ident = ident,
data_list = data_list,
site_env = site_env,
occ_centers = occ_centers,
occ_scales = occ_scales,
area_levels = area_levels,
cam_levels = cam_levels,
samp_mat = samp_mat,
SpeciesLookup = SpeciesLookup,
settings = list(occasion_days = occasion_days, n_occ = K,
min_det = min_det, min_plots = min_plots,
max_naive_occ = max_naive_occ,
mammals_only = mammals_only)
),
file.path(out_dir, "msPGOcc_results.rds")
)