First insights into the movement patterns of maned wolf with sarcoptic mange
Annotated code
Notes
This document includes the code used to conduct the analyses and is intended to support inspection of the analytical workflow. All files displayed in this annotated code are available for download from the project repository. The workflow follows this order:
- GPS data cleaning and processing with
ctmm - Home range estimation and home range overlap
- Time-Explicit Habitat Selection (TEHS) model, including data preparation and model fitting
- References
NoteHow to read this workflow
Each stage below is written as a short explanation followed by the corresponding code block for that specific task. You can see chunks of code by clicking on the little arrows you will find throughout the document.
Whenever a script depends on accessory functions or external model files, these are shown immediately after the main workflow step that calls them so the analytical logic can be followed in sequence.
1 GPS data cleaning and processing
The workflow begins with the gps_maned_wolf.rds dataset stored in data/raw/. The GPS collar data are from the two focal females used in the manuscript. In this section, we reformat the coordinates and timestamps, convert the data into a telemetry object, remove outliers, and export the processed GPS table used in the remaining analyses.
1.1 Data cleaning and processed GPS export
The complete cleaning script below reads the GPS data, standardizes the variables required by ctmm, applies the outlier-removal function, trims residual problematic fixes, and exports the processed GPS data used throughout the study.
01_data_cleaning.R: complete data cleaning workflow
# Clean GPS tracks and export the dataset to be used in all following analyses
# Load libraries
library(ctmm)
library(lubridate)
library(move2)
library(tidyverse)
# Two focal females used in the study
focal_animals <- c("704870A", "704871A")
# Create the output directory needed by the next script
dir.create(
"data/processed",
recursive = TRUE,
showWarnings = FALSE
)
# Load GPS data
gps_data <- read_rds("data/raw/gps_maned_wolf.rds")
# Standardize fields needed by ctmm and convert as telemetry objects
data_telemetry <- gps_data |>
mutate(
timestamp = ymd_hms(timestamp, quiet = TRUE),
timestamp_local = with_tz(timestamp, "America/Sao_Paulo"),
gps_latitude = as.numeric(gps_latitude),
gps_longitude = as.numeric(gps_longitude)
) |>
rename(
location.long = gps_longitude,
location.lat = gps_latitude,
individual.local.identifier = animal_id
) |>
filter(
individual.local.identifier %in% focal_animals,
!is.na(location.long),
!is.na(location.lat),
!is.na(timestamp),
year(timestamp_local) == 2024
) |>
select(-timestamp_local) |>
ctmm::as.telemetry()
# Animals' name
animals <- names(data_telemetry)
# Load the accessory function that applies sequential distance and speed filters
source("accessory_functions/01_function_ctmm_remove_outliers.R")
telemetry_no_out <- ctmm_remove_outliers(
x = data_telemetry,
individual_names = animals,
dist = 10000,
speed = 0.5
)
# Remove the very first days when residual points still fell inside the house and keep only the year for this study
telemetry_trimmed <- list()
for (animal in names(telemetry_no_out)) {
start_time <- min(mt_time(telemetry_no_out[[animal]]), na.rm = TRUE)
telemetry_trimmed[[animal]] <- telemetry_no_out[[animal]] |>
filter(
mt_time(telemetry_no_out[[animal]]) > (start_time + days(3))
)
}
# Keep only standard columns that will be reused by the next analysis scripts
class_target <- c(
"character",
"numeric",
"logical",
"factor",
"Date",
"integer64",
"POSIXct",
"POSIXlt"
)
clean_telemetry_final <- list()
for (animal in names(telemetry_trimmed)) {
clean_telemetry_final[[animal]] <- telemetry_trimmed[[animal]] |>
as.data.frame() |>
select(where(~ any(class(.) %in% class_target)))
}
# Save the single CSV consumed by the home-range and TEHS workflows
clean_telemetry_final |>
bind_rows() |>
write_csv("data/processed/gps_no_outliers_final.csv")1.2 Accessory function for outlier removal
The cleaning script depends on an accessory function that applies the same distance and speed checks to each individual.
This helper receives the telemetry objects, flags relocations that exceed the expected thresholds, and returns the cleaned move2 objects that are immediately reused by the main cleaning script.
01_function_ctmm_remove_outliers.R: helper function used in GPS cleaning
# Filter telemetry data using ctmm distance and speed diagnostics
# The function returns only the cleaned move2 objects needed by the next script
ctmm_remove_outliers <- function(
x,
individual_names = names(x),
dist = 3000,
speed = 1.5
) {
telemetry_mov2_no_out <- list()
for (animal in individual_names) {
cli::cli_h1("{animal} started")
# Read the telemetry object for one animal and diagnose extreme relocations
telemetry_data <- x[[animal]]
initial_outliers <- ctmm::outlie(telemetry_data)
# Apply the distance filter first to remove spatial spikes
telemetry_no_out_dist <- telemetry_data[initial_outliers$distance < dist, ]
# Recalculate outliers after the distance filter and then remove fast steps
distance_filtered_outliers <- ctmm::outlie(telemetry_no_out_dist)
telemetry_no_out_dist_speed <- telemetry_no_out_dist[
distance_filtered_outliers$speed < speed,
]
# Convert the filtered telemetry back to move2 so the main script can keep
# Using a tidy table-oriented workflow
telemetry_mov2_no_out[[animal]] <- move2::mt_as_move2(
telemetry_no_out_dist_speed
)
cli::cli_alert_success("Outlier filtering finished for {animal}")
}
telemetry_mov2_no_out
}2 Home range estimation and overlap
This section begins with the cleaned GPS dataset produced in stage 1. Using this file, we classified locations into wet and dry seasons, estimated annual utilization distributions, and quantified spatial overlap between individuals.
2.1 Home range models
The complete script below reads the cleaned GPS dataset, assigns the season-year combination for São Paulo state (Wet and Dry), and runs the ctmm workflow for each period. The outputs are the seasonal telemetry objects, fitted movement models, and AKDE estimates used in the overlap calculations.
02_home_range_function_ctmm.R: complete seasonal home range workflow
# Home range objects for females and save model outputs
# Load libraries
library(tidyverse)
# Source accessory function
source("accessory_functions/02_function_ctmm_seasonal.R")
# Load the cleaned dataset created in the previous step
data <- readr::read_csv(
"data/processed/gps_no_outliers_final.csv"
) |>
dplyr::mutate(
timestamp = lubridate::ymd_hms(timestamp, tz = "UTC"),
timestamp_local = lubridate::with_tz(timestamp, "America/Sao_Paulo")
) |>
dplyr::rename(individual_local_identifier = track) |>
dplyr::relocate(individual_local_identifier, .before = dplyr::everything()) |>
dplyr::filter(lubridate::year(timestamp_local) == 2024) |>
dplyr::mutate(
month_num = lubridate::month(timestamp_local),
calendar_year = lubridate::year(timestamp_local),
season_year = dplyr::if_else(
month_num <= 3,
calendar_year - 1,
calendar_year
),
season = as.factor(dplyr::if_else(
month_num >= 10 | month_num <= 3,
"Wet",
"Dry"
)),2.2 Accessory functions used for home range outputs
The home range script depends on an accessory function that executes the ctmm workflow period by period.
This accessory function loops over animal and period combinations, fits the movement model, computes AKDE home ranges, and writes the outputs that are later read by the overlap and plotting scripts.
02_function_ctmm_seasonal.R: helper functions for period-level ctmm analyses
# Process each period and save outputs
process_period_data <- function(data) {
periods <- unique(data$period)
for (period_name in periods) {
# Keep one period at a time and rename the coordinate columns for ctmm
data_filter <- data |>
filter(period == period_name) |>
mutate(timestamp = lubridate::as_datetime(timestamp, tz = "UTC")) |>
group_by(individual_local_identifier) |>
arrange(timestamp) |>
ungroup() |>
select(
individual.local.identifier = individual_local_identifier,
timestamp,
location.long = longitude,
location.lat = latitude
)
if (nrow(data_filter) < 10) {
message(sprintf(
"\nSkipping period %s due to insufficient data.\n",
period_name
))
next
}
message(sprintf("\nStarting analysis for period: %s\n", period_name))
# Convert the focal-animal data to telemetry objects for model fitting
data_telemetry <- ctmm::as.telemetry(
data_filter,
timezone = "America/Sao_Paulo"
)
if (inherits(data_telemetry, "telemetry")) {
animal_name <- unique(data_filter$individual.local.identifier)
data_telemetry <- list(data_telemetry)
names(data_telemetry) <- animal_name
}
data_identity <- names(data_telemetry)
# Save telemetry objects because they are used later
dir.create(
"data/processed/telemetry_object",
recursive = TRUE,
showWarnings = FALSE
)
readr::write_rds(
data_telemetry,
sprintf(
"data/processed/telemetry_object/data_telemetry_%s.rds",
period_name
)
)
# Fit one ctmm model per animal in the period
message(sprintf("\nFitting models for period: %s\n", period_name))
fit_list <- list()
for (animal in data_identity) {
my_tel <- data_telemetry[[animal]]
my_guess <- ctmm::ctmm.guess(my_tel, interactive = FALSE)
fit_list[[animal]] <- ctmm::ctmm.select(my_tel, my_guess, verbose = FALSE)
}
# Save fitted models because the final figure script can reuse them when it needs annual objects that are not yet on disk
dir.create(
"output/models/fit_models",
recursive = TRUE,
showWarnings = FALSE
)
readr::write_rds(
fit_list,
sprintf("output/models/fit_models/fit_list_%s.rds", period_name)
)
# Save AKDE objects because the overlap script and figure script both depend on them
message(sprintf("\nCalculating AKDE for period: %s\n", period_name))
akde_data <- ctmm::akde(data_telemetry, CTMM = fit_list, weights = TRUE)2.3 Home range overlap
Once the home range outputs exist, overlap can be computed directly from the saved AKDE estimates.
The script reads the saved AKDE objects for each period, calculates overlap metrics and their confidence intervals, and writes one consolidated table.
03_overlap_ctmm.R: complete home range overlap workflow
# Calculate seasonal home-range overlap between the female maned wolves.
# The exported table is used in Figure 1
# Load libraries
library(ctmm)
library(tidyverse)
# Load output data
akde_files <- list.files(
"output/models/fit_akde",
pattern = "list_akde_.*\\.rds$",
full.names = TRUE
)
final_overlap_results <- akde_files |>
map_dfr(
\(file_path) {
period_name <- stringr::str_extract(
basename(file_path),
"(?<=list_akde_).*(?=\\.rds)"
)
period_name <- stringr::str_replace(period_name, "Rainy", "Wet")
if (stringr::str_detect(period_name, "Total$")) {
return(tibble())
}
akde_list <- readRDS(file_path)
if (inherits(akde_list[[1]], "list")) {
akde_list <- purrr::list_flatten(akde_list)
}
# Skip periods with fewer than two focal animals because overlap cannot
# be calculated for a single home range
if (length(akde_list) < 2) {
return(tibble())
}
overlap_res <- ctmm::overlap(akde_list)
map_dfr(
c("low", "est", "high"),
\(ci_level) {
ci_matrix <- overlap_res$CI[,, ci_level]
upper_index <- which(
upper.tri(ci_matrix, diag = FALSE),
arr.ind = TRUE
)
tibble(
animal_1 = rownames(ci_matrix)[upper_index[, 1]],
animal_2 = colnames(ci_matrix)[upper_index[, 2]],
value = ci_matrix[upper_index],
ci_level = ci_level,
period = period_name
)
}
)
}
) |>
pivot_wider(names_from = ci_level, values_from = value) |>
mutate(
period = factor(
period,
levels = c(
"2023_Wet",
"2024_Dry",
"2024_Wet"
),
ordered = TRUE
)
) |>
arrange(period, animal_1, animal_2) |>
mutate(period = as.character(period))
dir.create("output/tables/overlap", recursive = TRUE, showWarnings = FALSE)
final_overlap_results |>
write_csv("output/tables/overlap/overlap_results_period_individuals.csv")Home-range estimates, seasonal overlap, and monthly distance moved for the two focal females are summarized in Figure 1. Their annual home-range map and the associated land-cover composition are presented in Figure 2.
3 Time-Explicit Habitat Selection (TEHS) model
This is the most model-intensive section of the workflow because it moves from cleaned relocation data to step construction, habitat extraction, alternative-step generation, Bayesian model fitting, and posterior consolidation.
To keep this section readable, the workflow is broken into the same order used during the analysis: first derive step-level habitat information, then prepare the objects needed by JAGS, then fit the two model components, and finally translate the posterior results into the manuscript panel.
3.1 Extract land-cover proportions for each movement step
The TEHS section starts from the same cleaned GPS file used previously in the project, but it transforms consecutive relocations into movement steps rather than home range objects. It requires the MapBiomas LULC raster for São Paulo (2023; 10-m pixels), which is not versioned on GitHub because of its size. Before running scripts 04 and 05, download the raster from MapBiomas and save it as data/raster/mapbiomas_10m_sp_2023-0000000000-0000000000.tif. The complete script below then extracts MapBiomas classes for each buffered movement segment and summarizes the step-level habitat information used in the TEHS workflow.
04_lulc_prop_TEHS.R: complete step-level LULC workflow
# Calculate step-level LULC proportions for maned wolf data
# Libraries ----
library(lubridate)
library(terra)
library(sf)
library(tidyverse)
options(scipen = 999) # avoid scientific notation
terra::terraOptions(progress = 0)
# Inputs and parameters ----
gps_path <- "data/processed/gps_no_outliers_final.csv"
raster_path <- "data/raster/mapbiomas_10m_sp_2023-0000000000-0000000000.tif"
output_dir <- "output/TEHS/mapbiomas"
raster_tag <- basename(raster_path) |> str_remove("\\.tif$")
focal_animals <- c("704870A", "704871A")
# Define parameters to be used in the code
buffer_m <- 30
speed_quantile <- 0.95
max_step_minutes <- 720
min_step_distance_m <- 1
min_lulc_use_prop <- 0.04
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
if (!file.exists(raster_path)) {
stop(
"Missing MapBiomas LULC raster: ", raster_path,
"\nDownload the 2023, 10-m LULC raster for São Paulo from ",
"https://brasil.mapbiomas.org/en/ and save it at this path before running this script."
)
}
# MapBiomas legend (Collection 10 - relevant classes) ----
lulc_legend <- tibble::tribble(
~class_code , ~class_name ,
3L , "Forest Formation" ,
4L , "Savanna Formation" ,
5L , "Mangrove" ,
6L , "Floodable Forest" ,
9L , "Forest Plantation" ,
11L , "Wetland" ,
12L , "Grassland" ,
13L , "Other Non Forest Formation" ,
15L , "Pasture" ,
18L , "Agriculture" ,
19L , "Temporary Crops" ,
20L , "Sugarcane" ,
21L , "Mosaic of Uses" ,
23L , "Beach, Dune and Sand Spot" ,
24L , "Urban Infrastructure" ,
25L , "Other Non Vegetated Area" ,
26L , "Water" ,
29L , "Rocky Outcrop" ,
30L , "Mining" ,
31L , "Aquaculture" ,
33L , "River, Lake and Ocean" ,
36L , "Perennial Crops" ,
49L , "Wooded Restinga"
)
# Read GPS data directly ----
# Read and standardize columns used in the workflow
gps <- read_csv(gps_path, show_col_types = FALSE) |>
rename(individual = track) |>
mutate(
timestamp = ymd_hms(timestamp, tz = "UTC", quiet = TRUE),
) |>
select(individual, timestamp, longitude, latitude) |>
filter(
individual %in% focal_animals,
!is.na(individual),
!is.na(timestamp),
!is.na(longitude),
!is.na(latitude)
)
# Read raster directly ----
# Use the confirmed MapBiomas raster file without tile selection logic
mapbiomas_raster <- rast(raster_path)
lulc_col <- names(mapbiomas_raster)[1]
# Convert GPS to projected coordinates for movement metrics ----
# Keep original lon/lat and add projected x/y in meters
gps_sf <- gps |>
st_as_sf(
coords = c("longitude", "latitude"),
crs = 4326,
remove = FALSE
) |>
st_transform(5880)
gps_coords <- st_coordinates(gps_sf)
# Add the columns 5880 to calulcate steps in meters
gps_5880 <- gps |>
mutate(
x_5880 = gps_coords[, 1],
y_5880 = gps_coords[, 2]
)
# Build movement steps ----
# Create start/end points and movement metrics for each individual
steps <- gps_5880 |>
arrange(individual, timestamp) |>
group_by(individual) |>
mutate(
x_sta = x_5880,
y_sta = y_5880,
x_end = lead(x_5880),
y_end = lead(y_5880),
time_sta = timestamp,
time_end = lead(timestamp),
time_min = as.numeric(difftime(time_end, time_sta, units = "mins")),
dist_m = sqrt((x_end - x_sta)^2 + (y_end - y_sta)^2),
speed_mps = dist_m / (time_min * 60),
month_num = month(time_sta),
calendar_year = year(time_sta),
season_year = if_else(
month_num <= 3,
calendar_year - 1L,
calendar_year
),
season = if_else(month_num >= 10 | month_num <= 3, "Wet", "Dry"),
period = paste(season_year, season, sep = "_")
) |>
ungroup() |>
filter(
!is.na(x_end),
!is.na(y_end),
!is.na(time_min),
time_min > 0,
!is.na(dist_m),
!is.na(speed_mps)
) |>
mutate(
dist_m = if_else(dist_m == 0, 1, dist_m)
)
# Filter movement outliers ----
# Remove very fast movements and very long temporal gaps
speed_limits <- steps |>
group_by(individual) |>
summarise(
speed_limit = quantile(speed_mps, speed_quantile, na.rm = TRUE),
.groups = "drop"
)
steps_filter <- steps |>
left_join(speed_limits, by = "individual") |>
filter(
speed_mps <= speed_limit,
time_min <= max_step_minutes,
dist_m > min_step_distance_m
) |>
select(-speed_limit)
if (nrow(steps_filter) == 0) {
stop("All movement steps were removed by filters. Review thresholds.")
}
# Periods by year and season
period_levels <- steps_filter |>
group_by(period) |>
summarise(start_date = min(time_sta), .groups = "drop") |>
arrange(start_date) |>
pull(period)
steps_filter <- steps_filter |>
mutate(
step_id = row_number(),
period = factor(period, levels = period_levels, ordered = TRUE)
)
# Create line buffers for each step ----
# Each step is converted to a 30 m buffer polygon
step_lines_geom <- lapply(seq_len(nrow(steps_filter)), function(i) {
st_linestring(
matrix(
c(
steps_filter$x_sta[i],
steps_filter$y_sta[i],
steps_filter$x_end[i],
steps_filter$y_end[i]
),
ncol = 2,
byrow = TRUE
)
)
})
# Create buffer
step_buffers <- st_sf(
step_id = steps_filter$step_id,
geometry = st_sfc(step_lines_geom, crs = 5880)
) |>
st_buffer(dist = buffer_m)
step_buffers_vect <- vect(step_buffers)
if (!same.crs(mapbiomas_raster, step_buffers_vect)) {
step_buffers_vect <- project(step_buffers_vect, mapbiomas_raster)
}
# Crop the raster once to the buffered step extent before extracting pixels
mapbiomas_crop <- terra::crop(
mapbiomas_raster,
terra::ext(step_buffers_vect),
snap = "out"
)
# Extract raster classes inside each step buffer ----
# Keep only non-missing class codes
extracted <- terra::extract(
mapbiomas_crop[[1]],
step_buffers_vect,
ID = TRUE
) |>
as_tibble() |>
rename(step_id = ID) |>
mutate(
class_code = as.integer(.data[[lulc_col]])
) |>
select(step_id, class_code) |>
filter(!is.na(class_code))
if (nrow(extracted) == 0) {
stop("No LULC pixels were extracted from step buffers.")
}
# Recode classes for the analysis setup ----
# Merge Temporary and Perennial Crops
step_lulc <- extracted |>
left_join(lulc_legend, by = "class_code") |>
mutate(
class_name = coalesce(class_name, paste0("Class ", class_code)),
class_group = case_when(
class_code %in% c(19L, 36L) ~ "Crops",
TRUE ~ class_name
)
)
# Calculate per-step proportions in long format ----
step_lulc_counts <- step_lulc |>
count(step_id, class_group, name = "n_pixels")
step_lulc_long_all <- step_lulc_counts |>
group_by(step_id) |>
mutate(prop_lulc = n_pixels / sum(n_pixels)) |>
ungroup()
# Calculate overall usage and apply threshold ----
# Keep only classes above the minimum usage threshold
lulc_usage_overall <- step_lulc_counts |>
group_by(class_group) |>
summarise(n_pixels = sum(n_pixels), .groups = "drop") |>
mutate(
prop_use = n_pixels / sum(n_pixels),
keep_class = prop_use >= min_lulc_use_prop
) |>
arrange(desc(prop_use))
kept_groups <- lulc_usage_overall |>
filter(keep_class) |>
pull(class_group)
if (length(kept_groups) == 0) {
stop(
"No LULC class reached the threshold of ",
scales::percent(min_lulc_use_prop, accuracy = 0.1),
"."
)
}
step_lulc_long <- step_lulc_long_all |>
filter(class_group %in% kept_groups)
# Build column names for wide output ----
# Convert class labels into safe snake_case column names
lulc_columns <- lulc_usage_overall |>
filter(keep_class) |>
mutate(
lulc_col = paste0(
"lulc_",
gsub(
"_+",
"_",
gsub(
"^_|_$",
"",
gsub(
"[^a-z0-9]+",
"_",
tolower(iconv(class_group, to = "ASCII//TRANSLIT"))
)
)
)
)
) |>
select(class_group, lulc_col)
# Create final step-level dataset (wide format) ----
# One row per step with one column per kept class
step_lulc_wide <- step_lulc_long |>
left_join(lulc_columns, by = "class_group") |>
select(step_id, lulc_col, prop_lulc) |>
pivot_wider(
names_from = lulc_col,
values_from = prop_lulc,
values_fill = 0
)
final_database <- steps_filter |>
left_join(step_lulc_wide, by = "step_id") |>
arrange(individual, time_sta)
lulc_col_names <- lulc_columns$lulc_col
if (length(lulc_col_names) > 0) {
final_database <- final_database |>
mutate(
across(all_of(lulc_col_names), ~ replace_na(.x, 0))
)
}
# Save only the files consumed by the next TEHS preparation step
write_csv(
final_database,
file.path(output_dir, str_glue("step_lulc_prop_common_{raster_tag}.csv"))
)
write_csv(
lulc_usage_overall,
file.path(
output_dir,
str_glue("step_lulc_overall_use_threshold_{raster_tag}.csv")
)
)
write_csv(
lulc_columns,
file.path(output_dir, str_glue("step_lulc_columns_kept_{raster_tag}.csv"))
)
cat(
"LULC step extraction finished.\n",
"Output dir: ",
output_dir,
"\n",
"Threshold: ",
scales::percent(min_lulc_use_prop, accuracy = 0.1),
"\n",
"Groups kept: ",
paste(kept_groups, collapse = ", "),
"\n",
sep = ""
)3.2 Prepare TEHS inputs and alternative potential steps
Once land-cover proportions have been calculated for each step, the next script organizes the data into the structure required by the TEHS models. It defines the habitat covariates, creates the temporal variables used by JAGS, and generates the observed step together with the alternative cardinal-direction steps for each animal.
This is the transition from descriptive movement processing to model-ready data preparation. The outputs written here are the exact tables consumed by the JAGS fitting script.
05_prep_data_to_run_TEHS.R: complete TEHS preparation workflow
# Time-Explicit Habitat Selection Model - data preparation
# Load libraries
library(terra)
library(sf)
library(tidyverse)
options(scipen = 999)
terra::terraOptions(progress = 0)
# Paths and parameters ----
raster_path <- "data/raster/mapbiomas_10m_sp_2023-0000000000-0000000000.tif"
raster_tag <- basename(raster_path) |> str_remove("\\.tif$")
tehs_dir <- "output/TEHS"
mapbiomas_dir <- file.path(tehs_dir, "mapbiomas")
prep_dir <- file.path(tehs_dir, "prep")
potential_steps_dir <- file.path(prep_dir, "potential_steps")
prep_complete_dir <- file.path(prep_dir, "prep_jags_complete")
buffer_m <- 30
miss_factor <- 1.5
cardinal_names <- c("west", "east", "north", "south")
walk(
c(prep_dir, potential_steps_dir, prep_complete_dir),
~ dir.create(.x, recursive = TRUE, showWarnings = FALSE)
)
if (!file.exists(raster_path)) {
stop(
"Missing MapBiomas LULC raster: ", raster_path,
"\nDownload the 2023, 10-m LULC raster for São Paulo from ",
"https://brasil.mapbiomas.org/en/ and save it at this path before running this script."
)
}
# Input files ----
steps_path <- file.path(
mapbiomas_dir,
str_glue("step_lulc_prop_common_{raster_tag}.csv")
)
lulc_columns_path <- file.path(
mapbiomas_dir,
str_glue("step_lulc_columns_kept_{raster_tag}.csv")
)
lulc_usage_path <- file.path(
mapbiomas_dir,
str_glue("step_lulc_overall_use_threshold_{raster_tag}.csv")
)
if (!file.exists(steps_path)) {
stop("Missing input file: ", steps_path)
}
if (!file.exists(lulc_columns_path)) {
stop("Missing input file: ", lulc_columns_path)
}
if (!file.exists(lulc_usage_path)) {
stop("Missing input file: ", lulc_usage_path)
}
# MapBiomas legend ----
lulc_legend <- tibble::tribble(
~class_code , ~class_name ,
3L , "Forest Formation" ,
4L , "Savanna Formation" ,
5L , "Mangrove" ,
6L , "Floodable Forest" ,
9L , "Forest Plantation" ,
11L , "Wetland" ,
12L , "Grassland" ,
13L , "Other Non Forest Formation" ,
15L , "Pasture" ,
18L , "Agriculture" ,
19L , "Temporary Crops" ,
20L , "Sugarcane" ,
21L , "Mosaic of Uses" ,
23L , "Beach, Dune and Sand Spot" ,
24L , "Urban Infrastructure" ,
25L , "Other Non Vegetated Area" ,
26L , "Water" ,
29L , "Rocky Outcrop" ,
30L , "Mining" ,
31L , "Aquaculture" ,
33L , "River, Lake and Ocean" ,
36L , "Perennial Crops" ,
49L , "Wooded Restinga"
)
# Read inputs ----
mapbiomas_raster <- rast(raster_path)
raster_lulc_col <- names(mapbiomas_raster)[1]
steps_common <- read_csv(steps_path, show_col_types = FALSE)
lulc_columns <- read_csv(lulc_columns_path, show_col_types = FALSE)
lulc_usage <- read_csv(lulc_usage_path, show_col_types = FALSE)
required_step_cols <- c(
"individual",
"step_id",
"time_min",
"dist_m",
"x_sta",
"y_sta",
"x_end",
"y_end"
)
missing_required <- setdiff(required_step_cols, names(steps_common))
if (length(missing_required) > 0) {
stop(
"Missing required columns in step input: ",
paste(missing_required, collapse = ", ")
)
}
lulc_cols_kept <- lulc_columns$lulc_col
if (!all(lulc_cols_kept %in% names(steps_common))) {
stop("Not all LULC columns from lulc_columns file exist in step input.")
}
# Define model covariates (baseline = most used class) ----
baseline_group <- lulc_usage |>
filter(keep_class) |>
arrange(desc(prop_use)) |>
slice(1) |>
pull(class_group)
baseline_col <- lulc_columns |>
filter(class_group == baseline_group) |>
pull(lulc_col)
model_covariates <- lulc_columns |>
filter(lulc_col != baseline_col) |>
pull(lulc_col)
covariate_config <- lulc_columns |>
mutate(
raster_tag = raster_tag,
baseline_group = baseline_group,
baseline_col = baseline_col,
use_in_model = lulc_col %in% model_covariates,
covariate_index = cumsum(use_in_model)
)
write_csv(
covariate_config,
file.path(prep_dir, str_glue("tehs_covariate_config_{raster_tag}.csv"))
)
# Build time variables used by JAGS time model ----
fix_interval_min <- median(steps_common$time_min, na.rm = TRUE)
miss_threshold_min <- fix_interval_min * miss_factor
steps_tehs <- steps_common |>
mutate(
time1 = time_min,
dist1 = dist_m,
Miss = if_else(time1 <= miss_threshold_min, 0, 1)
)
animals <- steps_tehs |>
distinct(individual) |>
pull(individual)
# Prepare per-animal files ----
for (animal in animals) {
cat("\\n--------------------------------------------\\n")
cat("Preparing TEHS files for animal: ", animal, "\\n", sep = "")
steps_animal <- steps_tehs |>
filter(individual == animal) |>
arrange(time_sta) |>
mutate(
x_end1 = x_sta - dist1,
y_end1 = y_sta,
x_end2 = x_sta + dist1,
y_end2 = y_sta,
x_end3 = x_sta,
y_end3 = y_sta + dist1,
x_end4 = x_sta,
y_end4 = y_sta - dist1
)
if (nrow(steps_animal) == 0) {
warning("No steps for animal: ", animal)
next
}
# Crop the raster once to the full extent covered by the observed and potential steps for this animal
animal_x <- c(
steps_animal$x_sta,
steps_animal$x_end,
steps_animal$x_end1,
steps_animal$x_end2,
steps_animal$x_end3,
steps_animal$x_end4
)
animal_y <- c(
steps_animal$y_sta,
steps_animal$y_end,
steps_animal$y_end1,
steps_animal$y_end2,
steps_animal$y_end3,
steps_animal$y_end4
)
animal_extent_sf <- sf::st_as_sf(
tibble(
geometry = sf::st_sfc(
sf::st_as_sfc(
sf::st_bbox(
c(
xmin = min(animal_x, na.rm = TRUE) - buffer_m,
xmax = max(animal_x, na.rm = TRUE) + buffer_m,
ymin = min(animal_y, na.rm = TRUE) - buffer_m,
ymax = max(animal_y, na.rm = TRUE) + buffer_m
),
crs = sf::st_crs(5880)
)
)
)
)
) |>
sf::st_transform(crs = terra::crs(mapbiomas_raster))
mapbiomas_crop_animal <- terra::crop(
mapbiomas_raster,
terra::ext(terra::vect(animal_extent_sf)),
snap = "out"
)
write_csv(
steps_animal,
file.path(
potential_steps_dir,
str_glue("potential_steps_coordinates_{animal}_{raster_tag}.csv")
)
)
# Build all four candidate directions together so the raster is extracted only once per animal
candidate_steps <- purrr::map_dfr(seq_along(cardinal_names), \(j) {
x_end_col <- str_glue("x_end{j}")
y_end_col <- str_glue("y_end{j}")
steps_animal |>
transmute(
step_id,
path_idx = j,
x_sta,
y_sta,
x_end = .data[[x_end_col]],
y_end = .data[[y_end_col]]
)
})
candidate_lines_geom <- lapply(seq_len(nrow(candidate_steps)), function(i) {
st_linestring(
matrix(
c(
candidate_steps$x_sta[i],
candidate_steps$y_sta[i],
candidate_steps$x_end[i],
candidate_steps$y_end[i]
),
ncol = 2,
byrow = TRUE
)
)
})
candidate_buffers <- st_sf(
row_id = seq_len(nrow(candidate_steps)),
step_id = candidate_steps$step_id,
path_idx = candidate_steps$path_idx,
geometry = st_sfc(candidate_lines_geom, crs = 5880)
) |>
st_buffer(dist = buffer_m)
candidate_buffers_vect <- vect(candidate_buffers)
if (!same.crs(mapbiomas_crop_animal, candidate_buffers_vect)) {
candidate_buffers_vect <- project(
candidate_buffers_vect,
mapbiomas_crop_animal
)
}
extracted <- terra::extract(
mapbiomas_crop_animal[[1]],
candidate_buffers_vect,
ID = TRUE
) |>
as_tibble() |>
rename(row_id = ID) |>
mutate(class_code = as.integer(.data[[raster_lulc_col]])) |>
select(row_id, class_code) |>
filter(!is.na(class_code)) |>
left_join(
candidate_buffers |>
st_drop_geometry() |>
select(row_id, step_id, path_idx),
by = "row_id"
)
if (nrow(extracted) == 0) {
stop("No LULC extracted for animal ", animal)
}
candidate_covariates <- extracted |>
left_join(lulc_legend, by = "class_code") |>
mutate(
class_name = coalesce(class_name, paste0("Class ", class_code)),
class_group = case_when(
class_code %in% c(19L, 36L) ~ "Crops",
TRUE ~ class_name
)
) |>
count(step_id, path_idx, class_group, name = "n_pixels") |>
group_by(step_id, path_idx) |>
mutate(prop_lulc = n_pixels / sum(n_pixels)) |>
ungroup() |>
left_join(lulc_columns, by = "class_group") |>
filter(!is.na(lulc_col)) |>
group_by(step_id, path_idx, lulc_col) |>
summarise(prop_lulc = sum(prop_lulc), .groups = "drop") |>
mutate(path_col = str_glue("{lulc_col}_{path_idx}")) |>
select(step_id, path_col, prop_lulc) |>
pivot_wider(
names_from = path_col,
values_from = prop_lulc,
values_fill = 0
)
missing_candidate_cols <- setdiff(
as.vector(outer(lulc_cols_kept, 1:4, paste, sep = "_")),
names(candidate_covariates)
)
if (length(missing_candidate_cols) > 0) {
candidate_covariates[missing_candidate_cols] <- 0
}
# Build full tSSF file (realized + 4 cardinal alternatives)
realized_covariates <- steps_animal |>
select(step_id, time1, dist1, Miss, all_of(lulc_cols_kept)) |>
rename_with(~ str_glue("{.x}_0"), all_of(lulc_cols_kept))
tssf_data_full <- realized_covariates |>
left_join(candidate_covariates, by = "step_id")3.3 Run the time model and habitat selection model in JAGS
After the preparation step, the workflow has everything needed to fit the Bayesian models. The modeling script first estimates the time component for each individual, saves the posterior draws, and then uses the prepared tSSF tables together with those draws to fit the habitat-selection component.
Separating the workflow into a time model and a selection model makes the fitting sequence easier to inspect and also keeps the saved outputs clean. In the published workflow, only the posterior draws needed by the next script are retained.
06_run_TEHS_models.R: complete TEHS model-fitting workflow
# Time-Explicit Habitat Selection Model - run time model and TEHS model
.libPaths(c("../r_libs", "r_libs", .libPaths()))
# Load libraries
library(jagsUI)
library(tidyverse)
options(scipen = 999)
# Paths and parameters ----
raster_path <- "data/raster/mapbiomas_10m_sp_2023-0000000000-0000000000.tif"
raster_tag <- basename(raster_path) |> str_remove("\\.tif$")
tehs_dir <- "output/TEHS"
prep_dir <- file.path(tehs_dir, "prep")
potential_steps_dir <- file.path(prep_dir, "potential_steps")
prep_complete_dir <- file.path(prep_dir, "prep_jags_complete")
models_dir <- file.path(tehs_dir, "models")
time_model_results_dir <- file.path(models_dir, "time_model_results")
selection_model_results_dir <- file.path(models_dir, "selection_model_results")
walk(
c(models_dir, time_model_results_dir, selection_model_results_dir),
~ dir.create(.x, recursive = TRUE, showWarnings = FALSE)
)
# Model files ----
time_model_file <- "accessory_functions/03_resist_avg_jags.R"
selection_model_file <- "accessory_functions/04_jags_tssf.R"
if (!file.exists(time_model_file)) {
stop("Missing JAGS model file: ", time_model_file)
}
if (!file.exists(selection_model_file)) {
stop("Missing JAGS model file: ", selection_model_file)
}
# Read covariate config ----
covariate_config_path <- file.path(
prep_dir,
str_glue("tehs_covariate_config_{raster_tag}.csv")
)
if (!file.exists(covariate_config_path)) {
stop("Missing covariate config: ", covariate_config_path)
}
covariate_config <- read_csv(covariate_config_path, show_col_types = FALSE)
model_covariates <- covariate_config |>
filter(use_in_model) |>
arrange(covariate_index) |>
pull(lulc_col)
if (length(model_covariates) == 0) {
stop("No model covariates available after baseline removal.")
}
# Discover animals from prepared files ----
animals <- list.files(
path = prep_complete_dir,
pattern = str_glue("^tSSF_data_.*_{raster_tag}[.]csv$")
) |>
discard(~ str_detect(.x, "^tSSF_data_full_")) |>
str_remove("^tSSF_data_") |>
str_remove(str_glue("_{raster_tag}[.]csv$"))
if (length(animals) == 0) {
stop("No prepared tSSF files found in: ", prep_complete_dir)
}
# 1) Time model ----
jags_time_results <- list()
for (animal in animals) {
cat("\\n--------------------------------------------\\n")
cat("Running TIME model for animal: ", animal, "\\n", sep = "")
potential_path <- file.path(
potential_steps_dir,
str_glue("potential_steps_coordinates_{animal}_{raster_tag}.csv")
)
time_output_path <- file.path(
time_model_results_dir,
str_glue("sims.list_{animal}_{raster_tag}.csv")
)
if (!file.exists(potential_path)) {
warning("Missing potential steps file for animal: ", animal)
next
}
if (file.exists(time_output_path)) {
message(
"Time model output already exists for ",
animal,
". Skipping refit."
)
next
}
time_data <- read_csv(potential_path, show_col_types = FALSE)
required_time_cols <- c("time1", "dist1", "Miss", model_covariates)
missing_time_cols <- setdiff(required_time_cols, names(time_data))
if (length(missing_time_cols) > 0) {
warning(
"Skipping animal ",
animal,
" due to missing columns: ",
paste(missing_time_cols, collapse = ", ")
)
next
}
dat_jags <- list(
delta.time = time_data$time1,
dist = time_data$dist1,
nobs = nrow(time_data),
ncov = length(model_covariates),
Miss = time_data$Miss,
xmat = data.matrix(time_data[, model_covariates])
)
set.seed(2)
params <- c("b0", "b", "g0", "g1")
n.iter <- 5000
n.thin <- 10
n.burnin <- n.iter / 2
n.chains <- 3
time_fit <- jags(
model.file = time_model_file,
parameters.to.save = params,
data = dat_jags,
n.chains = n.chains,
n.burnin = n.burnin,
n.iter = n.iter,
n.thin = n.thin,
DIC = FALSE
)
jags_time_results[[animal]] <- time_fit
# Save only the posterior draws consumed by the inspection script
sims_time_df <- as_tibble(as.data.frame(time_fit$sims.list))
write_csv(
sims_time_df,
time_output_path
)
}
# 2) TEHS model ----
jags_selection_results <- list()
for (animal in animals) {
cat("\\n--------------------------------------------\\n")
cat("Running TEHS model for animal: ", animal, "\\n", sep = "")
tssf_path <- file.path(
prep_complete_dir,
str_glue("tSSF_data_{animal}_{raster_tag}.csv")
)
sims_time_path <- file.path(
time_model_results_dir,
str_glue("sims.list_{animal}_{raster_tag}.csv")
)
selection_output_path <- file.path(
selection_model_results_dir,
str_glue("sims.list_{animal}_{raster_tag}.csv")
)
if (file.exists(selection_output_path)) {
message(
"Selection model output already exists for ",
animal,
". Skipping refit."
)
next
}
if (!file.exists(tssf_path) || !file.exists(sims_time_path)) {
warning(
"Skipping animal ",
animal,
" because required input file is missing."
)
next
}
tssf_data <- read_csv(tssf_path, show_col_types = FALSE)
sims_time <- read_csv(sims_time_path, show_col_types = FALSE)
b0_col <- names(sims_time) |>
keep(~ .x == "b0")
b_effect_cols <- names(sims_time) |>
keep(
~ str_detect(.x, "^b\\[[0-9]+\\]$") |
str_detect(.x, "^b[0-9]+$") |
str_detect(.x, "^b[.][0-9]+$")
) |>
discard(~ .x == "b0")
if (length(b0_col) == 0) {
warning("Skipping animal ", animal, ": missing b0 in time model output.")
next
}
if (length(b_effect_cols) > 0) {
b_effect_index <- case_when(
str_detect(b_effect_cols, "^b\\[[0-9]+\\]$") ~ as.numeric(str_extract(
b_effect_cols,
"[0-9]+"
)),
str_detect(b_effect_cols, "^b[.][0-9]+$") ~ as.numeric(str_extract(
b_effect_cols,
"[0-9]+"
)),
TRUE ~ as.numeric(str_remove(b_effect_cols, "^b"))
)
b_effect_cols <- b_effect_cols[order(b_effect_index)]
}
g_cols <- names(sims_time) |>
keep(~ str_detect(.x, "^g[0-9]+$")) |>
(\(x) x[order(as.numeric(str_remove(x, "^g")))])()
if (length(b_effect_cols) != length(model_covariates)) {
warning(
"Skipping animal ",
animal,
": number of b effects does not match model covariates."
)
next
}
if (length(g_cols) < 2) {
warning(
"Skipping animal ",
animal,
": missing g parameters in time model output."
)
next
}
required_tssf_cols <- c(
"time1",
"dist1",
"Miss",
str_glue("{model_covariates}_{0}"),
str_glue("{model_covariates}_{1}"),
str_glue("{model_covariates}_{2}"),
str_glue("{model_covariates}_{3}"),
str_glue("{model_covariates}_{4}")
)
missing_tssf_cols <- setdiff(required_tssf_cols, names(tssf_data))
if (length(missing_tssf_cols) > 0) {
warning(
"Skipping animal ",
animal,
" due to missing tSSF columns: ",
paste(missing_tssf_cols, collapse = ", ")
)
next
}
betas_mat <- cbind(
b0 = sims_time[[b0_col]],
sims_time |>
select(all_of(b_effect_cols)) |>
as.matrix()
)
gs_mat <- sims_time |>
select(all_of(g_cols[1:2])) |>
as.matrix()
prob <- matrix(NA_real_, nrow(tssf_data), 5)
for (i in seq_len(nrow(tssf_data))) {
for (path_idx in 0:4) {
covar_cols_idx <- str_glue("{model_covariates}_{path_idx}")
xmat <- c(1, as.numeric(tssf_data[i, covar_cols_idx]))
mean1 <- tssf_data$dist1[i] * exp(xmat %*% t(betas_mat))
b1 <- exp(gs_mat[, 1] + gs_mat[, 2] * tssf_data$Miss[i])
a1 <- b1 * mean1
prob[i, path_idx + 1] <- median(dgamma(tssf_data$time1[i], a1, b1))
}
}
colnames(prob) <- paste0("prob", 0:4)
dat <- bind_cols(tssf_data, as_tibble(prob))
covariate_cols_all_paths <- c(
str_glue("{model_covariates}_0"),
str_glue("{model_covariates}_1"),
str_glue("{model_covariates}_2"),
str_glue("{model_covariates}_3"),
str_glue("{model_covariates}_4")
)
prob_cols <- paste0("prob", 0:4)
dat <- dat |>
filter(if_all(all_of(c(covariate_cols_all_paths, prob_cols)), ~ !is.na(.x)))
if (nrow(dat) == 0) {
warning(
"Skipping animal ",
animal,
": no complete rows for TEHS model after NA filter."
)
next
}
dat_jags <- list(
nobs = nrow(dat),
ncov = length(model_covariates),
xmat0 = data.matrix(dat[, str_glue("{model_covariates}_0")]),
xmat1 = data.matrix(dat[, str_glue("{model_covariates}_1")]),
xmat2 = data.matrix(dat[, str_glue("{model_covariates}_2")]),
xmat3 = data.matrix(dat[, str_glue("{model_covariates}_3")]),
xmat4 = data.matrix(dat[, str_glue("{model_covariates}_4")]),
pmov0 = dat$prob0,
pmov1 = dat$prob1,3.4 Accessory JAGS models used by TEHS
The TEHS branch depends on two accessory model files: one for the time model and one for the habitat-selection component. These are not optional implementation details; they are part of the analysis itself because they define the statistical models that generated the manuscript results.
For that reason, both files are shown directly in this book so the published workflow contains not only the R scripts that call JAGS, but also the model code that JAGS executes.
03_resist_avg_jags.R: JAGS time model
# Time model for TEHS (dynamic number of LULC covariates)
# Expected data:
# - nobs: number of steps
# - ncov: number of LULC covariates
# - delta.time: observed travel time
# - dist: step distance
# - Miss: missing-fix indicator (0/1)
# - xmat: matrix [nobs x ncov] of LULC covariates
model{
for (i in 1:nobs) {
# Mean of gamma distribution
log_mu[i] <- b0 + inprod(xmat[i, 1:ncov], b[1:ncov])
mu[i] <- dist[i] * exp(log_mu[i])
# Dispersion terms
rate[i] <- exp(g0 + g1 * Miss[i])
shape[i] <- mu[i] * rate[i]
# Likelihood
delta.time[i] ~ dgamma(shape[i], rate[i])
}
# Priors
b0 ~ dnorm(0, 0.01)
for (k in 1:ncov) {
b[k] ~ dnorm(0, 1)
}
g0 ~ dnorm(0, 0.01)
g1 ~ dnorm(0, 0.01)
}
04_jags_tssf.R: JAGS habitat-selection model
# TEHS selection model
# Expected data:
# - nobs: number of observed steps
# - ncov: number of LULC covariates
# - xmat0..xmat4: covariate matrices for realized + 4 alternatives
# - pmov0..pmov4: movement probabilities from time model
# - y: ones vector (for one-trick)
model{
eps <- 1.0E-12
for (i in 1:nobs) {
# Relative probability of each alternative
p0[i] <- pmov0[i] * exp(inprod(xmat0[i, 1:ncov], betas[1:ncov]))
p1[i] <- pmov1[i] * exp(inprod(xmat1[i, 1:ncov], betas[1:ncov]))
p2[i] <- pmov2[i] * exp(inprod(xmat2[i, 1:ncov], betas[1:ncov]))
p3[i] <- pmov3[i] * exp(inprod(xmat3[i, 1:ncov], betas[1:ncov]))
p4[i] <- pmov4[i] * exp(inprod(xmat4[i, 1:ncov], betas[1:ncov]))
denom[i] <- p0[i] + p1[i] + p2[i] + p3[i] + p4[i] + eps
pi[i] <- p0[i] / denom[i]
# One-trick likelihood
y[i] ~ dbern(pi[i])
}
# Priors
for (k in 1:ncov) {
betas[k] ~ dnorm(0, 1)
}
}3.5 Consolidate posterior outputs into model-summary tables
After model fitting, the posterior draws still need to be reorganized into tidy summaries. This script performs that conversion by reading the saved posterior samples, standardizing them into a consistent tabular structure, and exporting the resulting tables.
07_inspect_TEHS_results.R: complete posterior consolidation workflow
# Consolidate TEHS posterior draws into the two files consumed directly by the final figure script
# Load library
library(tidyverse)
options(scipen = 999)
raster_path <- "data/raster/mapbiomas_10m_sp_2023-0000000000-0000000000.tif"
raster_tag <- basename(raster_path) |> str_remove("\\.tif$")
tehs_dir <- "output/TEHS"
prep_dir <- file.path(tehs_dir, "prep")
models_dir <- file.path(tehs_dir, "models")
time_model_results_dir <- file.path(models_dir, "time_model_results")
selection_model_results_dir <- file.path(models_dir, "selection_model_results")
results_dir <- file.path(tehs_dir, "results")
dir.create(results_dir, recursive = TRUE, showWarnings = FALSE)
# Read the model covariate order so the posterior columns stay aligned with the TEHS design matrix
covariate_config_path <- file.path(
prep_dir,
str_glue("tehs_covariate_config_{raster_tag}.csv")
)
if (!file.exists(covariate_config_path)) {
stop("Missing covariate config: ", covariate_config_path)
}
covariate_config <- read_csv(covariate_config_path, show_col_types = FALSE)
model_covariates <- covariate_config |>
filter(use_in_model) |>
arrange(covariate_index) |>
pull(lulc_col)
selection_files <- list.files(
path = selection_model_results_dir,
pattern = str_glue("^sims.list_.*_{raster_tag}[.]csv$"),
full.names = TRUE
)
if (length(selection_files) == 0) {
stop(
"No selection model output files found in: ",
selection_model_results_dir
)
}
animals <- basename(selection_files) |>
str_remove("^sims[.]list_") |>
str_remove(str_glue("_{raster_tag}[.]csv$"))
selection_sims_list <- list()
time_sims_list <- list()
for (animal in animals) {
selection_sims_path <- file.path(
selection_model_results_dir,
str_glue("sims.list_{animal}_{raster_tag}.csv")
)
time_sims_path <- file.path(
time_model_results_dir,
str_glue("sims.list_{animal}_{raster_tag}.csv")
)
if (file.exists(selection_sims_path)) {
selection_sims_list[[animal]] <- read_csv(
selection_sims_path,
show_col_types = FALSE
)
}
if (file.exists(time_sims_path)) {
time_sims_list[[animal]] <- read_csv(
time_sims_path,
show_col_types = FALSE
)
}
}
# Save one posterior table for the habitat-selection model
selection_posterior <- map2_dfr(
selection_sims_list,
names(selection_sims_list),
function(df, animal) {
beta_cols <- names(df) |>
keep(
~ str_detect(.x, "^betas\\[[0-9]+\\]$") |
str_detect(.x, "^betas[.][0-9]+$") |
str_detect(.x, "^betas[0-9]+$")
)
if (length(beta_cols) > 0) {
beta_index <- case_when(
str_detect(beta_cols, "^betas\\[[0-9]+\\]$") ~ as.numeric(str_extract(
beta_cols,
"[0-9]+"
)),
str_detect(beta_cols, "^betas[.][0-9]+$") ~ as.numeric(str_extract(
beta_cols,
"[0-9]+$"
)),
TRUE ~ as.numeric(str_remove(beta_cols, "^betas"))
)
beta_cols <- beta_cols[order(beta_index)]
}
df |>
select(all_of(beta_cols)) |>
set_names(model_covariates[seq_len(length(beta_cols))]) |>
mutate(individual = animal, .before = 1)
}
)
write_csv(
selection_posterior,
file.path(results_dir, str_glue("TEHS_selection_model_{raster_tag}.csv"))
)
# Save one posterior table for the time model
time_posterior <- map2_dfr(
time_sims_list,
names(time_sims_list),
function(df, animal) {
if (!"b0" %in% names(df)) {
return(tibble())
}
b_effect_cols <- names(df) |>
keep(
~ str_detect(.x, "^b\\[[0-9]+\\]$") |
str_detect(.x, "^b[0-9]+$") |
str_detect(.x, "^b[.][0-9]+$")
) |>
discard(~ .x == "b0")
if (length(b_effect_cols) > 0) {
b_effect_index <- case_when(
str_detect(b_effect_cols, "^b\\[[0-9]+\\]$") ~ as.numeric(str_extract(
b_effect_cols,
"[0-9]+"
)),
str_detect(b_effect_cols, "^b[.][0-9]+$") ~ as.numeric(str_extract(
b_effect_cols,
"[0-9]+"
)),
TRUE ~ as.numeric(str_remove(b_effect_cols, "^b"))
)
b_effect_cols <- b_effect_cols[order(b_effect_index)]
}
b_draws <- tibble(intercept = df$b0)
if (length(b_effect_cols) > 0) {
b_effect_df <- df |>
select(all_of(b_effect_cols))
names(b_effect_df) <- model_covariates[seq_len(ncol(b_effect_df))]
b_draws <- bind_cols(b_draws, b_effect_df)
}
b_draws |>
mutate(individual = animal, .before = 1)
}
)
write_csv(
time_posterior,
file.path(results_dir, str_glue("TEHS_time_model_{raster_tag}.csv"))Posterior estimates from the habitat-selection and time components of the TEHS model, shown separately for the two focal females, are summarized in Figure 3.
4 References
The workflow shown here relies mainly on two methodological branches: continuous-time home range estimation with ctmm and time-explicit habitat selection based on step-level movement covariates. Land use classifications are available on MapBiomas.
Fleming, C. H., and J. M. Calabrese. 2025.
ctmm: Continuous-Time Movement Modeling. R package version 1.3.0. https://CRAN.R-project.org/package=ctmm. DOI: https://doi.org/10.32614/CRAN.package.ctmmFleming, C. H., W. F. Fagan, T. Mueller, K. A. Olson, P. Leimgruber, and J. M. Calabrese. 2015. Rigorous home range estimation with movement data: a new autocorrelated kernel density estimator. Ecology 96:1182-1188. DOI: https://doi.org/10.1890/14-2010.1
Valle, D., Attias, N., Cullen, J. A., Hooten, M. B., Giroux, A., Oliveira-Santos, L. G. R., Desbiez, A. L. J. and Fletcher, R. J. 2024. Bridging the gap between movement data and connectivity analysis using the Time-Explicit Habitat Selection (TEHS) model. Movement Ecology 12:19. DOI: https://doi.org/10.1186/s40462-024-00461-1