diff --git a/.Rbuildignore b/.Rbuildignore index 52e21ae..9e49c05 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -12,3 +12,4 @@ ^LICENSE\.md$ ^README\.Rmd$ ^GeoTwilightR\.Rproj$ +^codemeta\.json$ diff --git a/AGENTS.md b/AGENTS.md index 60bcafb..221d33e 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -23,6 +23,7 @@ ## Performance +- Performance is a priority. - Prefer vectorised operations, matrix or array operations, `rowSums`, `colSums`, `sweep`, `%*%`, and `outer` where they fit. - Prefer efficient primitives over `apply()` when possible. Use `apply()` or `Map()` only when they remain clear and appropriate. - Avoid unnecessary copies, repeated computation, and hidden coercions. @@ -55,6 +56,20 @@ - Preserve numerical behaviour unless there is a clear risk such as `log(0)`, division by zero, or unstable normalisation. - When adding a numerical safeguard, keep it minimal and explain why it is needed. +## Error handling + +- Do not introduce `tryCatch()` or defensive error handling. +- Let errors occur naturally unless they are handled at a public entry point. + +## Silent failure risks + +Intervene only when there is a clear risk of unintended recycling, dimension dropping, implicit coercion, or invalid numerical operations. Otherwise, assume correct usage. + +## Rewriting vs patching + +- Default to a minimal patch. +- If code is clearly inconsistent, inefficient, or unclear, rewriting the function is allowed, but explicitly warn the user. + ## Output invariants Edits must preserve classes, dimensions, ordering, column names, and attributes. Do not silently change them. Explain any necessary trade-off before making it. @@ -70,6 +85,13 @@ Edits must preserve classes, dimensions, ordering, column names, and attributes. - Do not add or modify tests, roxygen documentation, README files, or vignettes unless requested. - Match validation to the change: source changed functions, run focused tests, and avoid long analyses unless needed. +## Generated files + +- Edit roxygen comments in `R/`, never `man/*.Rd` or `NAMESPACE`; run `roxygen2::roxygenise()` and commit the regenerated files. +- Edit `README.Rmd`, never `README.md`; run `devtools::build_readme()` and commit both. +- After editing the vignette or documentation, rebuild with `pkgdown::build_site()` or render the vignette to confirm it runs. +- Keep the software title, version, and DOI consistent across `DESCRIPTION`, `CITATION.cff`, `inst/CITATION`, `codemeta.json`, and `README.Rmd`. + ## Mandatory CI checks Before committing or pushing agent-made changes, run from the repository root: @@ -81,6 +103,46 @@ air format . --check Both commands must pass without warnings or errors introduced by the changes. If formatting fails, run `air format .`, review the resulting changes, and rerun both checks. Never defer failures to CI. +## Development, PR, and release workflow + +### Branch flow + +- Develop on `dev`: make focused commits, run the mandatory checks, and push to `origin/dev`. +- Open a pull request from `dev` to `main`. Merge `main` into `dev` and resolve conflicts before the PR is merged. +- Use a draft PR while release metadata or release checks are incomplete. Mark it ready only when the final package version, NEWS entry, and required GitHub checks are complete. + +### One canonical release block + +- For every release, write one Markdown release block first. It is the source of truth and must be copied **verbatim** to: + 1. the new top section of `NEWS.md`; + 2. the pull-request body; and + 3. the GitHub Release body created for the version tag. +- Do not shorten, paraphrase, reorder, or add items independently in any of those three places. +- Use the release version as the PR title (for example, `v0.2.0`) and as the top NEWS heading. +- Keep the heading, subsection headings, bullet text, Markdown links, and the full-changelog comparison link identical in all three copies. + +```md +# GeoPathSampleR vX.Y.Z + +## Main + +- [Describe the principal user-facing change](https://github.com/GeoPressure/GeoPathSampleR/commit/). + +## Minor + +- [Describe a smaller change or fix](https://github.com/GeoPressure/GeoPathSampleR/commit/). + +**Full Changelog**: +``` + +### Release checklist + +- Set `DESCRIPTION`, `CITATION.cff`, and `codemeta.json` to the final `X.Y.Z` version; do not merge a release with `.9000`. +- Run `cffr::cff_write()` after finalizing `DESCRIPTION`, then review and commit the generated `CITATION.cff` changes. Run `codemetar::write_codemeta()` to update `codemeta.json`. +- Add the canonical release block to `NEWS.md` before opening the PR, then paste that exact block into the PR body. +- Resolve all `R CMD check` warnings and release-relevant notes, and confirm the PR's GitHub Actions matrix is green. +- After merging to `main`, create tag `vX.Y.Z` and paste the unchanged canonical release block into the GitHub Release description. The Zenodo integration archives the release and takes its metadata from `CITATION.cff`. + ## Package constraints - Never use `:::` in package code. diff --git a/CITATION.cff b/CITATION.cff index 35d79cf..46fcc19 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -2,11 +2,11 @@ cff-version: 1.2.0 message: "To cite GeoPathSampleR in publications use:" type: software license: GPL-3.0-or-later -title: "GeoPathSampleR: Bayesian reconstruction of animal trajectories from geolocation likelihood maps" +title: "GeoPathSampleR: Posterior Sampling of Geolocator Paths" version: 0.1.0 date-released: 2026-08-15 doi: 10.5281/zenodo.21944194 -abstract: "Bayesian reconstruction of animal trajectories from geolocation likelihood maps." +abstract: "Reconstructs bird migration trajectories from stationary-period likelihood maps built with GeoPressureR. A Gibbs sampler combines residence-dependent departure, a ground-speed distribution conditional on movement, and a prior on cumulative route distance between long periods. Provides posterior stay summaries and sampling diagnostics." authors: - family-names: Nussbaumer given-names: Raphaël diff --git a/DESCRIPTION b/DESCRIPTION index 66e53e5..cd84e71 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,15 +1,16 @@ Package: GeoPathSampleR -Title: Bayesian Reconstruction of Animal Trajectories -Version: 0.1.0 +Title: Posterior Sampling of Geolocator Paths +Version: 0.1.0.9000 Authors@R: person("Raphaël", "Nussbaumer", , "rafnuss@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-8185-1020")) Author: Raphaël Nussbaumer [aut, cre] Maintainer: Raphaël Nussbaumer -Description: Bayesian reconstruction of animal trajectories from geolocation - likelihood maps using a Gibbs sampler with movement priors. Provides - prepared sampler inputs, posterior path summaries, and sampling - diagnostics for GeoPressureR tag objects. +Description: Reconstructs bird migration trajectories from stationary-period + likelihood maps built with GeoPressureR. A Gibbs sampler combines + residence-dependent departure, a ground-speed distribution conditional on + movement, and a prior on cumulative route distance between long periods. + Provides posterior stay summaries and sampling diagnostics. License: GPL (>= 3) + file LICENSE URL: https://github.com/GeoPressure/GeoPathSampleR, https://geopressure.org/GeoPathSampleR/ @@ -33,11 +34,14 @@ Suggests: leaflet, knitr, pkgdown, + plotly, quarto, rmarkdown, rnaturalearth, rnaturalearthdata, - testthat (>= 3.0.0) + scales, + testthat (>= 3.0.0), + tidyr Remotes: GeoPressure/GeoPressureR@dev Config/testthat/edition: 3 diff --git a/NAMESPACE b/NAMESPACE index d2d501d..f1728bd 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,6 +4,7 @@ export(path_collapse) export(path_summary) export(plot_movement) export(sampling_path) +export(sampling_path_default_movement) export(sampling_path_diagnostic) import(patchwork) importFrom(rlang,.data) diff --git a/NEWS.md b/NEWS.md index 81a0365..2b9f3d8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,20 @@ -# GeoPathSampleR 0.1.0 +# GeoPathSampleR v0.2.0 + +## Main + +- [Update the default ground-speed and residence-dependent departure parameters](https://github.com/GeoPressure/GeoPathSampleR/commit/d40999d) from the pooled multi-sensor calibration. Results sampled with the previous defaults should be recomputed before comparison. +- [Replace the route-detour model with a Gamma prior on route excess](https://github.com/GeoPressure/GeoPathSampleR/commit/7d4fb3f) whose scale depends on interval duration and direct endpoint distance ([3798978](https://github.com/GeoPressure/GeoPathSampleR/commit/3798978)). The prior applies to intervals with at least two movement days and endpoints at least 300 km apart, with predictors limited to the empirical calibration support. +- Export `sampling_path_default_movement()` and document the movement model, `route_detour` and the modular cut of `long_period_light_only` in `sampling_path()`. + +## Minor + +- [Fix the diagnostic report rendering and browser opening](https://github.com/GeoPressure/GeoPathSampleR/commit/fb0d61f), and [declare its optional report dependencies](https://github.com/GeoPressure/GeoPathSampleR/commit/3ff7a89). +- [Align the README and vignette with the movement model and sampler controls](https://github.com/GeoPressure/GeoPathSampleR/commit/ae3f0de), and add a provisional citation for the submitted methods manuscript alongside the software citation. +- [Format the elapsed time compactly in the sampling message](https://github.com/GeoPressure/GeoPathSampleR/commit/38fe161). + +**Full Changelog**: + +# GeoPathSampleR v0.1.0 * Initial public release of the stationary-path Gibbs sampler and its core summaries and diagnostics. diff --git a/R/path_collapse.R b/R/path_collapse.R index 69b5c03..fc00ae6 100644 --- a/R/path_collapse.R +++ b/R/path_collapse.R @@ -2,7 +2,9 @@ #' #' Summarise each realised path by merging consecutive stationary periods that #' were assigned to the same grid cell. Paths are identified by `j` and, when -#' present, `chain`. +#' present, `chain`. Each draw is collapsed separately, so the output keeps +#' one stay/move sequence per draw and all posterior variability; use +#' [path_summary()] to summarise across draws. #' #' @param paths A path data frame containing `stap_id`, `ind`, `lat`, and #' `lon`, and optionally `j`, `chain`, `start`, and `end`. diff --git a/R/path_summary.R b/R/path_summary.R index 764c7f9..9d075b6 100644 --- a/R/path_summary.R +++ b/R/path_summary.R @@ -1,9 +1,14 @@ #' Summarise an ensemble of paths #' -#' Compute either one posterior location per stationary period or one path -#' using the consensus stay/move structure across the ensemble. The function -#' accepts both [sampling_path()] output and simulation paths with `j` but no -#' `chain` column. +#' Summarise an ensemble of draws into a single path. With `by = "stap"`, the +#' positions of all draws are summarised separately for each stationary period. +#' With `by = "consensus_stay"`, the move probability of each boundary between +#' consecutive periods is the proportion of draws whose grid cell changes +#' there; boundaries above `move_threshold` split the track into consensus +#' stays, and positions are summarised over all draws and periods of each +#' stay. The function accepts both [sampling_path()] output and simulation +#' paths with `j` but no `chain` column. Use [path_collapse()] to keep each +#' draw instead. #' #' @param paths A path data frame containing `stap_id`, `ind`, `lat`, and #' `lon`, and optionally `j`, `chain`, `start`, and `end`. diff --git a/R/plot_movement.R b/R/plot_movement.R index f60839c..b8a3309 100644 --- a/R/plot_movement.R +++ b/R/plot_movement.R @@ -3,7 +3,8 @@ #' Plot the speed kernel and residence-dependent probability of movement for a #' `sampling_path()` movement specification. #' -#' @param movement movement list passed to `sampling_path()`. +#' @param movement movement list passed to `sampling_path()`, such as +#' [sampling_path_default_movement()]. #' @param speed speed values in km/h used to evaluate the movement kernel. #' @param stay_duration residence durations in days used to evaluate #' `movement$move_stay`. @@ -13,21 +14,7 @@ #' move/stay panel marks those parameters for the exponential-decay-to-plateau #' model. #' @examples -#' movement <- list( -#' method = "gamma", -#' shape = 1.26, -#' scale = 10.3, -#' low_speed_fix = 0.001, -#' zero_speed_ratio = 0, -#' move_stay_parameters = list( -#' p_0 = 0.607, -#' p_inf = 0.137, -#' tau = 0.932 -#' ), -#' move_stay = function(t) 0.137 + -#' (0.607 - 0.137) * exp(-t / 0.932) -#' ) -#' plot_movement(movement) +#' plot_movement(sampling_path_default_movement()) #' @import patchwork #' @importFrom rlang .data #' @family sampling_path diff --git a/R/sampling_path.R b/R/sampling_path.R index c814fbb..ebd3872 100644 --- a/R/sampling_path.R +++ b/R/sampling_path.R @@ -1,116 +1,52 @@ -#' Gibbs sampling for stationary paths +#' Default movement model for path sampling #' #' @description -#' Prepare the likelihood and movement inputs, run the Gibbs sampler for one or -#' more chains, and return the sampled path locations in long format. +#' Return the calibrated movement model used by default in [sampling_path()]. +#' Call `sampling_path_default_movement()` to inspect its values, plot it with +#' [plot_movement()], or modify some fields before passing it to +#' `sampling_path(movement = )`. #' -#' @param tag a GeoPressureR `tag` object containing the likelihood maps and -#' `tag_set_map` parameters. -#' @param movement one movement model applied across the complete track. The -#' parameter list must contain the fields passed to `speed2prob()` and a -#' `move_stay` function returning the probability of moving after a stay -#' duration. Users can optionally provide `approx_eps` to drop extremely small -#' movement weights from the precomputed neighbourhood (default `1e-6`). -#' @param iter number of Gibbs iterations. -#' @param chains number of chains to run. -#' @param warmup number of iterations to discard at the beginning. -#' @param thin thinning interval for saved samples. -#' @param refresh progress message update frequency, in iterations. -#' @param block_interval interval between block iterations. All remaining -#' iterations are site iterations that update each eligible period once. -#' Use `0` or `Inf` to disable block iterations. Values greater than one are -#' required so site iterations remain part of the chain. -#' @param long_period_light_only logical. `TRUE` (the default) draws unknown -#' `stap0` locations directly from their retained light likelihood and samples -#' the daily path conditionally without allowing movement to update them. -#' `FALSE` updates those locations with both light and movement information. -#' Tags without a `stap0` column have no long periods. -#' @param component_weights named non-negative weights for the light likelihood, -#' daily movement model, and interval route model. The default gives all three -#' components their calibrated weight. A zero removes that component's soft -#' contribution while retaining structural support constraints. -#' @param route_detour non-negative multiplier of the empirically predicted -#' excess route distance. `0` targets direct routes, `1` retains the empirical -#' prediction, and values above one favour more detoured routes. -#' @param route_model calibrated route-model parameter list. The default uses -#' the maintained duration-and-direct-distance calibration. Supply cross-fitted -#' parameters only for validation workflows. -#' @param likelihood field of the `tag` list containing the likelihood map. -#' @param thr_likelihood threshold of percentile to keep likely locations in -#' each stationary period. Lower values keep fewer nodes and speed up the -#' sampler, but exclude more of the light likelihood surface entirely. The -#' default `0.99` is a conservative state-space reduction; lower values such -#' as `0.95` are faster but more aggressive. -#' @param thr_gs maximum ground speed, in km/h, used as the hard radius of the -#' movement kernel. Larger values allow broader transitions and can improve -#' mixing when movements are uncertain, but they also increase computation. -#' Smaller values are faster but enforce a tighter movement radius. -#' @param workers number of parallel workers used for independent chains. -#' Use 1 for sequential execution. -#' @param seed optional RNG seed for reproducible chains. -#' @param quiet logical to hide messages. +#' @details +#' The movement model defines the daily transition between consecutive +#' stationary periods. It is a list with two parts. #' -#' @return A data.frame with columns `j`, `chain`, `stap_id`, `ind`, `lat`, and -#' `lon`. The result carries `type = "sampling"` and the scalar sampler -#' settings in a `sampling_parameters` attribute for reproducible diagnostics. +#' **Ground-speed kernel.** The fields `method`, `shape`, `scale`, +#' `low_speed_fix` and `zero_speed_ratio` are passed to +#' [GeoPressureR::speed2prob()] and follow the convention of +#' [GeoPressureR::graph_set_movement()]. Speed is the great-circle displacement +#' between two consecutive periods divided by 24 h, in km/h. With +#' `method = "gamma"`, the speed density is a Gamma distribution with shape +#' `shape` and scale `scale` (km/h). The probability of moving to a given cell +#' is this density divided by the distance and multiplied by `cos(latitude)`, +#' which accounts for the number and area of cells at that distance. The hard +#' movement radius is set separately by `sampling_path(thr_gs = )`. #' -#' @details -#' `sampling_path()` separates structural constraints from computational -#' approximations. Impossible geography, such as water after masking, is removed -#' from the state space. The arguments `thr_likelihood`, `thr_gs`, -#' and `movement$approx_eps` control how aggressively the sampler reduces the -#' fixed light support and movement kernel for speed. +#' **Residence-dependent departure.** `move_stay` is a function of the +#' residence `t`, the number of consecutive transitions spent at the current +#' location (`0` for the first transition after arrival), returning the +#' probability of departing at the next transition. The default is the +#' exponential decline to a plateau +#' `p_inf + (p_0 - p_inf) * exp(-t / tau)`. `move_stay_parameters` stores +#' `p_0`, `p_inf` and `tau` so that [plot_movement()] can annotate them; it is +#' not used by the sampler, so update `move_stay` when changing them. #' -#' A practical tuning strategy is: -#' - start with the defaults for final inference; -#' - use a smaller `thr_likelihood` only for explicitly approximate exploratory -#' runs; -#' - if no movement-feasible path exists, increase `thr_likelihood` or -#' `thr_gs`. +#' The optional field `approx_eps` (default `1e-6`) drops movement-kernel +#' weights below this value from the precomputed neighbourhood. #' -#' @examples -#' example_dir <- system.file("extdata", package = "GeoPathSampleR") -#' tag <- GeoPressureR::tag_create( -#' "14OI", -#' crop_start = "2015-07-17", -#' crop_end = "2016-07-11", -#' directory = file.path(example_dir, "data/raw-tag/14OI"), -#' assert_pressure = FALSE, -#' quiet = TRUE -#' ) -#' tag <- GeoPressureR::twilight_create(tag) -#' tag <- GeoPressureR::twilight_label_read( -#' tag, -#' file = file.path(example_dir, "data/twilight-label/14OI-labeled.csv") -#' ) -#' tag <- GeoPressureR::tag_stap_daily( -#' tag, -#' stap_long = file.path(example_dir, "data/stap-label/14OI.csv"), -#' movement_period = "day", -#' quiet = TRUE -#' ) -#' tag <- GeoPressureR::tag_set_map( -#' tag, -#' extent = c(-20, 40, -40, 60), -#' scale = 1 -#' ) -#' tag <- GeoPressureR::geolight_map( -#' tag, -#' twl_calib_adjust = 1, -#' fitted_location_duration = 30, -#' twl_llp = \(n) 1.5 * log(n) / n, -#' quiet = TRUE -#' ) +#' The default values were calibrated on multi-sensor reference trajectories +#' of land birds (Gamma shape `1.26` and scale `10.3` km/h; `p_0 = 0.607`, +#' `p_inf = 0.137`, `tau = 0.932`). They should not be assumed to transfer to +#' other taxa. #' -#' paths <- sampling_path( -#' tag, -#' iter = 3, -#' warmup = 1, -#' seed = 1, -#' quiet = TRUE -#' ) +#' @return A movement list with fields `method`, `shape`, `scale`, +#' `low_speed_fix`, `zero_speed_ratio`, `move_stay_parameters`, and +#' `move_stay`. +#' @examples +#' movement <- sampling_path_default_movement() +#' str(movement) +#' plot_movement(movement) #' @family sampling_path -#' @noRd +#' @export sampling_path_default_movement <- function() { list( method = "gamma", @@ -480,36 +416,153 @@ sampling_path_with_route_model <- function( #' #' @description #' Sample stationary trajectories from light likelihood maps and calibrated -#' movement components. +#' movement components. The posterior combines three components: the +#' stationary-period light likelihood, a daily movement model (see +#' [sampling_path_default_movement()]), and a route-distance prior on the +#' interval between consecutive long periods (`stap0`). #' -#' @param tag GeoPressureR tag object containing likelihood maps. -#' @param iter number of Gibbs iterations. +#' @param tag GeoPressureR tag object containing likelihood maps and +#' `tag_set_map` parameters. Long periods are identified by `tag$stap$stap0`, +#' as created by [GeoPressureR::tag_stap_daily()]. +#' @param iter number of Gibbs iterations per chain, including `warmup`. #' @param likelihood tag field containing the light likelihood map. -#' @param movement movement model containing a speed kernel and stay/move -#' probabilities. -#' @param component_weights named non-negative weights for the light, -#' movement, and route components. -#' @param route_detour non-negative multiplier of the empirical excess route -#' distance; `1` retains the calibrated default. -#' @param long_period_light_only whether unknown long periods are drawn from -#' their light likelihood without movement updates. +#' @param movement daily movement model: a ground-speed kernel passed to +#' [GeoPressureR::speed2prob()] and a `move_stay` function returning the +#' probability of departing after a residence time. The default, +#' [sampling_path_default_movement()], is the calibrated land-bird model; +#' call it to inspect the values and [plot_movement()] to plot them. +#' @param component_weights named non-negative weights for the `light`, +#' `movement`, and `route` components. Each weight multiplies that +#' component's log probability, so the default `1` uses the calibrated +#' components unchanged. A zero removes that component's soft contribution +#' while retaining hard support constraints (retained light support, water +#' mask and `thr_gs`). +#' @param route_detour non-negative multiplier of the scale of the Gamma prior +#' on route excess (see Details). `1` retains the calibrated prior, `0` +#' targets direct routes, and values above one favour more detoured routes. +#' @param long_period_light_only logical. `TRUE` applies a modular cut: at every +#' iteration, each unknown long-period location is redrawn from its light +#' likelihood alone, and only daily periods are updated with the movement and +#' route components. `FALSE` updates long periods with light, movement and +#' route information like daily periods. See Details. #' @param chains number of independent chains. #' @param warmup number of initial iterations to discard. #' @param thin interval between saved samples. -#' @param block_interval interval between block updates; use `0` to disable -#' them. -#' @param thr_likelihood retained light-likelihood percentile. -#' @param thr_gs maximum movement speed in km/h used for hard support. +#' @param block_interval interval between block iterations, which update each +#' current residence segment as a single block. All other iterations are +#' site iterations that update each eligible period once. Use `0` or `Inf` to +#' disable block iterations; values must otherwise be at least `2`. +#' @param thr_likelihood retained light-likelihood mass in each stationary +#' period. Lower values keep fewer cells and speed up the sampler, but exclude +#' more of the light likelihood surface entirely. The default `0.99` is a +#' conservative state-space reduction. +#' @param thr_gs maximum ground speed, in km/h, used as the hard radius of the +#' daily movement kernel (default `2000 / 24`, i.e. 2000 km per day). Speeds +#' below it are weighted by the `movement` kernel, so `thr_gs` is a +#' truncation of that kernel rather than a typical speed: set it above the +#' speeds the kernel supports. The default excludes less than 0.1% of the +#' default kernel. Larger values increase computation; smaller values make +#' faster transitions impossible rather than unlikely. #' @param refresh progress update interval in iterations. #' @param workers number of parallel workers for independent chains. -#' @param seed optional random seed. +#' Use `1` for sequential execution. +#' @param seed optional random seed for reproducible chains. #' @param quiet whether to suppress progress messages. #' @return A data.frame with columns `j`, `chain`, `stap_id`, `ind`, `lat`, and #' `lon`. The result carries `type = "sampling"` and sampler settings in a #' `sampling_parameters` attribute. +#' +#' @details +#' **Route-distance prior.** For each interval between consecutive long +#' periods \eqn{\ell}, the route ratio \eqn{R_\ell} is the sum of the daily +#' great-circle displacements divided by the direct distance \eqn{D_\ell} +#' between the two long-period locations. The sampler adds +#' `component_weights["route"]` times the log density +#' \deqn{\log \mathrm{Gamma}(R_\ell - 1 \mid \alpha_R,\ \mathrm{route\_detour} \times \theta_\ell),}{ +#' log dgamma(R - 1, shape = alpha_R, scale = route_detour * theta),} +#' with shape \eqn{\alpha_R = 1.22} and +#' \deqn{\log \theta_\ell = -1.88 + 0.552\,[\log(1 + T_\ell) - 3.47] - 0.318\,[\log D_\ell - 8.15],}{ +#' log(theta) = -1.88 + 0.552 * (log(1 + T) - 3.47) - 0.318 * (log(D) - 8.15),} +#' where \eqn{T_\ell} is the interval duration in days and \eqn{D_\ell} is in +#' km. Both covariates are clamped to the range of the calibration data. The +#' term is omitted when \eqn{D_\ell < 300} km or when the interval contains at +#' most one daily displacement. +#' +#' Because `route_detour` multiplies the Gamma scale, the prior mean and every +#' prior quantile of the route excess \eqn{R_\ell - 1} scale linearly with it. +#' For example, an interval of about 31 days with endpoints about 3,500 km +#' apart has \eqn{\theta_\ell \approx 0.15}, i.e. a prior mean route ratio of +#' about 1.19; `route_detour = 2` raises it to about 1.37 and `0.5` lowers it +#' to about 1.09. `route_detour = 0` collapses the scale, so any detour is +#' strongly penalised. The prior remains soft: the realised detour of the +#' posterior also depends on the light and movement components. +#' +#' **Long periods and modular cut.** Twilights accurately locate long periods +#' but leave daily-period locations ambiguous. With +#' `long_period_light_only = TRUE`, each unknown long-period location is +#' redrawn independently at every iteration from its retained light +#' likelihood (raised to `component_weights["light"]`), without using the +#' movement or route components. When a redraw violates hard movement support, +#' only the neighbouring daily locations needed to restore a feasible path are +#' resampled. Daily periods are then updated conditionally on these long-period +#' draws. This is a modular (cut) approximation to the joint posterior: +#' information flows from the long-period likelihoods to the daily path, but +#' the movement model cannot pull long-period locations away from their light +#' likelihood. Known-location periods are always fixed. +#' +#' **Computational approximations.** Impossible geography, such as water after +#' masking, is removed from the state space. The arguments `thr_likelihood`, +#' `thr_gs` and `movement$approx_eps` control how aggressively the sampler +#' reduces the fixed light support and movement kernel for speed. Start with +#' the defaults for final inference, use a smaller `thr_likelihood` only for +#' explicitly approximate exploratory runs, and increase `thr_likelihood` or +#' `thr_gs` if no movement-feasible path exists. +#' #' @examples #' \dontrun{ -#' paths <- sampling_path(tag, iter = 1000, chains = 4, seed = 1) +#' example_dir <- system.file("extdata", package = "GeoPathSampleR") +#' tag <- GeoPressureR::tag_create( +#' "14OI", +#' crop_start = "2015-07-17", +#' crop_end = "2016-07-11", +#' directory = file.path(example_dir, "data/raw-tag/14OI"), +#' assert_pressure = FALSE, +#' quiet = TRUE +#' ) +#' tag <- GeoPressureR::twilight_create(tag) +#' tag <- GeoPressureR::twilight_label_read( +#' tag, +#' file = file.path(example_dir, "data/twilight-label/14OI-labeled.csv") +#' ) +#' tag <- GeoPressureR::tag_stap_daily( +#' tag, +#' stap_long = file.path(example_dir, "data/stap-label/14OI.csv"), +#' movement_period = "day", +#' quiet = TRUE +#' ) +#' tag <- GeoPressureR::tag_set_map( +#' tag, +#' extent = c(-20, 40, -11, 60), +#' scale = 1 +#' ) +#' tag <- GeoPressureR::geolight_map( +#' tag, +#' twl_calib_adjust = 1, +#' fitted_location_duration = 30, +#' twl_llp = \(n) 1.5 * log(n) / n, +#' quiet = TRUE +#' ) +#' tag$map_light <- GeoPressureR::map_add_mask_water(tag$map_light) +#' +#' paths <- sampling_path( +#' tag, +#' movement = sampling_path_default_movement(), +#' iter = 500, +#' chains = 4, +#' warmup = 100, +#' seed = 1, +#' quiet = TRUE +#' ) #' } #' @family sampling_path #' @export @@ -570,10 +623,12 @@ sampling_path_route_model <- function() { min_direct_distance_km = 300, support = list( duration_log = seq(0.731, 4.907, length.out = 16), + # fmt: skip log_distance_lower = c( 5.749, 5.777, 5.786, 5.887, 6.064, 6.272, 6.359, 6.040, 6.266, 6.658, 6.723, 7.180, 7.109, 7.287, 7.285, 7.273 ), + # fmt: skip log_distance_upper = c( 7.315, 7.448, 7.634, 7.696, 8.028, 8.380, 8.417, 8.437, 8.513, 8.558, 8.674, 9.022, 9.143, 9.343, 9.419, 9.440 @@ -582,7 +637,6 @@ sampling_path_route_model <- function() { ) } -#' Prepare Span-Level Route Prior sampling_path_run_chain <- function( chain_id, lk, @@ -1342,7 +1396,7 @@ sampling_path_weight_log_prob <- function(log_prob, weight) { out } - +#' Prepare Span-Level Route Prior #' #' @return A route-prior list with sampler row indices, or `NULL`. #' @noRd diff --git a/R/sampling_path_diagnostic.R b/R/sampling_path_diagnostic.R index 74ccd36..3253cc9 100644 --- a/R/sampling_path_diagnostic.R +++ b/R/sampling_path_diagnostic.R @@ -488,17 +488,33 @@ sampling_path_diagnostic_render <- function( browse = TRUE, output_file = NULL ) { - if (is.null(output_file)) { - report_dir <- tempfile("sampling_path_diagnostic_") - dir.create(report_dir) - report_file <- file.path(report_dir, "sampling_path_diagnostic.html") - on.exit(unlink(report_file), add = TRUE) - } else { - report_file <- normalizePath(output_file, mustWork = FALSE) - dir.create(dirname(report_file), recursive = TRUE, showWarnings = FALSE) + quarto <- Sys.getenv("QUARTO_PATH", unset = unname(Sys.which("quarto"))) + if (!nzchar(quarto)) { + cli::cli_abort(c( + "Quarto is required to render the sampling-path diagnostic report.", + "i" = "Install Quarto ({.url https://quarto.org/docs/get-started/}) or set {.envvar QUARTO_PATH}." + )) } - diagnostic_file <- tempfile(fileext = ".rds") - params_file <- tempfile(fileext = ".json") + template <- system.file( + "report", + "sampling_path_diagnostic.qmd", + package = "GeoPathSampleR" + ) + if (!nzchar(template)) { + cli::cli_abort( + "The GeoPathSampleR diagnostic-report template is unavailable." + ) + } + + # Render in a private directory: the installed package directory is read-only + # on many systems and shared between concurrent renders. + render_dir <- tempfile("sampling_path_diagnostic_") + dir.create(render_dir) + on.exit(unlink(render_dir, recursive = TRUE), add = TRUE) + file.copy(template, file.path(render_dir, basename(template))) + diagnostic_file <- file.path(render_dir, "diagnostic.rds") + params_file <- file.path(render_dir, "params.json") + render_name <- "sampling_path_diagnostic.html" saveRDS( list( diagnostic = diagnostic, @@ -508,36 +524,22 @@ sampling_path_diagnostic_render <- function( ), diagnostic_file ) - on.exit(unlink(diagnostic_file), add = TRUE) - on.exit(unlink(params_file), add = TRUE) jsonlite::write_json( - list(diagnostic_file = diagnostic_file), + list( + diagnostic_file = normalizePath( + diagnostic_file, + winslash = "/", + mustWork = TRUE + ) + ), params_file, auto_unbox = TRUE ) - template <- system.file( - "report", - "sampling_path_diagnostic.qmd", - package = "GeoPathSampleR" - ) - if (!nzchar(template)) { - cli::cli_abort( - "The GeoPathSampleR diagnostic-report template is unavailable." - ) - } - render_name <- paste0("sampling_path_diagnostic_", Sys.getpid(), ".html") - rendered_file <- file.path(dirname(template), render_name) - resource_dir <- file.path( - dirname(template), - paste0(tools::file_path_sans_ext(render_name), "_files") - ) - on.exit(unlink(rendered_file), add = TRUE) - on.exit(unlink(resource_dir, recursive = TRUE), add = TRUE) render_output <- withr::with_dir( - dirname(template), + render_dir, system2( - "quarto", + quarto, c( "render", shQuote(basename(template)), @@ -546,7 +548,7 @@ sampling_path_diagnostic_render <- function( "--output", shQuote(render_name), "--execute-params", - shQuote(params_file) + shQuote(basename(params_file)) ), stdout = TRUE, stderr = TRUE @@ -559,13 +561,25 @@ sampling_path_diagnostic_render <- function( "x" = paste(utils::tail(render_output, 8), collapse = "\n") )) } - if (!file.copy(rendered_file, report_file, overwrite = TRUE)) { - cli::cli_abort( - "Failed to open the rendered sampling-path diagnostic report." + + # The report must outlive this call: browsers open it asynchronously. + report_file <- if (is.null(output_file)) { + tempfile("sampling_path_diagnostic_", fileext = ".html") + } else { + dir.create(dirname(output_file), recursive = TRUE, showWarnings = FALSE) + normalizePath(output_file, mustWork = FALSE) + } + if ( + !file.copy( + file.path(render_dir, render_name), + report_file, + overwrite = TRUE ) + ) { + cli::cli_abort("Failed to save the sampling-path diagnostic report.") } if (browse) { - utils::browseURL(report_file) + utils::browseURL(normalizePath(report_file, mustWork = TRUE)) } invisible(report_file) diff --git a/README.Rmd b/README.Rmd index 8121da0..33516e4 100644 --- a/README.Rmd +++ b/README.Rmd @@ -19,7 +19,12 @@ knitr::opts_chunk$set(collapse = TRUE, comment = "#>") [![format](https://github.com/GeoPressure/GeoPathSampleR/actions/workflows/format-check.yml/badge.svg)](https://github.com/GeoPressure/GeoPathSampleR/actions/workflows/format-check.yml) -GeoPathSampleR provides Bayesian reconstruction of animal trajectories from geolocation likelihood maps. It includes a Gibbs sampler, reusable prepared inputs, posterior stay summaries, and diagnostics for GeoPressureR tag objects. +GeoPathSampleR reconstructs bird migration trajectories from stationary-period +likelihood maps built with GeoPressureR. Its Gibbs sampler combines +residence-dependent departure, a ground-speed distribution conditional on +movement, and a prior on cumulative route distance between long periods. +Unknown long-period locations can be sampled from light evidence alone. The +package also provides posterior stay summaries and sampling diagnostics. ## 📦 Installation @@ -37,6 +42,11 @@ for inputs, movement settings, and interpretation. ## 📚 Citation -If you use GeoPathSampleR in your research, cite: +For the software, cite: -> Nussbaumer, R. (2026). *GeoPathSampleR: Bayesian reconstruction of animal trajectories from geolocation likelihood maps.* Zenodo. [https://doi.org/10.5281/zenodo.21944194](https://doi.org/10.5281/zenodo.21944194). Available at: [https://github.com/GeoPressure/GeoPathSampleR](https://github.com/GeoPressure/GeoPathSampleR) +> Nussbaumer, R. (2026). *GeoPathSampleR: Posterior Sampling of Geolocator Paths.* Zenodo. [https://doi.org/10.5281/zenodo.21944194](https://doi.org/10.5281/zenodo.21944194). Available at: [https://github.com/GeoPressure/GeoPathSampleR](https://github.com/GeoPressure/GeoPathSampleR) + +The methods and validation are described in the following submitted manuscript +(citation details will be updated after publication): + +> Nussbaumer, R., Gravey, M., Wong, J. B. & Benoit, L. (2026). Reconstructing bird migration trajectories from light: inspectable likelihoods and biologically informed movement priors. Manuscript submitted. diff --git a/README.md b/README.md index b57610c..a9ca657 100644 --- a/README.md +++ b/README.md @@ -13,10 +13,13 @@ [![format](https://github.com/GeoPressure/GeoPathSampleR/actions/workflows/format-check.yml/badge.svg)](https://github.com/GeoPressure/GeoPathSampleR/actions/workflows/format-check.yml) -GeoPathSampleR provides Bayesian reconstruction of animal trajectories -from geolocation likelihood maps. It includes a Gibbs sampler, reusable -prepared inputs, posterior stay summaries, and diagnostics for -GeoPressureR tag objects. +GeoPathSampleR reconstructs bird migration trajectories from +stationary-period likelihood maps built with GeoPressureR. Its Gibbs +sampler combines residence-dependent departure, a ground-speed +distribution conditional on movement, and a prior on cumulative route +distance between long periods. Unknown long-period locations can be +sampled from light evidence alone. The package also provides posterior +stay summaries and sampling diagnostics. ## 📦 Installation @@ -35,9 +38,16 @@ for inputs, movement settings, and interpretation. ## 📚 Citation -If you use GeoPathSampleR in your research, cite: +For the software, cite: -> Nussbaumer, R. (2026). *GeoPathSampleR: Bayesian reconstruction of -> animal trajectories from geolocation likelihood maps.* Zenodo. -> . Available at: -> +> Nussbaumer, R. (2026). *GeoPathSampleR: Posterior Sampling of +> Geolocator Paths.* Zenodo. . +> Available at: + +The methods and validation are described in the following submitted +manuscript (citation details will be updated after publication): + +> Nussbaumer, R., Gravey, M., Wong, J. B. & Benoit, L. (2026). +> Reconstructing bird migration trajectories from light: inspectable +> likelihoods and biologically informed movement priors. Manuscript +> submitted. diff --git a/_pkgdown.yml b/_pkgdown.yml index 0d988e7..178c7bd 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -19,6 +19,7 @@ reference: - title: Path sampling contents: - sampling_path + - sampling_path_default_movement - path_collapse - path_summary - sampling_path_diagnostic diff --git a/codemeta.json b/codemeta.json new file mode 100644 index 0000000..f0ce151 --- /dev/null +++ b/codemeta.json @@ -0,0 +1,341 @@ +{ + "@context": "https://doi.org/10.5063/schema/codemeta-2.0", + "@type": "SoftwareSourceCode", + "identifier": "GeoPathSampleR", + "description": "Reconstructs bird migration trajectories from stationary-period likelihood maps built with GeoPressureR. A Gibbs sampler combines residence-dependent departure, a ground-speed distribution conditional on movement, and a prior on cumulative route distance between long periods. Provides posterior stay summaries and sampling diagnostics.", + "name": "GeoPathSampleR: Posterior Sampling of Geolocator Paths", + "relatedLink": "https://geopressure.org/GeoPathSampleR/", + "codeRepository": "https://github.com/GeoPressure/GeoPathSampleR", + "issueTracker": "https://github.com/GeoPressure/GeoPathSampleR/issues", + "license": "https://spdx.org/licenses/GPL-3.0", + "version": "0.1.0.9000", + "programmingLanguage": { + "@type": "ComputerLanguage", + "name": "R", + "url": "https://r-project.org" + }, + "runtimePlatform": "R version 4.6.1 (2026-06-24)", + "author": [ + { + "@type": "Person", + "givenName": "Raphaël", + "familyName": "Nussbaumer", + "email": "rafnuss@gmail.com", + "@id": "https://orcid.org/0000-0002-8185-1020" + } + ], + "maintainer": [ + { + "@type": "Person", + "givenName": "Raphaël", + "familyName": "Nussbaumer", + "email": "rafnuss@gmail.com", + "@id": "https://orcid.org/0000-0002-8185-1020" + } + ], + "softwareSuggestions": [ + { + "@type": "SoftwareApplication", + "identifier": "dplyr", + "name": "dplyr", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=dplyr" + }, + { + "@type": "SoftwareApplication", + "identifier": "leaflet", + "name": "leaflet", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=leaflet" + }, + { + "@type": "SoftwareApplication", + "identifier": "knitr", + "name": "knitr", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=knitr" + }, + { + "@type": "SoftwareApplication", + "identifier": "pkgdown", + "name": "pkgdown", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=pkgdown" + }, + { + "@type": "SoftwareApplication", + "identifier": "plotly", + "name": "plotly", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=plotly" + }, + { + "@type": "SoftwareApplication", + "identifier": "quarto", + "name": "quarto", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=quarto" + }, + { + "@type": "SoftwareApplication", + "identifier": "rmarkdown", + "name": "rmarkdown", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=rmarkdown" + }, + { + "@type": "SoftwareApplication", + "identifier": "rnaturalearth", + "name": "rnaturalearth", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=rnaturalearth" + }, + { + "@type": "SoftwareApplication", + "identifier": "rnaturalearthdata", + "name": "rnaturalearthdata", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=rnaturalearthdata" + }, + { + "@type": "SoftwareApplication", + "identifier": "scales", + "name": "scales", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=scales" + }, + { + "@type": "SoftwareApplication", + "identifier": "testthat", + "name": "testthat", + "version": ">= 3.0.0", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=testthat" + }, + { + "@type": "SoftwareApplication", + "identifier": "tidyr", + "name": "tidyr", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=tidyr" + } + ], + "softwareRequirements": { + "1": { + "@type": "SoftwareApplication", + "identifier": "R", + "name": "R", + "version": ">= 4.1.0" + }, + "2": { + "@type": "SoftwareApplication", + "identifier": "assertthat", + "name": "assertthat", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=assertthat" + }, + "3": { + "@type": "SoftwareApplication", + "identifier": "cli", + "name": "cli", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=cli" + }, + "4": { + "@type": "SoftwareApplication", + "identifier": "GeoPressureR", + "name": "GeoPressureR", + "version": ">= 3.6.0" + }, + "5": { + "@type": "SoftwareApplication", + "identifier": "ggplot2", + "name": "ggplot2", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=ggplot2" + }, + "6": { + "@type": "SoftwareApplication", + "identifier": "glue", + "name": "glue", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=glue" + }, + "7": { + "@type": "SoftwareApplication", + "identifier": "jsonlite", + "name": "jsonlite", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=jsonlite" + }, + "8": { + "@type": "SoftwareApplication", + "identifier": "patchwork", + "name": "patchwork", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=patchwork" + }, + "9": { + "@type": "SoftwareApplication", + "identifier": "posterior", + "name": "posterior", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=posterior" + }, + "10": { + "@type": "SoftwareApplication", + "identifier": "rlang", + "name": "rlang", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=rlang" + }, + "11": { + "@type": "SoftwareApplication", + "identifier": "tibble", + "name": "tibble", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=tibble" + }, + "12": { + "@type": "SoftwareApplication", + "identifier": "withr", + "name": "withr", + "provider": { + "@id": "https://cran.r-project.org", + "@type": "Organization", + "name": "Comprehensive R Archive Network (CRAN)", + "url": "https://cran.r-project.org" + }, + "sameAs": "https://CRAN.R-project.org/package=withr" + }, + "SystemRequirements": null + }, + "fileSize": "52570.547KB", + "citation": [ + { + "@type": "SoftwareSourceCode", + "datePublished": "2026", + "author": [ + { + "@type": "Person", + "givenName": "Raphaël", + "familyName": "Nussbaumer", + "email": "rafnuss@gmail.com" + } + ], + "name": "GeoPathSampleR: Posterior Sampling of Geolocator Paths", + "identifier": "10.5281/zenodo.21944194", + "url": "https://github.com/GeoPressure/GeoPathSampleR", + "@id": "https://doi.org/10.5281/zenodo.21944194", + "sameAs": "https://doi.org/10.5281/zenodo.21944194" + } + ], + "releaseNotes": "https://github.com/GeoPressure/GeoPathSampleR/blob/main/NEWS.md", + "readme": "https://github.com/GeoPressure/GeoPathSampleR/blob/main/README.md", + "contIntegration": ["https://github.com/GeoPressure/GeoPathSampleR/actions/workflows/R-CMD-check.yaml", "https://github.com/GeoPressure/GeoPathSampleR/actions/workflows/pkgdown.yaml", "https://app.codecov.io/gh/GeoPressure/GeoPathSampleR", "https://github.com/GeoPressure/GeoPathSampleR/actions/workflows/jarl-check.yml", "https://github.com/GeoPressure/GeoPathSampleR/actions/workflows/format-check.yml"] +} diff --git a/inst/CITATION b/inst/CITATION index 644544b..3043bb7 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -2,7 +2,7 @@ citHeader("To cite GeoPathSampleR in publications:") citEntry( entry = "Manual", - title = "GeoPathSampleR: Bayesian reconstruction of animal trajectories from geolocation likelihood maps", + title = "GeoPathSampleR: Posterior Sampling of Geolocator Paths", author = person("Raphaël", "Nussbaumer", email = "rafnuss@gmail.com", role = c("aut", "cre")), year = "2026", doi = "10.5281/zenodo.21944194", diff --git a/man/path_collapse.Rd b/man/path_collapse.Rd index f352d85..6bf9bf6 100644 --- a/man/path_collapse.Rd +++ b/man/path_collapse.Rd @@ -25,7 +25,9 @@ A data frame with one row per collapsed stay and columns \description{ Summarise each realised path by merging consecutive stationary periods that were assigned to the same grid cell. Paths are identified by \code{j} and, when -present, \code{chain}. +present, \code{chain}. Each draw is collapsed separately, so the output keeps +one stay/move sequence per draw and all posterior variability; use +\code{\link[=path_summary]{path_summary()}} to summarise across draws. } \seealso{ Other path: diff --git a/man/path_summary.Rd b/man/path_summary.Rd index 57c6c12..ccec308 100644 --- a/man/path_summary.Rd +++ b/man/path_summary.Rd @@ -36,10 +36,15 @@ A data frame with one row per requested summary unit and columns probabilities at the start and end boundaries. } \description{ -Compute either one posterior location per stationary period or one path -using the consensus stay/move structure across the ensemble. The function -accepts both \code{\link[=sampling_path]{sampling_path()}} output and simulation paths with \code{j} but no -\code{chain} column. +Summarise an ensemble of draws into a single path. With \code{by = "stap"}, the +positions of all draws are summarised separately for each stationary period. +With \code{by = "consensus_stay"}, the move probability of each boundary between +consecutive periods is the proportion of draws whose grid cell changes +there; boundaries above \code{move_threshold} split the track into consensus +stays, and positions are summarised over all draws and periods of each +stay. The function accepts both \code{\link[=sampling_path]{sampling_path()}} output and simulation +paths with \code{j} but no \code{chain} column. Use \code{\link[=path_collapse]{path_collapse()}} to keep each +draw instead. } \seealso{ Other path: diff --git a/man/plot_movement.Rd b/man/plot_movement.Rd index 922942a..d6d62ec 100644 --- a/man/plot_movement.Rd +++ b/man/plot_movement.Rd @@ -7,7 +7,8 @@ plot_movement(movement, speed = seq(0, 120), stay_duration = seq(0, 20)) } \arguments{ -\item{movement}{movement list passed to \code{sampling_path()}.} +\item{movement}{movement list passed to \code{sampling_path()}, such as +\code{\link[=sampling_path_default_movement]{sampling_path_default_movement()}}.} \item{speed}{speed values in km/h used to evaluate the movement kernel.} @@ -25,25 +26,12 @@ Plot the speed kernel and residence-dependent probability of movement for a \code{sampling_path()} movement specification. } \examples{ -movement <- list( - method = "gamma", - shape = 1.26, - scale = 10.3, - low_speed_fix = 0.001, - zero_speed_ratio = 0, - move_stay_parameters = list( - p_0 = 0.607, - p_inf = 0.137, - tau = 0.932 - ), - move_stay = function(t) 0.137 + - (0.607 - 0.137) * exp(-t / 0.932) -) -plot_movement(movement) +plot_movement(sampling_path_default_movement()) } \seealso{ Other sampling_path: \code{\link[=sampling_path]{sampling_path()}}, +\code{\link[=sampling_path_default_movement]{sampling_path_default_movement()}}, \code{\link[=sampling_path_diagnostic]{sampling_path_diagnostic()}} } \concept{sampling_path} diff --git a/man/sampling_path.Rd b/man/sampling_path.Rd index 38f10f8..6a5dd22 100644 --- a/man/sampling_path.Rd +++ b/man/sampling_path.Rd @@ -25,23 +25,36 @@ sampling_path( ) } \arguments{ -\item{tag}{GeoPressureR tag object containing likelihood maps.} +\item{tag}{GeoPressureR tag object containing likelihood maps and +\code{tag_set_map} parameters. Long periods are identified by \code{tag$stap$stap0}, +as created by \code{\link[GeoPressureR:tag_stap_daily]{GeoPressureR::tag_stap_daily()}}.} -\item{iter}{number of Gibbs iterations.} +\item{iter}{number of Gibbs iterations per chain, including \code{warmup}.} \item{likelihood}{tag field containing the light likelihood map.} -\item{movement}{movement model containing a speed kernel and stay/move -probabilities.} - -\item{component_weights}{named non-negative weights for the light, -movement, and route components.} - -\item{route_detour}{non-negative multiplier of the empirical excess route -distance; \code{1} retains the calibrated default.} - -\item{long_period_light_only}{whether unknown long periods are drawn from -their light likelihood without movement updates.} +\item{movement}{daily movement model: a ground-speed kernel passed to +\code{\link[GeoPressureR:speed2prob]{GeoPressureR::speed2prob()}} and a \code{move_stay} function returning the +probability of departing after a residence time. The default, +\code{\link[=sampling_path_default_movement]{sampling_path_default_movement()}}, is the calibrated land-bird model; +call it to inspect the values and \code{\link[=plot_movement]{plot_movement()}} to plot them.} + +\item{component_weights}{named non-negative weights for the \code{light}, +\code{movement}, and \code{route} components. Each weight multiplies that +component's log probability, so the default \code{1} uses the calibrated +components unchanged. A zero removes that component's soft contribution +while retaining hard support constraints (retained light support, water +mask and \code{thr_gs}).} + +\item{route_detour}{non-negative multiplier of the scale of the Gamma prior +on route excess (see Details). \code{1} retains the calibrated prior, \code{0} +targets direct routes, and values above one favour more detoured routes.} + +\item{long_period_light_only}{logical. \code{TRUE} applies a modular cut: at every +iteration, each unknown long-period location is redrawn from its light +likelihood alone, and only daily periods are updated with the movement and +route components. \code{FALSE} updates long periods with light, movement and +route information like daily periods. See Details.} \item{chains}{number of independent chains.} @@ -49,18 +62,30 @@ their light likelihood without movement updates.} \item{thin}{interval between saved samples.} -\item{block_interval}{interval between block updates; use \code{0} to disable -them.} +\item{block_interval}{interval between block iterations, which update each +current residence segment as a single block. All other iterations are +site iterations that update each eligible period once. Use \code{0} or \code{Inf} to +disable block iterations; values must otherwise be at least \code{2}.} -\item{thr_likelihood}{retained light-likelihood percentile.} +\item{thr_likelihood}{retained light-likelihood mass in each stationary +period. Lower values keep fewer cells and speed up the sampler, but exclude +more of the light likelihood surface entirely. The default \code{0.99} is a +conservative state-space reduction.} -\item{thr_gs}{maximum movement speed in km/h used for hard support.} +\item{thr_gs}{maximum ground speed, in km/h, used as the hard radius of the +daily movement kernel (default \code{2000 / 24}, i.e. 2000 km per day). Speeds +below it are weighted by the \code{movement} kernel, so \code{thr_gs} is a +truncation of that kernel rather than a typical speed: set it above the +speeds the kernel supports. The default excludes less than 0.1\% of the +default kernel. Larger values increase computation; smaller values make +faster transitions impossible rather than unlikely.} \item{refresh}{progress update interval in iterations.} -\item{workers}{number of parallel workers for independent chains.} +\item{workers}{number of parallel workers for independent chains. +Use \code{1} for sequential execution.} -\item{seed}{optional random seed.} +\item{seed}{optional random seed for reproducible chains.} \item{quiet}{whether to suppress progress messages.} } @@ -71,16 +96,108 @@ A data.frame with columns \code{j}, \code{chain}, \code{stap_id}, \code{ind}, \c } \description{ Sample stationary trajectories from light likelihood maps and calibrated -movement components. +movement components. The posterior combines three components: the +stationary-period light likelihood, a daily movement model (see +\code{\link[=sampling_path_default_movement]{sampling_path_default_movement()}}), and a route-distance prior on the +interval between consecutive long periods (\code{stap0}). +} +\details{ +\strong{Route-distance prior.} For each interval between consecutive long +periods \eqn{\ell}, the route ratio \eqn{R_\ell} is the sum of the daily +great-circle displacements divided by the direct distance \eqn{D_\ell} +between the two long-period locations. The sampler adds +\code{component_weights["route"]} times the log density +\deqn{\log \mathrm{Gamma}(R_\ell - 1 \mid \alpha_R,\ \mathrm{route\_detour} \times \theta_\ell),}{ + log dgamma(R - 1, shape = alpha_R, scale = route_detour * theta),} +with shape \eqn{\alpha_R = 1.22} and +\deqn{\log \theta_\ell = -1.88 + 0.552\,[\log(1 + T_\ell) - 3.47] - 0.318\,[\log D_\ell - 8.15],}{ + log(theta) = -1.88 + 0.552 * (log(1 + T) - 3.47) - 0.318 * (log(D) - 8.15),} +where \eqn{T_\ell} is the interval duration in days and \eqn{D_\ell} is in +km. Both covariates are clamped to the range of the calibration data. The +term is omitted when \eqn{D_\ell < 300} km or when the interval contains at +most one daily displacement. + +Because \code{route_detour} multiplies the Gamma scale, the prior mean and every +prior quantile of the route excess \eqn{R_\ell - 1} scale linearly with it. +For example, an interval of about 31 days with endpoints about 3,500 km +apart has \eqn{\theta_\ell \approx 0.15}, i.e. a prior mean route ratio of +about 1.19; \code{route_detour = 2} raises it to about 1.37 and \code{0.5} lowers it +to about 1.09. \code{route_detour = 0} collapses the scale, so any detour is +strongly penalised. The prior remains soft: the realised detour of the +posterior also depends on the light and movement components. + +\strong{Long periods and modular cut.} Twilights accurately locate long periods +but leave daily-period locations ambiguous. With +\code{long_period_light_only = TRUE}, each unknown long-period location is +redrawn independently at every iteration from its retained light +likelihood (raised to \code{component_weights["light"]}), without using the +movement or route components. When a redraw violates hard movement support, +only the neighbouring daily locations needed to restore a feasible path are +resampled. Daily periods are then updated conditionally on these long-period +draws. This is a modular (cut) approximation to the joint posterior: +information flows from the long-period likelihoods to the daily path, but +the movement model cannot pull long-period locations away from their light +likelihood. Known-location periods are always fixed. + +\strong{Computational approximations.} Impossible geography, such as water after +masking, is removed from the state space. The arguments \code{thr_likelihood}, +\code{thr_gs} and \code{movement$approx_eps} control how aggressively the sampler +reduces the fixed light support and movement kernel for speed. Start with +the defaults for final inference, use a smaller \code{thr_likelihood} only for +explicitly approximate exploratory runs, and increase \code{thr_likelihood} or +\code{thr_gs} if no movement-feasible path exists. } \examples{ \dontrun{ -paths <- sampling_path(tag, iter = 1000, chains = 4, seed = 1) +example_dir <- system.file("extdata", package = "GeoPathSampleR") +tag <- GeoPressureR::tag_create( + "14OI", + crop_start = "2015-07-17", + crop_end = "2016-07-11", + directory = file.path(example_dir, "data/raw-tag/14OI"), + assert_pressure = FALSE, + quiet = TRUE +) +tag <- GeoPressureR::twilight_create(tag) +tag <- GeoPressureR::twilight_label_read( + tag, + file = file.path(example_dir, "data/twilight-label/14OI-labeled.csv") +) +tag <- GeoPressureR::tag_stap_daily( + tag, + stap_long = file.path(example_dir, "data/stap-label/14OI.csv"), + movement_period = "day", + quiet = TRUE +) +tag <- GeoPressureR::tag_set_map( + tag, + extent = c(-20, 40, -11, 60), + scale = 1 +) +tag <- GeoPressureR::geolight_map( + tag, + twl_calib_adjust = 1, + fitted_location_duration = 30, + twl_llp = \(n) 1.5 * log(n) / n, + quiet = TRUE +) +tag$map_light <- GeoPressureR::map_add_mask_water(tag$map_light) + +paths <- sampling_path( + tag, + movement = sampling_path_default_movement(), + iter = 500, + chains = 4, + warmup = 100, + seed = 1, + quiet = TRUE +) } } \seealso{ Other sampling_path: \code{\link[=plot_movement]{plot_movement()}}, +\code{\link[=sampling_path_default_movement]{sampling_path_default_movement()}}, \code{\link[=sampling_path_diagnostic]{sampling_path_diagnostic()}} } \concept{sampling_path} diff --git a/man/sampling_path_default_movement.Rd b/man/sampling_path_default_movement.Rd new file mode 100644 index 0000000..061253a --- /dev/null +++ b/man/sampling_path_default_movement.Rd @@ -0,0 +1,63 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/sampling_path.R +\name{sampling_path_default_movement} +\alias{sampling_path_default_movement} +\title{Default movement model for path sampling} +\usage{ +sampling_path_default_movement() +} +\value{ +A movement list with fields \code{method}, \code{shape}, \code{scale}, +\code{low_speed_fix}, \code{zero_speed_ratio}, \code{move_stay_parameters}, and +\code{move_stay}. +} +\description{ +Return the calibrated movement model used by default in \code{\link[=sampling_path]{sampling_path()}}. +Call \code{sampling_path_default_movement()} to inspect its values, plot it with +\code{\link[=plot_movement]{plot_movement()}}, or modify some fields before passing it to +\code{sampling_path(movement = )}. +} +\details{ +The movement model defines the daily transition between consecutive +stationary periods. It is a list with two parts. + +\strong{Ground-speed kernel.} The fields \code{method}, \code{shape}, \code{scale}, +\code{low_speed_fix} and \code{zero_speed_ratio} are passed to +\code{\link[GeoPressureR:speed2prob]{GeoPressureR::speed2prob()}} and follow the convention of +\code{\link[GeoPressureR:graph_set_movement]{GeoPressureR::graph_set_movement()}}. Speed is the great-circle displacement +between two consecutive periods divided by 24 h, in km/h. With +\code{method = "gamma"}, the speed density is a Gamma distribution with shape +\code{shape} and scale \code{scale} (km/h). The probability of moving to a given cell +is this density divided by the distance and multiplied by \code{cos(latitude)}, +which accounts for the number and area of cells at that distance. The hard +movement radius is set separately by \code{sampling_path(thr_gs = )}. + +\strong{Residence-dependent departure.} \code{move_stay} is a function of the +residence \code{t}, the number of consecutive transitions spent at the current +location (\code{0} for the first transition after arrival), returning the +probability of departing at the next transition. The default is the +exponential decline to a plateau +\code{p_inf + (p_0 - p_inf) * exp(-t / tau)}. \code{move_stay_parameters} stores +\code{p_0}, \code{p_inf} and \code{tau} so that \code{\link[=plot_movement]{plot_movement()}} can annotate them; it is +not used by the sampler, so update \code{move_stay} when changing them. + +The optional field \code{approx_eps} (default \code{1e-6}) drops movement-kernel +weights below this value from the precomputed neighbourhood. + +The default values were calibrated on multi-sensor reference trajectories +of land birds (Gamma shape \code{1.26} and scale \code{10.3} km/h; \code{p_0 = 0.607}, +\code{p_inf = 0.137}, \code{tau = 0.932}). They should not be assumed to transfer to +other taxa. +} +\examples{ +movement <- sampling_path_default_movement() +str(movement) +plot_movement(movement) +} +\seealso{ +Other sampling_path: +\code{\link[=plot_movement]{plot_movement()}}, +\code{\link[=sampling_path]{sampling_path()}}, +\code{\link[=sampling_path_diagnostic]{sampling_path_diagnostic()}} +} +\concept{sampling_path} diff --git a/man/sampling_path_diagnostic.Rd b/man/sampling_path_diagnostic.Rd index daf4908..2e84312 100644 --- a/man/sampling_path_diagnostic.Rd +++ b/man/sampling_path_diagnostic.Rd @@ -32,6 +32,7 @@ response. See the \href{https://geopressure.org/GeoPathSampleR/sampling-path-dia \seealso{ Other sampling_path: \code{\link[=plot_movement]{plot_movement()}}, -\code{\link[=sampling_path]{sampling_path()}} +\code{\link[=sampling_path]{sampling_path()}}, +\code{\link[=sampling_path_default_movement]{sampling_path_default_movement()}} } \concept{sampling_path} diff --git a/vignettes/figures/figure_1.png b/vignettes/figures/figure_1.png index 184922c..e32fb7d 100644 Binary files a/vignettes/figures/figure_1.png and b/vignettes/figures/figure_1.png differ diff --git a/vignettes/getting-started.qmd b/vignettes/getting-started.qmd index 3ff7dc9..0d368b9 100644 --- a/vignettes/getting-started.qmd +++ b/vignettes/getting-started.qmd @@ -11,9 +11,9 @@ execute: This vignette uses a light-only example tag bundled with GeoPathSampleR. It shows the complete workflow from raw geolocator data to sampled paths. The example is intentionally short and uses coarse settings so that it can be run while learning the workflow; it is not intended as a scientific analysis. -Figure @fig-workflow provides an overview of the steps that follow. Steps 0--3 cover the preparation of the likelihood map of light implemented in GeoPressureR and which is more extensively documented in the [GeoPressureManual](https://geopressure.org/GeoPressureManual/). We focus the explanation around key difference and provides links to the relevant section at the start of each subsection. Step 4 covers the sampling of the movement-model implemented in GeoPathSampleR. +@fig-workflow shows the two parts of the workflow. Part I uses [GeoPressureR](https://geopressure.org/GeoPressureR/) to turn curated twilights into stationary-period likelihood maps, with [GeoLightViz](https://geopressure.org/GeoPressureManual/geolightviz.html) supporting interactive review. Part II uses GeoPathSampleR to combine those maps with movement priors and sample posterior trajectories. The sections below follow this order; the [Light map chapter of the GeoPressureManual](https://geopressure.org/GeoPressureManual/light-map.html) documents Part I in more detail. -![GeoTwilight workflow. Steps 0--3 prepare, label, calibrate, and constrain the light likelihoods; Step 4 samples posterior trajectories with GeoPathSampleR. The subsection headings below follow this workflow.](figures/figure_1.png){#fig-workflow width="95%"} +![The two-part workflow. Part I constructs inspectable stationary-period likelihood maps from curated twilight observations using GeoPressureR and GeoLightViz. Part II combines these maps with movement priors in GeoPathSampleR to sample posterior trajectories.](figures/figure_1.png){#fig-workflow width="95%"} ## Getting ready @@ -93,6 +93,13 @@ tag <- tag_stap_daily( ) ``` +::: {.callout-tip} +## Revising long periods + +After `tag_stap_daily()`, `tag$stap` contains the daily periods as well as the long periods. To revise the long periods in GeoLightViz, start again with `tag_create()`, `twilight_create()`, and `twilight_label_read()` before opening `geolightviz(tag)`. Save the revised `data/stap-label/{id}.csv`, then rerun `tag_stap_daily()` and the following steps. + +::: + ### Define map [Defining a map](https://geopressure.org/GeoPressureManual/pressure-map.html#define-geographical-and-temporal-parameters-of-the-map) sets the spatial extent and grid resolution used by the likelihood maps. This step comes after stationary periods have been defined because the map uses them as its temporal dimension. The example does not contain known coordinates, so calibration is fitted automatically from the long stationary periods. @@ -136,7 +143,7 @@ tag <- geolight_map( ``` -::: {.callout-info} +::: {.callout-note} The arguments used here differ from the GeoPressureR defaults and follow the calibration described in the GeoTwilight analysis (see the relevant appendix): @@ -146,6 +153,14 @@ The arguments used here differ from the GeoPressureR defaults and follow the cal ::: +### Check the twilight calibration + +The calibration describes the distribution of twilight errors, expressed as solar zenith angles, used to build every likelihood map. Inspect it before continuing: the histogram shows the zenith angles of the calibration twilights, coloured by calibration period (here all fitted from long periods), and the black curve shows the fitted distribution used by the likelihood. A calibration period whose histogram departs strongly from the others can indicate mislabelled twilights or a period that is not stationary. + +```{r} +plot_twl_calib(tag) +``` + ### Add the water mask For landbird, a water mask removes cells from the set of locations that can be sampled. @@ -159,13 +174,21 @@ plot(tag, type = "map_light") ## Sampling trajectories The movement model implemented in `sampling_path()` combines three sources of information to define the trajectory: + - The light likelihood describes where each stationary period could be located. - The movement model describes plausible speeds and the probability of departing after a given residence time. - The route-length prior provides information about the expected detour of long intervals between stationary periods. ### Defining parameters -The default movement model is a Gamma speed distribution combined with a residence-dependent probability of movement. +#### Movement model + +The movement model defines the daily transition between consecutive stationary periods. It has two parts: + +- A ground-speed kernel (`method`, `shape`, `scale`, `low_speed_fix`, `zero_speed_ratio`), evaluated with `GeoPressureR::speed2prob()`. Speed is the displacement between two consecutive periods divided by 24 h. The default is a Gamma distribution with shape `1.26` and scale `10.3` km/h. +- A residence-dependent departure probability `move_stay(t)`, where `t` is the number of consecutive days already spent at the current location. The default declines exponentially from `p_0 = 0.607` just after arrival to a plateau `p_inf = 0.137` with time constant `tau = 0.932` days. `move_stay_parameters` only stores these values for `plot_movement()`; the sampler uses `move_stay`. + +The default `movement` argument of `sampling_path()` is `sampling_path_default_movement()`, calibrated on multi-sensor tracks of land birds. We write it out below to show its structure; `movement <- sampling_path_default_movement()` gives the same list. See `?sampling_path_default_movement` for details. ```{r} movement <- list( @@ -184,37 +207,40 @@ movement <- list( plot_movement(movement) ``` -The principal model controls are: +#### Model controls -- `component_weights = c(light = 1, movement = 1, route = 1)` balances the light, movement, and route components. Only relative values matter. Setting a weight to zero removes that component's soft contribution and should be used only for an explicit sensitivity analysis. -- `route_detour = 1` retains the empirical excess route-distance prediction. Values below one favour more direct routes; values above one favour more detoured routes. -- `long_period_light_only = TRUE` samples unknown long-period locations from their light likelihood without allowing the daily movement model to update them. This is useful when movement should not pull a long-period location toward another known or well-supported location. +The other arguments of `sampling_path()` that change the scientific model are: -The sampler uses Gibbs updates: it repeatedly updates one stationary-period location, or a block of locations, conditional on the current values of the other periods. The resulting draws are correlated, so the run must be long enough to explore the posterior and must be checked with multiple independent chains. - -The following arguments control computation and diagnostics rather than the scientific model: +- `component_weights = c(light = 1, movement = 1, route = 1)` uses each component without rescaling. Each weight scales that component's log probability, so changing even all three weights by the same factor changes the sampled distribution. Setting a weight to zero removes that component's soft contribution while retaining hard support constraints; use this only for an explicit sensitivity analysis. +- `route_detour = 1` keeps the calibrated route prior. The value multiplies the expected excess route distance (route length over direct distance, minus one) between consecutive long periods. For instance, for a prior mean route is 19% longer than direct; `2` doubles this excess, `0.5` halves it and `0` favours direct routes. +- `long_period_light_only = TRUE` applies a modular cut: unknown long-period locations are drawn from their light likelihood alone, and movement only informs the daily periods between them. `FALSE` lets movement also update long periods. -| Argument | Role and practical choice | -|---|---| -| `iter` | Number of saved iterations per chain after warmup. Increase it when effective sample sizes are too small. | -| `chains` | Number of independent chains. Use at least four for final inference. | -| `warmup` | Initial iterations discarded while the chain reaches its typical posterior region. | -| `thin` | Saving interval. Keep `1` unless storage is limiting; thinning does not improve mixing. | -| `block_interval` | Frequency of block updates. Keep the default unless explicitly comparing sampler efficiency. | -| `thr_likelihood` | Retained light-likelihood mass used to reduce the state space. The conservative default `0.99` is appropriate for final inference. | -| `thr_gs` | Maximum ground speed in km/h defining hard movement support. Use a biologically credible value. | -| `workers` | Number of parallel workers for independent chains. Use `1` while debugging and increase only when appropriate for the system. | -| `seed` | Random-number seed for reproducible chains. | +#### Computational controls +The following arguments control computation and diagnostics rather than the scientific model: +| Argument | Default | Role and practical choice | +|---------|-----|----------------------| +| `iter` | none | Total iterations per chain, including warmup. Increase it when effective sample sizes are too small. | +| `warmup` | `floor(iter / 4)` | Initial iterations discarded while the chain reaches its typical posterior region. | +| `thin` | `1` | Saving interval. The number saved is `floor((iter - warmup) / thin)`. Keep `1` unless storage is limiting; thinning does not improve mixing. | +| `chains` | `1` | Number of independent chains. Use at least four for final inference. | +| `block_interval` | `2` | Frequency of block updates, which move a whole stay at once. Keep the optimized default unless explicitly comparing sampler efficiency. | +| `thr_likelihood` | `0.99` | Retained light-likelihood mass used to reduce the state space. The conservative default is appropriate for final inference. | +| `thr_gs` | `2000 / 24` | Hard truncation of the speed kernel, in km/h (2000 km per day). Speeds are already weighted by the movement model, so set it above the speeds the kernel supports rather than at typical speeds; the default excludes less than 0.1% of the default kernel. Larger values cost computation; smaller values forbid transitions outright. | +| `workers` | `1` | Number of parallel workers for independent chains. Use `1` while debugging and increase only when appropriate for the system. | +| `seed` | `NULL` | Random-number seed for reproducible chains. | ### First run +The sampler uses Gibbs updates: it repeatedly updates one stationary-period location, or a block of locations, conditional on the current values of the other periods. The resulting draws are correlated, so the run must be long enough to explore the posterior and must be checked with multiple independent chains. + Begin with a short run to verify that the likelihood, movement support, and route prior are compatible. This run is for checking the workflow, not for inference. ```{r} paths <- sampling_path( tag, + movement = movement, iter = 500, chains = 4, warmup = 100, @@ -222,25 +248,26 @@ paths <- sampling_path( quiet = TRUE ) ``` - - -### Diagnose the sampler - -Always inspect convergence and computational-support diagnostics before interpreting a trajectory. The interactive report provides an overview of sampled paths, traces, R-hat, effective sample sizes, chain separation, and warnings about likelihood or movement boundaries. +Inspect convergence and computational-support diagnostics before interpreting a trajectory. The interactive report provides an overview of sampled paths, traces, R-hat, effective sample sizes, chain separation, and warnings about likelihood or movement boundaries. For interactive review, generate the report with: -An example of the resulting report is available on the [diagnostic example page](https://geopressure.org/GeoPathSampleR/sampling-path-diagnostic.html). - ```{r} #| eval: false diagnostic <- sampling_path_diagnostic(paths, tag, report = TRUE) ``` +An example of the resulting report is available on the [diagnostic example page](https://geopressure.org/GeoPathSampleR/sampling-path-diagnostic.html). + +### Final run + +Once the first run shows no problem with the likelihood, movement support or route prior, run longer chains for inference. Increase `iter` and `warmup` until R-hat is close to one and effective sample sizes are large enough for the quantities of interest, and check the diagnostics again before interpreting the paths. The settings below are illustrative only; a real analysis usually needs more iterations. + ```{r} #| eval: false paths <- sampling_path( tag, + movement = movement, iter = 1000, chains = 4, warmup = 500, @@ -252,43 +279,52 @@ sampling_path_diagnostic(paths, tag, report = FALSE) ### Summarise posterior trajectories -`path_collapse()` retains a stay summary for every sampled path. `path_summary()` returns either one posterior summary per stationary period or one posterior consensus trajectory by first deriving the modal stay/move structure. The consensus is useful for a compact summary, but it does not replace uncertainty assessment. +`paths` contains one row per posterior draw and stationary period: each draw is a complete trajectory, and daily periods spent at the same location appear as repeated rows. Two functions turn it into something easier to interpret, and they answer different questions: + +- `path_collapse()` works **within each draw**. It merges consecutive periods assigned to the same grid cell into a single stay, so each draw becomes a sequence of stays and moves. All draws are kept, and so is the posterior uncertainty. +- `path_summary()` works **across draws**. It reduces the ensemble to a single path, either one position per stationary period (`by = "stap"`) or one position per consensus stay (`by = "consensus_stay"`). + +| Summary | Call | One row per | Keeps uncertainty | +|------------|-------------------------------|------------------|--------| +| Collapsed draws | `path_collapse(paths, tag$stap)` | stay of each draw | yes | +| Mean path | `path_summary(paths, tag$stap, by = "stap")` | stationary period | no | +| Consensus path | `path_summary(paths, tag$stap, by = "consensus_stay")` | consensus stay | no | + +#### Individual draws + +We draw ten posterior trajectories at random and collapse each into its stays. + +```{r} +set.seed(1) +path_ids <- unique(paths[c("chain", "j")]) +path_ids <- path_ids[sample.int(nrow(path_ids), 10), ] +draws <- merge(paths, path_ids) |> + path_collapse(tag$stap) +``` + +#### Summary paths + +With `by = "stap"`, `path_summary()` summarises the positions of all draws separately for each stationary period, here with the mean. With `by = "consensus_stay"`, it first estimates the probability of a move at each boundary between consecutive periods as the proportion of draws whose location changes there. Boundaries with a probability above `move_threshold = 0.5` define the consensus stays, and positions are then summarised over all draws and periods of each stay. The columns `move_probability_start` and `move_probability_end` report how clear each boundary is. ```{r} -collapsed <- path_collapse(paths, tag$stap) +mean_path <- path_summary(paths, tag$stap, by = "stap", position = "mean") consensus <- path_summary(paths, tag$stap, by = "consensus_stay") +head(consensus[c("stap_id_start", "stap_id_end", "lat", "lon", "move_probability_end")]) ``` -To show both the central trajectory and posterior variability, we draw ten paths at random, calculate the posterior mean position for each stationary period with `path_summary()`, and overlay the consensus path. The random draws are deliberately shown as thin grey lines; the mean and consensus should remain visually prominent. +#### Compare the summaries + +@fig-path-summaries overlays the three summaries. ```{r} #| code-fold: true #| code-summary: "Show plotting code" #| label: fig-path-summaries -#| fig-cap: "Posterior trajectory summaries. Thin grey lines show ten randomly selected collapsed posterior draws. The green dashed line connects posterior mean positions for each stationary period, while the orange line shows the consensus stay/move path. The three summaries answer different questions and should not be interpreted as interchangeable estimates." +#| fig-cap: "Posterior trajectory summaries. Thin grey lines show ten randomly selected collapsed posterior draws. The green dashed line connects posterior mean positions for each stationary period, while the orange line shows the consensus stay/move path." -set.seed(1) -path_ids <- unique(paths[c("chain", "j")]) -path_ids <- path_ids[ - sample.int(nrow(path_ids), size = min(10, nrow(path_ids))), - , - drop = FALSE -] -draws <- paths[ - interaction(paths$chain, paths$j) %in% - interaction(path_ids$chain, path_ids$j), - , - drop = FALSE -] -draws <- path_collapse(draws, tag$stap) draws$j <- interaction(draws$chain, draws$j, drop = TRUE) - -mean_path <- path_summary( - paths, - tag$stap, - by = "stap", - position = "mean" -) +mean_path$j <- 1 +consensus$j <- 1 map_extent <- tag$param$tag_set_map$extent map <- leaflet::leaflet( @@ -310,7 +346,6 @@ map <- GeoPressureR::plot_path( polyline = list(color = "#8C8C8C", weight = 2, opacity = 0.45), circle = list(radius = 1, stroke = FALSE, fillOpacity = 0) ) -mean_path$j <- 1 map <- GeoPressureR::plot_path( mean_path, map = map, @@ -323,7 +358,6 @@ map <- GeoPressureR::plot_path( ), circle = list(radius = 1, stroke = FALSE, fillOpacity = 0) ) -consensus$j <- 1 GeoPressureR::plot_path( consensus, map = map, @@ -337,4 +371,8 @@ GeoPressureR::plot_path( ) ``` -The posterior mean is not necessarily a valid sampled trajectory: averaging locations can place a point between two distinct posterior modes. The consensus path preserves a representative stay/move structure, while the individual draws make multimodality and route uncertainty visible. +The three summaries should not be interpreted as interchangeable estimates: + +- The **individual draws** are valid trajectories and make multimodality and route uncertainty visible. Use them, or all collapsed draws, for any quantity derived from the trajectory, such as stopover durations or route length, and then summarise that quantity across draws. +- The **mean path** follows the stationary periods, but is not necessarily a valid trajectory: averaging locations can place a point between two distinct posterior modes, and it does not show when the bird stayed or moved. +- The **consensus path** preserves a representative stay/move structure, which gives a compact summary of the timing of movements. It does not replace an assessment of uncertainty.