Wildlife_Community_Composition_and_Stress_Runs
Alan Ivory
2026-10-06
What this file does
Part 1 Determine how many NMDS dimensions and what rarity level is needed to best fit the model (limit stress).
Part 2 Run the PERMANOVA analyses.
Data
Camera data
CameraRaw <- read_excel(cam_path, 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'CameraData <- CameraRaw %>%
mutate(
Plot = as.character(Plot),
Location = suppressWarnings(as.numeric(as.character(Location))),
N_Ind = suppressWarnings(as.numeric(as.character(Number.Of.Individuals))),
Common_Name = str_to_lower(str_trim(as.character(Common_Name))),
Date = as.Date(Date),
DateTime = as.POSIXct(Date, tz = "America/New_York") +
(as.numeric(Time) %% 86400),
Year = year(Date),
Area = str_extract(Plot, "^[A-Z]+")
)Filters
clock_errors <- CameraData %>%
filter(Year != 2026) %>% ## Only keep observations from 2026 (some camera were set to 2020-2023)
count(Plot, Year, name = "n_rows")
# Effort matrix: 1 = camera deployed that plot-day
effort_long <- read_excel(eff_path, sheet = "Sheet1") %>%
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)## New names:
## • `` -> `...1`CamFilt <- CameraData %>%
filter(Location == 1, Year == 2026,
!is.na(Common_Name), !is.na(DateTime)) %>%
semi_join(effort_long, by = c("Plot", "Date"))
# Only keep observations that occured on CWD
cat("Location codes in the raw review sheet:\n")## Location codes in the raw review sheet: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),
" + name/time present" = nrow(CameraData %>%
filter(Location == 1, Year == 2026, !is.na(Common_Name), !is.na(DateTime))),
" + camera deployed that day" = nrow(CamFilt))
cat("\nAttrition:\n"); print(data.frame(records = step_n))##
## Attrition:## records
## Raw records 5880
## On the CWD (Location == 1) 5127
## + 2026 only 4966
## + name/time present 4966
## + camera deployed that day 4411stopifnot(all(CamFilt$Location == 1), !any(is.na(CamFilt$Location)))
# Minimum 30 minute gap between records of the same species on the same CWD
gap_mins <- 30
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() %>%
mutate(Unique = 1L)
cat("\nIndependent on-CWD detections:", nrow(CamInd), "\n")##
## Independent on-CWD detections: 3530Sampling effort
effort <- effort_long %>%
group_by(Plot) %>%
summarise(first_day = min(Date),
last_day = max(Date),
camera_days = n(),
.groups = "drop") %>%
mutate(gap_days = as.numeric(last_day - first_day) + 1 - camera_days)
print(effort %>% arrange(camera_days), n = 50)## # A tibble: 50 × 5
## Plot first_day last_day camera_days gap_days
## <chr> <date> <date> <int> <dbl>
## 1 CW14 2026-04-22 2026-07-07 64 13
## 2 CW17 2026-04-15 2026-07-07 65 19
## 3 CW05 2026-04-06 2026-07-07 66 27
## 4 CW09 2026-04-06 2026-06-28 69 15
## 5 CE03 2026-04-08 2026-07-07 70 21
## 6 D03 2026-04-06 2026-07-07 70 23
## 7 D07 2026-04-07 2026-07-07 71 21
## 8 CE10 2026-04-08 2026-07-07 74 17
## 9 CW13 2026-04-08 2026-07-07 74 17
## 10 CE12 2026-04-08 2026-07-07 75 16
## 11 CE15 2026-04-14 2026-07-07 75 10
## 12 CW01 2026-04-06 2026-07-02 75 13
## 13 D02 2026-04-06 2026-07-07 78 15
## 14 D01 2026-04-06 2026-07-07 82 11
## 15 CW03 2026-04-06 2026-07-07 83 10
## 16 D15 2026-04-07 2026-07-07 83 9
## 17 CW18 2026-04-15 2026-07-07 84 0
## 18 CW19 2026-04-15 2026-07-07 84 0
## 19 CW20 2026-04-15 2026-07-07 84 0
## 20 D06 2026-04-07 2026-07-02 87 0
## 21 CE09 2026-04-08 2026-07-07 88 3
## 22 CW11 2026-04-08 2026-07-05 89 0
## 23 D08 2026-04-07 2026-07-07 89 3
## 24 CE08 2026-04-07 2026-07-07 90 2
## 25 CE01 2026-04-08 2026-07-07 91 0
## 26 CE02 2026-04-08 2026-07-07 91 0
## 27 CE05 2026-04-08 2026-07-07 91 0
## 28 CE06 2026-04-08 2026-07-07 91 0
## 29 CE07 2026-04-08 2026-07-07 91 0
## 30 CE11 2026-04-08 2026-07-07 91 0
## 31 CE13 2026-04-08 2026-07-07 91 0
## 32 CE14 2026-04-08 2026-07-07 91 0
## 33 CW12 2026-04-08 2026-07-07 91 0
## 34 CW15 2026-04-08 2026-07-07 91 0
## 35 CW16 2026-04-08 2026-07-07 91 0
## 36 CE04 2026-04-07 2026-07-07 92 0
## 37 CW07 2026-04-06 2026-07-07 92 1
## 38 D09 2026-04-07 2026-07-07 92 0
## 39 D10 2026-04-07 2026-07-07 92 0
## 40 D11 2026-04-07 2026-07-07 92 0
## 41 D12 2026-04-07 2026-07-07 92 0
## 42 D13 2026-04-07 2026-07-07 92 0
## 43 D14 2026-04-07 2026-07-07 92 0
## 44 CW02 2026-04-06 2026-07-07 93 0
## 45 CW04 2026-04-06 2026-07-07 93 0
## 46 CW06 2026-04-06 2026-07-07 93 0
## 47 CW08 2026-04-06 2026-07-07 93 0
## 48 CW10 2026-04-06 2026-07-07 93 0
## 49 D04 2026-04-06 2026-07-07 93 0
## 50 D05 2026-04-06 2026-07-07 93 0cat("\nTotal camera-days:", sum(effort$camera_days),
"| range:", min(effort$camera_days), "-", max(effort$camera_days), "\n")##
## Total camera-days: 4257 | range: 64 - 93Community matrix, as a function of the rarity threshold
det_long <- CamInd %>%
count(Plot, Common_Name, name = "detections") %>%
left_join(effort %>% dplyr::select(Plot, camera_days), by = "Plot") %>%
mutate(rate100 = detections / camera_days * 100)
sp_freq <- det_long %>% count(Common_Name, name = "n_plots") %>% arrange(n_plots)
wild_summary <- CamInd %>%
group_by(Plot) %>%
summarise(n_detections = n(), n_species = n_distinct(Common_Name),
.groups = "drop") %>%
left_join(effort %>% dplyr::select(Plot, camera_days), by = "Plot") %>%
mutate(det_rate100 = n_detections / camera_days * 100)
create_matrix_by_rarity <- function(rarity_threshold, det_df = det_long) {
wide <- det_df %>%
dplyr::select(Plot, Common_Name, rate100) %>%
pivot_wider(names_from = Common_Name, values_from = rate100,
values_fill = 0) %>%
arrange(Plot)
m <- as.matrix(wide[, -1])
rownames(m) <- wide$Plot
m <- m[, colSums(m) > 0, drop = FALSE]
# occurrence as a percentage of CWDs, computed on the full matrix
occ_pct <- colSums(m > 0) / nrow(m) * 100
m <- m[, occ_pct >= rarity_threshold, drop = FALSE]
empty <- rowSums(m) == 0
if (any(empty)) {
message(sprintf(" threshold %g%%: dropped %d log(s) with no detections left: %s",
rarity_threshold, sum(empty),
paste(rownames(m)[empty], collapse = ", ")))
m <- m[!empty, , drop = FALSE]
}
list(matrix = m, plots = rownames(m), n_species = ncol(m))
}
wild_matrix <- create_matrix_by_rarity(0)$matrix %>%
as.data.frame() %>% tibble::rownames_to_column("Plot") %>% as_tibble()
full <- create_matrix_by_rarity(0)
occ_spectrum <- tibble(
threshold_pct = c(0, 5, 10, 15, 20, 25, 30, 40, 50),
n_species = sapply(threshold_pct, function(t)
sum(colSums(full$matrix > 0) / nrow(full$matrix) * 100 >= t)))
cat("Species retained at each occupancy threshold:\n")## Species retained at each occupancy threshold:print(as.data.frame(occ_spectrum), row.names = FALSE)## threshold_pct n_species
## 0 54
## 5 35
## 10 26
## 15 21
## 20 19
## 25 15
## 30 14
## 40 12
## 50 9Site covariates
Daubenmire cover classes converted to midpoint percentages before averaging; tree density and basal area from the point-centred-quarter data by the Cottam and Curtis estimator.
Master <- read_excel(veg_path, sheet = "Master")
Point <- read_excel(veg_path, sheet = "Point")
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")
master_mid <- 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))
site_struct <- master_mid %>%
group_by(cwd) %>%
summarise(
n_quadrats = n_distinct(q),
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")
gc_matrix <- site_struct %>% dplyr::select(cwd, ends_with("_pct"))
site_struct$gc_shannon <- vegan::diversity(as.matrix(gc_matrix[, -1]),
index = "shannon")
pcq <- Point %>%
group_by(cwd) %>%
summarise(
mean_dist = mean(dist, na.rm = TRUE),
tree_dens = 10000 / (mean(dist, na.rm = TRUE)^2), # stems / ha
mean_dbh = mean(dbh, na.rm = TRUE),
basal_area = (10000 / (mean(dist, na.rm = TRUE)^2)) *
mean(pi * (dbh / 200)^2, na.rm = TRUE), # m2 / ha
tree_rich = n_distinct(species),
prop_pine = mean(str_detect(species, "^pin"), na.rm = TRUE),
.groups = "drop")
site_env <- site_struct %>%
left_join(pcq, by = "cwd") %>%
rename(Plot = cwd) %>%
mutate(Area = str_extract(Plot, "^[A-Z]+"),
decay_f = factor(decay, levels = 1:5)) %>%
left_join(wild_summary, by = "Plot")align_env <- function(plots) {
e <- site_env %>% filter(Plot %in% plots) %>% arrange(match(Plot, plots))
stopifnot(identical(e$Plot, plots))
e$decay_f <- factor(e$decay, levels = 1:5)
droplevels(e)
}Part 1 – How many dimensions, and which species?
rarity_thresholds <- c(0, 5, 10, 20, 30)
k_values <- 1:6
all_stress_results <- list()
ord_k2 <- list()
for (thresh in rarity_thresholds) {
cat("Rarity threshold:", thresh, "%\n")
mm <- create_matrix_by_rarity(thresh)
if (mm$n_species < 3) {
cat(" skipping: fewer than three species remain.\n"); next
}
for (k_val in k_values) {
if (k_val >= mm$n_species) next
set.seed(SEED)
nmds_run <- try(suppressWarnings(
metaMDS(mm$matrix, distance = "bray", k = k_val, trymax = 100,
autotransform = FALSE, trace = 0)), silent = TRUE)
if (inherits(nmds_run, "try-error")) {
cat(" k =", k_val, "failed to fit.\n"); next
}
sp <- stressplot(nmds_run)
all_stress_results[[length(all_stress_results) + 1]] <- data.frame(
rarity = thresh,
k = k_val,
stress = nmds_run$stress,
num_species = mm$n_species,
num_sites = nrow(mm$matrix),
converged = nmds_run$converged,
tries = nmds_run$tries,
nonmetric_R2 = 1 - nmds_run$stress^2,
linear_R2 = summary(lm(sp$y ~ sp$x))$r.squared)
if (k_val == 2) ord_k2[[as.character(thresh)]] <- nmds_run
}
}## Rarity threshold: 0 %
## Rarity threshold: 5 %
## Rarity threshold: 10 %
## Rarity threshold: 20 %
## Rarity threshold: 30 %
all_stress_results <- bind_rows(all_stress_results)
print(all_stress_results, row.names = FALSE)## rarity k stress num_species num_sites converged tries nonmetric_R2
## 0 1 0.27745191 54 50 0 100 0.9230204
## 0 2 0.14721968 54 50 1 21 0.9783264
## 0 3 0.11080584 54 50 1 20 0.9877221
## 0 4 0.08939230 54 50 4 20 0.9920090
## 0 5 0.07528838 54 50 2 20 0.9943317
## 0 6 0.06352285 54 50 1 20 0.9959648
## 5 1 0.27660483 35 50 0 100 0.9234898
## 5 2 0.14673253 35 50 1 21 0.9784696
## 5 3 0.11045086 35 50 1 20 0.9878006
## 5 4 0.08967853 35 50 1 20 0.9919578
## 5 5 0.07542908 35 50 1 20 0.9943105
## 5 6 0.06374058 35 50 1 20 0.9959371
## 10 1 0.26038629 26 50 0 100 0.9321990
## 10 2 0.14282093 26 50 3 20 0.9796022
## 10 3 0.10752605 26 50 1 20 0.9884381
## 10 4 0.08784644 26 50 2 20 0.9922830
## 10 5 0.07378341 26 50 1 20 0.9945560
## 10 6 0.06208730 26 50 1 20 0.9961452
## 20 1 0.25980296 19 50 0 100 0.9325024
## 20 2 0.14295950 19 50 1 20 0.9795626
## 20 3 0.10660754 19 50 4 20 0.9886348
## 20 4 0.08709031 19 50 1 20 0.9924153
## 20 5 0.07294574 19 50 2 20 0.9946789
## 20 6 0.06113806 19 50 1 20 0.9962621
## 30 1 0.25859286 14 50 1 23 0.9331297
## 30 2 0.13862362 14 50 1 38 0.9807835
## 30 3 0.10262514 14 50 6 20 0.9894681
## 30 4 0.08293283 14 50 3 20 0.9931221
## 30 5 0.06906161 14 50 1 20 0.9952305
## 30 6 0.05773091 14 50 1 59 0.9966671
## linear_R2
## 0.7186121
## 0.8581618
## 0.8937193
## 0.9092086
## 0.9273824
## 0.9480299
## 0.7165334
## 0.8590869
## 0.8941261
## 0.9099970
## 0.9299280
## 0.9489104
## 0.7339685
## 0.8651812
## 0.8981770
## 0.9104640
## 0.9296426
## 0.9500756
## 0.7373028
## 0.8668018
## 0.8990335
## 0.9107490
## 0.9325809
## 0.9557833
## 0.7356371
## 0.8703212
## 0.9042052
## 0.9164219
## 0.9434107
## 0.9604157write.csv(all_stress_results,
file.path(out_dir, "tables", "NMDS_stress_grid.csv"), row.names = FALSE)Elbow table
elbow <- all_stress_results %>%
arrange(rarity, k) %>%
group_by(rarity) %>%
mutate(drop_from_prev = lag(stress) - stress) %>%
ungroup()
print(elbow %>% filter(rarity == 0) %>%
dplyr::select(k, num_species, stress, drop_from_prev, converged) %>%
mutate(across(where(is.numeric), ~ round(.x, 4))), row.names = FALSE)## # A tibble: 6 × 5
## k num_species stress drop_from_prev converged
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 54 0.278 NA 0
## 2 2 54 0.147 0.130 1
## 3 3 54 0.111 0.0364 1
## 4 4 54 0.0894 0.0214 4
## 5 5 54 0.0753 0.0141 2
## 6 6 54 0.0635 0.0118 1NMDS stability across rarity thresholds
ref <- ord_k2[[as.character(min(rarity_thresholds))]]
procrustes_tab <- bind_rows(lapply(names(ord_k2), function(nm) {
o <- ord_k2[[nm]]
common <- intersect(rownames(scores(ref, "sites")),
rownames(scores(o, "sites")))
set.seed(SEED)
pt <- protest(scores(ref, "sites")[common, ],
scores(o, "sites")[common, ], permutations = 999)
data.frame(rarity = as.numeric(nm),
n_species = nrow(scores(o, "species")),
procrustes_r = pt$scale,
ss = pt$ss,
p = pt$signif,
n_shared = length(common))
}))
print(procrustes_tab %>% mutate(across(where(is.numeric), ~ round(.x, 4))),
row.names = FALSE)## rarity n_species procrustes_r ss p n_shared
## 0 54 1.0000 0.0000 0.001 50
## 5 35 0.9998 0.0005 0.001 50
## 10 26 0.9853 0.0292 0.001 50
## 20 19 0.9853 0.0293 0.001 50
## 30 14 0.9834 0.0329 0.001 50write.csv(procrustes_tab,
file.path(out_dir, "tables", "NMDS_procrustes_rarity.csv"),
row.names = FALSE)PERMANOVA sensitivity to rarity threshold
permanova_by_rarity <- bind_rows(lapply(rarity_thresholds, function(thresh) {
mm <- create_matrix_by_rarity(thresh)
if (mm$n_species < 3) return(NULL)
e <- align_env(mm$plots)
set.seed(SEED)
# same specification as the main model, with blocks rebuilt from whichever
# logs survived this threshold
a <- adonis2(vegdist(mm$matrix, "bray") ~ decay_f + canopy_cover +
shrubs_pct + basal_area,
data = e, permutations = how(blocks = factor(e$Area), nperm = 9999),
by = "margin")
as.data.frame(a) %>%
tibble::rownames_to_column("Term") %>%
filter(!Term %in% c("Residual", "Total")) %>%
mutate(rarity = thresh, n_species = mm$n_species, n_sites = nrow(mm$matrix))
}))
permanova_wide <- permanova_by_rarity %>%
dplyr::select(rarity, n_species, Term, R2, p = `Pr(>F)`) %>%
mutate(R2 = round(R2, 4), p = round(p, 4)) %>%
pivot_wider(names_from = Term, values_from = c(R2, p))
print(permanova_wide, row.names = FALSE, width = Inf)## # A tibble: 5 × 10
## rarity n_species R2_decay_f R2_canopy_cover R2_shrubs_pct R2_basal_area
## <dbl> <int> <dbl> <dbl> <dbl> <dbl>
## 1 0 54 0.0824 0.075 0.015 0.0316
## 2 5 35 0.0826 0.0762 0.0149 0.0317
## 3 10 26 0.0798 0.0808 0.0139 0.0321
## 4 20 19 0.0785 0.0827 0.0147 0.033
## 5 30 14 0.076 0.0842 0.0155 0.0337
## p_decay_f p_canopy_cover p_shrubs_pct p_basal_area
## <dbl> <dbl> <dbl> <dbl>
## 1 0.266 0.0008 0.666 0.113
## 2 0.261 0.0008 0.670 0.113
## 3 0.305 0.0008 0.703 0.117
## 4 0.324 0.0008 0.648 0.111
## 5 0.363 0.0009 0.596 0.114write.csv(permanova_by_rarity,
file.path(out_dir, "tables", "PERMANOVA_by_rarity.csv"), row.names = FALSE)
flip <- permanova_by_rarity %>%
group_by(Term) %>%
summarise(min_p = min(`Pr(>F)`), max_p = max(`Pr(>F)`),
crosses_05 = min(`Pr(>F)`) < 0.05 & max(`Pr(>F)`) >= 0.05,
.groups = "drop")
cat("\nTerms whose p-value crosses 0.05 somewhere on the rarity grid:\n")##
## Terms whose p-value crosses 0.05 somewhere on the rarity grid:print(as.data.frame(flip %>% mutate(across(where(is.numeric), ~ round(.x, 4)))),
row.names = FALSE)## Term min_p max_p crosses_05
## basal_area 0.1111 0.1167 FALSE
## canopy_cover 0.0008 0.0009 FALSE
## decay_f 0.2608 0.3633 FALSE
## shrubs_pct 0.5965 0.7029 FALSEPart 1 figures
p_scree <- ggplot(all_stress_results,
aes(x = k, y = stress, color = as.factor(rarity))) +
geom_hline(yintercept = 0.20, linetype = "dashed", color = "grey50") +
geom_hline(yintercept = 0.10, linetype = "dashed", color = "grey50") +
geom_line(linewidth = 0.9) +
geom_point(size = 2.6) +
annotate("text", x = max(k_values), y = 0.205, hjust = 1,
label = "Stress = 0.20 (fair)", size = 3, color = "grey35") +
annotate("text", x = max(k_values), y = 0.105, hjust = 1,
label = "Stress = 0.10 (good)", size = 3, color = "grey35") +
scale_color_viridis_d(name = "Rarity filter\n(% of logs)", end = 0.9) +
scale_x_continuous(breaks = k_values) +
scale_y_continuous(limits = c(0, NA)) +
labs(x = "Number of Dimensions (k)", y = "Stress") +
theme_bw(base_size = 12) +
theme(legend.position = "bottom")
print(p_scree)
ggsave(file.path(out_dir, "figures", "NMDS_scree_wildlife.png"), p_scree,
width = 7, height = 5, dpi = 300)p_elbow <- elbow %>%
filter(!is.na(drop_from_prev)) %>%
ggplot(aes(factor(k), drop_from_prev, fill = as.factor(rarity))) +
geom_hline(yintercept = 0.02, linetype = "dashed", color = "grey50") +
geom_col(position = position_dodge(width = 0.8), width = 0.75) +
scale_fill_viridis_d(name = "Rarity filter\n(% of logs)", end = 0.9) +
labs(x = "Dimension added (k-1 \u2192 k)", y = "Reduction in stress") +
theme_bw(base_size = 12) +
theme(legend.position = "bottom",
plot.caption = element_text(hjust = 0, color = "grey35"))
print(p_elbow)
ggsave(file.path(out_dir, "figures", "NMDS_elbow_wildlife.png"), p_elbow,
width = 7, height = 5, dpi = 300)term_cols <- c(decay_f = "#D55E00", canopy_cover = "#009E73",
shrubs_pct = "#0072B2", basal_area = "#CC79A7")
p_sens <- permanova_by_rarity %>%
mutate(Term = factor(Term, levels = names(term_cols))) %>%
ggplot(aes(rarity, `Pr(>F)`, color = Term, group = Term)) +
geom_hline(yintercept = 0.05, linetype = "dashed", color = "grey40") +
geom_line(linewidth = 0.9) + geom_point(size = 2.6) +
annotate("text", x = max(rarity_thresholds), y = 0.056, hjust = 1,
label = "p = 0.05", size = 3, color = "grey35") +
scale_color_manual(values = term_cols, name = "Term") +
scale_x_continuous(breaks = rarity_thresholds) +
labs(x = "Rarity filter (% of logs a species must occupy)",
y = "PERMANOVA p-value") +
theme_bw(base_size = 12)
p_r2 <- permanova_by_rarity %>%
mutate(Term = factor(Term, levels = names(term_cols))) %>%
ggplot(aes(rarity, R2, color = Term, group = Term)) +
geom_line(linewidth = 0.9) + geom_point(size = 2.6) +
scale_color_manual(values = term_cols, name = "Term") +
scale_x_continuous(breaks = rarity_thresholds) +
labs(x = "Rarity filter (% of logs a species must occupy)",
y = expression(PERMANOVA~R^2)) +
theme_bw(base_size = 12)
p_sensitivity <- (p_sens / p_r2) +
plot_layout(guides = "collect") +
plot_annotation(tag_levels = "A") &
theme(legend.position = "bottom")
print(p_sensitivity)
ggsave(file.path(out_dir, "figures", "PERMANOVA_rarity_sensitivity.png"),
p_sensitivity, width = 7, height = 8, dpi = 300)The settings Part 2 uses
chosen_k <- 4
chosen_rarity <- 0
chosen_row <- all_stress_results %>%
filter(k == chosen_k, rarity == chosen_rarity)
cat(sprintf("Selected: k = %d, rarity filter = %g%% (%d species, stress = %.3f)\n",
chosen_k, chosen_rarity, chosen_row$num_species, chosen_row$stress))## Selected: k = 4, rarity filter = 0% (54 species, stress = 0.089)p_justify <- p_scree +
geom_point(data = chosen_row, aes(x = k, y = stress),
color = "black", size = 6, shape = 1, stroke = 1,
inherit.aes = FALSE) +
labs(subtitle = sprintf("Selected: k = %d, no rarity filter (%d species, stress = %.3f)",
chosen_k, chosen_row$num_species, chosen_row$stress))
final_figure <- (p_justify / p_sens) +
plot_annotation(tag_levels = "A") &
theme(plot.tag = element_text(face = "bold"))
print(final_figure)
ggsave(file.path(out_dir, "figures", "Figure_S_NMDS_justification.png"),
final_figure, width = 7.5, height = 9, dpi = 300)Part 2 – The analysis at those settings
mm_main <- create_matrix_by_rarity(chosen_rarity)
comm_m <- mm_main$matrix
env_m <- align_env(mm_main$plots)
dist_w <- vegdist(comm_m, method = "bray")
cat("Analysis matrix:", nrow(comm_m), "logs x", ncol(comm_m), "species\n")## Analysis matrix: 50 logs x 54 speciesprint(env_m %>% dplyr::select(Plot, decay, n_detections, n_species,
camera_days, det_rate100) %>%
arrange(n_detections), n = 50)## # A tibble: 50 × 6
## Plot decay n_detections n_species camera_days det_rate100
## <chr> <dbl> <int> <int> <int> <dbl>
## 1 CE01 5 8 3 91 8.79
## 2 CW05 2 14 5 66 21.2
## 3 D08 1 14 8 89 15.7
## 4 D13 4 19 7 92 20.7
## 5 CE04 3 22 8 92 23.9
## 6 CW08 3 23 9 93 24.7
## 7 CE09 1 27 10 88 30.7
## 8 D12 5 28 9 92 30.4
## 9 D09 2 29 10 92 31.5
## 10 CW01 2 30 7 75 40
## 11 CW07 1 30 8 92 32.6
## 12 CE13 1 31 8 91 34.1
## 13 D14 4 33 11 92 35.9
## 14 D06 3 34 7 87 39.1
## 15 CE07 2 36 12 91 39.6
## 16 D05 3 38 8 93 40.9
## 17 CW10 3 39 8 93 41.9
## 18 D03 2 39 8 70 55.7
## 19 D15 2 39 11 83 47.0
## 20 CW19 3 41 9 84 48.8
## 21 CW20 5 43 12 84 51.2
## 22 CW02 4 45 9 93 48.4
## 23 D11 4 50 14 92 54.3
## 24 CE08 3 51 10 90 56.7
## 25 CW12 1 55 12 91 60.4
## 26 CE11 2 56 12 91 61.5
## 27 CE14 3 61 13 91 67.0
## 28 CE02 4 63 10 91 69.2
## 29 D02 3 69 12 78 88.5
## 30 CE12 5 70 12 75 93.3
## 31 CE15 3 72 19 75 96
## 32 CW13 4 76 18 74 103.
## 33 CW16 5 82 16 91 90.1
## 34 CW15 5 83 15 91 91.2
## 35 CW18 5 84 12 84 100
## 36 CW06 5 88 10 93 94.6
## 37 CE10 4 89 14 74 120.
## 38 CW03 2 89 13 83 107.
## 39 CW09 1 91 13 69 132.
## 40 CW17 4 92 17 65 142.
## 41 D10 3 104 12 92 113.
## 42 D04 3 110 21 93 118.
## 43 CW14 3 120 15 64 188.
## 44 CE06 2 129 19 91 142.
## 45 D07 3 132 15 71 186.
## 46 CW11 4 177 17 89 199.
## 47 CE05 4 179 15 91 197.
## 48 CW04 4 194 14 93 209.
## 49 CE03 5 195 13 70 279.
## 50 D01 5 207 16 82 252.The model
perm_formula <- dist_w ~ decay_f + canopy_cover + shrubs_pct + basal_area # variables of interest
na_check <- sapply(env_m[, c("canopy_cover", "shrubs_pct", "basal_area")],
function(x) sum(is.na(x)))
cat("Missing values in the model covariates:\n"); print(na_check)## Missing values in the model covariates:## canopy_cover shrubs_pct basal_area
## 0 0 0stopifnot(all(na_check == 0))
cat("\nCWD per stand x decay class (blocking needs within-block variation):\n")##
## CWD per stand x decay class (blocking needs within-block variation):print(table(env_m$Area, env_m$decay_f))##
## 1 2 3 4 5
## CE 2 3 4 3 3
## CW 3 3 4 5 5
## D 1 3 6 3 2Predictor collinearity
cand <- c("decay", "canopy_cover", "soil_moist", "gc_shannon",
"pl_pct", "dll_pct", "fw_pct", "lg_pct", "forbs_pct",
"shrubs_pct", "vines_pct", "bg_pct", "off_ground",
"tree_dens", "basal_area", "prop_pine")
cor_mat <- cor(env_m[, cand], use = "pairwise.complete.obs")
print(round(cor_mat, 2))## decay canopy_cover soil_moist gc_shannon pl_pct dll_pct fw_pct
## decay 1.00 0.27 0.00 -0.20 -0.35 0.36 -0.23
## canopy_cover 0.27 1.00 -0.14 -0.62 0.00 0.68 -0.27
## soil_moist 0.00 -0.14 1.00 0.08 0.11 -0.26 0.17
## gc_shannon -0.20 -0.62 0.08 1.00 0.04 -0.52 0.14
## pl_pct -0.35 0.00 0.11 0.04 1.00 -0.41 0.11
## dll_pct 0.36 0.68 -0.26 -0.52 -0.41 1.00 -0.47
## fw_pct -0.23 -0.27 0.17 0.14 0.11 -0.47 1.00
## lg_pct 0.00 -0.54 0.13 0.49 -0.05 -0.53 0.08
## forbs_pct -0.18 -0.60 0.11 0.57 -0.03 -0.51 0.02
## shrubs_pct 0.10 -0.05 -0.16 0.20 -0.25 0.07 0.23
## vines_pct -0.04 0.00 0.16 0.26 0.16 0.06 0.00
## bg_pct -0.21 -0.39 0.00 0.58 -0.22 -0.23 -0.02
## off_ground -0.37 -0.12 0.15 0.06 0.11 -0.16 0.30
## tree_dens 0.10 0.36 -0.07 -0.43 0.00 0.34 -0.26
## basal_area -0.07 0.41 0.00 -0.14 0.42 0.04 -0.21
## prop_pine -0.24 -0.71 0.23 0.57 0.32 -0.71 0.34
## lg_pct forbs_pct shrubs_pct vines_pct bg_pct off_ground tree_dens
## decay 0.00 -0.18 0.10 -0.04 -0.21 -0.37 0.10
## canopy_cover -0.54 -0.60 -0.05 0.00 -0.39 -0.12 0.36
## soil_moist 0.13 0.11 -0.16 0.16 0.00 0.15 -0.07
## gc_shannon 0.49 0.57 0.20 0.26 0.58 0.06 -0.43
## pl_pct -0.05 -0.03 -0.25 0.16 -0.22 0.11 0.00
## dll_pct -0.53 -0.51 0.07 0.06 -0.23 -0.16 0.34
## fw_pct 0.08 0.02 0.23 0.00 -0.02 0.30 -0.26
## lg_pct 1.00 0.44 -0.19 -0.11 0.33 -0.10 -0.17
## forbs_pct 0.44 1.00 -0.02 0.01 0.36 -0.01 -0.26
## shrubs_pct -0.19 -0.02 1.00 -0.04 0.03 0.04 -0.21
## vines_pct -0.11 0.01 -0.04 1.00 -0.07 0.18 -0.14
## bg_pct 0.33 0.36 0.03 -0.07 1.00 -0.06 -0.33
## off_ground -0.10 -0.01 0.04 0.18 -0.06 1.00 -0.01
## tree_dens -0.17 -0.26 -0.21 -0.14 -0.33 -0.01 1.00
## basal_area -0.12 -0.24 -0.33 0.18 -0.24 0.20 0.41
## prop_pine 0.35 0.43 0.10 0.17 0.31 0.13 -0.41
## basal_area prop_pine
## decay -0.07 -0.24
## canopy_cover 0.41 -0.71
## soil_moist 0.00 0.23
## gc_shannon -0.14 0.57
## pl_pct 0.42 0.32
## dll_pct 0.04 -0.71
## fw_pct -0.21 0.34
## lg_pct -0.12 0.35
## forbs_pct -0.24 0.43
## shrubs_pct -0.33 0.10
## vines_pct 0.18 0.17
## bg_pct -0.24 0.31
## off_ground 0.20 0.13
## tree_dens 0.41 -0.41
## basal_area 1.00 -0.07
## prop_pine -0.07 1.00write.csv(round(cor_mat, 3),
file.path(out_dir, "tables", "predictor_correlations.csv"))PERMANOVA
bd <- betadisper(dist_w, env_m$decay_f)
print(anova(bd))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 4 0.08713 0.021783 1.4287 0.24
## Residuals 45 0.68611 0.015247set.seed(SEED); print(permutest(bd, permutations = 9999))##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 9999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 4 0.08713 0.021783 1.4287 9999 0.2491
## Residuals 45 0.68611 0.015247set.seed(SEED)
h <- how(blocks = factor(env_m$Area), nperm = 9999) # Blocked by stand
perm_main <- adonis2(perm_formula, data = env_m, permutations = h,
by = "margin")
print(perm_main)## Permutation test for adonis under reduced model
## Marginal effects of terms
## Blocks: factor(env_m$Area)
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = perm_formula, data = env_m, permutations = h, by = "margin")
## Df SumOfSqs R2 F Pr(>F)
## decay_f 4 0.8639 0.08243 1.1161 0.2658
## canopy_cover 1 0.7859 0.07499 4.0616 0.0008 ***
## shrubs_pct 1 0.1571 0.01499 0.8119 0.6659
## basal_area 1 0.3312 0.03160 1.7113 0.1130
## Residual 42 8.1272 0.77545
## Total 49 10.4807 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tidy_adonis <- function(x, label) {
as.data.frame(x) %>%
tibble::rownames_to_column("Term") %>%
mutate(Model = label, across(where(is.numeric), ~ round(.x, 4)))
}
perm_tables <- tidy_adonis(
perm_main, "Marginal, decay as factor, blocked by stand (50 CWD)")
write.csv(perm_tables, file.path(out_dir, "tables", "PERMANOVA_50sites.csv"),
row.names = FALSE)H1: Does community variability increase with decay?
bd_h1 <- betadisper(dist_w, env_m$decay_f, type = "centroid")
bd_h1_med <- betadisper(dist_w, env_m$decay_f, type = "median")
set.seed(SEED)
pt_h1_med <- permutest(bd_h1_med, permutations = 9999)
sp_h1_med <- suppressWarnings(
cor.test(as.integer(as.character(env_m$decay_f)), bd_h1_med$distances,
method = "spearman"))
disp_by_class <- tibble(
decay = as.integer(as.character(env_m$decay_f)),
dist_cent = bd_h1$distances) %>%
group_by(decay) %>%
summarise(n_logs = n(), mean_dist = mean(dist_cent),
sd_dist = sd(dist_cent), .groups = "drop")
cat("Mean distance to decay-class centroid (higher = more variable community):\n")## Mean distance to decay-class centroid (higher = more variable community):print(as.data.frame(disp_by_class %>%
mutate(across(where(is.numeric), ~ round(.x, 4)))), row.names = FALSE)## decay n_logs mean_dist sd_dist
## 1 6 0.4459 0.0693
## 2 9 0.3550 0.0821
## 3 14 0.4015 0.1253
## 4 11 0.4419 0.0904
## 5 10 0.4798 0.1113set.seed(SEED)
pt_h1 <- permutest(bd_h1, permutations = 9999)
print(pt_h1)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 9999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 4 0.08738 0.021844 2.0679 9999 0.0982 .
## Residuals 45 0.47536 0.010563
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sp_h1 <- suppressWarnings(
cor.test(as.integer(as.character(env_m$decay_f)), bd_h1$distances,
method = "spearman"))
cat(sprintf("\nSpearman, decay class vs distance-to-centroid: rho = %.3f, p = %.4f\n",
sp_h1$estimate, sp_h1$p.value))##
## Spearman, decay class vs distance-to-centroid: rho = 0.222, p = 0.1214decay3 <- cut(as.integer(as.character(env_m$decay_f)),
breaks = c(0, 2, 3, 5),
labels = c("early (1-2)", "mid (3)", "late (4-5)"))
bd_h1b <- betadisper(dist_w, decay3, type = "centroid")
set.seed(SEED)
pt_h1b <- permutest(bd_h1b, permutations = 9999)
cat("\nMean distance to centroid, three-level grouping:\n")##
## Mean distance to centroid, three-level grouping:print(round(tapply(bd_h1b$distances, decay3, mean), 4))## early (1-2) mid (3) late (4-5)
## 0.4115 0.4015 0.4650print(pt_h1b)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 9999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.04216 0.021079 1.8061 9999 0.1806
## Residuals 47 0.54852 0.011671h1_table <- bind_rows(
disp_by_class %>%
transmute(Grouping = "Decay class (1-5)", Level = as.character(decay),
n_logs, mean_dist = round(mean_dist, 4),
sd_dist = round(sd_dist, 4)),
tibble(Grouping = "Decay grouped (3 levels)",
Level = levels(decay3),
n_logs = as.integer(table(decay3)),
mean_dist = round(as.numeric(tapply(bd_h1b$distances, decay3, mean)), 4),
sd_dist = round(as.numeric(tapply(bd_h1b$distances, decay3, sd)), 4)))
h1_tests <- tibble(
Test = c("permutest, 5 classes (centroid)",
"permutest, 5 classes (spatial median)",
"Spearman trend, 5 classes (centroid)",
"Spearman trend, 5 classes (spatial median)",
"permutest, 3 groups (centroid)"),
Statistic = c(round(pt_h1$tab$F[1], 3), round(pt_h1_med$tab$F[1], 3),
round(sp_h1$estimate, 3), round(sp_h1_med$estimate, 3),
round(pt_h1b$tab$F[1], 3)),
p = c(round(pt_h1$tab$`Pr(>F)`[1], 4), round(pt_h1_med$tab$`Pr(>F)`[1], 4),
round(sp_h1$p.value, 4), round(sp_h1_med$p.value, 4),
round(pt_h1b$tab$`Pr(>F)`[1], 4)))
cat("\nH1 test summary (all five tests, both dispersion definitions):\n")##
## H1 test summary (all five tests, both dispersion definitions):print(as.data.frame(h1_tests), row.names = FALSE)## Test Statistic p
## permutest, 5 classes (centroid) 2.068 0.0982
## permutest, 5 classes (spatial median) 1.429 0.2491
## Spearman trend, 5 classes (centroid) 0.222 0.1214
## Spearman trend, 5 classes (spatial median) 0.188 0.1922
## permutest, 3 groups (centroid) 1.806 0.1806if (any(h1_tests$p < 0.05) && any(h1_tests$p >= 0.05)) {
cat("\nWARNING: these tests disagree about H1 at the 0.05 level.",
"\n Report the range, not the favourable one.\n")
}
write.csv(h1_table, file.path(out_dir, "tables", "H1_dispersion_by_decay.csv"),
row.names = FALSE)
write.csv(h1_tests, file.path(out_dir, "tables", "H1_dispersion_tests.csv"),
row.names = FALSE)
cat("\nH1 is supported only if dispersion increases AND a test says so.",
"\nDirection without significance is a trend, not a result.\n")##
## H1 is supported only if dispersion increases AND a test says so.
## Direction without significance is a trend, not a result.h1_plot_df <- tibble(
decay = factor(as.integer(as.character(env_m$decay_f)), levels = 1:5),
dist_cent = bd_h1$distances)
p_h1 <- ggplot(h1_plot_df, aes(decay, dist_cent)) +
geom_boxplot(outlier.shape = NA, width = 0.55, fill = "grey92",
colour = "grey40") +
geom_jitter(width = 0.12, height = 0, size = 2, alpha = 0.6,
colour = "#0072B2") +
stat_summary(fun = mean, geom = "point", shape = 23, size = 3,
fill = "#D55E00", colour = "black") +
geom_smooth(aes(x = as.numeric(decay)), method = "lm", se = TRUE,
colour = "#D55E00", fill = "#D55E00", alpha = 0.12,
linewidth = 0.7) +
labs(x = "CWD decay class",
y = "Distance to decay-class centroid\n(Bray-Curtis)",
caption = sprintf(
"Diamonds are class means. Spearman rho = %.3f, p = %.3f; permutest F = %.2f, p = %.3f.",
sp_h1$estimate, sp_h1$p.value, pt_h1$tab$F[1], pt_h1$tab$`Pr(>F)`[1])) +
theme_classic(base_size = 12) +
theme(axis.title = element_text(face = "bold"),
plot.caption = element_text(hjust = 0, colour = "grey35"))
print(p_h1)## `geom_smooth()` using formula = 'y ~ x'
ggsave(file.path(out_dir, "figures", "H1_dispersion_by_decay.png"), p_h1,
width = 7, height = 5, dpi = 300, bg = "white")## `geom_smooth()` using formula = 'y ~ x'NMDS and fitted environmental vectors
set.seed(SEED)
nmds <- metaMDS(comm_m, distance = "bray", k = chosen_k, trymax = 500,
autotransform = FALSE)## Run 0 stress 0.09062785
## Run 1 stress 0.08948408
## ... New best solution
## ... Procrustes: rmse 0.05050491 max resid 0.1950346
## Run 2 stress 0.09057991
## Run 3 stress 0.09174185
## Run 4 stress 0.08939559
## ... New best solution
## ... Procrustes: rmse 0.006698504 max resid 0.028888
## Run 5 stress 0.0907999
## Run 6 stress 0.09030725
## Run 7 stress 0.09384273
## Run 8 stress 0.0893923
## ... New best solution
## ... Procrustes: rmse 0.00150335 max resid 0.004392659
## ... Similar to previous best
## Run 9 stress 0.0907962
## Run 10 stress 0.08939303
## ... Procrustes: rmse 0.0006807937 max resid 0.00242556
## ... Similar to previous best
## Run 11 stress 0.09082421
## Run 12 stress 0.08940388
## ... Procrustes: rmse 0.002236753 max resid 0.008106074
## ... Similar to previous best
## Run 13 stress 0.09079712
## Run 14 stress 0.09219705
## Run 15 stress 0.09047851
## Run 16 stress 0.08968206
## ... Procrustes: rmse 0.01471158 max resid 0.05220067
## Run 17 stress 0.09030929
## Run 18 stress 0.08963428
## ... Procrustes: rmse 0.01032478 max resid 0.04208863
## Run 19 stress 0.08939403
## ... Procrustes: rmse 0.001056475 max resid 0.003813016
## ... Similar to previous best
## Run 20 stress 0.08941592
## ... Procrustes: rmse 0.005413471 max resid 0.02837284
## *** Best solution repeated 4 timescat("NMDS stress (k =", chosen_k, "):", round(nmds$stress, 4),
"| convergent solutions:", nmds$converged, "\n")## NMDS stress (k = 4 ): 0.0894 | convergent solutions: 4png(file.path(out_dir, "figures", "NMDS_stressplot.png"), 900, 700)
stressplot(nmds); dev.off()## png
## 2set.seed(SEED)
env_fit <- envfit(nmds ~ decay + canopy_cover + soil_moist + gc_shannon +
shrubs_pct + fw_pct + dll_pct + lg_pct + basal_area +
tree_dens + off_ground,
data = env_m, permutations = 9999, na.rm = TRUE)
print(env_fit)##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## decay 0.71207 -0.70211 0.0568 0.2577
## canopy_cover 0.33335 0.94280 0.3386 0.0002 ***
## soil_moist -0.41180 -0.91127 0.0682 0.1939
## gc_shannon -0.05032 -0.99873 0.2390 0.0022 **
## shrubs_pct 0.94967 -0.31324 0.0498 0.2996
## fw_pct -0.57603 -0.81743 0.0234 0.5817
## dll_pct 0.35400 0.93524 0.2606 0.0009 ***
## lg_pct -0.32075 -0.94716 0.3544 0.0002 ***
## basal_area -0.28069 0.95980 0.1236 0.0494 *
## tree_dens -0.21294 0.97706 0.2395 0.0023 **
## off_ground 0.14719 0.98911 0.0028 0.9397
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 9999ef_tab <- data.frame(scores(env_fit, "vectors")) %>%
tibble::rownames_to_column("Variable") %>%
mutate(r2 = env_fit$vectors$r, p = env_fit$vectors$pvals) %>%
arrange(p)
write.csv(ef_tab, file.path(out_dir, "tables", "envfit_50sites.csv"),
row.names = FALSE)
nmds_df <- data.frame(scores(nmds, display = "sites")) %>%
tibble::rownames_to_column("Plot") %>%
left_join(env_m, by = "Plot")
sp_df <- data.frame(scores(nmds, display = "species")) %>%
tibble::rownames_to_column("Species")
vec_df <- ef_tab %>% filter(p < 0.05) %>%
mutate(NMDS1 = NMDS1 * ordiArrowMul(env_fit),
NMDS2 = NMDS2 * ordiArrowMul(env_fit))
p_nmds <- ggplot(nmds_df, aes(NMDS1, NMDS2)) +
geom_point(aes(colour = factor(decay), shape = Area), size = 3.2) +
stat_ellipse(aes(group = factor(decay), colour = factor(decay)),
type = "t", linewidth = 0.4, alpha = 0.6) +
geom_segment(data = vec_df, aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
arrow = arrow(length = unit(0.25, "cm")), colour = "grey20") +
geom_text_repel(data = vec_df, aes(NMDS1, NMDS2, label = Variable),
colour = "grey20", size = 3.2) +
scale_colour_viridis_d(name = "Decay class", option = "D", end = 0.9) +
annotate("text", x = Inf, y = -Inf, hjust = 1.1, vjust = -0.8,
label = paste0("Stress = ", round(nmds$stress, 3))) +
theme_classic(base_size = 12) +
labs(title = "Wildlife community on CWD by decay class")
ggsave(file.path(out_dir, "figures", "NMDS_decay_50sites.png"), p_nmds,
width = 8, height = 6, dpi = 300)
p_sp <- ggplot(nmds_df, aes(NMDS1, NMDS2)) +
geom_point(colour = "grey70", size = 2) +
geom_text_repel(data = sp_df, aes(NMDS1, NMDS2, label = Species),
size = 2.6, max.overlaps = 25) +
theme_classic(base_size = 12) +
labs(title = "Species scores")
ggsave(file.path(out_dir, "figures", "NMDS_species_50sites.png"), p_sp,
width = 8, height = 6, dpi = 300)Which species drive the decay-class signal?
env_m <- env_m %>%
mutate(decay_grp = cut(decay, breaks = c(0, 2, 3, 5),
labels = c("early (1-2)", "mid (3)", "late (4-5)")))
set.seed(SEED)
simp <- simper(comm_m, group = env_m$decay_grp, permutations = 999)
simp_sum <- summary(simp)
simp_tab <- bind_rows(lapply(names(simp_sum), function(nm) {
as.data.frame(simp_sum[[nm]]) %>%
tibble::rownames_to_column("Species") %>%
mutate(Contrast = nm) %>%
slice_head(n = 8)
}))
print(simp_tab %>% dplyr::select(Contrast, Species, average, cumsum, p))## Contrast Species average cumsum p
## 1 late (4-5)_mid (3) cotton_mouse 0.17043598 0.2666439 0.738
## 2 late (4-5)_mid (3) eastern_woodrat 0.08469374 0.3991456 0.300
## 3 late (4-5)_mid (3) nine_banded_armadillo 0.04500285 0.4695517 0.318
## 4 late (4-5)_mid (3) eastern_chipmunk 0.03928542 0.5310130 0.285
## 5 late (4-5)_mid (3) northern_cardinal 0.03427918 0.5846422 0.686
## 6 late (4-5)_mid (3) virginia_opossum 0.02958673 0.6309300 0.221
## 7 late (4-5)_mid (3) eastern_gray_squirrel 0.02546177 0.6707645 0.641
## 8 late (4-5)_mid (3) raccoon 0.02428011 0.7087503 0.550
## 9 late (4-5)_early (1-2) cotton_mouse 0.18183721 0.2768390 0.400
## 10 late (4-5)_early (1-2) eastern_woodrat 0.07346406 0.3886847 0.586
## 11 late (4-5)_early (1-2) nine_banded_armadillo 0.05230037 0.4683097 0.025
## 12 late (4-5)_early (1-2) northern_cardinal 0.03885471 0.5274643 0.350
## 13 late (4-5)_early (1-2) eastern_chipmunk 0.03772566 0.5848999 0.384
## 14 late (4-5)_early (1-2) eastern_gray_squirrel 0.02836911 0.6280906 0.447
## 15 late (4-5)_early (1-2) virginia_opossum 0.02592093 0.6675541 0.623
## 16 late (4-5)_early (1-2) hispid_cotton_rat 0.02425679 0.7044839 0.174
## 17 mid (3)_early (1-2) cotton_mouse 0.17854165 0.3024923 0.477
## 18 mid (3)_early (1-2) eastern_woodrat 0.06190704 0.4073776 0.774
## 19 mid (3)_early (1-2) northern_cardinal 0.03631250 0.4688997 0.545
## 20 mid (3)_early (1-2) raccoon 0.03059625 0.5207370 0.034
## 21 mid (3)_early (1-2) eastern_gray_squirrel 0.02782421 0.5678779 0.505
## 22 mid (3)_early (1-2) virginia_opossum 0.02552812 0.6111286 0.610
## 23 mid (3)_early (1-2) carolina_wren 0.02504546 0.6535616 0.075
## 24 mid (3)_early (1-2) eastern_chipmunk 0.02495538 0.6958420 0.923write.csv(simp_tab, file.path(out_dir, "tables", "SIMPER_decay_groups.csv"),
row.names = FALSE)
top_sp <- names(sort(colSums(comm_m), decreasing = TRUE))[1:10]
rate_by_decay <- as.data.frame(comm_m[, top_sp]) %>%
mutate(decay = env_m$decay) %>%
pivot_longer(-decay, names_to = "Species", values_to = "rate100") %>%
group_by(Species, decay) %>%
summarise(mean_rate = mean(rate100), .groups = "drop")
p_rate <- ggplot(rate_by_decay, aes(factor(decay), mean_rate,
fill = factor(decay))) +
geom_col() + facet_wrap(~ Species, scales = "free_y", ncol = 5) +
scale_fill_viridis_d(guide = "none", end = 0.9) +
theme_bw(base_size = 10) +
labs(x = "CWD decay class", y = "Mean detections / 100 camera-days",
title = "Ten most-detected species by decay class")
ggsave(file.path(out_dir, "figures", "species_by_decay.png"), p_rate,
width = 11, height = 5, dpi = 300)Sensitivity to distance measure and transformation
raw_matrix <- CamInd %>%
count(Plot, Common_Name, name = "detections") %>%
pivot_wider(names_from = Common_Name, values_from = detections,
values_fill = 0) %>%
arrange(Plot)
raw_m <- as.matrix(raw_matrix[, -1]); rownames(raw_m) <- raw_matrix$Plot
raw_m <- raw_m[env_m$Plot, colSums(raw_m) > 0, drop = FALSE]
sub25 <- env_m$n_detections >= 25
comm_25 <- comm_m[sub25, , drop = FALSE]
comm_25 <- comm_25[, colSums(comm_25) > 0, drop = FALSE]
sens <- list(
`rates, Bray-Curtis (main)` = list(d = vegdist(comm_m, "bray"), e = env_m),
`raw counts, Bray-Curtis` = list(d = vegdist(raw_m, "bray"), e = env_m),
`rates, Hellinger + Bray` = list(d = vegdist(decostand(comm_m, "hellinger"),
"bray"), e = env_m),
`presence/absence, Jaccard` = list(d = vegdist(comm_m, "jaccard", binary = TRUE),
e = env_m),
`logs with >= 25 detections` = list(d = vegdist(comm_25, "bray"),
e = env_m[sub25, ]))
sens_tab <- bind_rows(lapply(names(sens), function(nm) {
s <- sens[[nm]]
set.seed(SEED)
a <- adonis2(s$d ~ decay_f + canopy_cover + shrubs_pct + basal_area,
data = s$e,
permutations = how(blocks = factor(s$e$Area), nperm = 9999),
by = "margin")
as.data.frame(a) %>% tibble::rownames_to_column("Term") %>%
mutate(Version = nm, n_sites = nrow(s$e))
})) %>%
filter(!Term %in% c("Residual", "Total")) %>%
dplyr::select(Version, n_sites, Term, Df, R2, `F`, `Pr(>F)`) %>%
mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
as_tibble()
print(sens_tab, n = 40)## # A tibble: 20 × 7
## Version n_sites Term Df R2 F `Pr(>F)`
## <chr> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 rates, Bray-Curtis (main) 50 decay_f 4 0.0824 1.12 0.266
## 2 rates, Bray-Curtis (main) 50 canopy_cover 1 0.075 4.06 0.0008
## 3 rates, Bray-Curtis (main) 50 shrubs_pct 1 0.015 0.812 0.666
## 4 rates, Bray-Curtis (main) 50 basal_area 1 0.0316 1.71 0.113
## 5 raw counts, Bray-Curtis 50 decay_f 4 0.0832 1.14 0.24
## 6 raw counts, Bray-Curtis 50 canopy_cover 1 0.0738 4.03 0.0006
## 7 raw counts, Bray-Curtis 50 shrubs_pct 1 0.0175 0.953 0.520
## 8 raw counts, Bray-Curtis 50 basal_area 1 0.0302 1.65 0.134
## 9 rates, Hellinger + Bray 50 decay_f 4 0.104 1.45 0.0239
## 10 rates, Hellinger + Bray 50 canopy_cover 1 0.0814 4.56 0.0001
## 11 rates, Hellinger + Bray 50 shrubs_pct 1 0.016 0.898 0.621
## 12 rates, Hellinger + Bray 50 basal_area 1 0.0382 2.14 0.0526
## 13 presence/absence, Jaccard 50 decay_f 4 0.0973 1.30 0.0434
## 14 presence/absence, Jaccard 50 canopy_cover 1 0.07 3.73 0.0001
## 15 presence/absence, Jaccard 50 shrubs_pct 1 0.0194 1.04 0.491
## 16 presence/absence, Jaccard 50 basal_area 1 0.0225 1.20 0.387
## 17 logs with >= 25 detections 44 decay_f 4 0.085 1.05 0.357
## 18 logs with >= 25 detections 44 canopy_cover 1 0.0892 4.41 0.0004
## 19 logs with >= 25 detections 44 shrubs_pct 1 0.0249 1.23 0.307
## 20 logs with >= 25 detections 44 basal_area 1 0.0419 2.07 0.0434write.csv(sens_tab, file.path(out_dir, "tables", "sensitivity_PERMANOVA.csv"),
row.names = FALSE)Plant community on the 17 surveyed logs
Veg <- read_excel(veg_path, sheet = "Veg")
quad_n <- Veg %>% group_by(cwd) %>%
summarise(n_quad = n_distinct(q), .groups = "drop")
plant_cover <- Veg %>%
mutate(cover_pct = as.numeric(daub_mid(cover))) %>%
group_by(cwd, sp_name) %>%
summarise(total_cover = sum(cover_pct, na.rm = TRUE), .groups = "drop") %>%
left_join(quad_n, by = "cwd") %>%
mutate(avg_cover = total_cover / n_quad)
plant_matrix <- plant_cover %>%
dplyr::select(cwd, sp_name, avg_cover) %>%
pivot_wider(names_from = sp_name, values_from = avg_cover, values_fill = 0) %>%
arrange(cwd)
plant_div <- tibble(
Plot = plant_matrix$cwd,
plant_shannon = vegan::diversity(as.matrix(plant_matrix[, -1]), index = "shannon"),
plant_simpson = vegan::diversity(as.matrix(plant_matrix[, -1]), index = "simpson"),
plant_richness = specnumber(as.matrix(plant_matrix[, -1])),
plant_evenness = plant_shannon / log(plant_richness))
form_comp <- Veg %>%
mutate(cover_pct = as.numeric(daub_mid(cover))) %>%
group_by(cwd, type) %>%
summarise(cov = sum(cover_pct, na.rm = TRUE), .groups = "drop") %>%
group_by(cwd) %>% mutate(prop = cov / sum(cov)) %>% ungroup() %>%
filter(!is.na(type)) %>%
dplyr::select(cwd, type, prop) %>%
pivot_wider(names_from = type, values_from = prop, values_fill = 0,
names_prefix = "prop_") %>%
rename(Plot = cwd)
env17 <- env_m %>%
inner_join(plant_div, by = "Plot") %>%
left_join(form_comp, by = "Plot") %>%
arrange(Plot)
comm17 <- comm_m[rownames(comm_m) %in% env17$Plot, , drop = FALSE]
comm17 <- comm17[order(rownames(comm17)), , drop = FALSE]
comm17 <- comm17[, colSums(comm17) > 0, drop = FALSE]
stopifnot(identical(rownames(comm17), env17$Plot))
dist17 <- vegdist(comm17, method = "bray")
set.seed(SEED)
perm17 <- adonis2(dist17 ~ decay + plant_shannon + canopy_cover + soil_moist,
data = env17, permutations = 9999, by = "margin")
print(perm17)## Permutation test for adonis under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = dist17 ~ decay + plant_shannon + canopy_cover + soil_moist, data = env17, permutations = 9999, by = "margin")
## Df SumOfSqs R2 F Pr(>F)
## decay 1 0.2086 0.06478 1.1783 0.2879
## plant_shannon 1 0.1597 0.04959 0.9021 0.5152
## canopy_cover 1 0.3480 0.10807 1.9658 0.0434 *
## soil_moist 1 0.2720 0.08446 1.5364 0.1270
## Residual 12 2.1241 0.65970
## Total 16 3.2198 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(SEED)
perm17_seq <- adonis2(dist17 ~ decay + plant_shannon + plant_richness +
canopy_cover + prop_shrub,
data = env17, permutations = 9999, by = "terms")
print(perm17_seq)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = dist17 ~ decay + plant_shannon + plant_richness + canopy_cover + prop_shrub, data = env17, permutations = 9999, by = "terms")
## Df SumOfSqs R2 F Pr(>F)
## decay 1 0.2439 0.07576 1.4528 0.1491
## plant_shannon 1 0.1618 0.05025 0.9637 0.4644
## plant_richness 1 0.2365 0.07345 1.4085 0.1977
## canopy_cover 1 0.4675 0.14520 2.7845 0.0059 **
## prop_shrub 1 0.2632 0.08174 1.5676 0.1217
## Residual 11 1.8469 0.57360
## Total 16 3.2198 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1write.csv(bind_rows(tidy_adonis(perm17, "Marginal (17 CWD, plant data)"),
tidy_adonis(perm17_seq, "Sequential (17 CWD, plant data)")),
file.path(out_dir, "tables", "PERMANOVA_17sites.csv"), row.names = FALSE)
## Is the wildlife community tracking the plant community?
plant_m <- as.matrix(plant_matrix[plant_matrix$cwd %in% env17$Plot, -1])
rownames(plant_m) <- plant_matrix$cwd[plant_matrix$cwd %in% env17$Plot]
plant_m <- plant_m[order(rownames(plant_m)), , drop = FALSE]
dist_p <- vegdist(plant_m, method = "bray")
set.seed(SEED)
mant <- mantel(dist17, dist_p, method = "spearman", permutations = 9999)
print(mant)##
## Mantel statistic based on Spearman's rank correlation rho
##
## Call:
## mantel(xdis = dist17, ydis = dist_p, method = "spearman", permutations = 9999)
##
## Mantel statistic r: 0.2123
## Significance: 0.0389
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.155 0.198 0.234 0.275
## Permutation: free
## Number of permutations: 9999set.seed(SEED)
nmds17 <- metaMDS(comm17, distance = "bray", k = 2, trymax = 500,
autotransform = FALSE)## Run 0 stress 0.07418655
## Run 1 stress 0.09036692
## Run 2 stress 0.07418655
## ... Procrustes: rmse 1.464447e-05 max resid 3.403129e-05
## ... Similar to previous best
## Run 3 stress 0.1009792
## Run 4 stress 0.09103019
## Run 5 stress 0.07418656
## ... Procrustes: rmse 7.66671e-05 max resid 0.000173776
## ... Similar to previous best
## Run 6 stress 0.07418661
## ... Procrustes: rmse 0.0001563887 max resid 0.0003524506
## ... Similar to previous best
## Run 7 stress 0.09036692
## Run 8 stress 0.0741866
## ... Procrustes: rmse 0.0001505963 max resid 0.0003408466
## ... Similar to previous best
## Run 9 stress 0.07418655
## ... Procrustes: rmse 1.86911e-05 max resid 4.327929e-05
## ... Similar to previous best
## Run 10 stress 0.07418672
## ... Procrustes: rmse 0.0002346687 max resid 0.0005267354
## ... Similar to previous best
## Run 11 stress 0.2145652
## Run 12 stress 0.07418655
## ... Procrustes: rmse 5.630663e-05 max resid 0.0001275166
## ... Similar to previous best
## Run 13 stress 0.1636204
## Run 14 stress 0.07418658
## ... Procrustes: rmse 0.0001261246 max resid 0.0002842699
## ... Similar to previous best
## Run 15 stress 0.100979
## Run 16 stress 0.07418655
## ... Procrustes: rmse 1.877753e-05 max resid 3.653928e-05
## ... Similar to previous best
## Run 17 stress 0.07418663
## ... Procrustes: rmse 0.0001178893 max resid 0.0002683701
## ... Similar to previous best
## Run 18 stress 0.07418654
## ... New best solution
## ... Procrustes: rmse 1.95027e-05 max resid 4.20654e-05
## ... Similar to previous best
## Run 19 stress 0.07418655
## ... Procrustes: rmse 3.257045e-05 max resid 7.822157e-05
## ... Similar to previous best
## Run 20 stress 0.07418656
## ... Procrustes: rmse 1.930147e-05 max resid 4.679937e-05
## ... Similar to previous best
## *** Best solution repeated 3 timesset.seed(SEED)
nmdsP <- metaMDS(plant_m, distance = "bray", k = 2, trymax = 500,
autotransform = FALSE)## Run 0 stress 0.1404698
## Run 1 stress 0.1984314
## Run 2 stress 0.1404698
## ... New best solution
## ... Procrustes: rmse 1.313477e-06 max resid 3.203567e-06
## ... Similar to previous best
## Run 3 stress 0.1765051
## Run 4 stress 0.1404698
## ... Procrustes: rmse 1.416975e-06 max resid 3.683143e-06
## ... Similar to previous best
## Run 5 stress 0.1732316
## Run 6 stress 0.1404698
## ... Procrustes: rmse 4.701897e-06 max resid 1.111903e-05
## ... Similar to previous best
## Run 7 stress 0.1815212
## Run 8 stress 0.1494951
## Run 9 stress 0.229931
## Run 10 stress 0.1404698
## ... Procrustes: rmse 3.502418e-06 max resid 1.008224e-05
## ... Similar to previous best
## Run 11 stress 0.1494951
## Run 12 stress 0.1732316
## Run 13 stress 0.177116
## Run 14 stress 0.2530365
## Run 15 stress 0.1494951
## Run 16 stress 0.2123976
## Run 17 stress 0.1404698
## ... New best solution
## ... Procrustes: rmse 7.526193e-07 max resid 1.541177e-06
## ... Similar to previous best
## Run 18 stress 0.1494951
## Run 19 stress 0.1404698
## ... New best solution
## ... Procrustes: rmse 1.123276e-06 max resid 3.154582e-06
## ... Similar to previous best
## Run 20 stress 0.177116
## *** Best solution repeated 1 timesset.seed(SEED)
pro <- protest(nmds17, nmdsP, permutations = 9999)
print(pro)##
## Call:
## protest(X = nmds17, Y = nmdsP, permutations = 9999)
##
## Procrustes Sum of Squares (m12 squared): 0.7151
## Correlation in a symmetric Procrustes rotation: 0.5337
## Significance: 0.0113
##
## Permutation: free
## Number of permutations: 9999set.seed(SEED)
ef17 <- envfit(nmds17 ~ decay + plant_shannon + plant_richness + canopy_cover +
soil_moist + prop_shrub + prop_grass,
data = env17, permutations = 9999, na.rm = TRUE)
print(ef17)##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## decay 0.69922 0.71491 0.0687 0.6601
## plant_shannon -0.97210 -0.23459 0.0554 0.6792
## plant_richness -0.93338 0.35890 0.0818 0.5765
## canopy_cover -0.96628 0.25748 0.1761 0.2604
## soil_moist 0.60188 0.79859 0.5484 0.0023 **
## prop_shrub -0.83729 -0.54676 0.3172 0.0876 .
## prop_grass 0.95711 0.28974 0.4721 0.0062 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 9999nmds17_df <- data.frame(scores(nmds17, display = "sites")) %>%
tibble::rownames_to_column("Plot") %>% left_join(env17, by = "Plot")
ef17_tab <- data.frame(scores(ef17, "vectors")) %>%
tibble::rownames_to_column("Variable") %>%
mutate(r2 = ef17$vectors$r, p = ef17$vectors$pvals)
p17 <- ggplot(nmds17_df, aes(NMDS1, NMDS2)) +
geom_point(aes(colour = factor(decay), size = plant_shannon)) +
geom_text_repel(aes(label = Plot), size = 3) +
geom_segment(data = ef17_tab %>% filter(p < 0.1),
aes(x = 0, y = 0, xend = NMDS1 * ordiArrowMul(ef17),
yend = NMDS2 * ordiArrowMul(ef17)),
arrow = arrow(length = unit(0.25, "cm")), colour = "grey20") +
geom_text_repel(data = ef17_tab %>% filter(p < 0.1),
aes(NMDS1 * ordiArrowMul(ef17), NMDS2 * ordiArrowMul(ef17),
label = Variable), colour = "grey20", size = 3) +
scale_colour_viridis_d(name = "Decay class", end = 0.9) +
scale_size_continuous(name = "Plant H'") +
annotate("text", x = Inf, y = -Inf, hjust = 1.1, vjust = -0.8,
label = paste0("Stress = ", round(nmds17$stress, 3))) +
theme_classic(base_size = 12) +
labs(title = "Wildlife community vs. plant diversity (17 CWD with plant surveys)")
ggsave(file.path(out_dir, "figures", "NMDS_17sites_plants.png"), p17,
width = 8, height = 6, dpi = 300)Export
saveRDS(list(stress_grid = all_stress_results,
procrustes = procrustes_tab,
permanova_by_rarity = permanova_by_rarity,
final_nmds = nmds,
chosen_k = chosen_k,
chosen_rarity = chosen_rarity,
h1_tests = h1_tests,
perm_main = perm_tables,
envfit_tab = ef_tab,
seed = SEED),
file.path(out_dir, "tables", "community_results.rds"))
writexl::write_xlsx(
list(
independent_obs = CamInd %>%
dplyr::select(Plot, VideoNumber, DateTime, Common_Name, Class, Order,
Family, Genus, Species, Behavior, Location, N_Ind,
Camera.Type, Unique),
effort = effort,
wildlife_matrix = wild_matrix,
site_env = env_m,
site_env_17 = env17,
plant_matrix = plant_matrix,
clock_errors = clock_errors,
species_freq = sp_freq
),
file.path(proj, "03_Output", "CWD_analysis_clean_data.xlsx"))
cat("\nDone. Seed =", SEED, "| outputs in", out_dir, "\n")##
## Done. Seed = 97 | outputs in C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/13_CWD_Wildlife/03_Output/community