false zeros from imperfect detection
Hierarchical observation models on the tulpa inference engine: occupancy, N-mixture abundance, distance, removal, cover.
A site reads zero because the species is absent, or because you visited on the wrong day.
tulpaObs fits the models that hold those two apart: a latent ecological state under an
imperfect detection process, with the state integrated out exactly so occupancy and
detection come back as separate estimates. One fitter, one family object per model, and
the spatial fields, temporal structure, random effects, and latent factors written inside
the formula.
install.packages("pak")
# pak::pak("gcol33/tulpaObs")
install.packages("tulpaObs", repos = c('https://gcol33.r-universe.dev', 'https://cloud.r-project.org')) # latest from GitHub
pak::pak("gcol33/tulpaObs@v0.1.1") # a specific tagged release
Tagged releases are listed at https://github.com/gcol33/tulpaObs/releases.
pak resolves the dependency tree, pulling tulpa and tulpaMesh from the r-universe
repository DESCRIPTION declares in Additional_repositories, at or above the Imports
floor. A C++17 toolchain is needed (Rtools on Windows, Xcode CLI tools on macOS,
r-base-dev on Linux); both tulpa and tulpaObs compile their backends on first install.
library(tulpaObs)
fit <- tobs(
~ elevation + spatial(graph = adj), # occupancy: fixed effects + a spatial field
data = camera_data,
family = occu(), # MacKenzie et al. (2002) single-season
detection = ~ wind + time, # the detection process
y = y # N x J detection-history matrix
)
summary(fit)
predict(fit, newdata = grid) # marginalized per-site psi
tobs() is the only fitter; the family object chooses the model, and the response shape
follows from it.
tobs(~ forest, data = sites, family = abun(K.max = 50), detection = ~ effort, y = counts)
tobs(~ moisture, data = plots, family = cover(response = "beta"), y = pct_cover)
tobs(~ elev, data = sites, family = ms_occu(), detection = ~ effort, y = y_array,
species = species_id)
| Constructor | Model | Reference |
|---|---|---|
occu() |
Single-season occupancy | MacKenzie et al. (2002) |
dyn_occu() |
Dynamic occupancy (colonization / extinction) | MacKenzie et al. (2003) |
t_occu() |
Multi-season occupancy, AR1 year effect | spOccupancy tPGOcc
|
int_occu() |
Integrated multi-source occupancy | Miller et al. (2019) |
dyn_int_occu() |
Multi-season integrated occupancy | spOccupancy tIntPGOcc
|
fp_occu() |
False-positive occupancy | Miller et al. (2011) |
royle_nichols() |
Abundance-induced detection heterogeneity | Royle & Nichols (2003) |
occu_ttd() |
Time-to-detection occupancy | Garrard et al. (2008) |
occu_multi() |
Multi-species co-occurrence occupancy | Rota et al. (2016) |
jsdm() |
Joint species distribution (no detection) | |
abun() |
N-mixture abundance (Poisson / negbin) | Royle (2004) |
dyn_abun() |
Open (Dail-Madsen) N-mixture | Dail & Madsen (2011) |
count() |
Count / relative-abundance GLMM (no detection) | spAbundance abund
|
ms_count() |
Community relative-abundance GLMM (no detection) | spAbundance msAbund
|
distance() |
Binned distance sampling | Buckland et al. (2001) |
distsamp_open() |
Open-population distance sampling | Sollmann et al. (2015) |
gdistremoval() |
Joint distance + removal | Amundson et al. (2014) |
removal() |
Removal sampling | Dorazio et al. (2005) |
double_observer() |
Double-observer abundance | Nichols et al. (2000) |
cover() |
Cover hurdle (beta / lognormal / ordinal) | |
occu_categorical() |
Presence + nominal-class hurdle | |
occu_cover() |
Joint occupancy + cover hurdle | |
occu_multiscale_cover() |
Three-level (cell / plot / visit) + cover | Nichols et al. (2008) |
ms_occu() |
Community occupancy | spOccupancy msPGOcc
|
ms_dyn_occu() |
Community dynamic occupancy | |
ms_int_occu() |
Community integrated occupancy | |
ms_abun() |
Community N-mixture | spAbundance msNMix
|
ms_distance() |
Community distance sampling | Sollmann et al. (2016) |
ms_occu_cover() |
Community joint occupancy + cover |
Every family carries the same structured-effect vocabulary, so any of them composes with the engine's spatial, temporal, and latent structure.
Structured terms are written in the formula, the way lme4, mgcv, and INLA do. The
parser resolves bare symbols against data columns, so a spatial field is one term in the
formula:
~ elevation +
spatial(graph = adj) + # areal: icar / bym2 / car / car_proper
spatial(lon, lat, model = "spde") + # continuous: spde / gp / multiscale_gp
temporal(year, type = "ar1") + # temporal: ar1 / rw1 / rw2 / iid
(1 | observer) # random intercept (lme4 bar syntax)
spatial() and temporal() are the two front doors: model = picks the spatial
model ("icar" by default), type = the temporal one.
A coefficient that varies over space is a spatial field, so it is written under
spatial() and model = says whether space is a graph or a set of coordinates.
The areal form is a bar over graph nodes and has the broad family coverage:
~ spatial(~ 1 + year || cell, graph = adj) # areal SVC, method = "nested_laplace"
The continuous form names the varying coefficients directly and asks for a range
prior (P(range < r0) = alpha):
~ elevation + spatial(lon, lat, model = "svc",
coefficients = "elevation", prior.range = c(50, 0.05))
The continuous form fits on single-season occu() under laplace /
nested_laplace / nuts, and on removal(), distance(), fp_occu() and
dyn_abun() under laplace / nested_laplace; the surfaces come back as
fit$svc_field. Elsewhere the term errors and points at the areal bar.
svc(lon, lat, coefficients = ) is the direct constructor, as icar() is for
spatial(model = "icar").
A term enters whichever process it is written in, so detection = ~ (1 | observer) puts
the random effect on detection. share("id") shares one realization across both processes.
?tobs_terms documents every term and its arguments; ?svc and ?icar land there too.
Bar syntax follows lme4 throughout: (1 | site), (x | site) for a correlated random
slope, (x || site) uncorrelated, (1 | g:h) and (1 | g/h) for crossed and nested
factors.
Advanced usage: each spatial model also has a direct constructor, and spatial() forwards
its arguments to it, so the two forms build the same term.
icar(graph = adj) # spatial(graph = adj)
bym2(graph = adj) # spatial(graph = adj, model = "bym2")
spde(lon, lat) # spatial(lon, lat, model = "spde")
gp(lon, lat, cov = "matern", # spatial(lon, lat, model = "gp", cov = "matern",
prior.range = c(0.1, 0.05)) # prior.range = c(0.1, 0.05))
Model-specific arguments are checked against the model named in model =, so an argument
that model does not take is an error. gp() requires a PC prior on the range:
prior.range = c(r0, alpha) encodes P(range < r0) = alpha in the units of the
coordinates.
An occupancy site contributes a mixture over its latent state: occupied with probability
psi and detected on a Bernoulli schedule, or unoccupied and silent. tulpaObs evaluates
that marginalized likelihood directly, and the same holds for the N-mixture abundance N,
so psi and p come back as separate processes. Fitting the detection histories as one
stacked binomial identifies only their product psi * p; the
occupancy-vs-INLA vignette
runs that head-to-head on simulated data with known truth.
method names one fully specified route. Laplace is fast and the default; the stochastic
corrections debias it where it is biased; NUTS gives the full posterior on eighteen of the
families, the structured terms it samples listed on each family's help page.
tobs(..., method = "laplace") # EM + Laplace, deterministic
tobs(..., method = "laplace_sla") # skew-corrected marginals, adds a `skew` column
tobs(..., method = "laplace_gibbs") # post-EM Gibbs chain, Rubin-pooled
tobs(..., method = "nested_laplace") # multi-block: areal fields, joint cover hurdle
tobs(..., method = "nuts") # HMC, reports split-Rhat and bulk / tail ESS
Single-season occu(), t_occu() and the community families ms_occu(), ms_dyn_occu(),
ms_int_occu(), jsdm() and ms_count() also take method = "pg_gibbs", an exact
Polya-Gamma Gibbs sampler that recovers the community variance components the Laplace EM
attenuates. fit$method records
the resolved route, and an unsupported combination errors with a pointer to the methods
that family accepts.
Random-effect variance components carry the Laplace small-cluster bias for binary data (the
glmer nAGQ = 1 regime). An adaptive Gauss-Hermite refinement on the exact per-group
marginal runs by default (re.aghq), cutting the occupancy per-group-n = 8 sigma bias
from ~18% to ~4% and the detection-arm bias from ~70% to ~1%, both against NUTS.
Detection histories arrive as matrices, as long tables with one row per site-visit, or as
species arrays. Three helpers build the same tobs_data object from any of them:
tobs_format(y, occ.covs = sites, det.covs = list(wind = wind_mat), coords = xy)
tobs_data(visits_df, y = "detected", site = "site_id", visit = "visit")
tobs_format_ms(y_array, occ.covs = sites, species.names = spp)
occu_cover_inputs(plot_df, site = "plot", visit = "visit",
response = "present", y.pos = "cover")
summary(), plot(), and print() on a tobs_data report naive occupancy and detection,
per-visit rates, completeness, and a detection map when coordinates are present. Three
fixed-seed example datasets ship with a known truth attached: peatland_occu,
foray_counts, and meadow_cover.
Every family has a matching simulator that returns y, data, and the truth that
generated them, so a model can be checked against known parameters before it meets real
data. The package's own recovery tests are built on these:
sim <- simulate_occu(N = 200, J = 4,
beta.occ = c(-0.5, 1.2, 0.4),
beta.det = c(0.2, -0.8),
seed = 1)
fit <- tobs(~ occ_cov1 + occ_cov2, data = sim$data, family = occu(),
detection = ~ det_cov1, y = sim$y)
coef(fit)
sim$truth
simulate_abun(), simulate_cover(), simulate_dyn_occu(), simulate_ms_occu(),
simulate_distance(), and the rest cover the remaining families, including a
simulate_cover_joint() that shares a demeaned BYM2 field across both hurdle arms.
A fit is an ordinary R model object, whichever route produced it:
coef(fit); confint(fit); vcov(fit); logLik(fit)
fitted(fit); residuals(fit); predict(fit); simulate(fit)
ranef(fit); tidy(fit); glance(fit); update(fit)
converged(fit); convergence(fit)
fitted() returns list(psi, p, z) at the posterior mean, ranef() per-group BLUPs with
standard errors, and predict() works in-sample, on a design matrix (X.0=), or on terms
(terms=, which returns a tobs_prediction with a plot() method). Spatial fits
interpolate to new coordinates through tobs_predict_spatial(). $ accessors
(beta.samples, psi.samples, z.samples, run.time) follow spOccupancy naming.
check_model(fit, coords = xy) # roll-up: sampler, WAIC, PPC, zero-inflation, Moran's I
tobs_check_id(fit$model) # confounding, low detection, sparse data
waic(fit); dic(fit); cpo(fit)
ppc(fit, fit.stat = "chi-squared")
test_dispersion(fit); test_zero_inflation(fit); test_outliers(fit)
test_uniformity(pit_residuals(fit))
waic() and loo() are the loo generics, so loo::loo_compare() reads a tobs_fit
directly.
Moran's I, Durbin-Watson, variograms, and spatial_range / temporal_corr for the latent
structure come through from tulpa. tobs_stack() LOO-weights two or more fits on the same
observations into an ensemble whose predict() returns the combined predictive.
within_between() splits a covariate into its per-group mean and within-group deviation, so
~ x_btw + x_wtn separates the cross-site gradient from change inside a site across a
resurvey:
d <- within_between(plots, group = "plot", vars = c("temperature", "grazing"))
tobs(~ temperature_btw + temperature_wtn, data = d, family = cover(), y = d$cover)
For community fits, tobs_richness() returns posterior species richness and
tobs_associations() the residual species-association matrix from the latent factors.
Getting started
Model families
Structure
API.md is the full reference: every
argument of tobs(), the term registry, the method table, the control list, and the
complete export list.
The public API (tobs(), the family constructors, and the tobs_* S3 classes) is
stable. Internal entry points (.tobs_build_model(), .tobs_fit_model(),
.tobs_laplace()) and the C++ surface follow tulpa's ABI version and may change between
releases.
"Software is like sex: it's better when it's free."
Linus Torvalds
I'm a PhD student who builds R packages in my free time because I believe good tools should be free and open. I started these projects for my own work and figured others might find them useful too.
If this package saved you some time, buying me a coffee is a nice way to say thanks. It helps with my coffee addiction.
MIT (see the LICENSE file)
@software{tulpaObs,
author = {Colling, Gilles},
title = {tulpaObs: Hierarchical Latent-State Observation Models via tulpa},
year = {2026},
url = {https://github.com/gcol33/tulpaObs}
}