2. Data Import and Cleaning
Step 1 Load the data
read_excel() reads a named sheet from an
.xlsx file into a data frame — R’s
standard rectangular data structure (rows = observations, columns =
variables). Update file_path to match your own file
location.
file_path <- "/home/mwl04747/RTricks/A16_Camera_Traps/Camera Trap Data Entry.xlsx"
cam_raw <- read_excel(path = file_path, sheet = "Combined")
glimpse(cam_raw)
## Rows: 1,197
## Columns: 6
## $ Camera_ID <chr> "P70-03", "P70-03", "P70-03"…
## $ `DateTimeStamp (YYYY-MM-DD HH:MM:SS PST)` <chr> "46093.241377314815", "46093…
## $ Datetime_corrected <dttm> 2026-03-12 05:47:35, 2026-0…
## $ Taxa_code <chr> "C_latrans", "C_latrans", "C…
## $ Behavior_code <chr> "walk", "hop", "walk", "snif…
## $ Notes <chr> NA, NA, "exiting burrow", NA…
💡 glimpse() is your first-look friend.
It shows every column name, its data type (<chr> =
text, <dbl> = number, <dttm> =
date-time), and the first few values — far more useful than
head() for wide data frames.
Step 2 Standardize column names and
parse dates
Column names from Excel often contain spaces and mixed
capitalisation. clean_names() converts them to
snake_case (lower case, underscores instead of spaces) so
they are easier to type and less error-prone.
cam_raw$datetime <- cam_raw$Datetime_corrected
cam <- cam_raw %>%
clean_names() %>%
mutate(
date = as.Date(datetime), # strip the time component; keep date only
hour = hour(datetime) # extract hour of day as integer (0–23)
)
summary(cam$date)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## "2026-03-09" "2026-03-12" "2026-03-15" "2026-03-15" "2026-03-18" "2026-03-23"
## NA's
## "1"
Step 3 Remove ambiguous taxa
Some records contain notes like "bird",
"unknown", or combined entries
("squirrel/rat"). These are genuinely unidentifiable and
must be excluded before analysis — their presence would inflate richness
estimates and contaminate behavior summaries.
ambiguous_patterns <- c(
" and ", " or ", "/",
"unknown", "unidentified", "n/a", "^na$",
"bird$", "small brown bird"
)
cam_clean <- cam %>%
filter(
!str_detect(
str_to_lower(taxa_code),
paste(ambiguous_patterns, collapse = "|")
)
)
cat("Taxa retained after filtering:\n")
## Taxa retained after filtering:
unique(cam_clean$taxa_code)
## [1] "C_latrans" "A_californica" "O_beecheyi"
## [4] "S_bachmani" "M_crissalis" "C_cristata"
## [7] "S_niger" "D_virginiana" "M_musculus"
## [10] "Black-headed Grosbeak" "R_norvegicus" "M_mephitis"
📌 Regular expressions
str_detect() matches text against a
pattern. paste(…, collapse = “|”) joins
our list of ambiguous strings with the “or” operator (|),
so one call checks all of them simultaneously. The ! at the
front keeps only rows that do not match any ambiguous
pattern.
Step 4 Collapse behaviors into
functional types
Raw behavior codes ("walk", "run",
"hop") are too granular for most analyses.
case_when() groups them into five functional
categories — think of it as a multi-way if/else evaluated in
order.
cam <- cam_clean %>%
mutate(
behavior_function = case_when(
str_detect(behavior_code, "walk|run|hop|fly|climb|crawl") ~ "movement",
str_detect(behavior_code, "eat|dig|cache") ~ "foraging",
str_detect(behavior_code, "sniff") ~ "investigation",
str_detect(behavior_code, "stare") ~ "vigilance",
str_detect(behavior_code, "stop|stand|sit|lay") ~ "stationary",
TRUE ~ "other"
)
)
table(cam$behavior_function)
##
## foraging investigation movement other stationary
## 106 76 754 43 101
## vigilance
## 67
3. Defining Independent Detection Events
The pseudo-replication problem
A single animal may trigger the same camera dozens of times over a
few minutes. Treating each trigger as an independent observation would
pseudo-replicate the data — inflating sample sizes and
biasing every statistic downstream.
The 30-minute rule
The standard solution: any detections of the same species at
the same camera within 30 minutes of the previous detection are
treated as part of the same event. Only a detection after a ≥
30-minute gap counts as a new independent event.
cam_events <- cam %>%
arrange(camera_id, taxa_code, datetime) %>%
group_by(camera_id, taxa_code) %>%
mutate(
# Time gap in minutes between consecutive rows within each camera × species group
minutes_since_last = as.numeric(
difftime(datetime, lag(datetime), units = "mins")
),
# A new event begins if gap > 30 min OR this is the very first record
new_event = if_else(
is.na(minutes_since_last) | minutes_since_last > 30, 1L, 0L
),
# Cumulative sum gives a running event counter within each group
event_id = cumsum(new_event)
) %>%
ungroup() %>%
mutate(diel = if_else(hour >= 6 & hour < 18, "day", "night"))
# How many independent events per species?
cam_events %>%
distinct(camera_id, taxa_code, event_id) %>%
count(taxa_code, name = "independent_events") %>%
arrange(desc(independent_events))
## # A tibble: 12 × 2
## taxa_code independent_events
## <chr> <int>
## 1 O_beecheyi 179
## 2 S_bachmani 141
## 3 C_latrans 126
## 4 M_crissalis 32
## 5 S_niger 31
## 6 A_californica 21
## 7 M_musculus 15
## 8 C_cristata 3
## 9 Black-headed Grosbeak 1
## 10 D_virginiana 1
## 11 M_mephitis 1
## 12 R_norvegicus 1
📌 Why lag() is the right tool here
lag(datetime) returns the value from the previous
row. Because we first group_by(camera_id, taxa_code),
“previous” means the last detection of that species at that specific
camera — exactly what we need to measure within-burst gaps. Without
grouping, lag() would compare across species and cameras,
giving meaningless gaps.
⚠️ Choosing a time threshold
30 minutes is a widely used convention, but the right threshold
depends on your study organism’s typical bout length and movement speed.
For a chipmunk foraging in a small area, 30 minutes may be too long. For
a wide-ranging carnivore like a mountain lion, 60–120 minutes may be
more appropriate. Whatever threshold you choose, document it clearly in
your methods.
4. Descriptive Statistics
4.1 Detection event totals per taxon
event_summary <- cam_events %>%
distinct(camera_id, taxa_code, event_id) %>%
count(taxa_code, name = "detection_events") %>%
arrange(desc(detection_events))
print(event_summary)
## # A tibble: 12 × 2
## taxa_code detection_events
## <chr> <int>
## 1 O_beecheyi 179
## 2 S_bachmani 141
## 3 C_latrans 126
## 4 M_crissalis 32
## 5 S_niger 31
## 6 A_californica 21
## 7 M_musculus 15
## 8 C_cristata 3
## 9 Black-headed Grosbeak 1
## 10 D_virginiana 1
## 11 M_mephitis 1
## 12 R_norvegicus 1
ggplot(event_summary,
aes(x = reorder(taxa_code, detection_events),
y = detection_events,
fill = reorder(taxa_code, detection_events))) +
geom_col(color = "#1e3a2f", alpha = 0.88) +
geom_text(aes(label = detection_events), hjust = -0.2,
family = "mono", size = 3.4, color = "#1e3a2f") +
coord_flip() +
scale_fill_manual(
values = colorRampPalette(c("#d4e6da", "#1e3a2f"))(nrow(event_summary))
) +
scale_y_continuous(expand = expansion(mult = c(0, 0.12))) +
labs(title = "Independent detection events by taxon",
x = NULL, y = "Detection events") +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(family = "serif", face = "bold", color = "#1e3a2f"),
legend.position = "none",
panel.grid.major.y = element_blank()
)
4.2 Behavioral ethogram (stacked bar)
An ethogram documents the relative frequency of each
behavior. Proportional bars enable cross-species comparison even when
total event counts differ substantially between species.
ethogram <- cam_events %>%
count(taxa_code, behavior_function) %>%
group_by(taxa_code) %>%
mutate(prop = n / sum(n)) %>%
ungroup()
ggplot(ethogram, aes(x = taxa_code, y = prop, fill = behavior_function)) +
geom_col(color = "white", linewidth = 0.3) +
scale_fill_manual(
values = c(
"movement" = "#4a7c5e",
"foraging" = "#c8832a",
"investigation" = "#2e7d94",
"vigilance" = "#8aad96",
"stationary" = "#d4e6da",
"other" = "#7a7a72"
)
) +
scale_y_continuous(labels = scales::percent_format()) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(family = "serif", face = "bold", color = "#1e3a2f"),
axis.text.x = element_text(angle = 40, hjust = 1),
legend.title = element_blank()
) +
labs(title = "Behavioral composition by taxon",
x = NULL, y = "Proportion of observations")
4.3 Hourly activity heat map
A heat map encodes counts as color intensity. Each
cell is a taxon × hour combination — a fast way to spot nocturnal,
diurnal, or crepuscular patterns across all taxa simultaneously.
taxon_hour <- cam_events %>%
distinct(taxa_code, event_id, hour) %>%
count(taxa_code, hour, name = "events")
ggplot(taxon_hour, aes(x = hour, y = taxa_code, fill = events)) +
geom_tile(color = "white", linewidth = 0.4) +
scale_fill_gradientn(
colors = c(var_cream = "#f7f3ed", "#d4e6da", "#4a7c5e", "#1e3a2f"),
name = "Events"
) +
scale_x_continuous(breaks = seq(0, 23, 3),
labels = paste0(seq(0, 23, 3), ":00")) +
labs(title = "Hourly detection patterns by taxon",
x = "Hour of day", y = NULL) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(family = "serif", face = "bold", color = "#1e3a2f"),
axis.text.x = element_text(size = 9)
)
4.4 Camera-level activity profiles
Different cameras are placed in different micro-habitats. This plot
reveals whether behavioral composition varies by camera location — an
early signal of habitat-driven behavior differences.
camera_behavior <- cam_events %>%
distinct(camera_id, taxa_code, event_id, behavior_function) %>%
count(camera_id, behavior_function, name = "events")
ggplot(camera_behavior,
aes(x = camera_id, y = events, fill = behavior_function)) +
geom_col(color = "white", linewidth = 0.3) +
scale_fill_manual(
values = c(
"movement" = "#4a7c5e",
"foraging" = "#c8832a",
"investigation" = "#2e7d94",
"vigilance" = "#8aad96",
"stationary" = "#d4e6da",
"other" = "#7a7a72"
)
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(family = "serif", face = "bold", color = "#1e3a2f"),
axis.text.x = element_text(angle = 40, hjust = 1),
legend.title = element_blank()
) +
labs(title = "Camera activity profiles by behavior function",
x = "Camera ID", y = "Detection events")
5. Hypothesis Testing: Behavior × Diel Period
Rationale
Question: Does the proportion of behaviors
(e.g., foraging vs. vigilance) differ between day and night?
We use a chi-square test of independence — or
Fisher’s exact test when expected cell counts are small — to test
whether behavior type is statistically independent of diel period.
H₀: Behavior composition is the same during the day
and at night.
H₁: At least one behavior type is
more common at one time of day.
Run tests per taxon
test_behavior_diel <- function(df) {
# 1. Drop rows where diel is NA (records with hour == NA or bad datetime)
df <- df %>% filter(!is.na(diel), !is.na(behavior_function))
# 2. Guard: need observations in BOTH diel periods
if (length(unique(df$diel)) < 2)
return(tibble(test = "insufficient data — one diel period only",
stat_type = "—", statistic = NA_real_, p_value = NA_real_))
# 3. Build the contingency table
wide <- df %>%
count(diel, behavior_function) %>%
pivot_wider(names_from = behavior_function,
values_from = n,
values_fill = 0)
# 4. Guard: diel column must have no NAs before using as row names
wide <- wide %>% filter(!is.na(diel))
if (nrow(wide) < 2)
return(tibble(test = "insufficient data — pivot produced < 2 rows",
stat_type = "—", statistic = NA_real_, p_value = NA_real_))
tbl <- wide %>%
column_to_rownames("diel") %>%
as.matrix()
# 5. Guard: need at least 2 behavior columns
if (ncol(tbl) < 2)
return(tibble(test = "insufficient data — fewer than 2 behavior types",
stat_type = "—", statistic = NA_real_, p_value = NA_real_))
# 6. Choose chi-square vs. Fisher based on minimum expected count
exp_counts <- suppressWarnings(chisq.test(tbl)$expected)
if (any(exp_counts < 5)) {
# Fisher's exact test has no single test statistic (no chi-square value).
# For 2x2 tables it returns an odds ratio; for larger tables use NA.
res <- fisher.test(tbl, simulate.p.value = TRUE, B = 2000)
or_val <- if (!is.null(res$estimate)) round(as.numeric(res$estimate), 3)
else NA_real_
tibble(test = "Fisher",
statistic = or_val, # odds ratio when 2x2; NA for larger tables
stat_type = if (!is.na(or_val)) "odds ratio" else "—",
p_value = res$p.value)
} else {
res <- chisq.test(tbl)
tibble(test = "Chi-square",
statistic = round(as.numeric(res$statistic), 3),
stat_type = "X²",
p_value = res$p.value)
}
}
behavior_tests <- cam_events %>%
group_by(taxa_code) %>%
group_modify(~ test_behavior_diel(.x)) %>%
ungroup() %>%
mutate(
significant = ifelse(is.na(p_value), NA, p_value < 0.05),
p_value = round(p_value, 4)
)
# Re-order columns for a readable table
behavior_tests <- behavior_tests %>%
select(taxa_code, test, stat_type, statistic, p_value, significant)
knitr::kable(
behavior_tests,
col.names = c("Taxon", "Test used", "Statistic type",
"Value", "p-value", "p < 0.05?"),
caption = paste0(
"Diel behavior tests by taxon. ",
"For chi-square the value is X²; ",
"for Fisher it is the odds ratio (2×2 tables only; — for larger tables). ",
"Fisher\'s exact test has no chi-square statistic — ",
"NA in earlier versions was correct but misleadingly labeled."
)
)
Diel behavior tests by taxon. For chi-square the value is X²;
for Fisher it is the odds ratio (2×2 tables only; — for larger tables).
Fisher’s exact test has no chi-square statistic — NA in earlier versions
was correct but misleadingly labeled.
| A_californica |
Fisher |
— |
NA |
0.0480 |
TRUE |
| Black-headed Grosbeak |
insufficient data — one diel period only |
— |
NA |
NA |
NA |
| C_cristata |
insufficient data — one diel period only |
— |
NA |
NA |
NA |
| C_latrans |
Fisher |
— |
NA |
0.3653 |
FALSE |
| D_virginiana |
insufficient data — one diel period only |
— |
NA |
NA |
NA |
| M_crissalis |
Fisher |
— |
NA |
0.0240 |
TRUE |
| M_mephitis |
insufficient data — one diel period only |
— |
NA |
NA |
NA |
| M_musculus |
Fisher |
— |
NA |
0.1139 |
FALSE |
| O_beecheyi |
Fisher |
— |
NA |
0.0040 |
TRUE |
| R_norvegicus |
insufficient data — one diel period only |
— |
NA |
NA |
NA |
| S_bachmani |
Fisher |
— |
NA |
0.0190 |
TRUE |
| S_niger |
Fisher |
— |
NA |
0.0445 |
TRUE |
Visualize the diel behavior shift
behavior_diel <- cam_events %>%
# Remove records with no diel assignment or unclassified behavior
filter(!is.na(diel), !is.na(behavior_function)) %>%
count(taxa_code, diel, behavior_function) %>%
group_by(taxa_code, diel) %>%
mutate(prop = n / sum(n)) %>%
ungroup() %>%
# Convert diel to a factor with clean labels so NA never appears as a panel
mutate(
diel = factor(diel, levels = c("day", "night"),
labels = c("Day", "Night"))
)
# Only keep taxa that have observations in BOTH diel periods
both_periods <- behavior_diel %>%
distinct(taxa_code, diel) %>%
count(taxa_code) %>%
filter(n == 2) %>%
pull(taxa_code)
behavior_diel <- behavior_diel %>% filter(taxa_code %in% both_periods)
ggplot(behavior_diel,
aes(x = diel, y = prop, fill = behavior_function)) +
geom_col(color = "white", linewidth = 0.3) +
facet_wrap(~ taxa_code, ncol = 3) +
scale_fill_manual(
values = c(
"movement" = "#4a7c5e",
"foraging" = "#c8832a",
"investigation" = "#2e7d94",
"vigilance" = "#8aad96",
"stationary" = "#d4e6da",
"other" = "#7a7a72"
)
) +
scale_y_continuous(labels = scales::percent_format()) +
labs(title = "Behavioral composition: day vs. night",
x = NULL, y = "Proportion", fill = "Behavior") +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(family = "serif", face = "bold", color = "#1e3a2f"),
strip.text = element_text(family = "mono", face = "bold", color = "#2d5a3d"),
legend.title = element_blank(),
panel.spacing = unit(0.8, "lines")
)
💡 Interpreting a significant test A significant
p-value tells you that the behavior mix differs between day and
night, but not which behaviors drive the difference. Use the
faceted stacked bars above to identify which functional types shift
most. The chi-square standardized residuals (available
in the test object as chisq.test(tbl)$stdres) also pinpoint
which cells deviate most from independence.
6. Activity Overlap (Kernel Density)
Rationale
Two species active at the same hours are more likely to interact — as
competitors, predators, or prey. The overlap coefficient
Δ (Ridout & Linkie 2009) quantifies the similarity of two
diel activity curves on a 0–1 scale (0 = no temporal overlap, 1 =
identical schedules).
The overlap package fits smooth kernel
density curves to detection times treated as a circular
variable — because 23:00 and 01:00 are only 2 hours apart, not 22.
Why radians? Putting time on a circle
A standard 24-hour clock has an awkward edge: 23:45 and 00:15 are
only 30 minutes apart in reality, but subtracting their decimal hours
gives 23.5 — a huge apparent gap. The kernel density estimator would
treat those two detections as nearly opposite ends of the day.
The solution is to map the 24-hour cycle onto a
circle where midnight wraps seamlessly back to midnight. We do
this by converting each detection time to radians (0 to
2π), the natural unit for circular data.
The conversion has three steps:
- Decimal hours — combine the hour and fractional
minutes into a single number:
hour + minute/60. A detection
at 14:45 becomes 14.75.
- Fraction of the day — divide by 24 to get a value
between 0 and 1. 14.75 / 24 = 0.615.
- Scale to a full circle — multiply by 2π (≈ 6.283).
0.615 × 2π ≈ 3.86 radians.
cam_events <- cam_events %>%
mutate(
# Step 1: decimal hours → Step 2: fraction of day → Step 3: radians
time_rad = (hour(datetime) + minute(datetime) / 60) / 24 * 2 * pi
)
# Sanity check: values should span 0 to ~2π (6.28)
cat(sprintf("time_rad range: %.3f to %.3f (expected 0 to %.3f)\n",
min(cam_events$time_rad, na.rm = TRUE),
max(cam_events$time_rad, na.rm = TRUE),
2 * pi))
## time_rad range: 0.105 to 6.165 (expected 0 to 6.283)
Single pair: estimate Δ and bootstrap CI
# ── Minimum detections needed for bandwidth estimation ────────────────────────
MIN_N_OVERLAP <- 10 # raise further if bandwidth errors persist
# Helper: attempt overlap analysis, return NULL gracefully on failure
run_overlap <- function(ta, tb, sp1, sp2) {
if (length(ta) < MIN_N_OVERLAP || length(tb) < MIN_N_OVERLAP) {
message(sprintf("Skipping %s vs %s: too few detections (%d, %d).",
sp1, sp2, length(ta), length(tb)))
return(NULL)
}
if (length(unique(ta)) < 5 || length(unique(tb)) < 5) {
message(sprintf("Skipping %s vs %s: too few unique time values.", sp1, sp2))
return(NULL)
}
estimator <- if (min(length(ta), length(tb)) < 75) "Dhat4" else "Dhat1"
tryCatch({
delta_val <- overlapEst(ta, tb, type = estimator)[estimator]
set.seed(42)
bs_val <- bootstrap(ta, tb, nb = 999, type = estimator)
ci_val <- bootCI(delta_val, bs_val)
list(estimator = estimator, delta = delta_val,
ci_lower = ci_val["norm0", "lower"],
ci_upper = ci_val["norm0", "upper"])
}, error = function(e) {
message(sprintf("Overlap failed for %s vs %s: %s", sp1, sp2, conditionMessage(e)))
NULL
})
}
# ── Pick the two species with the most independent detections ─────────────────
top2 <- cam_events %>%
distinct(taxa_code, event_id) %>%
count(taxa_code, sort = TRUE) %>%
slice_head(n = 2) %>%
pull(taxa_code)
sp1_code <- top2[1]
sp2_code <- top2[2]
t1 <- cam_events %>% filter(taxa_code == sp1_code) %>% pull(time_rad)
t2 <- cam_events %>% filter(taxa_code == sp2_code) %>% pull(time_rad)
ov <- run_overlap(t1, t2, sp1_code, sp2_code)
if (!is.null(ov)) {
cat(sprintf("Activity overlap (%s): %s vs. %s \u2192 \u0394 = %.3f\n",
ov$estimator, sp1_code, sp2_code, ov$delta))
cat(sprintf("Bootstrap 95%% CI: [%.3f, %.3f]\n", ov$ci_lower, ov$ci_upper))
overlapPlot(
t1, t2,
rug = TRUE,
main = paste("Activity overlap:", sp1_code, "vs.", sp2_code),
xlab = "Time of day",
linecol = "#4a7c5e",
linecol2 = "#c8832a",
linewidth = 2,
linewidth2 = 2
)
legend("topright", legend = c(sp1_code, sp2_code),
col = c("#4a7c5e", "#c8832a"), lwd = 2, bty = "n")
} else {
# Fallback: plain hourly histograms when KDE cannot run
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(t1 / (2*pi) * 24, breaks = 24, col = adjustcolor("#4a7c5e", 0.7),
border = "white", xlab = "Hour of day", main = sp1_code, xlim = c(0,24))
hist(t2 / (2*pi) * 24, breaks = 24, col = adjustcolor("#c8832a", 0.7),
border = "white", xlab = "Hour of day", main = sp2_code, xlim = c(0,24))
par(mfrow = c(1, 1))
message("KDE overlap unavailable — showing hourly histograms instead.")
}
## Activity overlap (Dhat1): O_beecheyi vs. S_bachmani → Δ = 0.627
## Bootstrap 95% CI: [0.565, 0.690]
📌 Interpreting Δ
Δ > 0.75 = high overlap; 0.50–0.75 = moderate; < 0.50 = low. If
the bootstrap CI excludes a reference value of interest (e.g., Δ = 0.5),
the difference is statistically supported. See Ridout & Linkie
(2009) for formal guidance on estimator choice and CI
interpretation.
All pairwise overlaps (heat map)
taxa_list <- unique(cam_events$taxa_code)
n_sp <- length(taxa_list)
time_by_sp <- setNames(
lapply(taxa_list, function(sp)
cam_events %>% filter(taxa_code == sp) %>% pull(time_rad)),
taxa_list
)
# Helper: safe single-pair overlap (returns NA on any failure)
safe_overlap <- function(ta, tb) {
if (length(ta) < MIN_N_OVERLAP || length(tb) < MIN_N_OVERLAP) return(NA_real_)
if (length(unique(ta)) < 5 || length(unique(tb)) < 5) return(NA_real_)
est <- if (min(length(ta), length(tb)) < 75) "Dhat4" else "Dhat1"
tryCatch(
overlapEst(ta, tb, type = est)[est],
error = function(e) NA_real_
)
}
# Fill the symmetric overlap matrix, computing each pair only once
overlap_mat <- matrix(NA, n_sp, n_sp,
dimnames = list(taxa_list, taxa_list))
for (i in seq_len(n_sp)) {
for (j in seq_len(n_sp)) {
if (i == j) {
overlap_mat[i, j] <- 1
} else if (!is.na(overlap_mat[j, i])) {
overlap_mat[i, j] <- overlap_mat[j, i] # mirror already computed
} else {
overlap_mat[i, j] <- safe_overlap(time_by_sp[[i]], time_by_sp[[j]])
}
}
}
overlap_df <- as.data.frame(overlap_mat) %>%
rownames_to_column("sp1") %>%
pivot_longer(-sp1, names_to = "sp2", values_to = "delta")
ggplot(overlap_df, aes(x = sp1, y = sp2, fill = delta)) +
geom_tile(color = "white", linewidth = 0.5) +
geom_text(aes(label = round(delta, 2)), size = 3,
color = ifelse(overlap_df$delta > 0.6, "white", "#1a1a18")) +
scale_fill_gradientn(
colors = c("#f7f3ed", "#d4e6da", "#4a7c5e", "#1e3a2f"),
limits = c(0, 1), name = "Δ"
) +
labs(title = "Pairwise activity overlap (Δ)", x = NULL, y = NULL) +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(family = "serif", face = "bold", color = "#1e3a2f"),
axis.text.x = element_text(angle = 40, hjust = 1)
)
7. Species Diversity and Rarefaction
Rationale
Some cameras may be in richer habitats and detect more species. But
cameras also vary in effort — a busier camera
accumulates more detections by chance, inflating apparent diversity.
Rarefaction corrects for this by asking: how many
species would we expect if all cameras had the same total event
count?
Build the site × species matrix
# One row per camera, one column per taxon, cell = independent event count
site_sp <- cam_events %>%
distinct(camera_id, taxa_code, event_id) %>%
count(camera_id, taxa_code, name = "events") %>%
pivot_wider(names_from = taxa_code,
values_from = events,
values_fill = 0) %>%
column_to_rownames("camera_id")
print(site_sp)
## A_californica Black-headed Grosbeak C_latrans M_crissalis M_musculus
## P70-01 8 1 3 18 14
## P70-02 1 0 17 11 0
## P70-03 5 0 7 0 0
## P70-04 0 0 38 0 0
## T005 2 0 7 0 0
## T20-01 0 0 0 0 0
## T20-03 0 0 27 2 0
## T20-04 2 0 19 0 0
## TC22-01 1 0 1 0 1
## TC29-01 2 0 7 1 0
## O_beecheyi R_norvegicus S_bachmani S_niger M_mephitis C_cristata
## P70-01 11 1 24 2 0 0
## P70-02 59 0 52 0 0 0
## P70-03 1 0 1 0 0 0
## P70-04 1 0 0 2 0 0
## T005 5 0 5 0 0 0
## T20-01 20 0 0 0 1 0
## T20-03 4 0 26 0 0 3
## T20-04 27 0 12 5 0 0
## TC22-01 22 0 0 3 0 0
## TC29-01 29 0 21 19 0 0
## D_virginiana
## P70-01 0
## P70-02 0
## P70-03 0
## P70-04 0
## T005 0
## T20-01 0
## T20-03 0
## T20-04 0
## TC22-01 0
## TC29-01 1
Diversity indices
div_df <- tibble(
camera_id = rownames(site_sp),
richness = specnumber(site_sp),
shannon = diversity(site_sp, "shannon"),
simpson = diversity(site_sp, "simpson"),
total_events = rowSums(site_sp)
) %>%
arrange(desc(shannon))
knitr::kable(div_df,
digits = 3,
col.names = c("Camera", "Richness (S)", "Shannon (H′)",
"Simpson (1−D)", "Total events"),
caption = "Diversity indices by camera. Rarefied richness (below) controls for effort differences.")
Diversity indices by camera. Rarefied richness (below) controls
for effort differences.
| P70-01 |
9 |
1.810 |
0.807 |
82 |
| TC29-01 |
7 |
1.475 |
0.735 |
80 |
| T20-04 |
5 |
1.341 |
0.701 |
65 |
| T005 |
4 |
1.307 |
0.715 |
19 |
| P70-02 |
5 |
1.223 |
0.663 |
140 |
| T20-03 |
5 |
1.161 |
0.627 |
62 |
| P70-03 |
4 |
1.091 |
0.612 |
14 |
| TC22-01 |
5 |
0.786 |
0.367 |
28 |
| P70-04 |
3 |
0.308 |
0.138 |
41 |
| T20-01 |
2 |
0.191 |
0.091 |
21 |
div_df %>%
pivot_longer(c(richness, shannon, simpson),
names_to = "index", values_to = "value") %>%
ggplot(aes(x = camera_id, y = value, fill = camera_id)) +
geom_col(color = "#1e3a2f", alpha = 0.85, show.legend = FALSE) +
facet_wrap(~ index, scales = "free_y",
labeller = labeller(index = c(
richness = "Richness (S)",
shannon = "Shannon (H′)",
simpson = "Simpson (1−D)"
))) +
scale_fill_manual(
values = colorRampPalette(c("#d4e6da", "#1e3a2f"))(nrow(div_df))
) +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(family = "serif", face = "bold", color = "#1e3a2f"),
axis.text.x = element_text(angle = 40, hjust = 1),
strip.text = element_text(family = "mono", face = "bold", color = "#2d5a3d")
) +
labs(title = "Diversity indices by camera", x = "Camera ID", y = "Index value")
📌 Which diversity index to report?
Richness (S) — counts species, ignores abundance;
sensitive to sampling effort.
Shannon (H′) —
weights common and rare species equally on a log scale; the most widely
reported index in ecology.
Simpson (1−D) — weighted
toward dominant species; less sensitive to rare taxa.
Report all
three and discuss what differences reveal about community structure.
High S with low H′ indicates one highly dominant taxon with many rare
ones.
Rarefaction curves
The problem rarefaction solves. Camera A recorded 80
independent events and detected 6 species. Camera B recorded 12 events
and detected 4 species. Is Camera A genuinely richer, or did it simply
run longer and accumulate more detections by chance? Raw species counts
cannot answer this — they confound ecological richness with sampling
effort.
What a rarefaction curve does. For each camera, the
vegan function rarecurve() randomly
sub-samples its detection records 1, 2, 3 … n times and asks: on
average, how many species had we seen after exactly k detections?
The result is a curve of expected richness as a function of
effort. Plotting all cameras on the same axes — with the same x-scale —
makes a fair comparison possible: read off each curve at a common
x-value (the rarefaction target, set to the smallest
total across all cameras) to get effort-controlled richness.
How to read the plot. Each coloured line is one
camera. The vertical dashed line marks the rarefaction target — the
point where all cameras can be compared fairly. At that x-value, a
higher curve means genuinely greater local richness, not just more
deployment time.
- Steep and still rising at the target → the camera
is under-sampled; more events would likely reveal additional
species.
- Levelling off (asymptote) before the target → the
local species pool is well-characterized; additional effort would add
little.
- Short curves that barely reach the target → those
cameras had very few total events and their richness estimates are
highly uncertain.
min_n <- min(rowSums(site_sp))
cat(sprintf(
"Rarefaction target: %d detections (the minimum across %d cameras)\n",
min_n, nrow(site_sp)
))
## Rarefaction target: 14 detections (the minimum across 10 cameras)
cat(sprintf(
"Cameras with more than 2× the target: %d\n",
sum(rowSums(site_sp) > 2 * min_n)
))
## Cameras with more than 2× the target: 6
# Build a color palette with one color per camera
cam_cols <- colorRampPalette(
RColorBrewer::brewer.pal(8, "Dark2")
)(nrow(site_sp))
rarecurve(
site_sp,
step = 1,
col = cam_cols,
lwd = 1.8,
xlab = "Detection events (sampling effort)",
ylab = "Expected species richness",
main = "Rarefaction curves by camera",
label = TRUE,
cex = 0.75
)
abline(v = min_n, lty = 2, col = "#7a7a72", lwd = 1.5)
text(min_n, par("usr")[3] + 0.3,
paste0("Rarefaction target = ", min_n, " events"),
adj = c(-0.08, 0), cex = 0.8, col = "#7a7a72")
# Effort-controlled richness at the rarefaction target for each camera
rare_vals <- rarefy(site_sp, sample = min_n)
rarefied_df <- tibble(
Camera = names(rare_vals),
`Raw richness` = specnumber(site_sp),
`Total events` = rowSums(site_sp),
`Rarefied S` = round(rare_vals, 2),
`Rank (raw)` = rank(-specnumber(site_sp), ties.method = "min"),
`Rank (rarefied)` = rank(-rare_vals, ties.method = "min")
) %>% arrange(`Rank (rarefied)`)
knitr::kable(
rarefied_df,
caption = paste0(
"Effort-controlled species richness rarefied to ", min_n,
" detection events per camera. ",
"Cameras whose rarefied rank differs from their raw rank had ",
"apparent richness inflated (or deflated) by effort differences."
),
align = "lrrrrr"
)
Effort-controlled species richness rarefied to 14 detection
events per camera. Cameras whose rarefied rank differs from their raw
rank had apparent richness inflated (or deflated) by effort
differences.
| P70-01 |
9 |
82 |
5.69 |
1 |
1 |
| TC29-01 |
7 |
80 |
4.40 |
2 |
2 |
| T20-04 |
5 |
65 |
4.06 |
3 |
3 |
| P70-03 |
4 |
14 |
4.00 |
7 |
4 |
| T005 |
4 |
19 |
3.94 |
7 |
5 |
| P70-02 |
5 |
140 |
3.65 |
3 |
6 |
| T20-03 |
5 |
62 |
3.60 |
3 |
7 |
| TC22-01 |
5 |
28 |
3.39 |
3 |
8 |
| P70-04 |
3 |
41 |
1.91 |
9 |
9 |
| T20-01 |
2 |
21 |
1.67 |
10 |
10 |
💡 Rarefied S vs. raw richness — when do they
disagree?
If a camera’s raw rank is higher than its rarefied rank, its raw
species count was inflated by extra effort — it detected more species
simply because it ran longer, not because its habitat is richer. The
rarefied value is the fairer number to report in any cross-camera
comparison. Note that rarefaction always reduces or maintains
richness relative to the raw count — it can never add species that were
not observed.
8. Single-Season Occupancy Modeling
Rationale and key insight
Which cameras does species X truly occupy? The
complication is imperfect detection: even when a
species is present, a camera will not always fire. Treating “no
detection” as “truly absent” underestimates occupancy — sometimes
dramatically.
The MacKenzie et al. (2002) model simultaneously estimates:
- ψ (psi) — the probability a site is occupied.
- p — the probability of detecting the species on a
given survey occasion, given it is present.
We split the monitoring week into daily occasions. For each camera ×
day, the detection history is 1 (detected) or
0 (not detected).
Build detection histories
focal_sp <- unique(cam_events$taxa_code)[1]
cat("Focal species:", focal_sp, "\n")
## Focal species: A_californica
all_dates <- sort(unique(cam_events$date))
all_cameras <- sort(unique(cam_events$camera_id))
# 1 = species detected that camera × day; 0 = not detected
det_wide <- cam_events %>%
filter(taxa_code == focal_sp) %>%
distinct(camera_id, date) %>%
mutate(detected = 1) %>%
complete(camera_id = all_cameras, date = all_dates,
fill = list(detected = 0)) %>%
pivot_wider(names_from = date, values_from = detected) %>%
column_to_rownames("camera_id")
print(det_wide)
## 2026-03-09 2026-03-10 2026-03-11 2026-03-12 2026-03-13 2026-03-14
## P70-01 1 1 0 1 0 1
## P70-02 0 0 0 0 0 0
## P70-03 0 0 0 0 0 0
## P70-04 0 0 0 0 0 0
## T005 0 0 0 1 0 1
## T20-01 0 0 0 0 0 0
## T20-03 0 0 0 0 0 0
## T20-04 0 0 0 0 1 1
## TC22-01 0 0 0 0 0 0
## TC29-01 0 0 0 0 0 0
## 2026-03-15 2026-03-16 2026-03-17 2026-03-18 2026-03-19 2026-03-20
## P70-01 0 0 0 0 0 0
## P70-02 0 0 0 0 0 0
## P70-03 1 0 1 1 0 1
## P70-04 0 0 0 0 0 0
## T005 0 0 0 0 0 0
## T20-01 0 0 0 0 0 0
## T20-03 0 0 0 0 0 0
## T20-04 0 0 0 0 0 0
## TC22-01 1 0 0 0 0 0
## TC29-01 0 0 1 1 0 0
## 2026-03-21 2026-03-22 2026-03-23
## P70-01 0 0 0
## P70-02 0 1 0
## P70-03 0 0 0
## P70-04 0 0 0
## T005 0 0 0
## T20-01 0 0 0
## T20-03 0 0 0
## T20-04 0 0 0
## TC22-01 0 0 0
## TC29-01 0 0 0
Fit the null occupancy model
umf <- unmarkedFrameOccu(y = as.matrix(det_wide))
fit_null <- occu(~ 1 ~ 1, data = umf)
summary(fit_null)
##
## Call:
## occu(formula = ~1 ~ 1, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## 1.32 1.06 1.25 0.212
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## -1.86 0.321 -5.78 7.47e-09
##
## AIC: 104.4236
## Number of sites: 10
# Back-transform from logit scale to probability
occ_est <- backTransform(fit_null, type = "state")
det_est <- backTransform(fit_null, type = "det")
cat(sprintf(
"\nEstimated occupancy ψ = %.3f (95%% CI: %.3f – %.3f)\n",
occ_est@estimate,
confint(occ_est)[1], confint(occ_est)[2]
))
##
## Estimated occupancy ψ = 0.789 (95% CI: 0.320 – 0.968)
cat(sprintf(
"Estimated detection p = %.3f (95%% CI: %.3f – %.3f)\n",
det_est@estimate,
confint(det_est)[1], confint(det_est)[2]
))
## Estimated detection p = 0.135 (95% CI: 0.077 – 0.227)
Compare naïve vs. model-estimated occupancy
naive_occ <- mean(rowSums(det_wide) > 0)
cat(sprintf(
"Naïve occupancy (fraction of cameras with ≥1 detection): %.3f\n", naive_occ
))
## Naïve occupancy (fraction of cameras with ≥1 detection): 0.700
cat(sprintf(
"Model-estimated occupancy ψ: %.3f\n",
occ_est@estimate
))
## Model-estimated occupancy ψ: 0.789
cat(sprintf(
"Inflation due to imperfect detection: +%.3f\n",
occ_est@estimate - naive_occ
))
## Inflation due to imperfect detection: +0.089
Covariate model: does vegetation density affect occupancy?
If you have a site-level covariate (e.g., canopy cover measured at
each camera), you can test whether occupancy varies with habitat
quality. Replace the simulated veg_density below with your
real covariate.
# ── Replace this block with your real covariate data ────────────────────────
set.seed(7)
site_covs <- data.frame(
veg_density = scale(runif(nrow(det_wide), 0, 100))
)
# ────────────────────────────────────────────────────────────────────────────
umf_cov <- unmarkedFrameOccu(
y = as.matrix(det_wide),
siteCovs = site_covs
)
fit_cov <- occu(~ 1 ~ veg_density, data = umf_cov)
summary(fit_cov)
##
## Call:
## occu(formula = ~1 ~ veg_density, data = umf_cov)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 1.351 1.11 1.215 0.224
## veg_density 0.357 1.14 0.314 0.754
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## -1.86 0.321 -5.79 7.25e-09
##
## AIC: 106.3112
## Number of sites: 10
# AIC comparison
aic_df <- data.frame(
Model = c("Null (ψ(.), p(.))",
"Veg (ψ(veg), p(.))"),
AIC = round(c(fit_null@AIC, fit_cov@AIC), 2),
dAIC = round(c(0, fit_cov@AIC - fit_null@AIC), 2)
)
knitr::kable(
aic_df,
col.names = c("Model", "AIC", "ΔAIC"),
align = "lcc",
caption = paste0(
"AIC model comparison. The model with lower AIC is better supported. ",
"ΔAIC > 2 is the conventional threshold for meaningful improvement. ",
"ΔAIC < 2 means the simpler null model is preferred."
)
)
AIC model comparison. The model with lower AIC is better
supported. ΔAIC > 2 is the conventional threshold for meaningful
improvement. ΔAIC < 2 means the simpler null model is
preferred.
| Null (ψ(.), p(.)) |
104.42 |
0.00 |
| Veg (ψ(veg), p(.)) |
106.31 |
1.89 |
📌 Interpreting AIC
The model with the lower AIC has better predictive
support. ΔAIC > 2 is the conventional threshold for “meaningfully
better”. A positive coefficient on veg_density in the
occupancy sub-model means the species is more likely to occupy densely
vegetated sites. AIC penalizes complexity — a more parameterized model
must fit considerably better to justify its extra parameters.
Predicted occupancy across the vegetation gradient
The covariate model estimated a single coefficient for
veg_density on the logit (log-odds) scale.
That number is hard to interpret directly — a coefficient of +0.8 does
not mean “occupancy increases by 0.8.” The logit link is non-linear, so
the same unit increase in vegetation density produces a larger change in
ψ near 0.50 than it does near 0 or 1.
predict() translates the coefficient back into
probability space across the full observed range of the covariate,
producing the response curve — the most interpretable
summary of the model. The code below:
- Creates a prediction grid — 100 evenly spaced
values of
veg_density spanning the observed range (from the
most open to the densest site).
- Calls
predict(type = "state") —
"state" means the occupancy sub-model (ψ), as opposed to
"det" which would give the detection sub-model (p). The
function returns four columns: Predicted (ψ̂),
SE (standard error on the probability scale), and
lower/upper (95% confidence interval
bounds).
- Plots the curve and ribbon — the shaded band is the
95% CI, reflecting uncertainty in the estimated coefficient. Wide bands
indicate either a small sample of cameras or weak signal in the
data.
# Step 1: prediction grid across the full covariate range
newdat <- data.frame(
veg_density = seq(min(site_covs$veg_density),
max(site_covs$veg_density),
length.out = 100)
)
# Step 2: predict ψ at each grid point (back-transformed to probability scale)
pred <- predict(fit_cov, type = "state", newdata = newdat)
# Show what predict() returns — four columns
head(round(pred, 3))
## Predicted SE lower upper
## 1 0.722 0.284 0.140 0.976
## 2 0.724 0.278 0.146 0.976
## 3 0.726 0.273 0.152 0.975
## 4 0.728 0.268 0.158 0.974
## 5 0.730 0.263 0.164 0.974
## 6 0.731 0.259 0.171 0.973
# Step 3: plot the response curve
pred_df <- bind_cols(newdat, pred)
# Annotate where ψ = 0.5 (the inflection point of the logistic curve)
psi50_x <- pred_df$veg_density[which.min(abs(pred_df$Predicted - 0.5))]
psi50_y <- 0.5
ggplot(pred_df, aes(x = veg_density, y = Predicted)) +
# 95% CI ribbon
geom_ribbon(aes(ymin = lower, ymax = upper),
fill = "#d4e6da", alpha = 0.75) +
# Response curve
geom_line(color = "#1e3a2f", linewidth = 1.4) +
# ψ = 0.5 reference lines
geom_hline(yintercept = 0.5, linetype = "dashed",
color = "#c8832a", linewidth = 0.7) +
geom_vline(xintercept = psi50_x, linetype = "dashed",
color = "#c8832a", linewidth = 0.7) +
annotate("point", x = psi50_x, y = psi50_y,
shape = 21, size = 3.5, fill = "#c8832a", color = "white") +
annotate("text",
x = psi50_x + 0.08 * diff(range(pred_df$veg_density)),
y = 0.53,
label = sprintf("ψ = 0.50
at veg = %.2f", psi50_x),
color = "#c8832a", size = 3.2, hjust = 0, fontface = "italic") +
# Rug: actual camera covariate values
geom_rug(data = data.frame(veg_density = site_covs$veg_density),
aes(x = veg_density), y = NULL,
sides = "b", color = "#4a7c5e", alpha = 0.6, linewidth = 0.5) +
scale_y_continuous(limits = c(0, 1),
labels = scales::percent_format(),
breaks = seq(0, 1, 0.25)) +
labs(
title = paste("Predicted occupancy (ψ) for", focal_sp),
subtitle = "Shaded band = 95% CI | Rug marks = actual camera locations | Orange lines = ψ = 0.50",
x = "Vegetation density (standardized z-score)",
y = "Estimated occupancy probability (ψ)"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(family = "serif", face = "bold", color = "#1e3a2f"),
plot.subtitle = element_text(size = 9, color = "#7a7a72")
)
📌 Reading the response curve — four questions to
ask
1. Direction: Does ψ increase or decrease with the
covariate? A rising curve means the species is more likely to occupy
sites with more vegetation; a falling curve means the opposite.
2. Magnitude: How much does ψ change across the full
covariate range? If ψ only moves from 0.42 to 0.51, the effect is
statistically detectable but ecologically trivial. If it moves from 0.10
to 0.85, habitat selection is strong.
3. Uncertainty: How wide is the shaded band? Wide
CIs at the covariate extremes mean few cameras were sampled there —
interpret the predicted ψ at those values with caution.
4. Inflection point: The orange dot marks ψ = 0.50,
the covariate value at which occupancy and non-occupancy are equally
likely. This is the logistic regression analogue of an LC50 — a useful
threshold for management (e.g., “sites with vegetation density above X
are more likely than not to be occupied”).
⚠️ This example uses simulated covariate data
The veg_density values above were generated with
runif() and have no relationship to the actual cameras. The
curve shape and inflection point are therefore meaningless for your real
focal species. Replace the simulation block with measured field
covariates — canopy cover from a densiometer, distance to water,
trail-camera height class, or any site-level variable you recorded —
before drawing ecological conclusions.