Skip to content
BranchesTags

Repository files navigation

msomgom: Multi-Season Occupancy Model of the Gulf of Maine

A dynamic (multi-season, colonization/persistence) occupancy model — following the MacKenzie et al. (2003) formulation, itself an extension of the single-season detection model of MacKenzie et al. (2002) — for cetacean vessel-survey sightings, fit in JAGS (Plummer 2003). Originally built around Bay of Fundy North Atlantic right whale (RIWH) survey data from the NARWC Sightings Database (Kenney 2021), now generalized so file paths, the study-area polygon, date range, species, and MCMC settings are all driven by a YAML config rather than hardcoded.

See References at the bottom for full citations.

The model is currently fit one species per run (selected via config); the data-prep/gridding stages already loop over a configurable list of species and build one detection-history array per species, so extending this to a true hierarchical multi-species model later doesn't require redoing that part. See docs/refactor_plan.md for the full history of how this pipeline got to its current shape and why.

msomgom is a regular R package (devtools::check() passes with 0 errors/0 warnings/0 notes). The fastest way to see it work end-to-end is the vignette:

vignette("getting-started", package = "msomgom")

(Only found if msomgom was installed with build_vignettes = TRUE - see Install the package below; it's not the default.)

Contents: Quick start · Setup · Usage · Data · Repository layout · References

Quick start

No real survey data needed — this generates a config, a synthetic survey CSV, and fits the model against it, all in-memory in a temp directory:

# install.packages("devtools")
# devtools::install_github("chross22/msomgom", build_vignettes = TRUE, dependencies = TRUE)
library(msomgom)

dir.create("configs")
generate_config("mock_test", data_file = "data/mock_survey_data.csv",
                 n_chains = 2, n_adapt = 100, n_burn = 100, n_iter = 200)
generate_mock_data("configs/mock_test.yaml")

result <- run_occupancy_model("configs/mock_test.yaml")
result$evaluation$parameters # posterior summary: mean, sd, Rhat, effective size, ...

This needs rjags/dclone installed (plus a working JAGS binary) to actually fit the model — see Setup below for that. If you just want to see the whole thing run without installing JAGS yourself, vignette("getting-started", package = "msomgom") walks through this exact example already rendered.

Setup

You need R, JAGS, and a handful of R packages with native dependencies (spatial libraries for sf, JAGS bindings for rjags).

1. System dependencies (macOS / Homebrew)

brew install gdal geos proj udunits cmake

JAGS (Plummer 2003) needs to be 4.x, not 5.x. CRAN's rjags (Plummer 2025, as of this writing, 4-17) hard-rejects any JAGS version other than 4.x, but Homebrew's jags formula has moved on to 5.0.0. Install JAGS 4.3.2 from the last formula revision that built it, via a local tap:

brew tap-new local/jags4
curl -s "https://raw.githubusercontent.com/Homebrew/homebrew-core/cafb141891/Formula/j/jags.rb" \
  -o "$(brew --repository)/Library/Taps/local/homebrew-jags4/Formula/jags.rb"
sed -i '' '/no_autobump!/d' "$(brew --repository)/Library/Taps/local/homebrew-jags4/Formula/jags.rb"
brew install local/jags4/jags

If you hit fatal error: 'cmath' file not found while installing R packages below, your Xcode Command Line Tools are broken/incomplete. Fix with:

sudo rm -rf /Library/Developer/CommandLineTools
xcode-select --install

2. Install the package

# install.packages("devtools")
devtools::install_github("chross22/msomgom", build_vignettes = TRUE, dependencies = TRUE)

build_vignettes = TRUE is not the default - devtools/remotes skip building vignettes unless asked (remotes::install_github()'s build_opts actually includes --no-build-vignettes by default), so vignette("getting-started", package = "msomgom") will say the vignette doesn't exist if you install without it. dependencies = TRUE pulls in knitr/rmarkdown, needed to build it, since they're Suggests-only.

If you already installed without either flag, vignette() still won't find it after adding them later - reinstall with both flags rather than trying to patch an existing install. In the meantime, the same walkthrough is on GitHub: vignettes/getting-started.Rmd.

This pulls in the hard dependencies (dplyr, lubridate, parallel, readr, sf, sfheaders, stringr, tibble, yaml) automatically. A few packages are optional (in Suggests) and only needed for specific features - installed separately, and only if you use them:

  • rjags (Plummer 2025) and dclone (Sólymos 2010) - actually fitting the model (fit_occupancy_model()/run_occupancy_model()); also needs a working JAGS install (see step 1)
  • coda (Plummer et al. 2006) - convergence diagnostics (evaluate_occupancy_model())
  • terra (Hijmans et al. 2026) - reading local NetCDF covariate files (load_covariate_netcdf())
  • mapview, tmap, webshot - only if you set output.make_figs: true in a config
  • googledrive - only if you set paths.google_drive_filename
  • Microsoft365R - only if you set paths.onedrive_filename (see Fetching survey data below)
  • datamatch, derivoce - fetching/deriving environmental covariates (see Environmental covariates)
  • fancyfx - plotting a fitted covariate's effect (see Covariate effect plots)
  • knitr, rmarkdown - building the vignette
install.packages(c("rjags", "dclone", "coda", "terra"))

Spatial gridding uses sf (Pebesma 2018; Pebesma & Bivand 2023); config parsing uses yaml (Stephens & Simonov 2025); data cleaning uses dplyr (Wickham et al. 2026a), readr (Wickham et al. 2026b), stringr (Wickham 2025), lubridate (Grolemund & Wickham 2011), and tibble (Müller & Wickham 2026). This pipeline is written in R (R Core Team 2026); see References for full citations.

Usage

Configure a run

Generate a config programmatically rather than hand-editing YAML:

library(msomgom)
generate_config(
  "my_run",
  data_file = "data/my_survey_data.csv",
  beg_year = 2015, end_year = 2020,
  species_codes = c("RIWH", "HUWH"), active_species = "HUWH"
)
# -> writes configs/my_run.yaml (relative to the current working directory)

Two example configs ship with the package under system.file("extdata/configs", package = "msomgom"): bof_riwh.yaml reproduces the pipeline's original hardcoded Bay of Fundy / RIWH settings, and mock_test.yaml is a small/fast config for the mock-data smoke test below. See ?generate_config for what every field means (paths, survey platform/FILEID/on-effort filters, date range, season boundaries, study-area polygon, species, JAGS/MCMC settings, evaluation, cleanup).

The survey filters (platform_code, fileid_prefixes, on_effort_legtypes) are all config-driven rather than hardcoded to any one platform, so the pipeline isn't limited to vessel surveys - an aerial survey in the same NARWC handbook format works too, with its own platform code, FILEID prefix, and on-effort LEGTYPE codes (which differ from a vessel's; check the handbook for the codes that apply to your platform).

Run the pipeline

library(msomgom)
result <- run_occupancy_model("configs/my_run.yaml")
result$fit          # the fitted mcmc.list
result$evaluation    # convergence diagnostics + posterior summary, see below

or from the command line, via the CLI wrapper installed with the package:

# run once to find where it lives on your machine:
system.file("scripts/run_pipeline.R", package = "msomgom")
#> "/usr/local/lib/R/site-library/msomgom/scripts/run_pipeline.R"
Rscript /usr/local/lib/R/site-library/msomgom/scripts/run_pipeline.R configs/my_run.yaml

run_occupancy_model() calls, in order: prep_survey_data() (load + clean the survey CSV) -> build_detection_arrays() (build the spatial grid and detection-history arrays) -> fit_occupancy_model() (fit the model) -> evaluate_occupancy_model() (diagnostics, unless disabled in the config) -> cleanup_outputs() (only if cleanup.run_after_model: true).

Try it without real data

The real NARWC survey CSV isn't included in this repo (see Data below). To exercise the whole pipeline with synthetic data:

library(msomgom)
generate_config("mock_test", data_file = "data/mock_survey_data.csv",
                 output_dir = "output/mock_test",
                 n_chains = 2, n_adapt = 100, n_burn = 100, n_iter = 200)
generate_mock_data("configs/mock_test.yaml")
run_occupancy_model("configs/mock_test.yaml")

Mock data is grounded in the NARWC handbook's actual field codes (PLATFORM, FILEID, LEGTYPE, LEGSTAGE, VISIBLTY, BEAUFORT, IDREL, SPECCODE), including some decoy records the pipeline's filters should drop, so it's a real exercise of the filtering logic, not just a happy-path stub. See vignette("getting-started", package = "msomgom") for this same walkthrough, fully worked through and rendered.

Evaluate a saved run

evaluate_occupancy_model() also works standalone against a previously saved result, without re-fitting:

library(msomgom)
config <- load_config("configs/my_run.yaml")
evaluate_occupancy_model("output/RIWH.colext.RData", config = config)

It reports a posterior parameter summary table, Gelman-Rubin Rhat and effective sample size per parameter (with warnings for likely non-convergence or a noisy posterior), trace/density plots saved to <output_dir>/mcmc_diagnostics.pdf, and, when the model was run with jags.params: Z, a naive-vs-modeled occupancy comparison per year.

Debugging and assessing data/outputs

Before fitting, diagnose_pipeline() runs the data-prep and grid-building stages (never JAGS itself) and reports the most common reasons a run fails, hangs, or "succeeds" with meaningless output - a filter that leaves nothing, a study-area polygon that misses the survey tracks, a species with zero detections, a record that fell outside every configured season, or a covariate matrix with the wrong shape:

library(msomgom)
diagnose_pipeline("configs/my_run.yaml")
# with covariates, checked against config$covariates$psi/phi/gamma too:
diagnose_pipeline("configs/my_run.yaml", occ_covariates = list(sst = sst_avg$sst))

It prints a plain PASS/WARN/FAIL report and returns list(config, prep, arrays) invisibly (whichever were reached), so you can pick up investigating right where it stopped - e.g. plot_survey_coverage(result$arrays).

When diagnose_pipeline() isn't enough - a crash inside a pipeline stage itself, not just a suspicious-looking result - R's own debugging tools work on any msomgom function like they would on your own code:

  • traceback() - run immediately after an error, before anything else, to print the exact call stack that led to it.
  • prep_survey_data(config, verbose = TRUE) - reports how many records survive each filter step (platform, FILEID prefix, date range) as it runs, so you can see exactly where a filter zeroes things out.
  • debugonce(prep_survey_data) (or any other exported function) - runs that one call in R's interactive line-by-line debugger: n steps to the next line, c runs to completion, and typing a variable's name at the browser prompt shows its current value. Only fires once, so it doesn't need to be turned back off.
  • options(error = recover) - turned on for the rest of the session, this drops you into an interactive prompt at every uncaught error, letting you pick any frame in the call stack and inspect it - useful when you don't know in advance which function is going to fail. Turn it back off with options(error = NULL).
# see the exact call stack right after a crash
traceback()

# step through prep_survey_data() line by line
debugonce(prep_survey_data)
prep_survey_data(config, verbose = TRUE)

# drop into the debugger at any uncaught error, for the rest of the session
options(error = recover)
run_occupancy_model("configs/my_run.yaml")
options(error = NULL) # turn it back off when done

Three plotting functions, for when a fit looks wrong and you need to see why rather than just the summary numbers. Each returns the sf grid it plotted (with the plotted column added) invisibly, so you can inspect the underlying values or hand them to mapview() for an interactive version.

library(msomgom)
config <- load_config("configs/my_run.yaml")
prep <- prep_survey_data(config)
arrays <- build_detection_arrays(prep$tmpdat, prep$season_info, config)

# Where was the study area actually surveyed? Blank cells you expected
# coverage in usually mean a study-area-polygon or grid problem.
plot_survey_coverage(arrays)
plot_survey_coverage(arrays, season = 3)              # one season only

# Where are the sightings landing, relative to coverage above?
plot_sightings(arrays, "RIWH")
plot_sightings(arrays, "RIWH", season = 3)

# After fitting with jags_params = "Z": posterior occupancy probability
# per cell, spatially - the payoff visualization for an occupancy model.
fit <- fit_occupancy_model(arrays, config)             # config$jags$params: "Z"
plot_occupancy_map(fit, arrays)
plot_occupancy_map(fit, arrays, year = 2)

# A covariate that's NA everywhere, constant, or spatially implausible is
# much faster to catch here than after a fit quietly does nothing with it.
sst_avg <- average_covariates(env_dat, arrays$area_grid_sf, windows)
plot_covariate_map(sst_avg$sst, arrays, var_name = "sst")

plot_survey_coverage()/plot_sightings() are built directly from the same arrays the model fits on (arrays$reps and arrays$species_arrays), so what you see is exactly what went into the model, not a rederivation from the raw CSV. plot_occupancy_map() reports the same posterior-mean-Z numbers as compare_naive_vs_modeled_occupancy()'s table, just spatially instead of per-year. plot_covariate_map() takes any one covariate matrix from average_covariates() (see Environmental covariates below) and shows one window of it at a time (window, by position or by windows$label; defaults to the last one).

Environmental covariates

whale.mod (in jags.R) puts environmental covariates on detection probability (day-of-year, sea state, effort), and optionally on initial occupancy (psi), persistence (phi), and colonization (gamma) too - each process independently gets zero or more named covariates, chosen per config. A process with none configured stays intercept-only, exactly as if this didn't exist; a config with no covariates: section at all behaves identically to before this feature was added.

The two steps: prep the covariate data (average_covariates.R), then pass it into a run (covariates.psi/phi/gamma in the config + the occ_covariates argument).

1. Prep covariate data. average_covariates() works with any sf point object tagged with YEAR/MONTH/DAY, one column per variable - that shape is shared across a small family of sibling packages, so nothing needs reshaping between them:

  • datamatch (Ross, n.d.) fetches that shape live from the E.U. Copernicus Marine Service (Copernicus Marine Service, n.d.) via datamatch::accessEnvDat()
  • derivoce computes derived covariates (spatial/temporal gradients, distance to shore/front/isobath, lags, integrals, eddy kinetic energy, ...) from that same shape, returning it enriched with new columns - so it composes with either source below
  • load_covariate_netcdf() (in this package) builds the same shape directly from a folder of local daily NetCDF files, for when you already have files on disk instead of fetching live (e.g. from another source, or from datamatch run separately and saved)

Neither datamatch nor derivoce is a hard dependency - both are optional (Suggests), installed the same way as this package itself:

devtools::install_github("chross22/datamatch")
devtools::install_github("chross22/derivoce")

Whichever source you use, average_covariates() spatially averages onto this pipeline's hex grid and temporally averages over whatever windows you give it:

library(msomgom)
config <- load_config("configs/my_run.yaml")

# area_grid_sf must come from build_detection_arrays(), so the covariate grid
# matches the occupancy model's grid exactly, cell-for-cell
prep <- prep_survey_data(config)
arrays <- build_detection_arrays(prep$tmpdat, prep$season_info, config)

# tied to the occupancy model's own season/year structure
windows <- season_windows_from_config(config)
# or arbitrary fixed-interval windows, independent of any occupancy config
# windows <- regular_windows("2018-01-01", "2020-12-31", by = "1 month")

# from local files:
env_dat <- load_covariate_netcdf("data/covariates/sst", var_names = "sst")

# or live from Copernicus, optionally enriched with a derived covariate -
# both return the same env_dat shape average_covariates() expects, so this
# is a drop-in alternative to the load_covariate_netcdf() line above:
# env_dat <- datamatch::accessEnvDat(vars = "SST", years = 2018:2020, months = 1:12,
#                                     bounding_box = list(xmin = -70, xmax = -66, ymin = 41, ymax = 44))
# env_dat <- derivoce::horizontal_gradient(env_dat, "SST")  # adds an SST_grad column

sst_avg <- average_covariates(env_dat, arrays$area_grid_sf, windows)
# sst_avg$sst is a [num_cells x num_ssn] matrix, same num_ssn indexing as
# effort3d/jday3d/bft3d

2. Configure and run. List which covariate names (matching names in occ_covariates, e.g. sst_avg above) apply to which process, then pass the covariate list into run_occupancy_model():

library(msomgom)
generate_config("my_run", covariates_psi = c("sst"), covariates_phi = c("sst"))
# -> configs/my_run.yaml has covariates: { psi: [sst], phi: [sst], gamma: [] }

result <- run_occupancy_model("configs/my_run.yaml", occ_covariates = list(sst = sst_avg$sst))
result$evaluation$parameters # now includes mu.b.cov (psi) and mu.e.cov (phi)

Each process can take a different number of covariates (or none), including a different set of covariate names entirely - covariates.gamma above is empty, so colonization stays intercept-only while occupancy and persistence both get sst. Internally, since JAGS can't compile a covariate loop of length zero, the model's BUGS code is generated per run (build_whale_model_code() in jags.R) rather than being one static model, so each process's covariate block is included only when that process actually has covariates.

Extending this further - to a true hierarchical multi-species model that fits all configured species jointly with shared priors - is a deliberately separate, larger follow-up.

Covariate effect plots (via fancyfx)

A fit with covariates configured (above) is tagged msomgom_fit and carries the metadata needed to plot a covariate's fitted effect with fancyfx - the effect curve stacked above a rug of the raw covariate data, the same way fancyfx plots an mgcv::gam's partial effects or any other model's predictions. fancyfx is optional (Suggests); msomgom depends on it, not the other way around - effect_estimates.msomgom_fit() is what makes fancyfx work here, not anything fancyfx itself knows about this package:

# install.packages("devtools"); devtools::install_github("chross22/fancyfx")
library(msomgom)
library(fancyfx)

fit <- fit_occupancy_model(arrays, config, occ_covariates = list(sst = sst_avg$sst))

# the tidy estimate frame on its own:
effect_estimates(fit, "sst")

# or the full plot: effect curve + a rug of the raw covariate data above it
rug_dat <- data.frame(sst = as.vector(sst_avg$sst))
plotEffects(fit, rug_dat, "sst", xlab = "SST")

scale = "link" (the default) gives the covariate's own contribution to the linear predictor, centered at zero at the covariate's mean - the occupancy equivalent of a GAM's partial effect. scale = "response" gives the full predicted occupancy probability as the covariate varies, with any other covariates on the same process held at their mean. Since every msomgom fit is summarized from posterior draws, interval defaults to the credible interval ("ci") rather than a +/- 1 SE band, the same way fancyfx treats a brms/rstanarm fit. Requires a jags_params = "colext" fit (the default) - a "Z"-mode fit only tracks occupancy states, not the coefficients this reads.

Data

This pipeline expects a CSV in the NARWC Sightings Database format (Kenney 2021). data/ and output/ are gitignored — real survey data has its own data-sharing terms and shouldn't be committed here, and model outputs are regenerable from a run.

A real export's column names rarely match the NARWC schema exactly (Event instead of EVENTNO, Sp. Code instead of SPECCODE, ...). prep_survey_data() runs standardize_survey_columns() on the raw CSV first, which renames columns to their canonical name on a case-insensitive match or a small list of common aliases, so this usually doesn't need manual renaming first. ALT (altitude) additionally defaults to 750 if it's missing entirely, since this pipeline's own filtering doesn't use it. A column it genuinely can't match errors clearly, naming what's missing and what was actually found - see ?standardize_survey_columns for the alias list, or call it directly to check your own file before running the full pipeline: standardize_survey_columns(readr::read_csv("data/survey_data.csv")).

Fetching survey data (Google Drive / OneDrive)

If paths.data_file doesn't exist locally, prep_survey_data() can fetch it for you first - from Google Drive (paths.google_drive_filename, via the optional googledrive package) or OneDrive (paths.onedrive_filename, via the optional Microsoft365R package). Set at most one; if both are set, Google Drive takes priority.

generate_config("my_run", data_file = "data/survey_data.csv",
                 onedrive_filename = "Shared/survey_data.csv", # path *within* OneDrive
                 onedrive_type = "business")                   # "personal" (default) or "business"

onedrive_filename is the file's path within OneDrive, not a local path - Microsoft365R::get_personal_onedrive()/get_business_onedrive() resolve it via Microsoft Graph, prompting an interactive login on first use. Use onedrive_type = "business" for an organization/shared OneDrive.

If OneDrive is already synced to this machine (the usual case with the OneDrive desktop app), none of this is needed - just point data_file directly at wherever it's synced locally (e.g. "~/OneDrive - Org/.../survey_data.csv"), and skip onedrive_filename entirely.

Repository layout

msomgom is a standard R package; the layout follows the usual conventions (R CMD build/devtools::check() both pass clean):

DESCRIPTION, NAMESPACE            # package metadata + generated exports (roxygen2)
R/
  load_config.R                    # load_config(path) -> validated config list
  generate_config.R                # generate_config(name, ...) -> writes <name>.yaml
  generate_mock_data.R             # synthetic survey CSV for testing, no real data needed
  average_covariates.R             # spatial/temporal averaging of covariates onto the grid
  padstr0.R                        # zero-pads GMT time-of-day values
  makeSeasons.R                    # builds the season lookup table from config date ranges
  data_prep.R                      # prep_survey_data(config) -> cleaned survey data
  jagsPrep.R                       # build_detection_arrays(...) -> spatial grid + detection arrays
  jags.R                           # fit_occupancy_model(...) -> fits the JAGS model
  effect_estimates.R               # effect_estimates.msomgom_fit() -> fancyfx integration
  evaluate_model.R                 # evaluate_occupancy_model(...) -> convergence diagnostics
  plot_diagnostics.R               # plot_survey_coverage/plot_sightings/plot_occupancy_map/plot_covariate_map
  cleanup.R                        # cleanup_outputs(config) -> removes data file/figs
  main.R                           # run_occupancy_model(config_path) -> orchestrates all of the above
  msomgom-package.R                # package-level doc + centralized @import/@importFrom
man/                               # generated .Rd files (roxygen2, do not edit by hand)
vignettes/getting-started.Rmd      # worked end-to-end walkthrough on mock data
tests/testthat/                    # testthat suite (unit tests + end-to-end pipeline tests)
inst/
  extdata/configs/                 # example configs: bof_riwh.yaml, mock_test.yaml
  scripts/
    run_pipeline.R                 # CLI wrapper: Rscript run_pipeline.R config.yaml
tools/
  citations.csv                    # DOI/URL registry the shared citation engine checks
  citation-hooks.R                 # repo-specific citation checks (e.g. NARWC handbook version)
.github/workflows/check-citations.yaml  # calls chross22/distsamp's shared citation-check engine monthly
legacy/                            # archived pre-refactor scripts (gitignored, kept locally)
docs/refactor_plan.md              # full history of the generalization refactor

Citing msomgom

citation("msomgom")

That gives three entries: the package; the dynamic (multi-season, colonization/persistence) occupancy model of MacKenzie et al. (2003) that it implements; and the single-season detection model of MacKenzie et al. (2002) that the 2003 paper extends. Data products a run leans on carry their own citations separately - see References below. A methods section usually needs all of these.

References

Citations above reflect package versions installed at the time this was written (see Setup); run citation("pkgname") in R for the exact citation matching your installed version. A scheduled workflow runs monthly against the shared citation-check engine in chross22/distsamp, checking tools/citations.csv's registry against current CRAN/DOI metadata and running repo-specific checks from tools/citation-hooks.R (e.g. whether a newer NARWC handbook has been published); it opens an issue if anything needs review.

About

Fits a dynamic (multi-season) occupancy model for cetacean vessel-survey sightings in JAGS. Config-driven data cleaning, spatial gridding, environmental covariates, and diagnostics for any species or study area.

Topics

bayesian-statisticscetaceansecologygulf-of-mainejagsmarine-ecologymcmcoccupancy-modelingoceanographypopulation-ecologyrr-packageright-whalespecies-distribution-modelling

Resources

Readme
Activity

Stars

0 stars

Watchers

0 watching

Forks

0 forks
Report repository

Releases

Packages

Contributors

Languages