Skip to contents

Fit a Bayesian presence-only niche model

Source: R/fit_niche.R
fit_niche.Rd

Compiles and fits a Bayesian fundamental niche model using Stan via cmdstanr. Estimates the fundamental niche as a multivariate normal ellipsoid from occurrence data only.

Usage

fit_niche(
  env_occ,
  eta = 1,
  sigma_rate = 0.1,
  model = NULL,
  grainsize = NULL,
  chains = 4,
  iter_warmup = 1000,
  iter_sampling = 1000,
  seed = NULL,
  ...
)

Arguments

env_occ

Data frame or matrix of environmental values at presence points (rows = points, columns = environmental variables).

eta

Numeric, shape parameter for the LKJ prior on the correlation matrix. Default is 1 (uniform over correlation matrices). Larger values concentrate mass near the identity matrix.

sigma_rate

Numeric, rate parameter for the exponential prior on standard deviations. Default is 0.1 (weakly informative).

model

Optional; either NULL (default, compile on the fly), a character string with the model name ("presence_only"), or a pre-compiled CmdStanModel object. Passing a pre-compiled model avoids redundant compilation when running multiple fits in parallel (e.g., on HPC).

grainsize

Integer, grainsize for reduce_sum parallelization. If NULL or 0, Stan chooses automatically. Used only when the model is compiled with stan_threads = TRUE.

chains

Integer, number of MCMC chains (default 4).

iter_warmup

Integer, number of warmup iterations per chain (default 1000).

iter_sampling

Integer, number of sampling iterations per chain (default 1000).

seed

Integer, random seed for reproducibility.

...

Additional arguments passed to cmdstanr::CmdStanModel$sample().

Value

A CmdStanMCMC object with the fitted model.

Details

Requires the cmdstanr package and a working CmdStan installation. Install cmdstanr with:


install.packages("cmdstanr", repos = c(
  "https://stan-dev.r-universe.dev", getOption("repos")
))
cmdstanr::install_cmdstan()

For parallel HPC workflows, pre-compile the model once and pass the CmdStanModel object to avoid concurrent compilation:


# Single compilation (e.g., in a setup job) with thread support
mod <- compile_model(stan_threads = TRUE)

# Run with threading (e.g., 4 threads per chain)
fit_niche(occ, model = mod, grainsize = 10, chains = 2,
          threads = 4)  # passed via ... to mod$sample()

Examples

if (FALSE) { # \dontrun{
# Simulate 2D environmental data
set.seed(42)
occ <- data.frame(
  bio01 = rnorm(50, 20, 3),   # annual mean temperature (C)
  bio12 = rnorm(50, 800, 100) # annual precipitation (mm)
)

# Fit the presence-only model (compiles on the fly)
fit <- fit_niche(occ, chains = 2)

# Summarize posterior
fit$summary()

# Extract parameters
pars <- extract_pars(fit)
pars$mu     # centroid (optimal environment)
pars$sigma  # tolerances (breadth per dimension)
pars$R_corr # correlation (trade-offs)

# 95% most suitable region (Mahalanobis ellipsoid) in niche space:
# points x such that (x - mu)' Sigma^-1 (x - mu) <= qchisq(0.95, df=2)

# ---- Parallel HPC workflow ----
# Compile once with thread support
mod <- compile_model(stan_threads = TRUE)

# Run parallel fit (grainsize = 10, 4 threads per chain)
fit_par <- fit_niche(occ, model = mod, grainsize = 10,
                     chains = 2, threads = 4)
} # }