Package 'tulpaObs'

Description

Allows accessing spOccupancy-style fields (e.g., ⁠$beta.samples⁠, ⁠$psi.samples⁠) on tobs_fit objects. Since tobs stores actual posterior draws, this is a thin remapping layer.

Usage

## S3 method for class 'tobs_fit'
x$name

Arguments

Value

The requested field value.

Examples

sim <- simulate_occu(N = 100, J = 3, n.occ.covs = 1, n.det.covs = 1,
                     seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ det_cov1, y = sim$y, method = "laplace",
            control = list(verbose = FALSE, progress = FALSE))
dim(fit$beta.samples)
dim(fit$alpha.samples)

Description

Latent Poisson (or NB) abundance with binomial detection per visit (Royle 2004).

Usage

abun(K.max = NULL, mixture = c("poisson", "negbin", "zip", "zinb"))

Arguments

Details

Royle's (2004) N-mixture model: latent abundance N_i ~ Poisson(lambda_i) with ⁠log lambda_i = X_lambda beta_lambda⁠, and counts y_ij | N_i ~ Binomial(N_i, p_ij) with ⁠logit p_ij = X_p beta_p⁠. The abundance formula is the tobs() formula; the per-visit detection formula is detection. The marginal likelihood integrates N out exactly, so no route runs an EM over N.

Value

A tobs_family object.

Engines

• method = "laplace" – a direct Laplace approximation over the closed-form marginal with an observed-Fisher covariance; a grouped random effect on one arm is integrated by adaptive Gauss-Hermite quadrature.

• method = "nested_laplace" – an areal icar() / bym2() / car_proper() field on the abundance arm, its hyperparameters integrated on an outer grid (grid-integrated coefficient covariance).

• method = "nuts" – the in-tree C++ NUTS sampler over the same marginal, Poisson or negative binomial (log_r sampled), warm-started at the Laplace mode; with an areal term a fixed-hyper non-centered field on the abundance arm (icar() / car_proper() / bym2()); "zip" / "zinb" are "laplace"-only, and a field does not combine with a random effect here.

Examples

sim <- simulate_abun(N = 120, J = 4, n.abund.covs = 1, n.det.covs = 1, seed = 1)
fit <- tobs(~ abund_cov1, data = sim$data, family = abun(),
            detection = ~ det_cov1, y = sim$y, method = "laplace")
summary(fit)

Description

The tulpa::check_model() method for tobs_fit: the criteria, the posterior-predictive check and the goodness-of-fit tests in one call, printed as a report and drawn as the diagnostic panel.

Usage

## S3 method for class 'tobs_fit'
check_model(
  object,
  coords = NULL,
  n.samples = 250,
  nsim = NULL,
  seed = NULL,
  plot = TRUE,
  ...
)

Arguments

Details

The engine's default method builds its panel from fitted() and residuals() as numeric vectors, which a latent-state fit does not have – fitted() returns list(psi, p, z) and residuals() one series per process, so the response cannot be resolved. This method reads the same quantities through the family's own doors instead: pit_residuals(), ppc(), test_dispersion() and tulpa::moran_i().

Value

Invisibly, the diagnostic results: waic, dic, cpo, ppc, zero_inflation, dispersion, pit, uniformity, and moran when coords is supplied.

Examples

sim <- simulate_occu(N = 100, J = 3, n.occ.covs = 1, n.det.covs = 1,
                     seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ det_cov1, y = sim$y, method = "laplace",
            control = list(verbose = FALSE, progress = FALSE))
check_model(fit, seed = 1, plot = FALSE)

Description

Per-species coefficients from a batched occu_cover fit.

Usage

## S3 method for class 'tobs_batch'
coef(object, ...)

Arguments

Value

A named list of per-species coefficient vectors.

Examples

sims <- lapply(1:2, function(s) simulate_occu_cover(N = 40, J = 3, seed = s))
y <- lapply(sims, function(s) s$y)
y_pos <- lapply(sims, function(s) { p <- s$y_pos; p[is.na(p)] <- 0; p })
fit <- tobs(~ occ_cov1, data = sims[[1]]$data,
            family = occu_cover("lognormal"), detection = ~ 1,
            y = y, y.pos = y_pos, species = c("a", "b"), method = "laplace")
coef(fit)

Description

The fixed-effect estimates as one named numeric vector, named ⁠<arm>_<term>⁠ as stats::vcov() and stats::confint() name them: psi_(Intercept), p_det_cov1, lambda_elev, ... The visit-level detection coefficients (⁠p_visit_<cov>⁠) belong to the detection arm. A coordinate that enters no single arm, such as the log_r overdispersion, keeps its own name.

Usage

## S3 method for class 'tobs_fit'
coef(object, arm = NULL, ...)

Arguments

Value

A named numeric vector.

Examples

sim <- simulate_occu(N = 100, J = 3, n.occ.covs = 1, n.det.covs = 1,
                     seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ det_cov1, y = sim$y, method = "laplace",
            control = list(verbose = FALSE, progress = FALSE))
coef(fit)
coef(fit, arm = "psi")

Description

The public accessor for whether a fit converged, with one return shape across every tobs() family. Every family's optimiser / EM / sampler verdict is recorded under fit$convergence, and a Laplace-family fit that stops without meeting its criterion warns once, naming control$max.iter. convergence() returns the full record (converged, n_iter, sla_status when the simplified-Laplace marginals were used, and the per-group AGHQ solve status on the families fitted by that engine), and converged() returns the single logical.

Usage

convergence(object, ...)

## S3 method for class 'tobs_fit'
convergence(object, ...)

converged(object, ...)

## S3 method for class 'tobs_fit'
converged(object, ...)

Arguments

Value

convergence(): a list with converged (logical), n_iter (integer, NA for grid / closed-form fits with no iteration count), and sla_status (character, when present). A fit whose random effects were integrated by adaptive Gauss-Hermite quadrature also carries group_ok (logical, one entry per group: FALSE where that group's posterior solve failed and its estimates are NA), groups_failed (their indices) and groups_failed_names (their labels where the family names its groups – species, for the community models). converged is FALSE whenever any group failed, since the quantities that group contributes to are not estimates of anything; filter on these rather than on the text of the warning the engine raises. A single-species fit whose grouped random effects were integrated that way carries re_boundary: one record per random-effect SD (⁠sigma_<g>_<coef>⁠) and, on a correlated block, per correlation (⁠cor_<g>_<ci>_<cj>⁠), with the boundary Wald statistic, its critical value and distinguishable; FALSE means the SD cannot be told from zero or the correlation from +-1, and the fit says so in one message ("boundary (singular) fit"). converged(): a single TRUE / FALSE.

Examples

sim <- simulate_occu(N = 100, J = 3, n.occ.covs = 1, n.det.covs = 1,
                     seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ det_cov1, y = sim$y, method = "laplace",
            control = list(verbose = FALSE, progress = FALSE))
converged(fit)
convergence(fit)$n_iter

Description

A generalized linear model on an observed count (or continuous) response directly, with no detection sub-model and no latent abundance to marginalize. The abundance analogue of jsdm(): where abun() fits Royle's N-mixture (latent N plus imperfect detection), count() fits the observed response as a plain GLMM. This is the "relative abundance" model of spAbundance (abund): one value per site, log-link Poisson / negative-binomial counts or an identity-link Gaussian response.

Usage

count(response = c("poisson", "negbin", "gaussian", "binomial"))

Arguments

Details

The response is a numeric vector (one value per site), supplied via ⁠y =⁠ or on a two-sided formula left-hand side (count.value ~ predictors).

Value

A tobs_family object.

Engines

• method = "laplace" – one tulpa GLMM block; the negative-binomial size and the Gaussian residual variance are estimated by an outer dispersion loop.

• method = "nested_laplace" – a plain areal field (icar() / car_proper()) on the formula as a latent GMRF, integrated over its hyperparameters by the nested-Laplace EM; the areal term is required on this route.

See Also

abun() (N-mixture: latent abundance + detection), jsdm() (occurrence GLMM, no detection).

Examples

f <- count("poisson")
f

Description

Latent presence (Bernoulli) plus conditional positive cover (beta or lognormal). Does not share the replicate-detection assumption of the other families – see vignette("families") for the conceptual caveat.

Usage

cover(
  response = c("beta", "beta_oi", "lognormal", "lognormal_trunc", "ordinal", "gaussian"),
  breaks = NULL
)

Arguments

Value

A tobs_family object.

Engines

• method = "laplace" – two independent tulpa_laplace() fits, one per arm (presence, positive), non-spatial; the default.

• method = "laplace_sla" – the same two fits with simplified-Laplace skew-corrected marginals.

• method = "nested_laplace" – the joint nested-Laplace engine: shared or arm-specific areal fields, weighted trend fields, temporal() and re() blocks, every hyperparameter integrated on an outer grid; the route the "lognormal_trunc" and "ordinal" arms require.

• method = "nested_laplace_sla" – the nested route with skew-corrected marginals.

• method = "nuts" – the non-spatial sampler over the exact two-arm coefficient vector c(beta_presence, beta_positive, log_disp), giving calibrated intervals and a per-draw pointwise likelihood for WAIC / LOO; "beta", "lognormal" and "gaussian" arms only, and any structured term errors with a pointer to "nested_laplace".

Response on the formula left-hand side

The cover response is a single length-N vector, so it may sit on the top formula left-hand side and ⁠y =⁠ is dropped. The two calls are equivalent:

# response on the LHS (y = omitted)
tobs(cover.flat ~ time.sc + habitat +
       spatial(~ 1 + time.sc || cell_idx, graph = adj),
     data = dat, family = cover(response = "beta"),
     method = "nested_laplace")

# the same fit with a one-sided formula and an explicit y =
tobs(~ time.sc + habitat +
       spatial(~ 1 + time.sc || cell_idx, graph = adj),
     y = dat$cover.flat, data = dat, family = cover(response = "beta"),
     method = "nested_laplace")

Naming the response makes the per-arm spatial labels read naturally: cover() splits cover.flat into a presence arm and a positive arm, the arm names that per-arm formulas and share() address. The LHS is evaluated against data (then the calling environment), so it may be a bare column or an expression.

Joint nested-Laplace engine – spatial-prior parameterisation

When fitted with method = "nested_laplace" and an areal spatial term in the latent-presence formula (bym2(graph = adj), or car() / car_proper()), the engine identifies a single latent field z per region and parameterises the two arms as ⁠eta_occ = X beta_occ + sigma * z⁠ and ⁠eta_pos = X beta_pos + alpha * sigma * z⁠. Both BYM2 sub-blocks (phi, theta) are subject to hard sum-to-zero constraints (see .joint_inner_var() for the math): the cover-arm intercept beta_pos[1] is identified as the population mean of eta_pos under mean(z) = 0. Any simulator generating data for this engine must demean both phi_f and theta_f before scaling, otherwise the estimator targets beta_pos_0_truth + alpha * mean(w_s_sim) and coverage of the population truth collapses with alpha. See simulate_cover_joint() for a ready-made demeaned simulator.

Spatially varying trend (a weighted areal term)

A spatially varying trend is model structure, so it is declared in the formula as a SECOND, weighted areal term on the same graph as the intercept field, not through control:

~ time.sc +
  icar(graph = adj, group.var = "cell_idx") +
  icar(graph = adj, weight = time.sc, group.var = "cell_idx")

The unweighted term is the shared intercept field; the weighted term icar(..., weight = col) is the trend field, whose contribution to each arm's predictor is weight_i * z[cell_i], coupled onto the cover arm with its own scale (alpha_trend, reported in fit$alpha_trend / fit$sigma_trend) integrated over the outer grid. The umbrella spelling spatial(graph = adj, model = "icar", weight = col) resolves identically. The trend block inherits the intercept block's amplitude; give it its own with share(spatial(), terms = list(intercept = ..., trend = ...)). Requires method = "nested_laplace". A coupled trend cannot currently combine with temporal() / re() blocks in the same fit.

The axis the amplitude rides carries prior structure – a point mass at alpha = 0 ("no coupling") and a log-spaced slab above it – so it is set in one of two ways. share(alpha = grid(...)) STATES its nodes, and with them that structure. share(alpha = grid(n = 9)) states a RESOLUTION: the engine re-reads its own axis with that many slab nodes, point mass and slab bounds unchanged, so sharpening the axis never restates its structure. ⁠terms =⁠ gives either, per block. The resolution is the only way to raise this axis, because it does not densify when the donor control$sigma.grid does: on an informative data set the outer grid's quadrature effective sample size saturates on the copy amplitude while every other axis tracks the request (measured engine-side, NOTES_measurements.md). One block takes one of the two; giving both is an error. control$alpha.grid[.trend] and control$alpha.n[.trend] are the wire format these compile into and are refused as user input.

Set it when the fit's hyperparameter intervals are reported. On the coupled SBC fixture the declared resolution leaves the field SD miscalibrated once the data are informative – at 10 visits per site its uniformity p-value is 9.1e-05, against 0.17 at 3 visits – and control$alpha.n = 21 returns every scored quantity to nominal at both. Thirteen slab nodes do not (the field SD still reads 4.5e-03), and the price is a 2.3-3x longer fit, since the outer grid is a tensor. Point estimates are not affected the way the intervals are; the measurement is in NOTES_measurements.md.

Varying-coefficient spatial bar (the compact single-term form)

The intercept field plus its weighted trend field can also be written as one spatial() term carrying an lme4-style coefficient formula:

~ time.sc + habitat +
  spatial(~ 1 + time.sc || cell_idx, graph = adj)

The bar left-hand side spells the coefficient fields: the intercept column (1) is the unweighted field; each covariate column (time.sc) is a weight-scaled coefficient field (weight_i * z[cell_i]). The bar right-hand side (cell_idx) is the graph node index (the areal group.var); || requests independent intercept and slope fields, a single | makes them correlated. This desugars to exactly the two-term weighted-areal form above, so the two spellings give the same fit.

Choosing a field's arm

placement and share(): The cover hurdle's two arms are presence (the y > 0 Bernoulli arm) and positive (the y | y > 0 arm); summary() and the coefficient output print these same labels. A field is placed on an arm by writing it in that arm's per-arm formula (presence = ~ ..., positive = ~ ...), and shared across arms with share(). A field in the single shared formula reaches both arms.

Both arms (shared / copied). A field in the shared formula is one presence-anchored latent copied onto the positive arm with an estimated coupling per coefficient field (alpha for the intercept field, alpha_trend for the trend field), marginalized on the outer grid:

# one latent field, shared by both arms
tobs(cover.flat ~ time.sc + habitat +
       spatial(~ 1 + time.sc || cell_idx, graph = adj),
     data = dat, family = cover(response = "beta"), method = "nested_laplace")

presence: eta_presence = ... + u_cell + time.sc * s_cell
positive: eta_positive = ... + alpha * u_cell + alpha_trend * time.sc * s_cell

Per-arm formulas make the shared field explicit: place it on presence and copy it onto positive. share(spatial()) estimates the coupling on the default axis; share(spatial(), alpha = grid(c(...))) integrates it over supplied nodes, alpha = grid(n = 9) over the engine's own axis read at n slab nodes, and alpha = 0.5 fixes it. share(spatial(), prior = list(...)) regularizes the coefficient itself; a fit copying both the intercept and a weighted trend block is refused it, since one prior reaches the engine per fit.

tobs(presence = ~ time.sc + habitat +
                  spatial(~ 1 + time.sc || cell_idx, graph = adj),
     positive = ~ time.sc + habitat + share(spatial()),
     data = dat, family = cover(response = "beta"), method = "nested_laplace")

One arm only (free / separate). A field written in one arm's formula is a separate latent on that arm alone, with its own precision and no cross-arm coupling. A field in each arm's formula gives two independent latents:

tobs(presence = ~ time.sc + spatial(~ 1 + time.sc || cell_idx, graph = adj),
     positive = ~ time.sc + spatial(~ 1 + time.sc || cell_idx, graph = adj),
     data = dat, family = cover(response = "beta"), method = "nested_laplace")

presence: eta_presence = ... + u_presence_cell + time.sc * s_presence_cell
positive: eta_positive = ... + u_positive_cell + time.sc * s_positive_cell

The free fit reports sigma_armspecific. Arm-specific fields are their own spatial structure: they do not combine with a shared field, a weighted trend, or temporal() / re() in the same formula, and at most one targets each arm.

The || and | axis is separate from shared / free: || makes the intercept and slope fields independent, while a single | makes them correlated (a free cross-covariance, MCAR). A correlated | bar shared across both arms (in the shared formula, or on presence with share() on positive) is copied with one estimated amplitude; placed on one arm's formula it is a free-Sigma correlated field on that arm alone, no cross-arm copy.

field in the shared formula, no share()   one presence-arm latent, pinned at alpha = 0
field + share() (either spelling)         one shared / copied latent (presence anchor, coupling estimated)
a field in each arm's formula             separate / free latents, no coupling
||                                        independent intercept and slope coefficient fields
|                                         correlated (MCAR) coefficient fields, copy-only

Without a share() (and without control$alpha.grid) the presence field is NOT carried onto the cover arm: the amplitude is pinned at alpha = 0 and the field rides the presence arm alone. Structure nobody wrote is not in the model, the reading occu_cover() and occu_multiscale_cover() also have. Write + share(spatial()) in the shared formula, or place the field on presence and share(spatial()) on positive, to couple the two arms. A correlated | bar is the exception: that spelling is copy-only by definition, so the bar itself states the coupling.

A formula bar is a random effect, not a spatial field

A bare lme4 bar in the formula – (1 | cell), (1 + x | cell), (x || cell) – is a grouped random effect, not a spatial field, even when the grouping factor names the areal graph nodes. The engine's inline-MCAR call tulpa::spatial(graph, ~ 1 + x | cell) reads ⁠| cell⁠ as a separable spatial field, but the same spelling inside a cover() formula is parsed as an IID random effect on cell. For a spatial field on the cells write either the compact spatial() bar of the section above, spatial(~ 1 + x || cell, graph = adj), or the two-term weighted-areal form, ⁠icar(graph = adj, group.var = "cell") +⁠ icar(graph = adj, weight = x, group.var = "cell"). When a formula carries a bar whose grouping factor is also an areal term's group.var, cover() emits a one-time message noting the bar is fitted as a random effect; suppress it with base::suppressMessages().

Checkpoint / resume

A spatial cover-hurdle fit integrates over a large outer hyperparameter grid and can run for hours. control$checkpoint = list(path = "fit.ckpt", resume = TRUE) makes the joint engine append each completed grid cell to path; a resume = TRUE run loads the finished cells and solves only the rest, reproducing the from-scratch fit, so a killed or rebooted fit resumes instead of restarting. resume = FALSE starts a fresh file. Forwarded to tulpa::tulpa_nested_laplace_joint().

Examples

sim <- simulate_cover(N = 200, seed = 1)
fit <- tobs(~ x, data = sim$data, family = cover("beta"), y = sim$y)
summary(fit)

Description

Builds an opt-in fixed-effect prior for cover() models fit with a Laplace method. The cover hurdle has two arms with their own coefficients – an occurrence arm (Bernoulli presence) and a positive-cover arm (beta or lognormal) – each with an intercept and a slope bucket. The penalty is the quadratic ⁠+ sum((beta_j - mu_j)^2 / (2 sd_j^2))⁠ added to each arm's negative log-likelihood, the same mechanism occu_priors() uses.

Usage

cover_priors(
  occ.intercept = list(mean = 0, sd = 2),
  occ.slope = list(mean = 0, sd = 2.5),
  pos.intercept = list(mean = 0, sd = 3),
  pos.slope = list(mean = 0, sd = 2.5)
)

Arguments

Details

Cover priors are opt-in: without priors, a cover fit is unpenalised. The main use is regularising perfect separation in the occurrence arm at small N. Set any sd = Inf to drop that component.

Value

A cover_priors object, ready to pass to tobs(..., priors = ...).

Coverage

Applied on the separate-Laplace path (method = "laplace" / "laplace_sla") when the formula has no spatial term. Both arms are penalised: the occurrence and lognormal-positive arms through tulpa::tulpa_laplace(), and the beta-positive arm through tulpa::tulpa_laplace_beta(). The joint nested-Laplace path (method = "nested_laplace" / "nested_laplace_sla", used for spatial cover formulas) also threads the priors, as a per-arm beta_prior_mean / beta_prior_prec on the joint engine's occurrence and positive responses. The separate-Laplace path keeps rejecting a prior alongside a spatial term (that solver carries its own); add the spatial term through the nested-Laplace path to combine the two.

See Also

occu_priors()

Examples

# regularise the occurrence arm, leave the positive arm unpenalised
priors <- cover_priors(pos.intercept = list(mean = 0, sd = Inf),
                       pos.slope     = list(mean = 0, sd = Inf))
priors


sim <- simulate_cover(N = 150, seed = 1)
fit <- tobs(cover ~ x, data = sim$data,
            family = cover(response = "lognormal"),
            method = "laplace", priors = priors)
coef(fit)

Description

Latent abundance N_i ~ Poisson(lambda_i) (or negative binomial) in a covered region, observed through a half-normal or hazard-rate detection function over distance bins. With B bins the detected counts are multinomial over ⁠(bin 1, ..., bin B, undetected)⁠ with cell probabilities ⁠pi_b = integral_bin g(x; sigma) f(x) dx⁠ and ⁠1 - sum_b pi_b⁠, where f(x) is the distance density (uniform for a line transect, proportional to distance for a point transect). The latent N is summed out in closed form (truncation K.max), so the fit is a direct Laplace approximation (no EM), with a NUTS path over the same marginal; the per-bin detection integrals are evaluated by Gauss-Legendre quadrature.

Usage

distance(
  key = c("halfnorm", "hazard"),
  transect = c("line", "point"),
  cutpoints = NULL,
  K.max = NULL,
  mixture = c("poisson", "negbin")
)

Arguments

Details

The tobs() formula is the abundance (⁠log lambda⁠) model; detection is the site-level detection-scale (⁠log sigma⁠) model. The response y is an ⁠n_sites x n_bins⁠ integer matrix of per-bin detected counts. The bin edges and transect geometry travel with the family: distance(cutpoints = ...).

Value

A tobs_family object.

Engines

• method = "laplace" – a direct Laplace approximation over the closed-form marginal (Poisson or negative binomial); a grouped random effect on the abundance arm is integrated by adaptive Gauss-Hermite quadrature.

• method = "nested_laplace" – an areal icar() / car_proper() / bym2() field on the abundance arm, optionally composed with a temporal() field, through the shared areal-BFGS driver; a hazard-rate fit recovers the shape alongside the field.

• method = "nuts" – the in-tree C++ NUTS sampler over the same marginal, with a fixed-hyper areal or temporal field block on the abundance arm (one of the two, not both); the hazard-rate log-shape rides alongside.

References

Buckland, S. T., Anderson, D. R., Burnham, K. P., Laake, J. L., Borchers, D. L., Thomas, L. (2001). Introduction to Distance Sampling. Oxford. Royle, J. A., Dawson, D. K., Bates, S. (2004). Modeling abundance effects in distance sampling. Ecology 85, 1591-1597.

Examples

sim <- simulate_distance(N = 200, key = "halfnorm", transect = "line",
                         n.abund.covs = 1, n.sigma.covs = 1, seed = 1)
fit <- tobs(~ abund_cov1, data = sim$data,
            family = distance(key = "halfnorm", transect = "line",
                              cutpoints = sim$cutpoints),
            detection = ~ sigma_cov1, y = sim$y, method = "laplace")
summary(fit)

Description

The unmarked distsampOpen model: a Dail-Madsen open N-mixture (as dyn_abun()) observed by distance sampling at each primary period – the open-population counterpart of the single-season gdistremoval(). Initial abundance N_1 ~ Poisson(lambda), then ⁠N_t = Binomial(N_{t-1}, omega) + Poisson(gamma)⁠ (apparent survival omega + recruitment gamma); at each primary period the detected birds are distance-sampled into bins.

Usage

distsamp_open(
  transect = c("line", "point"),
  cutpoints = NULL,
  K.max = NULL,
  mixture = c("poisson", "negbin", "zip", "zinb"),
  dynamics = c("constant", "notrend", "trend", "autoreg", "ricker", "gompertz")
)

Arguments

Details

The distance-band allocation is conditional on the period total detected, so it factors out of the abundance HMM (the gdistremoval() trick): the marginal is the dyn_abun() forward recursion with the detection probability set to the overall distance detection pdist, plus the per-period band multinomials.

Four site-level arms: log abundance lambda (the tobs() formula), log distance scale sigma (detection), logit survival omega (omega = ~ ...), and log recruitment gamma (gamma = ~ ...), the last two intercept-only by default.

Value

A tobs_family object.

Inputs

y is a 3D array ⁠[n_sites x n_bins x n_seasons]⁠ of per-distance-band counts at each primary period (secondary occasions absorbed into the period total).

Engines

• method = "laplace" – the forward-HMM marginal maximised (BFGS) with an observed-information covariance; the only route (no field, no sampler).

References

Dail, D., Madsen, L. (2011). Models for estimating abundance from repeated counts of an open metapopulation. Biometrics 67, 577-587. Sollmann, R., Gardner, B., Chandler, R. B., Royle, J. A., Sillett, T. S. (2015). An open-population hierarchical distance sampling model. Ecology 96, 325-331.

See Also

gdistremoval() (single-season), dyn_abun() (open N-mixture), distance().

Examples

f <- distsamp_open(transect = "line", cutpoints = c(0, 10, 20, 30, 40),
                   dynamics = "notrend")
f

Description

Abundance from a double-observer protocol (unmarked multinomPois with a double-observer pi-function). Site abundance N ~ Poisson(lambda) is surveyed by two observers with detection p1 / p2. The state formula models ⁠log lambda⁠; detection the shared site-level per-observer detection design (observers 1 and 2 carry separate coefficients). By Poisson-multinomial thinning the observable cell counts are independent Poissons, so the marginal is closed form with no latent-abundance summation.

Usage

double_observer(type = c("independent", "dependent"))

Arguments

Value

A tobs_family object for tobs().

Independent vs dependent protocol

type = "independent" (the default): the two observers detect independently and each individual is recorded by observer 1 only, observer 2 only, or both, so y is an ⁠N x 3⁠ matrix of cell counts in that column order and the three cells are ⁠pi = (p1(1 - p2), (1 - p1)p2, p1 p2)⁠.

type = "dependent": a removal-style protocol where a primary observer records what it detects and a secondary observer records only what the primary missed, so y is an ⁠N x 2⁠ matrix ⁠(primary-detected, secondary-only)⁠ with cells ⁠pi = (p_pri, (1 - p_pri)p_sec)⁠. A single fixed primary observer gives only two cells for three parameters (lambda, p1, p2), a ridge that does not separate the two observer detections; observer role-swapping identifies them – pass ⁠primary =⁠ (a length-N vector in ⁠{1, 2}⁠ naming each site's primary observer), and with observer 1 primary at some sites and observer 2 at others ⁠(p_pri, p_sec)⁠ alternates between ⁠(p1, p2)⁠ and ⁠(p2, p1)⁠, giving the four cell means that recover ⁠(lambda, p1, p2)⁠. With p1 = p2 (an interchangeable pair) the dependent protocol reduces to a two-pass removal() model.

Engines

• method = "laplace" – the closed-form Poisson-multinomial marginal maximised (BFGS) with an observed-information covariance; the only route (no field, no sampler).

Examples

sim <- simulate_double_observer(N = 200, seed = 1)
fit <- tobs(~ abund_cov1, data = sim$data, family = double_observer(),
            detection = ~ det_cov1, y = sim$y, control = list(verbose = FALSE))
coef(fit)

# Dependent (role-swapping) protocol:
sd <- simulate_double_observer(N = 300, type = "dependent", seed = 1)
fd <- tobs(~ abund_cov1, data = sd$data, family = double_observer("dependent"),
           detection = ~ det_cov1, y = sd$y, primary = sd$primary,
           control = list(verbose = FALSE))
coef(fd)

Description

Latent abundance evolves across primary seasons via apparent survival and recruitment (Dail & Madsen 2011): N_1 ~ Poisson(lambda); for t >= 2, N_t = S_t + G_t with survivors ⁠S_t ~ Binomial(N_{t-1}, omega)⁠ and recruits G_t ~ Poisson(gamma); observed via Binomial(N_t, p) over secondary visits. The latent abundance sequence is summed out by an exact HMM forward recursion (it is not closed form, unlike the static abun()); analytic gradients come from forward-mode differentiation of the scaled forward algorithm, so the fit is a direct maximum-likelihood / Laplace fit with a NUTS path over the same marginal.

Usage

dyn_abun(K.max = NULL, mixture = c("poisson", "negbin", "zip", "zinb"))

Arguments

Details

Four arms: initial abundance lambda (the tobs() formula) and detection p (detection) are site-level; apparent survival omega (the omega argument, default ~ 1) and recruitment gamma (the gamma argument, default ~ 1) span the T - 1 transition intervals. A constant omega / gamma is shared across a site's seasons; supplying a season-varying covariate (a ⁠[n_sites x (T - 1)]⁠ matrix column of data, one column per transition interval) on omega / gamma gives interval-specific vital rates. The response y is a 3D array ⁠[n_sites x max_visits x n_seasons]⁠ (or a list of per-season count matrices); missing visits are NA.

Value

A tobs_family object.

Engines

• method = "laplace" – analytic-gradient BFGS over the forward marginal with an observed-information covariance; a grouped random effect on the initial-abundance arm is integrated by adaptive Gauss-Hermite quadrature.

• method = "nested_laplace" – an areal icar() / car_proper() field on the initial-abundance arm (BFGS over the forward marginal plus the CAR prior), optionally composed with a temporal() field, through the shared areal-BFGS driver.

• method = "nuts" – the in-tree C++ NUTS sampler over the same forward marginal, with a fixed-hyper areal field block on the initial-abundance arm; a temporal() term is not sampled.

References

Dail, D., Madsen, L. (2011). Models for estimating abundance from repeated counts of an open metapopulation. Biometrics 67, 577-587.

Examples

f <- dyn_abun(mixture = "negbin")
f

Description

Dynamic (multi-season) occupancy observed by several detection sources (spOccupancy tIntPGOcc): the product of a dynamic occupancy HMM (season-1 occupancy psi1, colonization gamma, extinction eps) and integrated occupancy (a per-season emission that pools multiple detection sources). Pooling sources across seasons is how colonization / extinction estimates come out of individually-sparse opportunistic data. The latent occupancy sequence integrates out by the two-state HMM forward recursion; the exact marginal is maximised with an observed-information vcov.

Usage

dyn_int_occu()

Details

y is a length-S list of ⁠[sites x visits x seasons]⁠ detection arrays (0/1/NA). The state formula models ⁠logit psi1⁠; colonization = ~ ... and extinction = ~ ... the site-level transitions (required, as in dyn_occu()); detection the shared per-source detection design (each source carries its own coefficients). v1: every source covers all sites and the same season grid, constant transitions, site-level detection.

Value

A tobs_family object for tobs().

Engines

• method = "laplace" – the non-spatial fit: the HMM-forward marginal maximised with analytic forward-backward gradients and an observed-information covariance.

• method = "nested_laplace" – a shared icar() field on the first-season occupancy formula (spOccupancy stIntPGOcc) through the shared areal-BFGS driver, the field hyperparameters integrated on an outer grid; that field is the only structured term accepted and is required on this route.

Examples

sim <- simulate_dyn_int_occu(N = 200, T.seasons = 4, S = 2, seed = 1)
fit <- tobs(~ 1, data = sim$data, family = dyn_int_occu(),
            detection = ~ 1, colonization = ~ 1, extinction = ~ 1,
            y = sim$y, sources = sim$sources, control = list(verbose = FALSE))
coef(fit)

Description

Bernoulli occupancy state with colonisation + extinction transitions across seasons (the MacKenzie et al. dynamic model). Colonisation / extinction are constant across a site's seasons by default; a season-varying rate is supplied as a ⁠[n_sites x (T - 1)]⁠ matrix column of data (one column per transition interval) on colonization / extinction. Detection is site-level by default; a ⁠[n_sites x T]⁠ matrix column of data (one column per primary season) on the detection formula gives per-season detection. All season-varying paths run under method = "laplace" (the exact HMM-forward marginal refine calibrates the coefficients); they are gated under "nuts" with a pointer.

Usage

dyn_occu()

Value

A tobs_family object.

Engines

• method = "laplace" – EM + Laplace over the exact HMM-forward marginal, the default and the route every season-varying path runs on.

• method = "laplace_sla" – the Laplace fit with simplified-Laplace skew-corrected marginals.

• method = "laplace_gibbs" / method = "laplace_mi" – the Laplace fit with a post-EM Gibbs / multiple-imputation variance correction.

• method = "nested_laplace" – an areal field on the occupancy arm through the nested-Laplace EM, its hyperparameters integrated on an outer grid.

• method = "nuts" – the unified C++ NUTS sampler over the HMM marginal; season-varying colonisation / extinction / detection error with a pointer.

Examples

f <- dyn_occu()
f

Description

Fitted values (occupancy and detection probabilities)

Usage

## S3 method for class 'tobs_fit'
fitted(object, ...)

Arguments

Value

A list with psi (occupancy probabilities), p (detection probabilities), and z (posterior state probability at posterior-mean parameters). For an int_occu() fit p is a named list, one detection probability vector per source (matching y's names), since each source carries its own detection design and, when present, its own detection-arm spatial field. For single-season, int_occu() and community models z is the per-site Bayes posterior P(z=1 | y) (an int_occu() site is occupied with certainty if ANY source detected it, and otherwise pools every source's non-detection evidence); for dynamic models z is the forward-backward (HMM smoothing) state posterior ⁠P(z_t=1 | y_{1:T})⁠ as an ⁠[n_sites x n_seasons]⁠ matrix.

Examples

sim <- simulate_occu(N = 100, J = 3, n.occ.covs = 1, n.det.covs = 1,
                     seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ det_cov1, y = sim$y, method = "laplace",
            control = list(verbose = FALSE, progress = FALSE))
str(fitted(fit))

Description

A binomial N-mixture dataset (Royle 2004) of repeated counts at 100 sites over 3 visits, drawn from a known generative model. Latent abundance increases with shrub cover and decreases with elevation; per-visit detection increases with survey effort. Synthetic, not field data; the generative coefficients and the realised latent abundances are in truth.

Usage

foray_counts

Format

A list with components:

y

100 x 3 integer matrix of counts.

occ.covs

Data frame of site covariates: elevation, shrub (both standardised).

det.covs

Named list with one 100 x 3 visit covariate matrix, effort (standardised).

truth

List of generative quantities: beta_lambda, beta_p, the realised latent abundance N, and per-site lambda.

See Also

abun(), tobs(), simulate_abun()

Examples

data(foray_counts)
summary(as.vector(foray_counts$y))

fit <- tobs(~ elevation + shrub, data = foray_counts$occ.covs,
            family = abun(), detection = ~ effort,
            y = foray_counts$y, visits = foray_counts$det.covs)
coef(fit)

Description

Occupancy with both false negatives and false positives in the detection process (e.g. acoustic classifiers, citizen-science misidentification), using the Miller et al. (2011) confirmed-detection design that makes the model robustly identifiable. Each visit yields a detection state ⁠y in {0, 1, 2}⁠: 0 = no detection, 1 = ambiguous detection (a true detection OR a false positive), 2 = certain / confirmed detection (only possible when the site is truly occupied). The latent occupancy z marginalises in closed form (two states), so the fit maximises the exact marginal likelihood directly (analytic-gradient BFGS, observed-information covariance) with a NUTS path over the same marginal.

Usage

fp_occu()

Details

Four site-level logit arms: occupancy psi (the tobs() formula), true detection p11 (detection), false-positive rate p10, and the probability a true detection is certain b. The p10 and b predictors default to intercept-only and are set with the p10 = ~ ... and certainty = ~ ... arguments to tobs() (certainty is the b arm). The response y is an ⁠n_sites x J⁠ integer matrix in ⁠{0, 1, 2}⁠ (NA visits dropped).

Value

A tobs_family object.

Engines

• method = "laplace" – analytic-gradient BFGS over the exact marginal with an observed-information covariance; a grouped random effect on the psi or p11 arm is integrated by adaptive Gauss-Hermite quadrature.

• method = "nested_laplace" – an areal icar() / car_proper() field on the occupancy arm (BFGS over the marginal plus the CAR prior), optionally composed with a temporal() field, through the shared areal-BFGS driver.

• method = "nuts" – the in-tree C++ NUTS sampler over the same marginal, with a fixed-hyper areal or temporal field block on the occupancy arm (one of the two, not both).

References

Miller, D. A. W., Nichols, J. D., McClintock, B. T., Grant, E. H. C., Bailey, L. L., Weir, L. A. (2011). Improving occupancy estimation when two types of observational error occur. Ecology 92, 1422-1428. Royle, J. A., Link, W. A. (2006). Generalized site occupancy models allowing for false positive and false negative errors. Ecology 87, 835-841.

Examples

f <- fp_occu()
f

Description

The unmarked gdistremoval model (Amundson et al. 2014): a single-season point-count design in which the detected individuals are recorded two ways at once – their distance band (distance sampling) and the removal period of first detection (removal sampling). This is NOT the open-population distance model (distsampOpen); the abundance is static.

Usage

gdistremoval(transect = c("line", "point"), cutpoints = NULL)

Arguments

Details

Site abundance N_i ~ Poisson(lambda_i); the detected birds are cross-classified by a distance band and a removal period. Writing pdist_i for the overall distance detection, prem_i for the overall removal detection, the total detected is a binomial thinning of N_i, and Poisson is closed under binomial thinning, so the marginal is closed-form:

$\mathrm{ysum}_i \sim \mathrm{Poisson}(\lambda_i\, p^{dist}_i\, p^{rem}_i),$

with the distance-band counts and removal-period counts as two conditional multinomials – the double_observer() Poisson-multinomial pattern, here with a distance multinomial (half-normal key band integrals) and a depleting-removal multinomial pi_k = r (1 - r)^{k-1}.

Three site-level arms: log abundance lambda (the tobs() formula), log distance scale sigma (detection), and logit per-period removal capture r (removal = ~ ..., default intercept-only).

Value

A tobs_family object.

Inputs

y is an ⁠n_sites x n_bins⁠ integer matrix of per-distance-band counts; y.rem an ⁠n_sites x n_periods⁠ integer matrix of per-removal-period counts. The per-site row totals must match (the same detected birds cross-classified).

Engines

• method = "laplace" – the closed-form marginal maximised (BFGS) with an observed-information covariance; the only route (no field, no sampler).

References

Amundson, C. L., Royle, J. A., Handel, C. M. (2014). A hierarchical model combining distance sampling and time removal to estimate detection probability during avian point counts. The Auk 131, 476-494.

See Also

distance(), removal(), double_observer().

Examples

f <- gdistremoval(transect = "point", cutpoints = c(0, 10, 20, 30, 40))
f

Description

Every family, and every method within a family, glances to one column set: nobs, df, logLik, n_fixed, n_samples, n_divergent, mean_accept and converged, the sampler columns NA on a fit with no draws. The joint nested-Laplace fits add the outer grid placement and, when requested, the outer Pareto-k diagnostic. The joint-coupled families (occu_cover(), cover(), occu() spatial, occu_multiscale_cover()) carry both at the fit top level.

Usage

## S3 method for class 'tobs_fit'
glance(x, ...)

Arguments

Value

The one-row glance data frame. A joint-coupled fit adds two outer-grid placement columns, reported whether or not the grid moved so an inert auto-recenter is visible in a batch summary:

outer_grid_placement

"fixed" (the grid the fit was given) or "auto_recentered" (the engine moved it onto the hyperparameter mode).

outer_grid_recenter_declined

On a "fixed" placement, why the recenter did not apply – for instance a user-pinned axis, or "auto_recenter_disabled". NA when the grid was recentered.

A fit that also requested the Pareto-k diagnostic (control$diagnose.k = TRUE, off by default) adds three more:

pareto_k

The outer importance-sampling $\hat{k}$ for the hyperparameter Gaussian summary; ⁠< 0.7⁠ indicates a reliable summary, NA when the diagnostic did not run or the proposal was degenerate.

pareto_k_is_ess

The importance-sampling effective sample size on the PSIS-smoothed weights (numeric); pareto_k_is_ess / control$k.samples is the relative IS efficiency.

pareto_k_proposal_source

How the importance proposal was built: "mode_hessian" from the Laplace curvature at the hyperparameter mode – curvature-backed, so the $\hat{k}$ stays trustworthy even when a sharp posterior collapses the integration grid to ~1 cell; "grid_moment" from the grid-weighted covariance of the integration nodes (with "moment_matched" its refinement); or "grid_mixture", the local-bump-per-cell mixture matching what the engine samples on a spread grid.

Examples

sim <- simulate_occu(N = 100, J = 3, n.occ.covs = 1, n.det.covs = 1,
                     seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ det_cov1, y = sim$y, method = "laplace",
            control = list(verbose = FALSE, progress = FALSE))
glance(fit)

Description

Multiple data sources informing a shared latent occupancy state, each with its own detection process.

Usage

int_occu()

Details

A continuous-mesh spde() term on the detection formula loads a spatially-structured detection field, broadcast to every source, under method = "laplace". The per-source fields are reported as a named list in fit$spatial_field_det. All sources must share one detection formula; areal terms belong on the state arm.

Value

A tobs_family object.

Engines

• method = "laplace" – EM + Laplace over the multi-source two-state marginal, the default.

• method = "laplace_sla" – the Laplace fit with simplified-Laplace skew-corrected marginals.

• method = "laplace_gibbs" / method = "laplace_mi" – the Laplace fit with a post-EM Gibbs / multiple-imputation variance correction.

• method = "nested_laplace" – an areal field on the state arm through the nested-Laplace EM, its hyperparameters integrated on an outer grid.

• method = "nuts" – the unified C++ NUTS sampler over the same marginal.

Examples

f <- int_occu()
f

Description

Multivariate occurrence with shared latent factors. No observation replication – treats observed presence/absence as the response.

Usage

jsdm()

Value

A tobs_family object.

Engines

• method = "laplace" – the community Laplace-EM shared with ms_count() (a logit-link community GLMM); latent() factors alone (lfJSDM) run here through the block-coordinate driver.

• method = "nuts" – samples the exact joint community posterior (community means, per-species deviations, community covariance) over the Bernoulli response, warm-started at the Laplace-EM mode; non-spatial, an areal or latent() term errors with a pointer.

• method = "pg_gibbs" – a per-species Polya-Gamma conjugate Gibbs with community mean and Inverse-Gamma variance updates, giving a calibrated community-variance posterior; non-spatial.

• method = "nested_laplace" – a shared areal field (icar() / car_proper() / bym2()), with or without latent() factors (sfJSDM), by block coordinate ascent; the areal term is required on this route.

Examples

f <- jsdm()
f

Description

The value and degrees of freedom are the engine's (the tulpa_fit method); the "nobs" attribute is nobs.tobs_fit() on the same fit, so stats::BIC() and anything else reading the attribute counts the observations nobs() reports.

Usage

## S3 method for class 'tobs_fit'
logLik(object, ...)

Arguments

Value

A logLik object.

Examples

sim <- simulate_occu(N = 100, J = 3, n.occ.covs = 1, n.det.covs = 1,
                     seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ det_cov1, y = sim$y, method = "laplace",
            control = list(verbose = FALSE, progress = FALSE))
logLik(fit)

Description

A vegetation cover-hurdle dataset for 150 grassland plots surveyed across three years, drawn from a known generative model: the species occurs with a probability driven by moisture, and where it occurs, its conditional cover (a Beta-distributed proportion in ⁠(0, 1)⁠) increases with moisture and decreases with grazing. Roughly half the plots record zero cover. year_c is the centred survey year, so within_between() can split the cross-plot baseline from the within-plot temporal deviation. Synthetic, not field data; the generative coefficients are attached as attr(meadow_cover, "truth").

Usage

meadow_cover

Format

A data frame with 150 rows (one per plot) and columns:

plot

Plot identifier (factor).

year

Survey year (2018-2020).

year_c

Survey year centred on its mean.

moisture

Standardised soil-moisture index.

grazing

Standardised grazing-pressure index.

cover

Recorded cover as a proportion in ⁠[0, 1]⁠ (0 where absent).

See Also

cover(), tobs(), within_between(), simulate_cover()

Examples

data(meadow_cover)
mean(meadow_cover$cover > 0)

fit <- tobs(cover ~ moisture, data = meadow_cover,
            family = cover(response = "beta"))
coef(fit)

Description

Per-species N-mixture with shared community-level hyperparameters.

Usage

ms_abun(K.max = NULL, mixture = c("poisson", "negbin", "zip", "zinb"))

Arguments

Value

A tobs_family object.

Engines

• method = "laplace" – the in-tree C++ community Laplace-EM (per-species coefficient random effects, closed-form community covariance M-step, AGHQ variance debias), the default; latent() factors run here through the block-coordinate driver.

• method = "nested_laplace" – a shared areal field (icar() / bym2() / car_proper()) on the abundance arm (the spAbundance sfMsNMix model), Poisson or grid-integrated negative-binomial size; with latent() factors beside the field the spatial-factor model by block coordinate ascent.

• method = "nuts" – the non-spatial community sampler over the closed-form per-(species, site) marginal, sampling the community means, per-species deviations and community covariances jointly (Poisson, or per-species negative-binomial log_r_s), warm-started at the Laplace-EM mode; with a shared car_proper() field (Poisson) a fixed-hyper non-centered field joins the sampler, while icar() / bym2(), temporal() and re() terms stay on "nested_laplace".

Reading the negbin community dispersion

Under mixture = "negbin" the community mean mu_log_r reaches coef() / vcov() / confint() with a marginal Wald SE, and the variance component sigma_log_r is reported on fit$ms_dispersion. That interval is calibrated conditional on the variance component being recovered, and the two are not independent: sigma_log_r is a scalar variance over species and at few species it can settle near its lower boundary, which shifts mu_log_r and narrows its SE at the same time.

Measured on simulated data (39 laplace fits at 8 and 36 species, sigma.logr = 0.5): where sigma_log_r came back at least 0.30, the nominal 95% interval covered 33 of 34 with a sqrt(mean(z^2)) of 0.88; where it came back below, it covered 2 of 5, the point estimate was 2.2x further from the truth and the SE 28% narrower. So read fit$ms_dispersion$sigma_log_r before quoting a mu_log_r interval: a value near zero, well under what the data should support, means the interval is not trustworthy even though the fit converged and reported no warning. Those thresholds are the ones this fixture separated on and are not a general rule; what transfers is the conditioning, not the number.

Away from that boundary the interval is calibrated, and calibrated at every group count measured. Over 97 fits at 8 / 18 / 36 species the scale sd(mu_log_r) / mean(se) reads 1.066 / 1.265 / 0.981, and the middle figure is the seed block rather than the estimator: mu_log_r is a POPULATION mean, each seed draws S log-dispersions around it, and that draw supplies about two thirds of the across-seed spread. The 18-species blocks drew theirs 18-21% wider than sigma.logr / sqrt(S). Put the draw at its expectation and the scale is 1.077 / 1.101 / 0.977; rebuild the SE at the simulated sigma as well and it is 0.990 / 1.035 / 0.963 (NOTES_measurements.md). The residual there is the attenuation of sigma_log_r itself (0.418 / 0.448 / 0.487 against 0.5), which shrinks as species are added and which the SE inherits, since se^2 = (sigma_log_r^2 + c) / S.

A simulation study measuring this needs simulate_ms_abun()'s truth$mu_log_r_real – the mean of the log-dispersions that seed actually drew – beside truth$mu_log_r. Score coverage of a community mean against the constant; score point recovery and interval SCALE against the realized mean, or the draw is charged to the estimator.

control$logr.sigma.prior puts a Penalized-Complexity prior on that variance, which adds curvature at sigma_log_r -> 0. Measured on the same fixture at c(1, 0.05), it is worth reaching for when a fit reports a sigma_log_r near zero or fails outright with a singular marginal Hessian – one seed that errored under pure maximum likelihood converged and covered under the prior, and variance components near 0.01-0.06 lifted by an order of magnitude. It does not repair coverage: on 19 paired seeds the nominal 95% interval covered 17 either way, with the same two seeds missing, because the prior bites near the boundary while those misses sit at sigma_log_r around 0.2. It is off by default because the fits that were already calibrated pick up a small systematic shift under it (mu_log_r by -0.006, p = 0.001).

Examples

f <- ms_abun(mixture = "negbin")
f

Description

The community analogue of count(): a per-species GLMM on an observed count or continuous response with Gaussian community hyperpriors on the per-species coefficients (the spAbundance msAbund model). No detection and no latent state – the abundance counterpart of ms_occu() without the occupancy layer. Poisson / negative-binomial (log link) or Gaussian (identity). y is an ⁠n_sites x n_species⁠ matrix (or a named list of n_species count vectors), one value per site per species; species names the columns.

Usage

ms_count(response = c("poisson", "negbin", "gaussian", "binomial"))

Arguments

Details

As with the other community Laplace-EM families (ms_occu(), ms_occu_cover()), the community means are recovered essentially unbiased for the Gaussian (identity link, exact Laplace) and Poisson responses; the negative-binomial slope carries a mild first-order-Laplace (PQL) attenuation of order the community variance (a few percent, shrinking with the number of species and observations per species). The community-mean Wald intervals are calibrated to the package rubric (pooled coverage at least 0.85).

method = "nuts" samples the exact joint posterior (community means, per-species deviations, and the community covariance) for all three responses – the negative binomial carrying a per-species dispersion random effect, the Gaussian a per-species free residual variance – which removes the Laplace-EM's negative-binomial attenuation and returns calibrated, non-Gaussian community intervals. A shared areal field or latent factors are available through method = "nested_laplace" / the latent() term; see the package overview for the spatial and factor variants.

Value

A tobs_family object.

Engines

• method = "laplace" – the shared community Laplace-EM, the default; latent() factors alone run here through the block-coordinate driver.

• method = "nested_laplace" – a shared areal field icar() on the abundance formula (the spAbundance sfMsAbund model), with or without latent() factors, by block coordinate ascent; the field is required on this route.

• method = "nuts" – the exact joint community posterior sampled by the in-tree C++ target, warm-started at the Laplace-EM mode, as described above; non-spatial.

• method = "pg_gibbs" – a per-species Polya-Gamma conjugate Gibbs with community mean and Inverse-Gamma variance updates, for the logistic "binomial" response only (another response errors with a pointer); non-spatial.

See Also

count() (single species), ms_occu(), ms_abun()

Examples

sim <- simulate_ms_count(N = 120, n.species = 8, seed = 1)
fit <- tobs(~ x, data = sim$data, family = ms_count(), y = sim$y,
            species = colnames(sim$y), method = "laplace")
summary(fit)

Description

The community analogue of distance(): per-species binned distance sampling with Gaussian community hyperpriors on the per-species abundance and detection-scale coefficients, so rare species borrow strength from common ones through the shared community means and covariances – the same pooling that stabilises ms_abun(), on the distance-sampling response.

Usage

ms_distance(
  key = c("halfnorm", "hazard"),
  transect = c("line", "point"),
  cutpoints = NULL,
  K.max = NULL,
  mixture = c("poisson", "negbin")
)

Arguments

Details

Per species s, site i, distance bin b:

N_{s,i}        ~ Poisson(lambda_{s,i})
y_{s,i,.} | N  ~ Multinomial(N_{s,i}; pi_{s,i,1..B}, 1 - p_det)
log lambda_{s,i} = X_i     . (mu_lambda + b_lambda_s)
log sigma_{s,i}  = X_sig_i . (mu_sigma  + b_sigma_s)
b_lambda_s ~ N(0, Sigma_lambda),  b_sigma_s ~ N(0, Sigma_sigma)

with pi_b the integral of the detection function over bin b. The latent ⁠N_{s,i}⁠ integrates out per species-site in closed form, exactly as for distance(); the per-species coefficient deviations are the random effects, integrated by the shared community Laplace-EM.

Adding latent(n) to the abundance formula gives residual species co-occurrence through n per-site latent factors with per-species loadings; adding an areal term (icar() / car_proper() / bym2()) or spde() alongside gives a shared spatial field. Both route through the same block-coordinate engine as the other community families.

Value

A tobs_family object.

Engines

• method = "laplace" – the shared community Laplace-EM over the existing distance kernel, plain or with latent() factors alone (lfMsDS).

• method = "nested_laplace" – a shared field (icar() / car_proper() / bym2() / spde()) on the abundance formula, with or without factors (sfMsDS), by block coordinate ascent; the field is required on this route.

Poisson only, and no sampler route.

Inputs

y is a 3D array ⁠[n_sites x n_bins x n_species]⁠ or a named list of ⁠n_sites x n_bins⁠ per-bin count matrices. formula is the abundance (⁠log lambda⁠) predictor, detection the detection-scale (⁠log sigma⁠) predictor, and species the species labels.

See Also

distance() (single species), ms_abun() (community N-mixture).

Examples

f <- ms_distance(key = "halfnorm", transect = "line",
                 cutpoints = c(0, 25, 50, 75, 100))
f

Description

Per-species dynamic (HMM) occupancy with Gaussian community hyperpriors on the per-species first-season occupancy and detection coefficients, and shared community-wide colonisation / extinction transition coefficients.

Usage

ms_dyn_occu()

Value

A tobs_family object.

Engines

• method = "laplace" – the shared community Laplace-EM, the default: the shared colonisation / extinction dynamics and the per-species first-season occupancy / detection components recover across seeds (see tests/testthat/test-ms-dyn-occu.R, community-mean 95% CI coverage measured ~0.98).

• method = "pg_gibbs" – the community Polya-Gamma machinery of msPGOcc with a two-state HMM forward-filtering backward-sampling latent step (spOccupancy tMsPGOcc), giving a calibrated community-variance posterior; constant transitions, site-level detection, no structured term.

• method = "nested_laplace" – a shared areal field on the first-season occupancy formula (~ 1 + icar(graph = adj)) fits the spOccupancy stMsPGOcc model: the field is shared across species, and because the first-season occupancy psi1 only sets the initial mixing weight of each species' HMM, the block-coordinate driver alternates the community EM (field as a psi1 offset) with an areal field Newton, the field recovering cleanly (cor ~0.94). A spatially-varying-coefficient bar (~ spatial(~ 1 + w || cell, graph = adj)) adds a shared covariate-weighted field alongside the intercept field, fitting the svcTMsPGOcc model through the same K-field weighted-ICAR solve as the community count SVC (both fields recover, cor ~0.90 / ~0.89). icar() only; bym2() / car_proper() fields are a follow-up. The field is required on this route.

• method = "nuts" – the non-spatial community sampler over the exact per-species HMM-forward marginal (in-tree C++), sampling the community means, per-species first-season / detection deviations, both per-arm community covariances and the shared colonisation / extinction coefficients jointly, warm-started at the Laplace-EM mode; an areal term errors with a pointer to "nested_laplace".

See Also

dyn_occu(), ms_occu()

Examples

f <- ms_dyn_occu()
f

Description

Multiple detection data sources share one latent occupancy state per species, with Gaussian community hyperpriors on the per-species coefficients of the occupancy arm and of every per-source detection arm. Fit by a shared community Laplace-EM.

Usage

ms_int_occu()

Value

A tobs_family object.

Engines

• method = "laplace" – the shared community Laplace-EM, the default: the shared occupancy mean and the per-source detection community components recover near nominal across seeds and more than one source (see tests/testthat/test-ms-int-occu.R). The binary community-mean intervals carry the mild Laplace under-dispersion typical of occupancy data (measured 95% CI coverage ~0.89 at small per-species n).

• method = "pg_gibbs" – the community integrated Polya-Gamma Gibbs sampler (one latent state per species-site observed by every source), giving a calibrated community-variance posterior; site-level per-source detection.

• method = "nuts" – the community integrated sampler over the exact multi-source two-state per-(species, site) marginal (in-tree C++), sampling the community means, per-species deviations and the D + 1 per-arm community covariances jointly, warm-started at the Laplace-EM mode.

No route takes a structured term: a field or random effect on either formula errors with a pointer, and there is no nested-Laplace path for this family.

See Also

int_occu(), ms_occu()

Examples

f <- ms_int_occu()
f

Description

Per-species occupancy and detection with Gaussian community hyperpriors on the per-species coefficients of each arm (the spOccupancy msPGOcc model). The occupancy and detection per-species deviations are independent, each with its own community covariance, fit by a community Laplace-EM.

Usage

ms_occu()

Value

A tobs_family object.

Engines

• method = "laplace" – the shared community Laplace-EM, the default; latent() factors alone (lfMsPGOcc) run here through the block-coordinate driver.

• method = "nuts" – the non-spatial community sampler over the closed-form per-(species, site) two-state marginal, sampling the community means, per-species deviations and both per-arm community covariances jointly, warm-started at the Laplace-EM mode; an areal term errors with a pointer to "nested_laplace".

• method = "pg_gibbs" – the hierarchical Polya-Gamma Gibbs sampler (spOccupancy msPGOcc), giving a calibrated community-variance posterior; non-spatial.

• method = "nested_laplace" – a shared areal field (icar() / bym2() / car_proper()) on the occupancy arm: a plain intercept field runs the in-tree C++ community-spatial nested Laplace-EM, a varying-coefficient bar (svcMsPGOcc) or latent() factors beside the field (sfMsPGOcc) run the block-coordinate fitter; the field is required on this route.

See Also

occu(), ms_dyn_occu(), ms_int_occu()

Examples

sim <- simulate_ms_occu(N = 80, J = 4, n.species = 8, seed = 1)
fit <- tobs(~ x, data = sim$data, family = ms_occu(), detection = ~ 1,
            y = sim$y, species = paste0("sp", seq_len(8)), method = "laplace")
summary(fit)

Description

The community version of occu_cover(): a per-species joint occupancy-cover model with Gaussian community hyperpriors on the per-species coefficients of all three arms (occupancy psi, detection p, and positive cover). Rare species borrow strength from common ones through the shared community means and covariances, the same pooling that stabilises ms_occu() and ms_abun() but on the joint occupancy + vegetation-cover response.

Usage

ms_occu_cover(response = c("beta", "lognormal", "gaussian"))

Arguments

Details

Per species s, cell i, visit j:

z_{s,i}        ~ Bernoulli(psi_{s,i})
y_{s,i,j} | z  ~ Bernoulli(p_{s,i,j})
c_{s,i,j} | y  ~ f_pos(eta_pos_{s,i,j}, disp)
logit psi_{s,i}  = X_occ_i . (mu_occ + b_occ_s)
logit p_{s,i,j}  = X_p_{ij} . (mu_p   + b_p_s)
g(cover)         = X_pos_{ij} . (mu_pos + b_pos_s)
b_occ_s ~ N(0, Sigma_occ), b_p_s ~ N(0, Sigma_p),
b_pos_s ~ N(0, Sigma_pos)

The latent presence z integrates out per species-cell in closed form (the same two-state mixture as occu_cover()); the per-species coefficient deviations are the random effects, integrated by a Laplace-EM. The positive-arm dispersion is a shared community parameter.

Value

A tobs_family object.

Inputs

y and y.pos are 3D arrays ⁠[n_sites x max_visits x n_species]⁠ (or named lists of ⁠n_sites x max_visits⁠ matrices, one per species); ⁠species =⁠ is required. The occupancy formula, the detection formula, and the cover positive formula carry community covariates shared across species. coef() returns the community means; ranef() the per-species coefficient deviations.

Engines

• method = "laplace" – the community Laplace-EM (per-species deviations integrated, closed-form community covariance M-step), non-spatial or with the reduced-rank spatial factors below.

• method = "nuts" – the community sampler over the exact per-(species, cell) two-state marginal, warm-started at the Laplace-EM mode: non-spatial it samples the community means, per-species deviations, the three per-arm community covariances and the shared log-dispersion jointly; with a shared areal field it samples the spatial-factor model below.

Reduced-rank spatial factors

A single areal field term (icar(graph = adj), car_proper(graph = adj), or bym2(graph = adj)) on the occupancy formula fits the reduced-rank (HMSC / spatial-gllvm) community model: K shared latent fields with per-species loadings on the occupancy predictor, ⁠logit psi_sc = X_c (mu_occ + b_occ_s) + sum_k L_sk w_kc⁠. A single shared field with only a species intercept can shift each map's level but not its shape; the loadings give each species its own spatial shape as a combination of the shared factors, borrowing strength so rare taxa get a calibrated map. The number of factors is set by control = list(n.factors = K), or chosen automatically by control = list(n.factors = "auto", n.factors.max = M), which fits the identified (lower-triangular loading) ladder and picks the rank that maximises the empirical-Bayes Laplace marginal likelihood (the field is integrated out, so its prior supplies the Occam penalty – latent-level criteria such as WAIC track the field's effective dimension, not the rank, and under-select). The per-rank evidence is returned in fit$spatial$K_selection. sd.load sets the loading prior scale.

Adding the SAME icar(graph = adj) term to the cover positive formula shares the latent fields across the two processes: the fields then also load on the cover predictor through a free loading matrix, ⁠g(cover)_sc = X_c (mu_pos + b_pos_s) + sum_k Lpos_sk w_kc⁠, so a species' spatial occupancy pattern and its spatial cover pattern are linked through one set of factors. The cover loadings are returned in fit$spatial$loadings_cover. The field is shared, so the cover-arm term must match the occupancy arm (same type and graph), and a cover-arm field without a matching occupancy field errors.

The field structure is set by the term: icar() (improper intrinsic CAR, the default), car_proper(graph = adj) (a proper CAR whose per-factor correlation rho is returned in fit$spatial$rho_w), or bym2(graph = adj) (the Riebler 2016 convolution whose per-factor spatial-variance fraction phi is returned in fit$spatial$phi_w); the field hyperparameter is estimated by EM. All three share one engine; fit$spatial$field_type records the choice. Structured terms on the detection arm, or an unsupported field term, error from the dispatcher.

The shared fields imply a residual species-association matrix (the spatial-JSDM / HMSC output): tobs_associations() returns the occupancy association ⁠corr(L L')⁠ and, with a cover-arm factor, the cover association and the cross-arm occupancy-vs-cover association, each marginalised over the loading posterior so an interval accompanies the estimate. predict() returns the other JSDM output – the per-species per-cell occupancy maps (psi with psi_lower / psi_upper), marginalised over the loading + field posterior so a rare species borrows strength across the shared factors for a calibrated map.

Scope

The non-spatial fit is Laplace-EM. A community spatial occu_cover – a shared latent field coupled across the occupancy and cover arms with per-species RE on all three arms – is the reduced-rank spatial-factor fit: an icar() / car_proper() / bym2() term on the occupancy formula (and, with a matching term on the cover formula, a cover-arm factor) fits per-species loadings on the shared field via Laplace-EM and NUTS, with the per-species community covariances Sigma_occ / Sigma_p / Sigma_pos on top. The free per-species loading form generalises a single common field amplitude, so it subsumes the common-amplitude coupling of occu_cover()'s joint-coupled engine. A community model whose per-arm variance components are integrated on an outer grid (the engine route of occu_cover()'s joint-coupled path) is not used: the joint nested-Laplace engine integrates every variance component on its outer grid, so per-arm community RE variances plus the field hyperparameters exceed any tractable grid size – the closed-form covariance M-step of the Laplace-EM is the scaling route for community variance components. Structured terms beyond the shared field error from the dispatcher rather than being silently dropped.

Community variance debias

The community-MEAN estimates (coef(), vcov(), confint()) are unbiased. The community-VARIANCE components – the per-arm covariance matrices in fit$ms_community$Sigma_occ / Sigma_p / Sigma_pos and their ⁠sd_*⁠, which set the spread of the per-species deviations – pick up Laplace small-cluster attenuation at small per-species n under the EM M-step. They are debiased by default with an adaptive Gauss-Hermite quadrature of the exact per-species RE posterior (control$re.aghq = TRUE, the same correction the single-arm AGHQ path applies; disable with re.aghq = FALSE, set nodes with control$n.quad). The cover hurdle ties a species' psi / p / cover coefficients through the data, so the per-species RE posterior is not separable across arms and the AGHQ is a tensor product over the joint RE vector – n.quad^P nodes per species in the total RE dimension P. That cost is exponential in P, so the debias is a hard scope limit: it runs only up to a small total RE dimension (control$re.aghq.maxdim, default 4, e.g. intercept-only on all three arms), and larger RE designs keep the EM variance. The EM (Laplace) variance is a documented lower bound on the true component – attenuated toward zero at small per-species n, monotonically less attenuated as n grows – not a bias of unknown sign. fit$ms_community$var_attenuation$debias records "aghq" or "none", ⁠$affects⁠ names the lower-bounded components, and print() flags the EM case.

See Also

occu_cover() (single species), ms_occu() (community occupancy, no cover), ms_abun() (community N-mixture).

Examples

f <- ms_occu_cover(response = "lognormal")
f

Description

The count of observed response entries the fit's likelihood reads, which is the sample size stats::AIC() and stats::BIC() form an information criterion at. What one entry is follows the family's own response layout: a surveyed (site, visit) cell for the occupancy families, a (site, band) count for distance sampling, a (site, visit, source) cell for the integrated ones.

Usage

## S3 method for class 'tobs_fit'
nobs(object, ...)

Arguments

Value

Integer count of the fit's observed response entries. A model type with no registered count is an error naming that type.

Examples

sim <- simulate_occu(N = 100, J = 3, n.occ.covs = 1, n.det.covs = 1,
                     seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ det_cov1, y = sim$y, method = "laplace",
            control = list(verbose = FALSE, progress = FALSE))
nobs(fit)

Description

Latent Bernoulli occupancy state with binomial detection per visit. The MacKenzie et al. (2002) single-season model.

Usage

occu()

Value

A tobs_family object.

Engines

• method = "laplace" – penalized EM + Laplace with Gaussian marginals, the default.

• method = "laplace_sla" – the Laplace fit with simplified-Laplace skew-corrected marginals.

• method = "laplace_gibbs" / method = "laplace_mi" – the Laplace fit with a post-EM Gibbs / multiple-imputation variance correction (Rubin-pooled refits, not an MCMC chain).

• method = "pg_gibbs" – a Polya-Gamma Gibbs sampler over the exact posterior (the spOccupancy PGOcc engine): site-level detection, an areal icar() field is accepted; a temporal() / re() / latent() / svc() term errors with a pointer.

• method = "nested_laplace" – multi-block nested Laplace: areal, temporal and random-effect hyperparameters are integrated on an outer grid, and predict() interpolates all-NA detection sites from the field.

• method = "nested_laplace_sla" – the nested route with skew-corrected marginals.

• method = "nuts" – the unified C++ NUTS sampler over the structured terms (spatial, temporal, re(), svc(), latent()), reporting Rhat / ESS.

Examples

f <- occu()
f

Description

Scans candidate spatial cell sizes and yearly clusterings of long-format, typically single-visit plot data and scores each by how well a single-season occupancy model separates occupancy (psi) from detection (p). Replication is manufactured by pooling plots into a spatial cell (assumes plots in a cell share psi) and years into a contiguous block (assumes occupancy is closed across the block); an occasion is any record in a (cell, block) bucket, so single-visit plots, repeat-visited plots, and the no-pooling limit all flow through one code path.

Usage

occu_aggregation_scan(
  data,
  response,
  coords,
  year,
  plot = NULL,
  cell.sizes = NULL,
  block.lengths = NULL,
  score = c("info", "count"),
  family = c("occupancy"),
  control = list()
)

Arguments

Details

Two scoring modes:

• "info" fits the constant (intercept-only) occupancy model per candidate and reports the posterior SEs of psi/p and the smallest eigenvalue and condition number of the 2x2 ⁠(logit psi, logit p)⁠ observed information. A confounding ridge shows up directly as a near-zero eigenvalue.

• "count" checks only the structural necessary conditions (units with two or more occasions, detected units, a non-degenerate naive detection rate).

The recommended candidate is the least-pooling identifiable one (smallest mean occasions per unit, ties broken by largest information eigenvalue): pool only as much as identifiability requires. When no candidate is identifiable the data cannot support a standard occupancy model and a Royle-Nichols or count model is the alternative.

Value

A tobs_aggregation_scan object: candidates (scored data.frame), recommended (chosen row or NULL), segmentation (auto changepoint blocks, if used), and scan metadata.

Examples

set.seed(1)
d <- data.frame(east = runif(200), north = runif(200),
                year = sample(2001:2006, 200, replace = TRUE),
                detected = rbinom(200, 1, 0.3))
scan <- occu_aggregation_scan(d, response = "detected",
                              coords = c("east", "north"), year = "year",
                              cell.sizes = c(0.2, 0.4),
                              block.lengths = c(2, 3))
scan

Description

A hurdle for a response that is either absent or, when present, one of K nominal (unordered) classes – a colour morph, microhabitat / substrate use, or a cryptic-species / classifier label observed given the organism is present. The presence/absence part is a Bernoulli arm; the class given present is a baseline-category multinomial logit (the last class is the baseline). The two arms factorise the likelihood exactly:

$P(y = 0) = 1 - \psi, \qquad P(y = k) = \psi \, p_k,$

with $p = \mathrm{softmax}(X \beta_{class})$.

Usage

occu_categorical(classes = NULL)

Arguments

Details

This is the categorical counterpart of cover() (presence + a magnitude): here the positive part is an unordered class, so it uses a multinomial logit rather than beta / lognormal. For an ordered class response (Braun-Blanquet cover bands) use cover(response = "ordinal"), which exploits the ordering; this family is for classes with no ordering.

Value

A tobs_family object.

Response

y is a length-N integer vector in ⁠0..K⁠: 0 is absent, k is class k. It may sit on the formula left-hand side (y ~ ...), so ⁠y =⁠ can be dropped. The formula predictors are shared by both arms.

Engines

Non-spatial Laplace (method = "laplace"), the only route. The multinomial math is the FD-validated tulpa kernel (multinomial_logit.h); the non-spatial fit is the vectorised R Newton over the same closed forms. Spatial fields / NUTS (the native multi-process likelihood) and the latent-class misclassification variant (the K-class generalisation of fp_occu(), a confusion matrix on the observed label) are documented follow-ups.

See Also

cover() (presence + magnitude), fp_occu() (two-state false-positive detection).

Examples

f <- occu_categorical(classes = c("red", "green", "blue"))
f

Description

Cell-level latent presence psi, per-visit binomial detection p, and per-visit positive cover f_pos (beta or lognormal) on a third linear predictor. The N-mixture analogue for vegetation cover: where abun() couples occupancy to counts via Royle's binomial-N marginal, occu_cover() couples it to a positive-cover observation at each detected visit. The latent presence z marginalises out in closed form (two states), so the marginal log-likelihood is exact.

Usage

occu_cover(
  response = c("beta", "lognormal", "gaussian"),
  cover.aggregate = NULL
)

Arguments

Details

Per-cell likelihood:

any_det_i : L_i = psi_i * prod_j h_ij
no_det_i  : L_i = psi_i * prod_j (1 - p_ij) + (1 - psi_i)

h_ij = (1 - p_ij) * 1{y_ij = 0}
     + p_ij       * f_pos(y_pos_ij; eta_pos_ij, ...) * 1{y_ij = 1}

A detected visit (y_ij = 1) with a missing cover (y_pos_ij = NA) keeps its detection term but drops the f_pos factor: cover is taken missing-at-random, so the cover likelihood runs over the detected visits with an observed cover.

Reduces to occu() when the cover arm is degenerate, and to the plot-level cover hurdle (cover()) when J = 1 and detection is perfect.

Value

A tobs_family object.

Engines

The non-spatial fit is a direct Laplace approximation (method = "laplace") or a NUTS sampler over the same exact two-state marginal (method = "nuts", beta or lognormal cover), giving calibrated intervals and a per-draw pointwise likelihood for WAIC / LOO. A shared areal field across the occupancy and cover arms is the method = "nested_laplace" path (a structured icar() / bym2() term on the psi formula); a structured term on a non-spatial route errors from the dispatcher with a pointer to it.

method = "nuts" also samples the coupled areal field(s) (icar() / bym2() / car_proper()), written either as areal terms or as the bar spatial(~ 1 || cell, graph = adj), TOGETHER WITH their hyperparameters: the field SD, the mixing (bym2) or spatial-correlation (car_proper) parameter, and - where the formula asks for one - the cover-arm copy amplitude. Each is bounded to the span the nested_laplace path's outer quadrature integrates that axis over, in that grid's own coordinate, and the same control$sigma.grid / alpha.grid / rho.car.grid knobs set it on both routes - so the sampler is an independent reference for the hyperparameter layer rather than a fit conditioned on that layer's point estimate. The measure over that span is the one the outer grid declares for the axis: flat in the grid's coordinate for the field SD and the mixing parameter, and for the copy amplitude the penalized-complexity Exponential slab, with control$copy.slab = "flat" asking for the flat alternative the grid also accepts. The grid weighs that slab against a point mass at alpha = 0 ("no coupling"), which a gradient sampler cannot visit, so the sampled copy amplitude is the posterior conditional on a coupled field; the point mass's share of the prior is read on the nested_laplace route. fit$nuts$sampled_hyper and fit$nuts$fixed_hyper name which is which per fit (icar pins rho at 1: the intrinsic precision has no mixing parameter; a field with no share() pins the copy amplitude at 0), fit$hyper_draws carries their posterior alongside field_sd, the geometric-mean marginal SD the field implies. control$fixed.hyper = TRUE conditions on the warm nested-Laplace estimate instead.

The grid-integrated nested_laplace path reports field_sd too (fit$spatial$field_sd_mean / field_sd_sd, alongside sigma_mean), in the SAME geo-mean-marginal-SD convention as the NUTS path's field_sd. It is a fixed multiple of sigma (not an independent grid axis), so it is NOT in fit$means / fit$vcov – folding it in there would make the joint parameter vcov exactly singular. Do not compare fits, or a fit against simulate_occu_cover()'s sigma, by reading the raw sigma on this path: it is the field's amplitude against the unscaled intrinsic ICAR precision Q = D - W, which differs from field_sd by sqrt(scale_q), a graph-size-dependent factor (about 2.1 for a 30-node chain graph). field_sd is the number comparable across fits and to a simulation truth; sigma is the raw amplitude the engine's grid axis is spelled in.

A spatially-varying coefficient - a weighted areal term beside the unweighted intercept field, or a bar column spatial(~ 1 + year || cell, graph = adj) - is sampled as a SECOND field block, with its own whitened surface and its own (SD, mixing, copy amplitude) coordinates: two fields share no hyperparameter. The surfaces are fit$spatial_field and fit$trend_field / fit$trend_fields (named by the weight column), and each block's hypers carry that block's suffix in fit$hyper_draws and in fit$nuts$sampled_hyper (sigma_trend, alpha_trend, indexed when there are several). Each field's amplitude is addressed on its own through share(spatial(), terms = list(...)) – with stated nodes, or grid(n = ) for the engine's axis read more finely – or through control$alpha.grid[.trend] / control$alpha.n[.trend], the same knobs the grid-integrated route reads. A second field multiplies the warm fit's outer grid, so a defaulted axis is thinned to three nodes over the same span - the sampled prior is unchanged, the warm fit stays affordable. A correlated bar (|, one free-Sigma MCAR block across the fields) stays on the grid-integrated nested_laplace path; the community spatial-factor route is ms_occu_cover().

On the shared-field nested_laplace path the field-coupled occupancy slope Wald interval is mildly anti-conservative at small N (the grid-integrated Laplace under-disperses that coefficient; the pooled coverage across the three arms stays near nominal). The non-spatial nuts path gives fully calibrated occupancy intervals.

What a bare areal term means for the cover arm

A term's process is the formula it sits in. An areal term written in the occurrence formula therefore loads on the OCCURRENCE arm alone: the cover arm sees no field, and the copy amplitude alpha is 0. This holds under method = "nested_laplace" and method = "nuts" alike, so the same input fits the same model on either engine.

To carry the occurrence field onto the cover arm as well, write the copy in the positive formula:

positive = ~ pos_cov + share(spatial())                      # amplitude estimated
positive = ~ pos_cov + share(spatial(), alpha = grid(c(...))) # over given nodes
positive = ~ pos_cov + share(spatial(), alpha = grid(n = 9))  # engine axis, finer
positive = ~ pos_cov + share(spatial(), alpha = 0.5)          # fixed
positive = ~ pos_cov + share(spatial(), prior = list("pc.prec", c(4, 0.01)))

The amplitude is integrated over that axis on the nested_laplace path and sampled over its span on the nuts path; a scalar ⁠alpha =⁠ pins it on both. ⁠prior =⁠ regularizes the coefficient itself, in the joint driver's ⁠list(<family>, <params>)⁠ shape (see prior_alpha in tulpa::tulpa_nested_laplace_joint()); one reaches the engine per fit, so a fit copying both an intercept and a trend block is refused it rather than given it on the first.

control$alpha.grid[.trend] and control$alpha.n[.trend] are the WIRE FORMAT these compile into, not user surface: a fit that sets one is refused and told which formula form to write. control$prior.alpha is still the lower-level spelling of ⁠prior =⁠; set it in one place, not both.

Without a share() the occurrence field is NOT carried onto the cover arm – there is no implicit coupling – so the amplitude is not a parameter of the model, and neither engine reports one. Under nested_laplace alpha is absent from coef(), vcov() and the n_params count; under nuts it is absent from fit$nuts$sampled_hyper, fit$nuts$fixed_hyper and the fit$hyper_draws columns. predict() returns a cover arm the occurrence field does not enter. Add share(spatial()) to couple the two arms.

The axis the amplitude rides carries prior structure – a point mass at alpha = 0 ("no coupling") and a log-spaced slab above it – so it is set in one of two ways. share(alpha = grid(...)) STATES its nodes, and with them that structure. share(alpha = grid(n = 9)) states a RESOLUTION: the engine re-reads its own axis with that many slab nodes, point mass and slab bounds unchanged, so sharpening the axis never restates its structure. ⁠terms =⁠ gives either, per block. The resolution is the only way to raise this axis, because it does not densify when the donor control$sigma.grid does: on an informative data set the outer grid's quadrature effective sample size saturates on the copy amplitude while every other axis tracks the request (measured engine-side, NOTES_measurements.md). One block takes one of the two; giving both is an error. control$alpha.grid[.trend] and control$alpha.n[.trend] are the wire format these compile into and are refused as user input.

Set it when the fit's hyperparameter intervals are reported. On the coupled SBC fixture the declared resolution leaves the field SD miscalibrated once the data are informative – at 10 visits per site its uniformity p-value is 9.1e-05, against 0.17 at 3 visits – and control$alpha.n = 21 returns every scored quantity to nominal at both. Thirteen slab nodes do not (the field SD still reads 4.5e-03), and the price is a 2.3-3x longer fit, since the outer grid is a tensor. Point estimates are not affected the way the intervals are; the measurement is in NOTES_measurements.md.

Coupled fields and spatially-varying trends

The spatial engines (method = "nested_laplace", default control$engine = "joint", and method = "nuts") share one areal field (the cell intercept, an unweighted icar() / bym2() term) across the occupancy and cover arms, when the positive formula copies it. ADDITIONAL coupled fields - spatially-varying coefficients, e.g. a temporal trend - are added as weighted areal terms in the formula:

~ elev + icar(graph = adj) + icar(graph = adj, weight = year)

Each weighted term icar(graph, weight = col) is a second field whose contribution to a predictor row is weight_i * z[cell_i]. All such fields share the same graph; a share() carries one onto the cover arm with its own scale (alpha for the intercept field, alpha_trend for a trend field), integrated over the outer grid, and a field the positive formula does not copy stays on occurrence. share(spatial(), terms = list(...)) gives a per-block amplitude and must address every block. The intercept field is reported in fit$spatial_field, trend fields in fit$trend_field / fit$trend_fields. The trend coupling grid defaults to the intercept block's amplitude; give the trend its own with share(spatial(), terms = list(intercept = ..., trend = ...)). Under method = "nuts" each field is its own sampled block, with the scales sampled over those same spans rather than integrated over the grid.

A single trend field may also be requested out-of-band with control = list(trend = list(weight = "<col>")), naming a numeric per-cell covariate in the cell data; this is the equivalent of one weighted formula term. Specify the trend field one way or the other, not both.

Sites larger than cells (group.var)

By default each site (one row of y / data, one latent occupancy state) is its own field node, so the graph must have one node per site. Passing group.var = "<col>" to the icar() / bym2() term maps each site to a field node named by that integer column, so several sites can share one node. The field then stays length n_cells (the graph) while occupancy, detection, and cover run over n_sites. The motivating layout is a site = cell x time-period: plots in a cell-period are the detection replicates, occupancy is per cell-period, and a per-site time weight on a coupled trend field (icar(graph, weight = time, group.var = "cell")) gives a detection-corrected occupancy trend on a shared cell field.

Per-group random intercept on the shared-field path (re() / (1 | g))

On the spatial nested_laplace path a single random INTERCEPT term on the psi formula – re(g) or the (1 | g) bar – layers a per-group occupancy offset on top of the shared field. It joins the joint fit as one iid prior block whose latent rides the occupancy arm only; its variance integrates on the outer grid alongside the field sigma / alpha and is reported as the sigma_re hyperparameter, with the per-group BLUPs in fit$re and via ranef(). A default re.sigma.grid (log-spaced) sets the grid; pass control$re.sigma.grid to override. Scope: one random-intercept term (a scalar per group maps onto the one iid block); a random slope, a correlated multi-coefficient block, or an RE without a shared field is rejected (the joint engine integrates every variance component on its grid, so multiple variance components do not scale – richer RE is the non-spatial cover-hurdle EM's route). This is the single-species analogue of the community spatial occu_cover (ms_occu_cover()) and the tulpaObs consumer of tulpa's field + per-group RE engine composition.

A bar is therefore always a random effect, never a spatial field, even when the grouping factor names the areal graph nodes. The engine's inline-MCAR call tulpa::spatial(graph, ~ 1 + x | cell) reads ⁠| cell⁠ as a separable spatial field, but the same spelling in an occu_cover() formula is an IID random effect on cell. For a spatial field on the cells use an areal term – icar(graph = adj, group.var = "cell") (plus icar(graph = adj, weight = x, group.var = "cell") for a spatially-varying trend) – not a bar. When a formula carries a bar whose grouping factor is also an areal term's group.var, occu_cover() emits a one-time message noting the bar is fitted as a random effect; suppress it with base::suppressMessages().

Independent field on the cover arm (placement)

A share(spatial()) in the positive formula shares the occupancy field: the cover arm sees it as ⁠alpha * (occupancy field)⁠, the coregionalization copy on the outer ⁠(sigma, alpha)⁠ grid. When the cover trend is spatially structured but is not a scalar multiple of the occupancy field, alpha collapses toward 0 and the cover arm inherits no field, so per-cell conditional cover (and its change over time, delta_cover_cond) comes out flat. A spatial() bar placed in the positive formula adds an INDEPENDENT, non-copied areal field on the cover arm alone:

tobs(occurrence = ~ x + icar(graph = adj, group.var = "cell"),
     detection  = ~ 1,
     positive   = ~ time + spatial(~ 1 + time || cell, graph = adj),
     data = cell_dat, y = y, y.pos = y_pos,
     family = occu_cover("lognormal"), method = "nested_laplace")

The occupancy field still drives psi and, via the alpha copy, delta_cover_exp; the independent cover field carries a cover-specific structure the alpha copy cannot express. Its intercept + coefficient columns become separate ICAR blocks on the cover arm, each with its own precision integrated on the outer grid and reported as sigma_pos_field (intercept) / ⁠sigma_pos_field_<col>⁠ (a covariate column). The field SD grid defaults to control$sigma.grid; set control$sigma.grid.pos.field to integrate it over its own (usually coarser) grid, which keeps the added axis from multiplying the outer-grid cost. The per-cell field posterior is in fit$pos_field / fit$pos_field_table (and fit$pos_field_tables per column). Scope: it composes with the shared occupancy field but uses per-visit cover (cover.aggregate = "none") and does not combine with the correlated | MCAR field, the latent cover RE, or the batched fused path; like the shared field it is fitted as ICAR (bym2/car read as ICAR). A static intercept cover field is only weakly identified against the occupancy field's alpha copy (both are per-cell cover offsets); a time-weighted trend cover field is identified separately and is what makes delta_cover_cond spatially varying.

The same placement works on the detection arm: a spatial-field term in the detection formula (detection = ~ 1 + spatial(~ 0 + time || cell, graph = adj)) fits an independent field on the detection predictor, for a spatially-structured detection probability. Its SD is reported as sigma_p_field / ⁠sigma_p_field_<col>⁠.

Checkpoint / resume

A full-field occu_cover() fit integrates over a large outer hyperparameter grid and can run for hours, so a reboot or OOM kill otherwise loses the whole run. control$checkpoint = list(path = "fit.ckpt", resume = TRUE) makes the joint engine append each completed grid cell to path; a resume = TRUE run loads the finished cells and solves only the remaining ones, reproducing the from-scratch fit. resume = FALSE starts a fresh file (any stale checkpoint at path is removed first). A torn tail from a killed write is discarded and re-solved, and a checkpoint written for different data or settings is rejected rather than resumed onto. Forwarded to tulpa::tulpa_nested_laplace_joint().

Cell-aggregated cover (cover.aggregate)

On the shared-field spatial path the cover arm has one observation per detected visit, so a cell with many detected plots informs the shared field far more than the single occupancy observation for that cell; the field is then driven almost entirely by the cover arm and the detection-corrected occupancy surface flattens. cover.aggregate collapses the cover arm to a single response per occupancy unit (cell-period) so the two arms inform the field with comparable weight:

• "mean" – model the mean cover over the unit's detected visits;

• "median" – model the median cover;

• "latent" – model a per-unit cover random effect $u_i \sim N(0, \sigma_u^2)$ shared across the unit's detected visits and integrated out per unit, so the cover arm contributes one marginal observation per unit while keeping every detected visit (the principled counterpart of mean / median aggregation). The lognormal arm integrates in closed form (compound-symmetry); the beta arm uses adaptive Gauss-Hermite quadrature (control$n.quad, default 15). The within-unit dispersion is pre-fit from the within-unit spread and held fixed; $\sigma_u$ is integrated on the outer grid (control$sigma.u.grid).

• "none" – the per-visit cover arm (one cover observation per detected visit).

NULL (the default) selects "mean" on the shared-field spatial nested_laplace path and "none" on the non-spatial laplace path (which has no shared field to over-weight). Aggregation needs a cell-level positive design: a visit-level covariate in the positive formula cannot be collapsed to one value per cell and errors. Aggregation is currently wired on the spatial path only; requesting "mean" / "median" on the non-spatial laplace fit errors rather than silently using per-visit cover.

Random effects on the observation arms

Under method = "nested_laplace" (shared-field spatial) and method = "nuts" (non-spatial, which samples the group SD), random intercepts may be written on the detection or positive-cover formula with the usual lme4 bar or re() spelling, e.g. detection = ~ effort + (1 | habitat). The grouping is per visit – one code per ⁠(site, visit)⁠ entry, distinct from the per-site occupancy-arm grouping – so a many-level categorical visit covariate (an EUNIS habitat class, an observer) enters as a partially pooled random intercept. Crossed ((1 | habitat) + (1 | observer)) and nested ((1 | region/site), which expands to re(region) + re(region:site)) groupings are supported: each term joins the joint fit as its own iid latent block. Each block's variance integrates on the outer grid (sigma_re_p / sigma_re_pos for a lone term on the detection / cover arm, suffixed by the grouping variable – sigma_re_p_habitat – when several terms share an arm; tune with control$re.sigma.grid.p / re.sigma.grid.pos). The per-group BLUPs are reported in fit$re (one entry per term, keyed by arm or "<arm>:<var>") and via ranef(). predict() sums every term's BLUP offset on the predicted arm when newdata carries the grouping columns and type = "detection" (or "occurrence" / "cover_cond"); an unseen / held-out level shrinks that term to the population mean. Several crossed terms multiply the outer grid, so set control$integration = "ccd" (and/or coarsen the RE grids) at scale.

Random slopes are also supported (tulpa engine >= 0.0.39): an uncorrelated slope ((x || g), (0 + x | g)) is one per-row weighted iid block per coefficient, and a correlated slope ((1 + x | g)) one multivariate free-Sigma miid block. A slope term's fit$re entry carries an ⁠[n_groups x n_coefs]⁠ BLUP matrix, a per-coefficient sigma, and (correlated) a cor matrix; the hyperparameters are ⁠sigma_re_p_<coef>⁠ and ⁠cor_re_p_<ci>_<cj>⁠, marginalized over the grid. predict() weights each coefficient by its covariate column in newdata (intercept = 1). Each slope covariate is standardized to unit SD before fitting (mirroring the fixed-effect design), so the variance grid is scale-invariant – the reported slope BLUPs and sigma are back-transformed to the covariate's natural units (correlation is scale-free). A correlated slope adds p(p+1)/2 Sigma axes to the outer grid, so its free-Sigma grid uses a compact principled default (symmetric correlation nodes including 0); widen it with control$re.logchol.grid.p / re.logchol.grid.pos.

The positive-cover RE needs per-visit cover (cover.aggregate = "none"). As with the occupancy-arm RE, each grid-integrated variance carries the binary / small-cluster inner-Laplace attenuation and is a lower bound on the truth; the BLUPs recover the per-group structure.

method = "nuts" samples an observation-arm random intercept instead of integrating it on a grid: each grouping factor is a non-centered block (⁠b_g = sigma_re z_g⁠, z_g ~ N(0, 1)) whose group SD is a coordinate of the sampled vector under a N(0, 1.5^2) prior on ⁠log sigma_re⁠ – the same weakly informative width the other observation-family samplers use, and the reason the sampled SD carries no small-cluster attenuation to correct for. Crossed and nested groupings are simply several blocks, on either or both observation arms. fit$re keys and ranef() match the nested_laplace path; each term additionally reports sigma_median and the per-draw sigma_draws, the summaries to quote for a right-skewed variance component, with per-SD convergence in fit$nuts$re_sigma_rhat / re_sigma_ess. The sampled coefficient surface (coef(), vcov(), fit$draws) is unchanged. Still on nested_laplace only: a random slope, and an observation-arm RE composed with the coupled areal field. waic() / cpo() on a sampled fit read the coefficient draws, so they score the RE at zero – the same limitation the sampled coupled field has.

See Also

occu() (no cover), cover() (plot-level hurdle, no detection), abun() (counts not cover).

Examples

N <- 120; J <- 4
sim <- simulate_occu_cover(N = N, J = J, n.occ.covs = 1, n.det.covs = 1,
                           n.pos.covs = 1, positive = "beta", seed = 1)
long <- data.frame(site_id = rep(seq_len(N), each = J),
                   visit    = rep(seq_len(J), times = N),
                   y        = as.vector(t(sim$y)),
                   det_cov1 = sim$visit_data$det_cov1,
                   pos_cov1 = sim$visit_data$pos_cov1)
od  <- tobs_data(long, y = "y", site = "site_id", visit = "visit",
                 det.covs = c("det_cov1", "pos_cov1"))
cell_dat <- cbind(data.frame(site_id = seq_len(N)), sim$data)
y_pos <- sim$y_pos; y_pos[is.na(y_pos)] <- 0
fit <- tobs(~ occ_cov1, data = cell_dat, family = occu_cover("beta"),
            detection = ~ det_cov1, positive = ~ pos_cov1,
            y = od$y, y.pos = y_pos, visits = od$det.covs, method = "laplace")
summary(fit)

Description

Assemble the paired occurrence / cover response arms, the visit-level covariate grid, and the site-level design frame that an occu_cover() fit consumes, from one long (one row per site-visit / plot) data frame. This is the same construction tobs() runs internally when a single occu_cover() fit is handed a long frame; call it directly when you want to inspect the arms before fitting, or pass them on as ⁠y = ⁠, ⁠y.pos = ⁠, ⁠visits = ⁠, ⁠data = ⁠ to tobs().

The occurrence and cover arms are pivoted onto one shared ⁠(sites, visits)⁠ grid, so they align row-for-row, and the site-level design is the first row per site. The cover arm is positive-only (meaningful where the species is detected): in the dense grid absent / unsampled cover cells are set to 0, and in the compact (ragged) carrier they are NA in the values and never read by the joint engine.

Usage

occu_cover_inputs(
  data,
  site,
  visit,
  response,
  y.pos,
  occ.covs = NULL,
  det.covs = NULL,
  coords = NULL,
  compact = TRUE,
  positive = "beta"
)

Arguments

Value

A list with y, y_pos, visits, site_data, coords, n_sites, max_visits, n_visits, and compact.

See Also

tobs(), tobs_data().

Examples

sim <- simulate_occu_cover(N = 20, J = 3, seed = 1)
long <- data.frame(site = rep(1:20, each = 3), visit = rep(1:3, 20),
                   occur = as.vector(t(sim$y)),
                   cover = as.vector(t(ifelse(is.na(sim$y_pos), 0,
                                              sim$y_pos))),
                   det_cov1 = sim$visit_data$det_cov1,
                   occ_cov1 = rep(sim$data$occ_cov1, each = 3))
inp <- occu_cover_inputs(long, site = "site", visit = "visit",
                         response = "occur", y.pos = "cover",
                         det.covs = "det_cov1", compact = FALSE,
                         positive = "lognormal")
dim(inp$y)

Description

Joint occupancy of S species where the occupancy state ⁠z in {0,1}^S⁠ follows a log-linear model with first-order (per species) and second-order (per species pair) natural parameters (Rota et al. 2016; unmarked occuMulti). The second-order parameters are the species interactions (positive = co-occur more than independent; negative = avoid). The state formula is the shared occupancy covariate design (each natural parameter carries its own coefficients); detection the shared site-level per-species detection design. The latent state integrates out by enumerating the 2^S states, so the exact marginal is maximised with an observed-information vcov.

Usage

occu_multi()

Details

y is a length-S list of ⁠N x J⁠ 0/1/NA detection matrices (or a 3D ⁠[sites x visits x species]⁠ array); species names the arms.

Value

A tobs_family object for tobs().

Engines

• method = "laplace" – the enumerated 2^S-state marginal maximised (BFGS) with an observed-information covariance; the only route (no field, no sampler).

Examples

sim <- simulate_occu_multi(S = 2, N = 300, seed = 1)
fit <- tobs(~ scov1, data = sim$data, family = occu_multi(),
            detection = ~ 1, y = sim$y, species = sim$species,
            control = list(verbose = FALSE))
coef(fit)

Description

A cell-level occupancy gate, a plot-level availability gate, per-visit detection, and the cover hurdle, for data where a site's "visits" are spatially distinct plots aggregated into a ⁠(cell, period)⁠ rather than temporal revisits (the EVA / MOTIVATE vegetation layout; Nichols et al. 2008; Mordecai et al. 2011). occu_cover() treats plots as detection replicates of one occupancy state, which on spatial subunits conflates within-cell prevalence into the detection arm (Kendall & White 2009); this family separates them with an explicit middle level:

Usage

occu_multiscale_cover(response = c("beta", "lognormal", "gaussian"))

Arguments

Details

z_c        ~ Bernoulli(psi_c)                 # cell / range occupancy
a_cj | z=1 ~ Bernoulli(theta_cj)              # plot availability / use
y_cjv|a=1  ~ Bernoulli(p_cjv)                 # detection
cover|y=1  ~ f_pos(.; eta_pos, disp)          # cover hurdle (beta / lognormal)

Both z (over cells) and a (over plots) marginalize in closed form (two states each), so the joint marginal log-likelihood is exact and reuses the same nested-Laplace cell-coupling machinery as occu_cover().

Value

A tobs_family object.

Inputs

y / y.pos are ⁠[n_plots x max_visits]⁠ matrices (plots are the rows, the availability units; visits the columns). The state-process formula is the cell-level occupancy predictor and MUST carry an areal field naming the per-plot cell column, icar(graph = adj, group.var = "cell"). availability = ~ ... is the plot-level theta predictor (default ~ 1); detection the per-visit p predictor; positive = ~ ... the cover predictor (default the detection formula). y.pos is read only where y == 1.

A detected visit (y_cjv = 1) with a missing cover (y_pos_cjv = NA) keeps its detection term but drops the f_pos factor: cover is taken missing-at-random, so the cover likelihood runs over the detected visits with an observed cover. This is the rule occu_cover() applies, and it holds on every engine below.

Identifiability

The availability (theta) and detection (p) levels separate only with replication WITHIN a plot. Single releves supply none, so a plain fit identifies psi (cell) and the product theta * p (plot) – it reduces to occu_cover() with p := theta * p. Within-plot temporal replication (e.g. a resurvey of the same plot in a later period) makes the third level estimable.

Engines

Three engines. method = "nested_laplace" carries a single shared areal field coupled across the occupancy (sigma) and cover (alpha * sigma) arms, integrated over the outer ⁠(sigma, alpha)⁠ grid. method = "laplace" and method = "nuts" are the non-spatial path (iid cells, no field): the exact three-level marginal (z over cells, a over plots both summed in closed form) optimised directly ("laplace", a Gaussian observed-Fisher posterior) or sampled ("nuts", the exact coefficient posterior with calibrated intervals and WAIC / LOO). Cells are declared the same way on every path, via an icar(graph = adj, group.var = "<cell>") term (the graph drives the field under "nested_laplace" and supplies only the plot -> cell map under "laplace" / "nuts"). On the "nested_laplace" path additional weighted areal terms in the psi formula (⁠icar(graph = adj, group.var = "<cell>", weight = <cell covariate>)⁠) add spatially-varying-coefficient trend fields, each coupled onto the cover arm with its own alpha_trend; the fitted fields are in fit$trend_field / fit$trend_fields. The coupled / trend field is not sampled, so method = "nuts" takes a single cell-declaring areal term only.

The coupling is written as share(spatial()) in the positive formula, the same spelling occu_cover() takes; it needs a field to copy, so it is accepted under method = "nested_laplace" and refused on the two non-spatial engines, which fix the field at 0. Omitting it decouples the arms: the block is pinned at alpha = 0 and the field rides occupancy alone, the same meaning the absence has on occu_cover(). There is no implicit default coupling.

The copy amplitude's axis carries prior structure – a point mass at alpha = 0 and a log-spaced slab above it – so share(alpha = grid(...)) STATES its nodes, while share(alpha = grid(n = )) states a RESOLUTION: the engine re-reads its own axis with that many slab nodes, point mass and bounds unchanged. ⁠terms =⁠ gives either, per block; one block takes one of the two. control$alpha.grid[.trend] / control$alpha.n[.trend] are the wire format these compile into and are refused as user input.

See Also

occu_cover() (two-level), cover() (plot hurdle, no detection).

Examples

f <- occu_multiscale_cover(response = "beta")
f

Description

Constructs a prior specification consumed by tobs() when a Laplace method (method = "laplace" etc.) is used for the occupancy family. The penalty is applied as a quadratic term + sum((beta_j - mu_j)^2 / (2 * sd_j^2)) added to the negative-log-posterior in each M-step submodel block by tulpa::tulpa_laplace() (via the per-block beta_prior).

Usage

occu_priors(
  p.intercept = list(mean = 0, sd = 1.5),
  p.slope = list(mean = 0, sd = 2.5),
  beta.occ.intercept = list(mean = 0, sd = 2),
  beta.occ.slope = list(mean = 0, sd = 5)
)

Arguments

Details

The defaults are weakly informative: they pull the detection intercept toward p = 0.5 strongly enough to break the psi-p identifiability ridge at small J, but stay wide enough not to bias estimates at N >= 600 with informative covariates. Setting any sd to Inf disables that component of the prior (same penalised objective, just 1/Inf^2 = 0).

Value

An occu_priors object, ready to pass to tobs(..., priors = ...).

Examples

# disable the detection-slope penalty
priors <- occu_priors(p.slope = list(mean = 0, sd = Inf))
priors


sim <- simulate_occu(N = 100, J = 3, n.occ.covs = 1, n.det.covs = 1,
                     seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ det_cov1, y = sim$y, method = "laplace",
            priors = occu_priors(),
            control = list(verbose = FALSE, progress = FALSE))
coef(fit)

Description

Occupancy where a survey records the TIME to first detection rather than a 0/1 outcome (Garrard et al. 2008; unmarked occuTTD). At an occupied site the time-to-detection is exponential with rate lambda (constant hazard); a survey of length survey.length that reaches its end without a detection is censored. An unoccupied site never detects. The state formula models ⁠logit psi⁠ (occupancy); detection models ⁠log lambda⁠ (the site-level detection rate). The latent occupancy state integrates out in closed form (two states), so the fit maximises the exact marginal with an observed-information vcov.

Usage

occu_ttd(survey.length = 1)

Arguments

Details

y is an ⁠N x J⁠ matrix of detection times: a value in ⁠(0, survey.length)⁠ is a detection; a value ⁠>= survey.length⁠ is a non-detection (censored); NA is a survey not conducted.

Value

A tobs_family object for tobs().

Engines

• method = "laplace" – the exact marginal maximised (BFGS) with an observed-information covariance; the only route (no field, no sampler).

Examples

sim <- simulate_occu_ttd(N = 200, J = 4, Tmax = 3, seed = 1)
fit <- tobs(~ psi_cov1, data = sim$data, family = occu_ttd(survey.length = 3),
            detection = ~ rate_cov1, y = sim$y, control = list(verbose = FALSE))
coef(fit)

Description

A single-season occupancy detection-history dataset for 120 wetland sites surveyed on 4 visits, drawn from a known generative model (MacKenzie et al. 2002). Occupancy increases with site wetness and decreases with elevation; per-visit detection increases with survey effort. Twelve detection cells are set to NA to mimic incomplete histories. Synthetic, not field data; the generative coefficients are in truth.

Usage

peatland_occu

Format

A list with components:

y

120 x 4 integer matrix of detections (0/1/NA).

occ.covs

Data frame of site covariates: elevation, wetness (both standardised).

det.covs

Named list with one 120 x 4 visit covariate matrix, effort (standardised).

coords

120 x 2 matrix of x / y coordinates in the unit square.

truth

List of generative quantities: beta_psi, beta_p, the realised latent occupancy z, and per-site psi.

See Also

occu(), tobs(), simulate_occu()

Examples

data(peatland_occu)
str(peatland_occu$y)

fit <- tobs(~ elevation + wetness, data = peatland_occu$occ.covs,
            family = occu(), detection = ~ effort,
            y = peatland_occu$y, visits = peatland_occu$det.covs)
coef(fit)

Description

The tulpa::pit_residuals() method for tobs_fit: the family's own predictive CDF, with the randomisation step a discrete response needs.

Usage

## S3 method for class 'tobs_fit'
pit_residuals(
  object,
  n.samples = 250,
  nsim = NULL,
  seed = NULL,
  observed = NULL,
  ...
)

Arguments

Value

Numeric vector of PIT residuals.

See Also

tulpa::test_uniformity() for the uniformity test on the result.

Examples

sim <- simulate_occu(N = 100, J = 3, n.occ.covs = 1, n.det.covs = 1,
                     seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ det_cov1, y = sim$y, method = "laplace",
            control = list(verbose = FALSE, progress = FALSE))
pit <- pit_residuals(fit, seed = 1)
hist(pit, breaks = 10)

Description

Heatmap of the identifiability score over the scanned (cell size x temporal block) grid: information smallest-eigenvalue for score = "info", fraction of replicated units for score = "count". Identifiable cells are outlined.

Usage

## S3 method for class 'tobs_aggregation_scan'
plot(x, ...)

Arguments

Value

Invisible NULL.

Examples

set.seed(1)
d <- data.frame(east = runif(200), north = runif(200),
                year = sample(2001:2006, 200, replace = TRUE),
                detected = rbinom(200, 1, 0.3))
scan <- occu_aggregation_scan(d, response = "detected",
                              coords = c("east", "north"), year = "year",
                              cell.sizes = c(0.2, 0.4),
                              block.lengths = c(2, 3))
plot(scan)

Description

2x2 panel: detection frequency histogram, per-visit detection rates, visit completeness, and spatial detection map (if coordinates available).

Usage

## S3 method for class 'tobs_data'
plot(x, ...)

Arguments

Value

Invisible NULL.

Examples

sim <- simulate_occu(N = 30, J = 3, seed = 1)
long <- data.frame(site = rep(seq_len(30), times = 3),
                   visit = rep(1:3, each = 30),
                   det = as.vector(sim$y),
                   occ_cov1 = sim$data$occ_cov1,
                   effort = runif(90))
dat <- tobs_data(long, y = "det", site = "site", visit = "visit",
                 occ.covs = "occ_cov1", det.covs = "effort")
plot(dat)

Description

The discrepancy-statistic check (the spOccupancy ppcOcc construction): latent occupancy is drawn from its full conditional, a detection replicate from the fitted model, and the two are compared through a discrepancy statistic to give a Bayesian p-value. This is a test, not a plot, so it is its own generic rather than a method on tulpa::pp_check(), which draws the graphical check.

Usage

ppc(object, ...)

## S3 method for class 'tobs_fit'
ppc(object, fit.stat = c("freeman-tukey", "chi-squared"), n.samples = 500, ...)

Arguments

Value

A list with fit.y, fit.y.rep, and bayesian.p.

Examples

sim <- simulate_occu(N = 60, J = 3, seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ det_cov1, y = sim$y, method = "laplace",
            control = list(verbose = FALSE))
ppc(fit, n.samples = 100)$bayesian.p

Description

Occurrence probability is always p = plogis(X * beta_occ). The conditional positive cover mu depends on the positive-part family:

Usage

## S3 method for class 'cover_fit'
predict(
  object,
  newdata = NULL,
  type = NULL,
  include.RE = FALSE,
  times = NULL,
  time.col = NULL,
  level = 0.95,
  nsim = 1000L,
  draws = TRUE,
  weights = NULL,
  aggregate = FALSE,
  mc.tol = 0.05,
  nsim.max = 10000L,
  mc.floor = 0.001,
  ...
)

Arguments

Details

• positive = "lognormal": mu = exp(eta_pos + sigma_pos^2 / 2) (lognormal back-transform on log-cover).

• positive = "beta": mu = plogis(eta_pos) (mean of the beta on the natural cover scale with logit link).

Expected cover is E[y] = p * mu under both positive parts.

The separate-Laplace fit (method = "laplace") returns a fixed-effects-only numeric vector. The nested-Laplace shared-field fit (method = "nested_laplace") instead projects the shared occupancy-cover field and returns a tobs_prediction of posterior draws – the same tidy / change contract as predict.tobs_fit() for occu_cover(): pass type = "change" with times = c(t1, t2) and time.col for a per-cell delta map. Each prediction unit is a row of newdata (or a cell column indexing the field cells), and every quantity is marginalized per draw over the grid-integrated joint posterior.

Value

Separate-Laplace fit: a numeric vector. Nested-Laplace fit: a tobs_prediction.

Examples

fit <- tobs(cover ~ moisture, data = meadow_cover,
            family = cover(response = "beta"))
nd <- data.frame(moisture = c(-1, 0, 1))
predict(fit, newdata = nd, type = "occupancy")
predict(fit, newdata = nd, type = "conditional")
predict(fit, newdata = nd, type = "expected")

Description

Returns, per row of newdata, the presence probability psi, the conditional class probabilities p_1..p_K (given present), and the unconditional class probabilities psi * p_k plus the absent probability 1 - psi.

Usage

## S3 method for class 'occu_categorical_fit'
predict(object, newdata, ...)

Arguments

Value

A list with psi (length n), cond (n x K conditional class probabilities), and joint (n x (K+1): absent, then class 1..K).

Examples

sim <- simulate_occu_categorical(N = 300, seed = 1)
fit <- tobs(~ x, data = sim$data, family = occu_categorical(), y = sim$y)
pr <- predict(fit, newdata = data.frame(x = c(-1, 0, 1)))
round(pr$joint, 3)

Description

Four modes:

• In-sample: predict(fit) returns fitted values.

• State posterior / NA-response: predict(fit, type = "state") returns the marginalised per-site occupancy posterior of a "nested_laplace" fit: the per-row eta posterior is a Gaussian mixture over the hyperparameter grid (per-cell fitted linear predictor and predictive variance), so psi is its Gauss-Hermite mean and psi_lower / psi_upper are the mixture-CDF quantiles – a calibrated 95% credible interval, exact for every latent prior. Rows with heldout = TRUE are the INLA-style NA-response prediction targets (single-season sites whose detection history was all-missing), where occupancy is interpolated from the latent field rather than informed by detections.

• Design-matrix: predict(fit, X.0 = ...) predicts at new covariate values.

• Terms-based: predict(fit, terms = "elev") varies one covariate, others at mean.

• Joint occu_cover: for a occu_cover joint-coupled fit, predict(fit, newdata, type = "occurrence" | "cover_cond" | "cover_exp" | "change") samples the joint latent from the grid-integrated posterior (the outer-grid mixture via tulpa::tulpa_posterior_draws()) and marginalises every derived quantity per draw. type = "change" with times = c(t1, ..., tK) (K >= 2) returns a per-cell change table. Each quantity – occupancy psi, conditional cover cover_cond, expected cover cover_exp – has a level column per step (psi_T1..psi_TK, ...) carrying the posterior mean with .sd / .lwr / .upr at level, and a delta per step against t1 (delta_psi, delta_cover_cond, delta_cover_exp) with .lwr / .upr and .prob_pos, the posterior probability P(delta > 0) per cell. The expected-cover change is split exactly into its occupancy- and abundance-driven parts (delta_cover_from_occ, delta_cover_from_ab), reported the same way. With two times the delta columns are unsuffixed; with K > 2 each carries ⁠_T<k>⁠. Every step is evaluated on ONE draw set, so levels and deltas share a posterior and are jointly valid. A standalone occu() joint fit returns the same layout with occupancy only. The result is a tobs_prediction table (one row per cell) carrying per-unit ⁠[cell x nsim]⁠ draw matrices in attr(, "draws"); map it yourself, e.g. left_join(cents, pr, by = "cell") then geom_tile(aes(x, y, fill = delta_psi)) (or geom_sf() on polygon cells). type = "trajectory" with times = c(t1, ..., tK) predicts occupancy psi, conditional cover cover_cond and expected cover cover_exp at each time, for any times the user chooses, as a long table with the posterior mean, median and level interval of each (one row per cell and time); aggregate = TRUE averages over the cells per draw and gives the overall trajectory, one row per time. All times share one draw set.

• Community families: ms_occu(), ms_dyn_occu(), ms_int_occu(), ms_abun(), ms_distance() and ms_occu_cover() predict per-species matrices ⁠[rows x species]⁠, in-sample from fitted() and at newdata from the per-species coefficients. Each family answers to its own response types: "occupancy" / "detection" (ms_occu, ms_dyn_occu, ms_int_occu), "abundance" / "detection" (ms_abun), "lambda" / "sigma" (ms_distance), and "occupancy" / "detection" / "cover_cond" / "cover_exp" (ms_occu_cover). A fit carrying an areal field or a latent-factor loading is in-sample only, the field having no value at an unseen site.

• Spatial-factor community: for a reduced-rank spatial-factor ms_occu_cover() fit, predict(fit, type = "occupancy" | "cover_cond" | "cover_exp") returns the per-species per-cell posterior as a long table (one row per cell x species): "occupancy" gives psi, "cover_cond" the conditional cover mean E[cover | present], "cover_exp" the unconditional expected cover psi * E[cover | present], each with a ⁠_median⁠ and a ⁠_lower⁠ / ⁠_upper⁠ interval. The maps are marginalised over the loading + field posterior, so a rare species borrows strength across the shared factors for a calibrated map. The latent fields are tied to the cell graph, so X.0 / newdata prediction is not supported.

Usage

## S3 method for class 'tobs_fit'
predict(
  object,
  X.0 = NULL,
  type = c("occupancy", "detection", "both", "state"),
  quantiles = c(0.025, 0.5, 0.975),
  terms = NULL,
  n.points = 50L,
  newdata = NULL,
  times = NULL,
  level = 0.95,
  nsim = 1000L,
  draws = TRUE,
  time.col = NULL,
  X.det.0 = NULL,
  weights = NULL,
  aggregate = FALSE,
  mc.tol = 0.05,
  nsim.max = 10000L,
  mc.floor = 0.001,
  ...
)

Arguments

Value

Depends on mode. In-sample: fitted() result. "state": a data.frame with row, psi (marginalised posterior mean), psi_lower / psi_upper (equal-tailed 95% credible interval; NA when the engine did not return per-cell predictive variance), and heldout. Design-matrix/ terms: data.frame with estimate and CIs; type = "detection" on an int_occu() fit returns a named list of such data.frames (one per source), and type = "both" returns list(occupancy = , detection = ). occu_cover: a tobs_prediction table (one row per cell) with per-unit draw matrices in attr(, "draws"); type = "trajectory" a long table (cell, time, then ⁠<q>_mean⁠ / ⁠_median⁠ / ⁠_lwr⁠ / ⁠_upr⁠ for psi, cover_cond and cover_exp), one row per cell and time, or per time with aggregate = TRUE. A trajectory states its Monte Carlo precision in attr(, "mc_se_max") (the largest MC standard error of any reported median or bound, in posterior SDs) and attr(, "nsim_used"); under nsim = "auto" also attr(, "mc_tol") and attr(, "nsim_capped") (TRUE when nsim.max stopped the draws before mc.tol was met). attr(, "mc_binding") names the bound that sets mc_se_max: a one-row data frame with its cell (NA with aggregate = TRUE), time, quantity, bound ("median", "lwr" or "upr") and mc_se.

Examples

sim <- simulate_occu(N = 100, J = 3, n.occ.covs = 1, n.det.covs = 1,
                     seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ det_cov1, y = sim$y, method = "laplace",
            control = list(verbose = FALSE, progress = FALSE))
pr <- predict(fit, type = "occupancy")
summary(pr$psi)

Description

Print a batched multi-response occu_cover fit.

Usage

## S3 method for class 'tobs_batch'
print(x, ...)

Arguments

Value

x, invisibly.

Examples

sims <- lapply(1:2, function(s) simulate_occu_cover(N = 40, J = 3, seed = s))
y <- lapply(sims, function(s) s$y)
y_pos <- lapply(sims, function(s) { p <- s$y_pos; p[is.na(p)] <- 0; p })
fit <- tobs(~ occ_cov1, data = sims[[1]]$data,
            family = occu_cover("lognormal"), detection = ~ 1,
            y = y, y.pos = y_pos, species = c("a", "b"), method = "laplace")
print(fit)

Description

Print method for tobs_family

Usage

## S3 method for class 'tobs_family'
print(x, ...)

Arguments

Value

x, invisibly.

Examples

occu()
abun(mixture = "negbin")

Description

Print method for tobs_fit

Usage

## S3 method for class 'tobs_fit'
print(x, ...)

Arguments

Value

x, invisibly.

Examples

sim <- simulate_occu(N = 100, J = 3, n.occ.covs = 1, n.det.covs = 1,
                     seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ det_cov1, y = sim$y, method = "laplace",
            control = list(verbose = FALSE, progress = FALSE))
fit

Description

Returns the per-group random-effect posterior summaries. Under NUTS the non-centred draws are reconstructed to the natural BLUP scale (⁠b_{g,c} = sigma_c * (L z_g)_c⁠) and summarised; the deterministic Laplace path returns the variance-component EM modes and their Schur-complement standard errors. A community fit reports its per-species deviations from the community mean. Every family and engine returns the same columns.

Usage

## S3 method for class 'tobs_fit'
ranef(object, ...)

Arguments

Value

A data frame with one row per random-effect level and coefficient and the columns arm (the linear predictor the effect enters: psi, p, lambda, ...; psi+p for an effect shared across arms; NA when the fit does not distinguish arms), group (the grouping variable; "species" for the per-species deviations of a community fit), level (the group level, or the species), term (the coefficient the effect shifts, "(Intercept)" for a random intercept), estimate and std.error (NA where the fitter reports no posterior SD). Zero rows when the fit carries no random effects.

On a community fit, std.error is the per-species posterior SD of the deviation, the square root of the matching diagonal entry of the species' posterior covariance fit$ms_community$Cinv[[s]] (the analogue of lme4's conditional variance). On the Laplace routes that covariance is conditional on the community means; on a spatial fit integrated over an outer hyperparameter grid it combines the per-node covariances and the spread of the per-node BLUPs (law of total covariance); on the NUTS and Polya-Gamma Gibbs routes it is the covariance of the sampled deviations. The ms_abun() log_r (negative binomial) and logit_omega (zero-inflated) deviations carry an SD as well: every route that fits them (joint AGHQ, and NUTS for log_r) keeps them in the covariance. std.error is NA for a species whose per-species posterior solve failed.

Examples

sim <- simulate_occu(N = 150, J = 4, n.occ.covs = 1, n.det.covs = 1,
                     n.visit.groups = 8, sigma.visit = 1, seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ (1 | visit_group), visits = sim$visits,
            y = sim$y, method = "laplace",
            control = list(verbose = FALSE, progress = FALSE))
ranef(fit)

Description

Sequential-depletion removal sampling: latent abundance N_i ~ Poisson(lambda_i) (or negative binomial) observed through K ordered removal passes, where pass k removes ⁠Binomial(N_i - sum_{l<k} y_{il}, p_{ik})⁠ of the individuals still present. The declining catch sequence identifies detection p and abundance N; the latent N is summed out in closed form (truncation K.max), so the fit is a direct Laplace approximation (no EM), with a NUTS path over the same marginal.

Usage

removal(K.max = NULL, mixture = c("poisson", "negbin"))

Arguments

Details

The tobs() formula is the abundance (⁠log lambda⁠) model; detection is the per-pass detection (⁠logit p⁠) model. The response y is an ⁠n_sites x K⁠ integer matrix of per-pass removals with the passes in column order; complete pass sequences are required (no NA).

Value

A tobs_family object.

Engines

• method = "laplace" – a direct Laplace approximation over the closed-form marginal (Poisson or negative binomial); a grouped random effect is integrated by adaptive Gauss-Hermite quadrature.

• method = "nested_laplace" – an areal icar() / car_proper() field on the abundance arm, optionally composed with a temporal() field, through the shared areal-BFGS driver.

• method = "nuts" – the in-tree C++ NUTS sampler over the same marginal, with a fixed-hyper areal or temporal field block on the abundance arm (one of the two, not both).

References

Royle, J. A. (2004). N-mixture models for estimating population size from spatially replicated counts. Biometrics 60, 108-115. Dorazio, R. M., Jelks, H. L., Jordan, F. (2005). Improving removal-based estimates of abundance. Biometrics 61, 1093-1101.

Examples

f <- removal(K.max = 100)
f

Description

One contract for every family: a list with occ, the STATE-level residual series (per site, or per site and season for a multi-season family), and det, the OBSERVATION-level series (per visit / pass / distance bin). A family carrying only one of the two levels fills that one and leaves the other NULL; nothing returns a bare matrix, so a caller reading ⁠$occ⁠ gets the state residual or NULL, never an error.

Usage

## S3 method for class 'tobs_fit'
residuals(object, type = c("deviance", "pearson", "response"), ...)

Arguments

Value

A list with occ (state-level) and det (observation-level) residuals, either of which may be NULL for a family that has no such level. A model type with no registered handler is an error naming it.

Examples

sim <- simulate_occu(N = 100, J = 3, n.occ.covs = 1, n.det.covs = 1,
                     seed = 1)
fit <- tobs(~ occ_cov1, data = sim$data, family = occu(),
            detection = ~ det_cov1, y = sim$y, method = "laplace",
            control = list(verbose = FALSE, progress = FALSE))
res <- residuals(fit, type = "pearson")
str(res)

Description

Occupancy with abundance-induced detection heterogeneity (Royle & Nichols 2003; unmarked occuRN). Latent abundance N_i ~ Poisson(lambda_i) drives the per-visit detection probability 1 - (1 - r_ij)^{N_i}, where r_ij is the per-individual detection probability. The state formula models ⁠log lambda⁠ (abundance); detection models ⁠logit r⁠. Detection is site-level by default; passing visits makes it visit-varying (⁠logit r_ij⁠ gains the visit-level covariates), exactly as for the occupancy / N-mixture front doors. The latent N marginalises in closed form (a Poisson sum to K.max), so the fit maximises the exact marginal with an observed-information vcov.

Usage

royle_nichols(K.max = NULL)

Arguments

Value

A tobs_family object for tobs().

Engines

• method = "laplace" – the exact marginal maximised (BFGS) with an observed-information covariance; the only route (no field, no sampler).

Examples

sim <- simulate_royle_nichols(N = 150, J = 5, seed = 1)
fit <- tobs(~ x, data = sim$data, family = royle_nichols(),
            detection = ~ 1, y = sim$y, control = list(verbose = FALSE))
coef(fit)

# Visit-varying detection via `visits`:
sv <- simulate_royle_nichols(N = 150, J = 5, beta.r.visit = 0.8, seed = 1)
fv <- tobs(~ x, data = sv$data, family = royle_nichols(),
           detection = ~ w, y = sv$y, visits = sv$visits,
           control = list(verbose = FALSE))
coef(fv)

Description

Runs the POSTERIOR simulation-based calibration experiment (Talts et al. 2018; Sailynoja, Schmitt, Buerkner and Vehtari 2026, Algorithm 2) on a fitted tobs_fit, through the engine's tulpa::sbc() front door. The truth is drawn from the fit's own posterior at the observed data, a replicate is simulated at that truth on FRESH cells, the model is refitted on the two data sets together, and the truth's rank under that augmented posterior is recorded. Under exact inference those ranks are Uniform(0, 1); the departure from uniformity is a verdict on the approximation.

Usage

## S3 method for class 'tobs_fit'
sbc(
  object,
  n.sim = 100L,
  n.draws = 1000L,
  n.ref = 200L,
  quantities = NULL,
  controls = character(),
  bad.factor = 1.25,
  level = 0.95,
  seed = 0L,
  model.only = FALSE,
  fit.control = list(),
  control = list(),
  ...
)

Arguments

Details

Registered families: occu_cover (the coupled occupancy + cover hurdle on the joint nested-Laplace engine), and occu, abun, count, removal, distance, fp_occu, royle_nichols, occu_ttd and double_observer, whose site marginals multiply. For those the replicate is the family's own simulate() kernel at the drawn theta, the rank arm is the family's exact marginal log-likelihood, and two kinds of fit are refused rather than approximated: a structured term (its field is a latent quantity shared across sites that theta does not hold, which is what the coupled occu_cover route handles by drawing on fresh cells) and a visit-level observation design (the replicate kernel assembles the observation arm from the site-level design, so the two sides would carry different models).

dyn_occu (constant-rate only) pools on the SITE axis and leaves the season axis alone, since its response is 3D ⁠[n_sites x max_visits x n_seasons]⁠: a bespoke replicate generator forward-simulates the colonization/extinction HMM, but the posterior draws and the rank statistic are the same shared machinery every other registered family uses.

int_occu (full source-site overlap only) pools on the SITE axis and leaves the per-source axis alone, since its response is a list of one detection matrix per source; the rank statistic is the shared machinery unchanged. Registering it surfaced a real bug (since fixed): the detection intercept was silently left in standardized-covariate units on every fit with an autoscaled detection covariate, independent of SBC.

gdistremoval pools two response matrices (yDist band counts, yRem period counts) together on the SITE axis, reusing the same list-of-matrices pool int_occu uses; its replicate generator is its own simulate() handler.

dyn_abun (constant-rate only) shares dyn_occu's 3D ⁠[n_sites x max_visits x n_seasons]⁠ response and site-axis pooling, but unlike dyn_occu it already has a working simulate() handler (the Dail-Madsen open N-mixture forward is not a two-state HMM that needs a bespoke generator), so its replicate is the same family-simulate() route as occu/abun/etc.

occu_categorical is a tobs_multiarm_fit with two independent Laplace-Gaussian blocks (presence, class) rather than a joint MVN draw matrix; draws() samples the two blocks independently and loglik_many() scores the two-arm likelihood directly, since this family has no .tobs_pointwise_loglik dispatch to reuse.

distsamp_open (constant-dynamics, Poisson only) shares dyn_abun's 3D response and site-axis pooling, and its fit$means/fit$draws are the standard .tobs_bfgs_marginal_fit() shape, so draws/simulate/ loglik_many are all the shared generic ones.

occu_multi's response is a list of S per-species matrices, the same shape int_occu() pools; its simulate() is custom, since species share one joint multi-species state rather than being independently observed.

dyn_int_occu is the product of the multi-season and multi-source shapes: a named list of S per-source 3D arrays, pooled by a new .tobs_sbc_pool_named_3d that composes both existing pooling rules.

t_occu is a pg_gibbs family whose fit$draws is already the real pooled posterior sample; it has no simulate() handler and no .tobs_pointwise_loglik dispatch, so both are custom – simulate draws a fresh AR1 year effect at the theta's own ⁠(sigma, rho)⁠, and loglik_many is a Laplace approximation to that year effect's marginal (FD-validated, brute-force cross-checked at T = 2).

cover is a tobs_multiarm_fit with the same two-independent-block shape as occu_categorical; positive = "lognormal" only for v1, with dispersion held fixed (no SE is reported for it anywhere in the package).

occu_multiscale_cover is a standard single-block fit (unlike cover's two-block one); the exchangeable unit is the CELL rather than the plot, so pooling and site labeling track cell indices.

A family is registered by adding one entry to the internal registry – a replicate generator, a refit call and optionally a joint statistic; the pooling, the grouping labels, the arm construction and the controls are shared.

Value

An object of class sbc from tulpa, carrying the per-quantity rank ECDFs, the exact simultaneous band, and the uniformity tests. When model.only = TRUE, the callback list instead.

What is scored

The arm fixed effects, the areal field SD on each arm, the copy scale alpha, and the dispersion when the fit estimates one. The per-cell field itself is not scored – it is integrated out by the fit and carries no truth theta could hold. The cover-arm field SD and alpha are read per draw off the outer hyperparameter grid, whose cells are sampled by their own normalized weight, so what is ranked is the grid-marginalized posterior of each rather than a function of component modes.

A quantity whose posterior has no spread – a dispersion the joint engine holds fixed in its cell-coupling spec, for instance – is dropped rather than scored, and named in attr(x, "fixed").

An extra log_lik quantity ranks a joint statistic of the whole parameter vector, which catches a posterior getting each marginal right while getting the dependence between them wrong.

Fresh cells

The replicate is drawn on new cells whose areal graph is block-diagonal against the observed one. occu_cover couples its occurrence and cover arms through one shared field, and theta carries no per-cell field value, so a replicate re-observing the observed cells would depend on the observed data given theta and break the factorization posterior SBC rests on. group_ids is supplied to tulpa::sbc(), which verifies the observable half of that (disjoint group labels) rather than assuming it.

Controls

controls = c("wide", "narrow") adds arms reporting the same draws rescaled about their mean by bad.factor and its reciprocal. They are deliberately mis-scaled posteriors and are expected to leave the band: a calibration read that nothing can fail is not evidence.

References

Talts, S., Betancourt, M., Simpson, D., Vehtari, A. and Gelman, A. (2018). Validating Bayesian inference algorithms with simulation-based calibration. arXiv:1804.06788.

Sailynoja, T., Schmitt, M., Buerkner, P.-C. and Vehtari, A. (2026). Posterior SBC: simulation-based calibration checking conditional on data. Statistics and Computing 36:78.

See Also

tulpa::sbc(), waic()

Examples

N <- 20L; J <- 3L
adj <- matrix(0L, N, N)
for (s in seq_len(N)) {
  if (s > 1L) adj[s, s - 1L] <- 1L
  if (s < N)  adj[s, s + 1L] <- 1L
}
sim <- simulate_occu_cover(N = N, J = J, positive = "lognormal",
                           adj = adj, sigma = 0.7, alpha = 1, seed = 1L)
long <- data.frame(site_id = rep(seq_len(N), each = J),
                   visit = rep(seq_len(J), times = N),
                   y = as.vector(t(sim$y)),
                   det_cov1 = sim$visit_data$det_cov1,
                   pos_cov1 = sim$visit_data$pos_cov1)
od <- tobs_data(long, y = "y", site = "site_id", visit = "visit",
                det.covs = c("det_cov1", "pos_cov1"))
y_pos <- sim$y_pos; y_pos[is.na(y_pos)] <- 0
# The dispersion goes on the outer grid so it is estimated and scored, and
# the field-SD grid is pinned so both stages integrate the same support.
ctl <- list(engine = "joint", verbose = FALSE,
            sigma.grid   = exp(seq(log(0.15), log(2.0), length.out = 9)),
            phi.grid.pos = exp(seq(log(0.20), log(0.90), length.out = 7)))
fit <- tobs(~ occ_cov1 + icar(graph = adj),
            data = cbind(data.frame(site_id = seq_len(N)), sim$data),
            family = occu_cover("lognormal"),
            detection = ~ det_cov1,
            positive = ~ pos_cov1 + share(spatial()),
            y = od$y, y.pos = y_pos, visits = od$det.covs,
            method = "nested_laplace", control = ctl)
sbc(fit, n.sim = 20L, controls = "narrow", fit.control = ctl)