Package 'tulpa'

Description

Build the symmetric adjacency matrix that spatial() and spatial_car() consume, from a common spatial layout, instead of hand-coding it. The model still receives an explicit graph (spatial(graph = g$adjacency, ...)), so the graph stays inspectable before fitting – the engine never guesses connectivity from coordinates silently.

adjacency() is a single front-door verb that dispatches on the kind of input:

a data.frame

cell centroids on a regular grid – queen (edge or corner) or rook (edge only) contiguity over the lattice. This covers both plain coordinate grids and rasterised cells passed as a table of centres.

an sf object

polygon contiguity from shared boundaries (queen = share a point or edge; rook = share an edge). Requires sf.

a SpatRaster (terra) or stars object

raster cells – the lattice of (non-NA) cell centres, queen or rook. Requires terra or stars respectively.

The result is a tulpa_adjacency object holding the matrix plus the cell identifier for each node, so a model can pass g$adjacency while the observation data is remapped to 1-based node indices with node_index().

Usage

adjacency(x, ...)

## Default S3 method:
adjacency(x, ...)

## S3 method for class 'data.frame'
adjacency(
  x,
  type = c("queen", "rook"),
  x_coord = "x",
  y_coord = "y",
  id = NULL,
  order = 1L,
  tolerance = 1.5,
  offsets = NULL,
  ...
)

## S3 method for class 'sf'
adjacency(x, type = c("queen", "rook", "touches"), id = NULL, ...)

## S3 method for class 'SpatRaster'
adjacency(
  x,
  type = c("queen", "rook"),
  order = 1L,
  tolerance = 1.5,
  offsets = NULL,
  na_rm = TRUE,
  ...
)

## S3 method for class 'stars'
adjacency(
  x,
  type = c("queen", "rook"),
  order = 1L,
  tolerance = 1.5,
  offsets = NULL,
  na_rm = TRUE,
  ...
)

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

Arguments

Value

A tulpa_adjacency object: a list with

adjacency

the ⁠[n x n]⁠ symmetric sparse adjacency matrix (dgCMatrix, 0/1, zero diagonal) to pass as graph / adjacency.

ids

the cell identifier for each node, in node order (length n).

n

the number of nodes.

cellsize

the inferred cell size c(x, y) for grid / raster layouts, or NA for polygons.

type

the contiguity rule used.

See Also

node_index() to map cell identifiers to node indices, check_adjacency() to validate a hand-built matrix, spatial() and spatial_car() which consume the graph.

Examples

# A 3 x 3 regular grid of cell centres
grid <- expand.grid(x = 1:3, y = 1:3)
grid$cell <- paste0("c", seq_len(nrow(grid)))

g <- adjacency(grid, x_coord = "x", y_coord = "y", id = "cell",
               type = "queen")
g
g$adjacency

# Second-order (24-neighbour) queen contiguity: any order is settable
g2 <- adjacency(grid, id = "cell", order = 2)

# Advanced: a custom stencil (symmetrized for the undirected field)
g3 <- adjacency(grid, id = "cell", offsets = list(c(1, 0), c(0, 1)))

# Use it in a model: graph stays explicit and inspectable
spatial(graph = g$adjacency, formula = ~ 1 || cell_idx)

# Remap observation data (original cell ids -> 1:n node indices) by key
obs <- data.frame(cell = c("c5", "c1", "c5", "c9"))
obs$cell_idx <- node_index(g, obs$cell)
obs

Description

Marginal-likelihood approximation for a generalised linear mixed model with one cluster-level intercept random effect. Adaptive Gauss-Hermite quadrature (AGQ) generalises Laplace by replacing the single-point Gaussian integral around the cluster's posterior mode with n_quad-point quadrature. AGQ at n_quad = 1 recovers the Laplace approximation; higher n_quad reduces approximation error, especially for clusters with few observations or non-Gaussian likelihoods.

This is an engine block, not a front door: its tuning knobs (max_iter, tol, n_quad) sit in the signature rather than in a control list.

Scoped to the lme4::glmer(..., nAGQ = N) use case: one intercept- only RE term, families binomial, poisson, or gaussian. For multi-RE or random-slope models, use tulpa_laplace() (Laplace at the joint mode) or HMC.

Usage

agq_fit(
  y,
  X,
  group,
  n_groups = max(group),
  family = c("binomial", "poisson", "gaussian"),
  n_trials = NULL,
  sigma_eps = 1,
  n_quad = 7L,
  offset = NULL,
  beta_init = NULL,
  sigma_init = 1,
  max_iter = 200L,
  tol = 1e-06,
  verbose = FALSE
)

Arguments

Value

A list with class tulpa_fit carrying:

• means: named vector c(beta, log_sigma_re).

• mode: same as means (Gaussian fit at the optimum).

• cov: covariance matrix of ⁠(beta, log_sigma_re)⁠ from inverse observed information.

• sigma_re: estimated RE SD.

• log_marginal: AGQ-approximated log marginal likelihood at the fitted values.

• n_quad: quadrature order used.

• n_iter, converged: optimiser diagnostics.

• inference_mode: "structured".

• inference_tier: 2L.

• backend: "agq".

An optimum at which the AGHQ objective is undefined (some group's solve fails there, so its value is the failure sentinel rather than a marginal likelihood) is an error naming the groups, not a returned fit.

Tier

Tier 2 (Structured). AGQ is exact in the limit n_quad -> infinity but at finite n_quad it is a controlled approximation – same epistemic class as Laplace.

References

Pinheiro, J. C., & Bates, D. M. (1995). Approximations to the log-likelihood function in the nonlinear mixed-effects model. Journal of Computational and Graphical Statistics, 4(1), 12-35.

See Also

tulpa_laplace() for the joint-mode Laplace path (handles multi-RE and spatial structure).

Examples

set.seed(1)
n_g <- 20; n_per <- 5; n <- n_g * n_per
group <- rep(seq_len(n_g), each = n_per)
x <- rnorm(n)
X <- cbind(1, x)
u <- rnorm(n_g, 0, 0.8)
eta <- 0.3 + 0.7 * x + u[group]
y <- rbinom(n, 1, plogis(eta))

# Compare Laplace (n_quad = 1) and AGQ-7.
fit_lap <- agq_fit(y, X, group, family = "binomial", n_quad = 1)
fit_agq <- agq_fit(y, X, group, family = "binomial", n_quad = 7)
c(fit_lap$log_marginal, fit_agq$log_marginal)

Description

AIC and BIC penalise a maximised log-likelihood. logLik.tulpa_fit() reports one only for a fit that maximised over every parameter it reports (quantity = "log_likelihood", as agq_fit() does); a mean log posterior, a log evidence and a conditional log marginal likelihood are not, and both criteria refuse them rather than return a number. Compare those fits by their evidence, or by compare_models(criterion = "waic") / "loo".

Usage

## S3 method for class 'tulpa_fit'
AIC(object, ..., k = 2)

## S3 method for class 'tulpa_fit'
BIC(object, ...)

Arguments

Value

As stats::AIC() / stats::BIC() for maximised log-likelihoods; errors otherwise.

Description

Convert a fit's posterior draws to a posterior draws object. as_draws() returns the draws_array shape; as_draws_array(), as_draws_matrix(), as_draws_df() and as_draws_rvars() return theirs. When the posterior package is installed these are also registered against its generics, so posterior::as_draws(fit) works.

Usage

as_draws(x, ...)

## S3 method for class 'tulpa_fit'
as_draws(x, n_draws = NULL, seed = NULL, ...)

as_draws_array(x, ...)

## S3 method for class 'tulpa_fit'
as_draws_array(x, n_draws = NULL, seed = NULL, ...)

as_draws_matrix(x, ...)

## S3 method for class 'tulpa_fit'
as_draws_matrix(x, n_draws = NULL, seed = NULL, ...)

as_draws_df(x, ...)

## S3 method for class 'tulpa_fit'
as_draws_df(x, n_draws = NULL, seed = NULL, ...)

as_draws_rvars(x, ...)

## S3 method for class 'tulpa_fit'
as_draws_rvars(x, n_draws = NULL, seed = NULL, ...)

Arguments

Details

Fits differ in whether they carry draws at all. Sampler and nested-Laplace fits do, and convert directly. A Gaussian-approximation fit (mode = "laplace", mode = "eb") carries a mode and a precision instead, and converting it means drawing from the approximation – which is a modelling decision, not a format change, because every downstream posterior summary would then treat a normal approximation as a posterior sample. So it is opt-in: pass n_draws to synthesize that many draws from N(coef(object), vcov(object)), or get an error naming the alternative. Synthesized draws form a single chain and cover the fixed effects only.

Value

A posterior draws object of the requested shape.

See Also

tulpa_draws_array() for the base R array without the dependency, posterior_sample() for the raw matrix.

Examples

set.seed(1)
df <- data.frame(x = rnorm(80))
df$y <- rpois(80, exp(0.5 + 0.3 * df$x))
fit <- tulpa(y ~ x, data = df, family = "poisson")
if (requireNamespace("posterior", quietly = TRUE)) {
  posterior::summarise_draws(as_draws(fit))
}

Description

Declares that a setting shaping the outer hyperparameter grid carries a default the caller computed, not a choice the user made. The auto-recenter pass (outer_grid_placement) leaves a user-pinned setting exactly as given, and re-centres (or, for a prior, engages its own regularizer over) a marked one when the fit rails against its ceiling.

Four kinds of setting take the mark:

• a grid axis on a nested-Laplace prior block (sigma_grid, tau_grid, ...) – a numeric vector of nodes, or the ⁠[n_cells x k]⁠ matrix of pre-paired coordinates the families whose axis is a matrix take (mcar / miid's logchol_grid, tgmrf's theta_grid_built);

• an entry of tulpa_nested_laplace_joint()'s phi_grid – one arm's dispersion axis;

• a scalar grid-construction knob in control, for a driver that builds its axes rather than taking them (fit_st_nested()'s n_grid_spatial, tau_upper, ...);

• a prior_sigma hyperprior specification – a list, e.g. list("pc.prec", c(U = 3, alpha = 0.01)).

Wrapper packages are the intended caller: one that builds a default of its own – because it derives a second axis from it, hands the same vector to several blocks, or exposes its own argument with a default – would otherwise be indistinguishable from a user who pinned that setting deliberately. Mark it and the rescue stays live. A setting whose value is exactly the engine's own default is recognised without a mark; anything else needs one – and a dispersion axis ALWAYS needs one, since the engine has no default of its own to recognise it against.

The mark is an attribute, so it is dropped by sort(), [, c() and as.numeric(): build the value first, mark it last.

Usage

auto_grid(x, place = TRUE)

Arguments

Value

x carrying the marker attribute. Numeric input is coerced to double IN PLACE, so everything else it carries – dim() and dimnames() above all – survives the mark; a list is returned unchanged apart from the attribute.

Declaring a default the engine must integrate as written

place = FALSE separates the two questions the mark otherwise answers at once. Provenance – whose choice these nodes are – and placement policy – whether the pass may move them – are different questions, and a package that measured its own default as the one to integrate has an answer to the second that is not "the user pinned it". Such an axis is treated exactly as a pin everywhere the engine ACTS on it (the placement pass leaves it, refinement densifies within its span rather than following the posterior past the end nodes, and no curvature is computed for it), and differs only in what the fit REPORTS: "default_axis_pinned" rather than "axis_pinned", so a reader is not told they pinned an axis they never wrote down. Leaving the mark off instead buys the same integration and says the wrong thing about it.

place is read wherever the engine would otherwise RE-PLACE the setting: a grid axis, and a scalar grid-construction knob. A prior_sigma specification is not re-placed but replaced, by the engine's own regularizer, and carries no place.

See Also

is_auto_grid(), auto_grid_place(), tulpa_nested_laplace_joint(), fit_st_nested()

Examples

prior <- list(type = "icar", sigma_grid = auto_grid(c(0.1, 0.5, 1, 2, 3)))
is_auto_grid(prior$sigma_grid)
is_auto_grid(auto_grid(list("pc.prec", c(U = 3, alpha = 0.01))))
auto_grid_place(auto_grid(c(0.5, 1, 2), place = FALSE))

Description

Reads back what auto_grid()'s place argument recorded, so a wrapper that rebuilds a value (as.numeric() drops every attribute) can re-apply both halves of the declaration rather than only the provenance half.

Usage

auto_grid_place(x)

Arguments

Value

FALSE when x was marked auto_grid(place = FALSE), TRUE otherwise – including for a value carrying no mark at all, which the engine holds because it reads as a pin rather than because it asked to be held.

See Also

auto_grid(), is_auto_grid()

Examples

auto_grid_place(auto_grid(c(0.5, 1, 2)))
auto_grid_place(auto_grid(c(0.5, 1, 2), place = FALSE))

Description

Per posterior draw s, R2_s = Var_i(mu_si) / (Var_i(mu_si) + Var_res_s), where mu_si are the response-scale fitted means and Var_res_s is the family's residual variance averaged over observations (Gelman et al. 2019). The linear predictor is rebuilt per draw exactly as in posterior_predict() (fixed effects + random effects + offset at the training data).

Usage

bayes_R2(object, ...)

## S3 method for class 'tulpa_fit'
bayes_R2(
  object,
  ndraws = NULL,
  summary = TRUE,
  probs = c(0.025, 0.975),
  seed = NULL,
  ...
)

Arguments

Value

With summary = TRUE (default) a one-row data frame with estimate (posterior median), std.error, and the probs quantiles; with summary = FALSE the vector of per-draw R^2 values.

References

Gelman, Goodrich, Gabry & Vehtari (2019). R-squared for Bayesian regression models. The American Statistician 73(3):307-309.

Examples

set.seed(1)
d <- data.frame(x = rnorm(200))
d$y <- rnorm(200, 2 * d$x, 1)
fit <- tulpa(y ~ x, data = d, family = "gaussian", mode = "laplace", phi = 1)
bayes_R2(fit)

Description

Estimate the log marginal likelihood $\log Z = \log p(y)$ from posterior draws via the Meng-Wong / Gronau bridge-sampling identity. Useful for Bayes factors, model comparison, and as a sanity check on Laplace-approximated marginal likelihoods.

The classical iterative scheme (Meng & Wong 1996; Gronau et al. 2017) is used with a multivariate-normal proposal fit to half of the posterior draws – the other half is used for the bridge ratio so the proposal is independent of the samples that score it. All computations are done in log-space via logsumexp for numerical stability.

Usage

bridge_sampling(
  draws,
  log_posterior,
  n_proposal = nrow(draws),
  max_iter = 1000L,
  tol = 1e-10,
  split = TRUE,
  verbose = FALSE
)

Arguments

Value

A list with:

• log_marginal: the bridge-sampling estimate of $\log Z$.

• n_iter: bridge iterations to convergence.

• converged: logical.

• re_sd: relative MSE estimate (Fruehwirth-Schnatter 2004); a rough quality check, smaller is better.

• proposal: the fitted proposal list(mean, cov).

Tier

Bridge sampling is a post-hoc marginal-likelihood estimator, not a sampling backend, so it does not appear in the INFERENCE_TIERS registry. It operates on draws from any Tier-1 (Exact) backend.

References

Meng, X.-L., & Wong, W. H. (1996). Simulating ratios of normalizing constants via a simple identity: a theoretical exploration. Statistica Sinica, 6, 831-860.

Gronau, Q. F., Sarafoglou, A., Matzke, D., Ly, A., Boehm, U., Marsman, M., ... & Steingroever, H. (2017). A tutorial on bridge sampling. Journal of Mathematical Psychology, 81, 80-97.

Examples

# Toy: marginal likelihood of N(theta | 0, 1) under Y = theta + eps,
# eps ~ N(0, 1), prior theta ~ N(0, 10). Closed-form available.
y <- 1.5
log_post <- function(theta) {
  dnorm(y, theta, 1, log = TRUE) + dnorm(theta, 0, 10, log = TRUE)
}
draws <- matrix(rnorm(2000, mean = y * 100 / 101, sd = sqrt(100 / 101)),
                ncol = 1)
bs <- bridge_sampling(draws, log_post)
bs$log_marginal  # should match dnorm(y, 0, sqrt(101), log = TRUE)

Description

stats::fitted(), stats::residuals(), stats::predict() and posterior_predict() for a categorical response: a tulpa_multinomial() fit (family = "multinomial") or a tulpa_ordinal() fit (family = "ordinal"). The response is one of K classes, so the response-scale quantities are N x K matrices with one column per class.

• fitted() returns the class probabilities at the posterior mode.

• residuals() returns, per class, the indicator of the observed class minus its fitted probability ("response"), or that difference divided by the indicator's standard deviation sqrt(p (1 - p)) ("pearson"). A categorical response has no single scalar residual; each class column is the Bernoulli residual of that class's indicator.

• predict() returns the link scale – the K - 1 baseline-category logits of the multinomial model (N x (K - 1)), or the ordinal location predictor ⁠x' beta⁠ (length N) – or the class probabilities (type = "response").

• posterior_predict() draws one class per observation for each posterior draw and returns the class indices, with the response levels attached; simulate.tulpa_fit() turns them into factor columns.

The pointwise log-likelihood behind cpo(), dic(), loo::waic() and loo::loo() is the log probability of each observed class at each stored posterior draw; dic() evaluates it at the posterior-mean parameters.

Usage

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

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

## S3 method for class 'tulpa_categorical'
predict(
  object,
  newdata = NULL,
  type = c("link", "response"),
  se.fit = FALSE,
  level = 0.95,
  ...
)

## S3 method for class 'tulpa_categorical'
posterior_predict(object, newdata = NULL, ndraws = NULL, seed = NULL, ...)

## S3 method for class 'tulpa_categorical'
pp_check(object, ndraws = 50, ...)

Arguments

Value

fitted(): an N x K probability matrix. residuals(): an N x K matrix. predict(): a matrix (or, for the ordinal link scale, a vector); with se.fit = TRUE a list of fit, se.fit, lower, upper of that shape. posterior_predict(): an ⁠ndraws x N⁠ integer matrix of class indices with attributes levels and ordered.

See Also

tulpa_multinomial(), tulpa_ordinal()

Description

Check that a hand-built adjacency matrix is a well-formed spatial graph before passing it to spatial() / spatial_car(): square, finite, symmetric, zero on the diagonal, 0/1 valued, free of negative weights and of isolated nodes, and connected. adjacency() runs the same checks on the graphs it constructs. The spatial constructors refuse a graph with a missing, non-finite or negative entry.

Usage

check_adjacency(adjacency, ids = NULL)

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

Arguments

Value

Invisibly, a tulpa_adjacency_check list with the per-check results (square, finite, symmetric, zero_diag, binary, nonneg, isolated-node indices, n_components and component_sizes, edge count) and an overall ok flag. Issues are reported via warning() and printed; the function does not stop, so every problem surfaces in one pass.

See Also

adjacency() to construct a graph, node_index().

Examples

adj <- matrix(0, 4, 4)
for (i in 1:3) adj[i, i + 1] <- adj[i + 1, i] <- 1
check_adjacency(adj)

Description

Runs the core convergence diagnostics on a fit and reports whether Rhat, bulk-ESS, and divergence thresholds are all met. A terse companion to diagnostic_summary() for use in scripts and tests.

Usage

check_diagnostics(
  fit,
  rhat_threshold = 1.01,
  ess_threshold = 400,
  quiet = FALSE
)

Arguments

Value

Invisibly, TRUE if all checks pass, FALSE if any fail, or NA for a non-chain (approximation) fit where Rhat/ESS do not apply.

See Also

diagnostic_summary(), diagnostics(), n_divergent()

Examples

set.seed(123)
df <- data.frame(x = rnorm(60))
df$y <- rpois(60, exp(0.5 + 0.3 * df$x))
fit <- tulpa(y ~ x, data = df, family = "poisson", mode = "hmc",
             control = list(n_iter = 500L, warmup = 250L, n_chains = 2L,
                            seed = 1L))
check_diagnostics(fit)

Description

Produces a 2x2 or 1x3 panel of diagnostic plots:

1. QQ plot of PIT residuals vs Uniform (with KS p-value)

2. Residuals vs fitted (with lowess smoother)

3. Dispersion histogram (observed variance vs simulated)

4. Spatial correlogram via Moran's I (if coords provided)

Usage

check_model(object, ...)

## Default S3 method:
check_model(object, coords = NULL, nsim = 250L, seed = 123L, ...)

Arguments

Value

Invisible list with ks_p, disp_ratio, moran (if spatial)

Examples

set.seed(1)
d <- data.frame(x = rnorm(80))
d$y <- rpois(80, exp(0.3 + 0.5 * d$x))
fit <- tulpa(y ~ x, data = d, family = "poisson", mode = "laplace")
check_model(fit, nsim = 50)

Description

Fixed-effect coefficients

Usage

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

Arguments

Value

Named numeric vector of fixed-effect posterior means (the Laplace mode for the Laplace tier). Random effects come from ranef().

Description

Rank fitted models best-first by an information criterion. "waic" and "loo" use the native pointwise log-likelihood layer (tulpa_criteria() / tulpa_psis()) – WAIC and PSIS-LOO respectively, computed from each fit's ⁠[n_draws x n_obs]⁠ pointwise log-likelihood (fit$draws$log_lik), with no loo package dependency. "loglik" returns the (integrated) joint log-likelihood with the parameter count. A fit carrying no pointwise log-likelihood (a deterministic / point approximation) yields NA criterion columns rather than an error, so the table always has one row per model.

"wbic" ranks by wbic(), an estimate of each fit's negative log evidence, and needs every fit to carry its tempered draws (mode = "smc", control = list(wbic = TRUE)); a fit without them gets an NA row.

Usage

compare_models(..., criterion = c("waic", "loo", "wbic", "loglik"))

Arguments

Value

A data frame. For "wbic" (ranked best-first): model, wbic, delta (WBIC gap to the best model, >= 0), and weight (exp(-delta) normalized: the posterior model probability under equal prior model weights, with each log evidence read off its WBIC). For "loglik": model, n_params, logLik, quantity and conditioned_on (which log-scale quantity logLik() reports, and the hyperparameters it holds fixed; fits that disagree on either are refused rather than ranked). For "waic" / "loo" (ranked best-first): model, elpd, se_elpd, p_eff, ic (-2 * elpd), delta (elpd gap to the best model), se_diff (SE of that pointwise elpd difference), and weight (the Akaike-style weight on the criterion).

See Also

model_average() for model-averaged predictions, tulpa_criteria() and tulpa_psis() for the native criteria layer.

Examples

set.seed(1)
df <- data.frame(x = rnorm(120))
df$y <- rpois(120, exp(0.4 + 0.5 * df$x))
f1 <- tulpa(y ~ x, data = df, family = "poisson")
f2 <- tulpa(y ~ 1, data = df, family = "poisson")
compare_models(full = f1, null = f2, criterion = "waic")

Description

Credible intervals for the fixed effects

Usage

## S3 method for class 'tulpa_fit'
confint(object, parm = NULL, level = 0.95, ...)

Arguments

Value

Matrix with lower and upper columns, labelled as stats::confint.default() labels them ("2.5 %" and "97.5 %" at the default level). A nested-Laplace fit carries interval_source / interval_declined (which read produced the bounds – by default the grid's Gaussian-mixture CDF), retained_mass (the share of the grid weight the bounds are conditional on), and skew_applied, one logical per reported coefficient saying whether its bounds are the inner-Laplace skew-corrected quantiles or not. See summary.tulpa_fit().

Description

Generic front doors onto the two criteria tulpa_criteria() computes that the loo package owns no generic for. WAIC and PSIS-LOO have theirs (loo::waic(), loo::loo()), so a model package registers methods on those rather than on new names that would mask them.

Usage

pointwise_loglik(object, ...)

## Default S3 method:
pointwise_loglik(object, ...)

## S3 method for class 'tulpa_fit'
pointwise_loglik(object, ndraws = NULL, ...)

dic(object, ...)

## Default S3 method:
dic(object, loglik_at_mean = NULL, ...)

## S3 method for class 'tulpa_fit'
dic(object, ...)

cpo(object, ...)

## Default S3 method:
cpo(object, ...)

## S3 method for class 'tulpa_fit'
cpo(object, ...)

## S3 method for class 'tulpa_fit'
waic(x, ...)

## S3 method for class 'tulpa_fit'
loo(x, ...)

Arguments

Details

pointwise_loglik() is the one door onto the ⁠[n_draws x n_obs]⁠ matrix itself, the input every criterion above is computed from. compare_models() and model_average() call it (through an internal wrapper) rather than assuming the engine's own fit layout, so a model package that registers ⁠pointwise_loglik.<its fit class>()⁠ – alongside its own waic() / loo() / dic() / cpo() methods – reaches model comparison and averaging too.

The default methods take a draws x observations pointwise log-likelihood matrix, the same input tulpa_criteria() takes. A model package registers a method taking its own fit object, builds the matrix from the posterior, and delegates here.

The tulpa_fit methods build that matrix from the fit itself, from the same source compare_models(criterion = "waic") reads: a log_lik the backend stored with its draws, else the family density evaluated at posterior draws of the in-sample linear predictor (sampler draws, the outer-grid mixture of a nested-Laplace fit, or Gaussian draws at the fixed-effect mode and covariance). The draws are pinned to a fixed internal seed, so repeated calls return the same numbers and leave the session RNG untouched. dic() plugs in the posterior mean of the linear predictor, so its p_dic counts effective parameters in the linear predictor's parameterization; where the fit stored its log_lik and carries no linear predictor draws, the DIC fields are NA and dbar is still reported. A fit with no pointwise log-likelihood (no stored response, a family that is not one built-in family name, or a linear predictor the fit cannot reproduce) is refused with an error naming the fit's class and the reason.

loo::waic() and loo::loo() dispatch to the tulpa_fit methods once loo is loaded, and return loo's own waic / psis_loo objects, so loo::loo_compare() reads them. On an MCMC chain fit loo() passes relative effective sample sizes computed over the fit's chains; on an i.i.d. or approximation fit the draws are independent and r_eff is 1.

Value

dic() and cpo() return a tulpa_criteria object. waic() returns a loo waic object and loo() a loo psis_loo object.

See Also

tulpa_criteria() for every criterion at once and for what the LOO unit means; compare_models() to rank several fits.

Examples

set.seed(1)
y  <- rnorm(40)
mu <- matrix(rnorm(200 * 40, sd = 0.2), 200, 40)
ll <- dnorm(matrix(y, 200, 40, byrow = TRUE), mean = mu, log = TRUE)
cpo(ll)

d <- data.frame(x = rnorm(120))
d$y <- rpois(120, exp(0.4 + 0.6 * d$x))
fit <- tulpa(y ~ x, data = d, family = "poisson", mode = "laplace")
dic(fit)
cpo(fit)$lpml
if (requireNamespace("loo", quietly = TRUE)) loo::waic(fit)

Description

Provides a comprehensive diagnostic report for a tulpa model, combining convergence metrics, divergence information, and actionable recommendations.

Usage

diagnostic_summary(fit, quiet = FALSE)

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

Arguments

Value

A list with class tulpa_diagnostic_summary containing:

status

Overall status: "PASS", "WARN", or "FAIL"

n_divergent

Number of divergent transitions

divergent_pct

Percentage of divergent transitions

worst_rhat

Data frame of parameters with worst Rhat

worst_ess

Data frame of parameters with worst ESS

e_bfmi

E-BFMI value (HMC only)

pareto_k, quad_ess

approximation fits only: the outer PSIS k-hat, or the grid quadrature ESS when no k-hat was produced

hyper_share_min, hyper_share_max

approximation fits only: the range over coefficients of the share of the fixed-effect marginal variance the hyperparameter integration contributed, between / (within + between) under the law of total variance over the outer grid. Near zero means the reported interval is what conditioning on a single hyperparameter point would have given. Reported without a threshold and never scored: a small share is not by itself a miscalibrated interval. NA when the fit retains no per-cell pieces to decompose

pareto_k_declined

approximation fits only, and only when there is no k-hat: WHY – "not_requested" and "unguessable_axis: <axis>" are benign or permanent, "degenerate_proposal" and "grid_too_small" are signals about the fit, and "internal_inconsistency" is an engine bug and raises the status to "WARN"

inner_skew_max, inner_skew_declined

approximation fits only: the largest scored inner-Laplace ⁠|gamma_3|⁠, or why nothing was scored ("coupled_arm" marks arms the inner layer could score neither per observation nor through the cell tensor)

axis_fields_dropped

data frame of grid axes the fit's own resolved path could not read and dropped as engine defaults: one row per field, with the block, family, path and the axis that path integrated instead. Absent whenever every supplied axis was used

recommendations

Character vector of recommendations

See Also

check_diagnostics(), diagnostics(), plot_diagnostics()

Examples

set.seed(123)
df <- data.frame(x = rnorm(60))
df$y <- rpois(60, exp(0.5 + 0.3 * df$x))
fit <- tulpa(y ~ x, data = df, family = "poisson", mode = "hmc",
             control = list(n_iter = 500L, warmup = 250L, n_chains = 2L,
                            seed = 1L))
ds <- diagnostic_summary(fit)
print(ds)

Description

The diagnostic front door for any tulpa fit. What a fit's posterior draws can be asked depends on how they were produced, so this reads the draws provenance and returns the diagnostic that applies:

MCMC chain draws

improved Rhat (the maximum of rank-normalized split-Rhat and folded split-Rhat), bulk / tail / mean / sd / quantile effective sample size, and Monte Carlo standard errors, following Vehtari et al. (2021). Split-Rhat is defined for a single chain, so a result is produced for any number of chains.

i.i.d. approximation draws

the approximation-reliability table – the PSIS tail-shape pareto_k scored against the exact inner-Laplace marginal and the outer-grid quadrature effective sample size (the OUTER hyperparameter-grid integration layer), the inner-Laplace skewness diagnostic gamma_3 when computed (the INNER Gaussian approximation to the latent field, a separate layer pareto_k does not cover), a combined whole-fit verdict naming which layer degrades when one does, and a per-parameter posterior summary. See the "Approximation reliability" sections below for the full description of this table and its attributes.

point summaries

no sample to diagnose; returns NULL with a message naming the backend.

Provenance is read from ⁠$draws_kind⁠ (stamped by tulpa_dispatch), falling back to the backend registry's emits property and then to an inner ⁠$joint_fit⁠. A fit that predates the tag is treated as a chain, so an older fit is never silently refused.

Usage

diagnostics(fit, ...)

## Default S3 method:
diagnostics(
  fit,
  pars = NULL,
  measures = c("rhat", "ess_bulk", "ess_tail"),
  probs = c(0.05, 0.95),
  sbc = NULL,
  ...
)

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

## S3 method for class 'sbc'
diagnostics(fit, ...)

Arguments

Value

For a chain fit, a data frame with a parameter column followed by one column per requested measure; entries are NA for parameters that are constant or have too few draws. For an i.i.d. approximation fit, the laplace_diagnostics table (see the "Approximation-reliability table columns" section below for its attributes). For a point fit, NULL.

Extending

diagnostics() is an S3 generic so a model package can answer for its own fit class (diagnostics.ratiod_fit(), say) rather than shadowing this export with a same-named function. The default method does the provenance routing described above and is what a method should delegate to once it has assembled a draws array.

Calibration alongside reliability

The tables above score ONE fit's own internal reliability. Whether the backend's posterior is CALIBRATED is a different question, answered over many simulated data sets by sbc(), and the two disagree in both directions – a fit whose reliability band is clean on both layers can still fail calibration, and one whose outer k-hat is well past the escalation threshold can pass. So the band is a screen, not a verdict. Pass an sbc result as ⁠sbc =⁠ and the calibration verdict is attached to the table and printed underneath it; diagnostics() called on the sbc result itself returns its report table.

Approximation reliability (i.i.d. draws)

Per-parameter reliability diagnostics for a fit whose posterior draws are i.i.d. samples from a deterministic approximation (the nested-Laplace grid-mixture posterior ⁠sum_k w_k N(mode_k, V_k)⁠), where the between-chain Gelman-Rubin Rhat reported for a chain fit does not apply. This is the table that plays Rhat's role for the deterministic engine: it answers "did the approximation work", not "did the chains mix".

The headline is a Pareto-smoothed importance-sampling (PSIS) reliability diagnostic for the OUTER hyperparameter-grid integration. The nested integrator scores its hyperparameter grid against the exact inner-Laplace marginal posterior with a generalized-Pareto fit to the upper tail of the importance ratios ⁠log p_target(theta) - log q_proposal(theta)⁠ (see tulpa_psis()); the resulting tail-shape pareto_k is the "did the outer integration work" number – k-hat < 0.5 good, 0.5-0.7 usable, ⁠>= 0.7⁠ unreliable (Vehtari et al. 2024; Yao et al. 2018). It is computed at fit time and read back here; a fit that did not run the diagnostic, or whose grid proposal degenerated, reports it as NA and is assessed on the grid quadrature reliability instead.

pareto_k scores the outer integration only. A high pareto_k on a fit whose grid quadrature is healthy (ess_grid well above 1, largest cell weight modest) flags outer-integration (CI-width) calibration in the right-skewed hyperparameter tail and does not by itself invalidate the point estimates, which the grid quadrature governs.

outer_regime qualifies what a high pareto_k means, and is the reason a bare threshold on pareto_k is not a reliability verdict. A sharp hyperparameter posterior collapses the grid onto ~1 cell (ess_grid near 1); the outer integration has then degenerated to a point evaluation at the modal hyperparameter, so pareto_k is scoring how well a Gaussian at that mode stands in for the hyperparameter marginal, not how well a grid integrated it. Where the dominant cell is INTERIOR to the grid the collapse is benign – the grid bracketed the mode, the estimate is empirical Bayes there, and only integrated hyperparameter uncertainty is missing. Where it sits at a grid BOUNDARY the grid may simply be too narrow: grid_edge_axes / grid_edge_sides name the axes to widen. On a fit whose pareto_k cleared the good band, the outer diagnostic also fits a skew-normal proposal and reports the marginal's estimated skewness as outer_skew_max, so an inflated k-hat that was purely the symmetric proposal's mismatch with a skewed variance-component marginal is both corrected and explained. A skew-normal has Gaussian tails, so this can never mask a genuinely heavy-tailed target.

The grid quadrature reliability – the effective sample size ess_grid = 1 / sum(w_k^2) of the outer integration weights and the largest single cell weight – is always computed from the stored grid: a grid that collapses onto one cell (ess_grid near 1) integrates no hyperparameter uncertainty, while a spread grid does.

A SEPARATE layer – the inner Gaussian Laplace approximation to the latent-field conditional posterior pi(x | theta, y), which pareto_k does not cover – is scored by inner_skew: the leading-order Edgeworth skewness estimate gamma_3 (Rue, Martino & Chopin 2009 Sec 3.2.3) at the fitted MAP grid cell, computed when control$diagnose_skew = TRUE (the default) on the fitting call. Reading a high pareto_k alone as "the fit is broken" conflates the two layers: an occu_cover batch flagged 42/78 species "unreliable" on outer k-hat alone when their point estimates, governed by the healthy inner layer, were fine – the reliability attribute is the combined verdict that names which layer degrades, if either does.

Each parameter row also carries the rank-normalized split-Rhat and bulk / tail effective sample size of the draws (Vehtari et al. 2021). On i.i.d. draws these sit at ~1.00 and ~n_draws by construction; they are reported, clearly as i.i.d.-draw Monte-Carlo diagnostics and not chain mixing, to document that the reported posterior summaries are not Monte-Carlo-limited.

A posterior sample is what those per-parameter rows are computed from, and nothing else here needs one: every reliability quantity above is read off the fit. So a fit that carries no draws – a default single-block nested-Laplace fit, whose posterior is the retained outer-grid mixture rather than a sample – reports the full band with an empty per-parameter body and n_draws = NA, and records why in the param_table_declined attribute. tulpa_posterior_draws() samples that mixture where the rows are wanted. The one case that still returns NULL is a fit with neither draws nor any reliability quantity, such as a plain Laplace fit with no outer grid to score.

Approximation-reliability scope

The PSIS pareto_k diagnoses the OUTER (hyperparameter) integration: whether the Gaussian-proposal-over-grid approximation of the marginal hyperparameter posterior p(theta | data) can be importance-corrected to the exact inner marginal. This is the dominant approximation in nested Laplace and the one with an exactly evaluable target. A full latent-space PSIS against the exact joint posterior pi(x) is not computed: the latent prior marginal ⁠p(x) = integral p(x | theta) p(theta) dtheta⁠ has no closed form, and for the marginalized-occupancy / cover-hurdle likelihoods the exact joint density is evaluable only inside the C++ kernel, so a stored fit cannot reconstruct it. The grid quadrature reliability is the complementary stored-fit number.

The inner layer carries a SECOND score, inner_pareto_k, which needs no likelihood derivative at all and therefore answers where inner_skew declines. The inner Gaussian at the fitted hyperparameter is an importance proposal for the exact conditional posterior, and the joint density is the target, so PSIS on that ratio scores the inner approximation directly. It is computed on the same probed indices along the same conditional-mean curve, one dimension per index, since importance sampling degrades with dimension and a k-hat over the whole latent field would report n_x rather than the approximation. A Pareto shape index is scale-free, so it is banded only on indices whose realized importance efficiency shows a correction worth describing; inner_pareto_k_uniform records that none did, which is what a well-approximated inner layer looks like.

inner_skew diagnoses the INNER (latent-field) Laplace: whether the Gaussian approximation to pi(x_i | theta, y) is itself a good fit, at each scored latent index i. gamma_3 is exact for a gaussian-family coefficient (the log-likelihood is exactly quadratic in eta) and declines to NA – never a silently-wrong 0 ("perfectly Gaussian") – whenever no per-observation third-derivative oracle is available: a coupled multi-process likelihood (e.g. tulpaObs's occu_cover, whose arms combine non-separably through a CellCouplingSpec) or a family with no registered third derivative. Only the requested latent indices are scored (every arm's fixed-effects coefficients by default – see control$skew_idx), since each index costs one extra linear solve; the full latent field is not scored by default on a large spatial field.

Approximation-reliability table columns

The table has one row per parameter – parameter, mean, sd, n_draws, mcse_mean, moment_source – carrying attributes. There is deliberately no rhat / ⁠ess_*⁠ column: those are what the draws-provenance gate withholds on a non-chain fit, and this table is where it dispatches instead. Where the draws were sampled from an outer-grid mixture the fit carries (attr(draws, "grid_mixture"), stamped by tulpa_posterior_draws() on a single-block fit, by tulpa_re_cov_nested() and by model packages), mean / sd are that mixture's exact moments, mcse_mean is 0 and moment_source is "mixture_moments"; otherwise they are the draw moments, mcse_mean is sd / sqrt(n_draws) and moment_source is "draws".

pareto_k

the outer PSIS reliability k-hat (NA if not computed).

pareto_k_band

"good" / "ok" / "unreliable" / NA.

pareto_k_declined, pareto_k_declined_note

when pareto_k is NA, WHY: "not_requested", "not_applicable", "unguessable_axis" (naming the axis – a permanent limitation of that family, so read ess_grid instead), "draws_too_few", "grid_too_small", "no_varying_axis", "degenerate_proposal", or "internal_inconsistency" (an engine bug worth reporting), plus a one-line reading of it.

pareto_k_is_ess

importance-sampling ESS on the smoothed weights.

pareto_k_ci_low, pareto_k_ci_high, pareto_k_se_boot, pareto_k_band_confident, diagnose_draws

the bootstrap interval of the outer k-hat, its standard error, whether the interval sits inside one band, and the draw count it used; present when the fit ran the bootstrap.

k_quality_requested, k_quality_reached, k_quality_best, k_quality_reason, k_quality_rounds

the outcome of a control$k_quality request; present when one was made.

inner_skew, inner_skew_idx

the per-index inner-Laplace skewness gamma_3 and the latent indices it was probed at; present when the probe ran.

ess_grid, n_grid, rel_ess_grid, max_weight

grid quadrature reliability. n_grid counts the cells SOLVED, so a posterior sharp enough to underflow all but one of them reads ⁠ess_grid = 1.00 of 124 cells⁠ rather than of one, and rel_ess_grid is the share of the solved grid the quadrature actually uses.

outer_regime

"spread" / "collapsed_interior" / "collapsed_edge" – whether the outer grid integrated hyperparameter uncertainty at all, and if not whether its dominant cell is interior (empirical Bayes at the mode: point estimates sound, hyperparameter uncertainty not integrated) or against a grid boundary (widen it).

grid_edge_axes, grid_edge_sides

for an edge collapse, the axes the dominant cell sits against and on which side.

grid_edge_mass_axes

axes holding material weight on one of their own boundary nodes, as axis:side – that axis truncates its own marginal there whether or not its mode sits on the node, which is the weaker and more common statement grid_railed_axes cannot make. Read in every regime, a spread grid included.

outer_skew_max

largest estimated |skewness| of the hyperparameter marginal, computed only when the k-hat triggered the skew-normal proposal rescue (NA means the Gaussian proposal already fit, not "symmetric and unchecked").

outer_regime_note

a one-line reading of a collapsed regime or of boundary mass on an axis, absent on a spread grid carrying neither.

grid_railed_axes

outer axes whose OWN marginal is maximal at one of their own endpoints, as axis:side – the span does not contain that axis's posterior mode, so its marginal is a truncated tail at any spacing. Reported whether or not the engine was allowed to, able to, or built to move the axis.

grid_placement, grid_recentred_axes, grid_placement_declined, grid_placement_note

whether the outer grid was re-centred, on which axes, and – when it was not – why.

scope

the outer diagnostic's scope string.

scope_backend, scope_title, scope_layers, scope_text

what the table describes: the backend it was read for, the header title, the reliability layers that backend has (c("outer", "inner") for a nested-Laplace fit, empty for a backend with neither), and the scope sentence – for a nested fit what its outer layer integrates over, for any other backend what its draws are and which reliability quantity it computes. A backend without the two layers has no whole-fit verdict, and its reliability is NA.

inner_skew_max

the largest ⁠|gamma_3|⁠ among the scored latent indices (NA if control$diagnose_skew = FALSE or nothing scored).

inner_skew_band

"good" / "ok" / "unreliable" / NA, banded on inner_skew_max by the general skewness-magnitude convention (Bulmer 1979) – not a Rue-Martino-Chopin-specific cutoff.

inner_skew_scored, inner_skew_probed

how many of the probed latent indices returned a finite gamma_3 vs how many were probed.

inner_skew_declined, inner_skew_arms_declined, inner_skew_declined_note

when nothing was scored, WHY: "coupled_arm" (STRUCTURAL – the coupled arms have neither a per-observation sum nor a cell third-derivative tensor to read, so the outer k-hat is the only reliability number this fit has), "curvature3_unavailable", "no_finite_contribution", "no_probe_indices", "not_requested", "backend_unsupported", "solve_failed", or "not_converged" (the probe re-solve stopped short of a mode, so neither inner score has a point to read); the arms (1-based) a joint fit had no oracle for, which is also set on a PARTIALLY scored fit; and a one-line reading.

inner_pareto_k, inner_pareto_k_band

the inner-Laplace importance k-hat over the probed subspace, and its band on the same convention as the outer k-hat. Available wherever a mode was found, including a fit gamma_3 cannot score.

inner_pareto_k_rel_ess, inner_pareto_k_is_ess

the smallest realized importance efficiency and effective sample size across the probed indices – how much correcting the inner Gaussian actually needs, which is what makes the scale-free shape above readable.

inner_pareto_k_uniform

TRUE when no probed index carried a material correction, i.e. the inner Gaussian reproduces the conditional posterior over the sampled region.

inner_pareto_k_scored, inner_pareto_k_probed

how many probed indices returned a finite k-hat vs how many were probed.

inner_pareto_k_declined, inner_pareto_k_declined_note

when it is NA, WHY, from the same closed vocabulary the outer k-hat uses.

reliability

the combined whole-fit verdict: "reliable" only when both layers are good; otherwise names which layer is scoped or flags both as unreliable. The inner layer enters through the worse of its two scores, so a fit whose cubic term declined is still assessed.

and a trailing summary attribute (a one-row data frame of the headline numbers) for printing.

References

Vehtari, Gelman, Simpson, Carpenter & Burkner (2021). Rank-normalization, folding, and localization: an improved Rhat for assessing convergence of MCMC. Bayesian Analysis 16(2):667-718.

Vehtari, Simpson, Gelman, Yao & Gabry (2024). Pareto smoothed importance sampling. JMLR 25(72):1-58.

Yao, Vehtari, Simpson & Gelman (2018). Yes, but did it work?: Evaluating variational inference. ICML, PMLR 80:5581-5590.

Rue, Martino & Chopin (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. JRSS-B 71(2):319-392.

See Also

sbc() for calibration, the approximation-reliability sections above for the i.i.d.-fit table in full, tulpa_draws_array(), tulpa_posterior_draws(), tulpa_psis(), plot_rhat(), plot_ess(), diagnostic_summary(), check_diagnostics()

Examples

set.seed(1)
df <- data.frame(x = rnorm(60))
df$y <- rpois(60, exp(0.5 + 0.3 * df$x))

# chain fit -> Rhat / ESS
hmc <- tulpa(y ~ x, data = df, family = "poisson", mode = "hmc",
             control = list(n_iter = 500L, warmup = 250L, n_chains = 2L,
                            seed = 1L))
diagnostics(hmc)

# deterministic fit -> PSIS approximation reliability
smc <- tulpa(y ~ x, data = df, family = "poisson", mode = "smc")
diagnostics(smc)

Description

Tests first-order autocorrelation in temporally-ordered residuals.

Usage

durbin_watson(
  object,
  alternative = c("two.sided", "greater", "less"),
  time = NULL,
  resid_type = "pearson"
)

Arguments

Value

An htest object with DW statistic, lag-1 r, and p-value

Examples

set.seed(1)
e <- as.numeric(arima.sim(list(ar = 0.6), n = 100))
durbin_watson(e)
# Residuals recorded out of time order:
tt <- sample(100)
durbin_watson(e[tt], time = tt)

Description

Computes E-BFMI (Betancourt 2016) from the per-iteration Hamiltonian energies an HMC or NUTS sampler records: the variance of the energy transition between successive iterations over the variance of the marginal energy. Each chain is diffed separately, since a transition across a chain boundary compares unrelated trajectories, and the lowest chain value is reported: one badly exploring chain is enough to distrust a fit. Values below 0.3 indicate that momentum resampling moves the sampler too little between energy level sets.

Usage

energy_bfmi(energy, chain = NULL)

Arguments

Value

A list with value (the lowest per-chain E-BFMI, NA when no chain has two finite, non-constant energies), per_chain, energy_centred (energies centred within each chain) and energy_diff (within-chain transitions), the last two pooled across chains for plotting.

Examples

set.seed(1)
energy <- c(cumsum(rnorm(200, sd = 0.3)), rnorm(200, sd = 1))
energy_bfmi(energy, chain = rep(1:2, each = 200))$per_chain

Description

The response families ⁠family = ⁠ accepts on the engine's front doors, in their canonical (default-link) spelling. A family's ⁠*_<link>⁠ variants (e.g. "gamma_inverse") are listed separately by linked_family_names().

Usage

family_names()

Value

A character vector of family names.

Examples

family_names()

Description

The posterior of a fit's spatial field, one row per spatial unit: an areal region (a row of the adjacency matrix) or an SPDE mesh node. The value reported is the field's contribution to the linear predictor at that unit, the quantity a covariate is read against.

Usage

field_effects(object)

Arguments

Details

Where the posterior comes from follows what the fit computes:

• A ModelData sampler fit (mode = "hmc" and its siblings) with an ICAR, BYM2 or proper-CAR field: every draw's field contribution is evaluated through the engine's own linear-predictor assembly, so an intrinsic field is centred as its density centres it and a BYM2 field is the combined sigma * (sqrt(rho) * s * phi + sqrt(1 - rho) * theta) at that draw's sigma and rho (s the scaling factor). A unit that no observation falls in contributes nothing to the linear predictor and is reported as NA.

• A Polya-Gamma Gibbs fit (mode = "gibbs") stores the field it samples (phi_spatial[k], plus theta_spatial[k] for BYM2, summed); the rows summarize those draws.

• A nested-Laplace fit with an ICAR or proper-CAR field: each outer cell contributes its conditional mode, and its marginal variance where the fit retained one, weighted by the cell's posterior weight. With the variances the summary is the exact mean, SD and quantiles of that Gaussian mixture (source = "mixture"); without them the weighted mean of the cell modes, with no spread (source = "mode"). A BYM2 field is reported where the fit stores its grid-weighted combined field (fit_st_nested()); a tulpa() nested BYM2 fit carries its two components with cell-dependent mixing and errors.

• A conditional-Laplace areal fit (mode = "laplace" / "eb") reports the field at its mode, with no spread. For BYM2 this is the combined field at the sigma = 1, rho = 0.5 the conditional kernel fixes.

• An SPDE fit (fit_spde()) reports the mesh-node field: its draws on the sampler path, its mode (the posterior-mean hyperparameters' fit on the nested path) otherwise.

A fit whose field cannot be located this way (a continuous GP / NNGP / HSGP field, a restricted field, an inline spatial() varying-coefficient field) errors naming the field and the backend, as does a fit with no spatial field. A temporal field is read by temporal().

Value

A data frame with one row per spatial unit: unit (the node index), label (the adjacency row name where it has one, else the index), estimate, sd, the 2.5% / 97.5% bounds conf.low / conf.high, and source: "sampled" for a posterior draw summary, "mixture" for a nested-Laplace Gaussian mixture, "mode" for a point at the mode (sd and the bounds NA). The field type is attached as attribute "field".

See Also

smooth_effects() for an s(x) term, ranef() for random effects, temporal() for a temporal field.

Examples

set.seed(1)
K <- 12
W <- matrix(0, K, K)
for (i in 1:K) { j <- if (i < K) i + 1 else 1; W[i, j] <- W[j, i] <- 1 }
region <- sample(1:K, 200, replace = TRUE)
d <- data.frame(region = region, x = rnorm(200))
d$y <- rpois(200, exp(0.3 + 0.5 * d$x + sin(2 * pi * region / K)))
fit <- tulpa(y ~ x + spatial(region), data = d, family = "poisson",
             spatial = list(type = "icar", adjacency = W), mode = "laplace")
fe <- field_effects(fit)
plot(fe$unit, fe$estimate, type = "b")

Description

Fits a GLM with a Matern spatial field via the SPDE approach. Uses CHOLMOD sparse solver with optional nested Laplace for hyperparameter integration.

Usage

fit_spde(
  y,
  X,
  spatial,
  family = "binomial",
  n_trials = NULL,
  range = NULL,
  sigma = NULL,
  nested_laplace = is.null(range) || is.null(sigma),
  phi = 1,
  offset = NULL,
  re_idx = NULL,
  n_re_groups = 0L,
  sigma_re = 1,
  mode = c("laplace", "nuts"),
  hyperprior = c("proper", "flat"),
  control = list()
)

Arguments

Value

A list with:

• mode: mode of the latent field (beta + mesh node effects)

• log_marginal: log marginal likelihood

• converged: convergence flag

• spatial: the spatial specification (for prediction)

• pareto_k, pareto_k_is_ess: outer Pareto-$\hat{k}$ and its importance-sampling ESS (NA when diagnose_k = FALSE)

• pareto_k_proposal_source, pareto_k_first_pass: which proposal family the reported k-hat came from (several candidates are scored and the best kept), and the k-hat of the first pass – the proposal exactly as placed, before refinement. A large gap says the placement is poor even where the verdict is fine.

• nested: nested Laplace results (if used)

References

Lindgren, Rue & Lindstrom (2011). An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. JRSS-B 73(4):423-498. Rue, Martino & Chopin (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. JRSS-B 71(2):319-392.

Examples

if (requireNamespace("fmesher", quietly = TRUE)) {
  set.seed(1)
  n <- 200L
  coords <- cbind(runif(n), runif(n))
  mesh <- fmesher::fm_mesh_2d(loc = coords, max.edge = c(0.15, 0.4), cutoff = 0.05)
  fem  <- fmesher::fm_fem(mesh)
  A    <- as(fmesher::fm_basis(mesh, loc = coords), "CsparseMatrix")
  spec <- spatial_spde_custom(C = fem$c0, G = fem$g1, A = A, nu = 1,
                              prior_range = c(0.3, 0.5), prior_sigma = c(0.6, 0.05))
  w <- as.numeric(rnorm(spec$n_mesh, 0, 0.6)); w <- w - mean(w)
  x <- rnorm(n)
  y <- rpois(n, exp(2.0 + 0.5 * x + as.numeric(spec$A %*% w)))
  fit <- fit_spde(y = y, X = cbind(1, x), spatial = spec, family = "poisson")
  fit$nested$range_mean
}

Description

Fits ⁠y ~ X beta + u_spatial[s] + v_temporal[t]⁠ with a spatial field – areal (icar / bym2 / car_proper) or continuous (hsgp / nngp) – and a temporal field (rw1 / rw2 / ar1), integrating the spatial hyperparameter(s), temporal precision, and (for ar1) the temporal autocorrelation over a hyperparameter grid via the ⁠cpp_nested_laplace_st_*⁠ kernels. The fixed-effect posterior is the grid-marginalised mixture; the spatial and temporal field posterior means are the grid-weighted latent modes.

Usage

fit_st_nested(
  y,
  X,
  spatial_idx,
  adjacency,
  temporal_idx,
  n_times,
  spatial_type = c("icar", "bym2", "car_proper", "hsgp", "nngp"),
  temporal_type = c("ar1", "rw1", "rw2"),
  family = "binomial",
  n_trials = NULL,
  phi = 1,
  cyclic = FALSE,
  re_idx = NULL,
  n_re_groups = 0L,
  sigma_re = 1,
  hyperprior = c("proper", "flat"),
  control = list(),
  coords = NULL,
  nn = 10L,
  cov_type = 2L,
  hsgp_m = 6L,
  hsgp_c = 1.5
)

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

Arguments

Value

A tulpa_fit (subclass tulpa_nested_laplace) carrying the fixed-effect posterior (draws via the grid mixture), spatial_effects, temporal_effects, log_marginal, weights, theta_grid over ⁠(tau_spatial, tau_temporal, rho)⁠ for an areal spatial field or ⁠(sigma2, lengthscale, tau_temporal, rho)⁠ / ⁠(sigma2, phi_gp, tau_temporal, rho)⁠ for "hsgp" / "nngp", and family / n_trials / phi (the R-level convention), which is what fitted(), residuals(), posterior_predict() and simulate() read the response family from. Also carries pareto_k_regime ("spread" / "collapsed_interior" / "collapsed_edge", see tulpa_nested_laplace_joint()'s return docs for the definition) and outer_grid_placement ("fixed" or "auto_recentered") plus, on a "fixed" placement, outer_grid_recenter_declined ("grid_knobs_overridden" / "grid_not_collapsed" / "no_usable_curvature" / "refit_failed" / "sd_ceiling_unresolved" / "sd_floor_unresolved"). A recentred fit also carries outer_grid_pinned_axes, the axes whose knobs were pinned and whose nodes were therefore kept, and outer_grid_recenter_sd_clamp / ⁠_sd_raw⁠ / ⁠_sd_used⁠ – per moved axis, which mode-SD bound the placement hit, the SD the stencil measured, and the SD the axis was laid from. A bound-decline is PER AXIS here: the axes the mode-find did resolve are still re-placed, and outer_grid_recenter_sd_declined names the ones that kept their incoming nodes and on which bound, so a partially re-placed grid is not read as a fully re-placed one. With every free axis declined the pass reports the grid as the fixed one it still is.

See Also

tulpa() (front door), tulpa_nested_laplace() (single field).

Examples

set.seed(1)
n_s <- 16L; n_t <- 8L; N <- 400L
adj <- matrix(0, n_s, n_s)
for (i in 1:(n_s - 1)) adj[i, i + 1] <- adj[i + 1, i] <- 1
s <- sample(n_s, N, TRUE); tt <- sample(n_t, N, TRUE)
us <- as.numeric(scale(cumsum(rnorm(n_s)))); vt <- as.numeric(scale(cumsum(rnorm(n_t))))
x <- rnorm(N)
y <- rbinom(N, 1, plogis(0.2 + 0.5 * x + 0.7 * us[s] + 0.6 * vt[tt]))
fit <- fit_st_nested(y, cbind(1, x), s, adj, tt, n_t, family = "binomial")

Description

In-sample mean response from the fixed effects and the observation offset (⁠E[y] = g^{-1}(X beta + offset)⁠, trial-scaled for binomial). Random effects are held at their prior mean of zero; group-level effects are in ranef(). y - fitted(object) equals residuals(object, type = "response"). On a zero-inflated fit (ziformula) the mean is the mixture's, ⁠(1 - pi) E[y | eta]⁠ with ⁠pi = plogis(X_zi beta_zi)⁠; over a zero-truncated family that is the hurdle mean.

Usage

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

Arguments

Value

Numeric vector of fitted mean responses, length nobs.

Description

The fixed-effect point estimates, equivalent to coef() on a tulpa_fit. Provided under the fixef name for code written against the lme4 / nlme interface: lme4::fixef(fit) and nlme::fixef(fit) dispatch here too when either package is installed, so a tulpa_fit can be dropped into an lme4-shaped workflow.

Usage

fixef(object, ...)

## S3 method for class 'tulpa_fit'
fixef(object, ...)

Arguments

Details

Note the deliberate difference from lme4::coef.merMod, which returns per-group sums of the fixed and random effects. On a tulpa_fit, coef() returns the fixed effects alone and fixef() is its synonym; the random effects are ranef().

Value

Named numeric vector of fixed-effect estimates.

See Also

coef.tulpa_fit(), ranef()

Examples

set.seed(1)
df <- data.frame(x = rnorm(80))
df$y <- rpois(80, exp(0.5 + 0.3 * df$x))
fit <- tulpa(y ~ x, data = df, family = "poisson")
fixef(fit)

Description

Performs Geweke's convergence diagnostic, comparing the mean of the first portion of a chain to the last portion. It is a within-chain diagnostic, so on a multi-chain fit it is computed for each chain separately.

Usage

geweke_test(fit, frac1 = 0.1, frac2 = 0.5, pars = NULL)

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

Arguments

Details

The Geweke test computes a z-score comparing the means of early and late portions of a chain. Large z-scores (|z| > 2) indicate the chain has not converged. The chains are never pooled: concatenated, the "first" window would come from the first chain and the "last" from another, which compares two chains rather than a chain's start with its end (use diagnostics() for the between-chain R-hat). A parameter that is constant within a chain has an undefined z-score, reported as NA.

Value

A data frame of class tulpa_geweke with columns chain, parameter, z_score and p_value, one row per chain and parameter.

See Also

diagnostics(), check_diagnostics()

Examples

set.seed(123)
df <- data.frame(x = rnorm(60))
df$y <- rpois(60, exp(0.5 + 0.3 * df$x))
fit <- tulpa(y ~ x, data = df, family = "poisson", mode = "hmc",
             control = list(n_iter = 500L, warmup = 250L, n_chains = 2L,
                            seed = 1L))
geweke_test(fit)

Description

Model-level summary statistics (broom-compatible)

Usage

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

Arguments

Value

Single-row data frame.

Examples

set.seed(1)
df <- data.frame(x = rnorm(100))
df$y <- rpois(100, exp(0.3 * df$x))
fit <- tulpa(y ~ x, data = df, family = "poisson")
glance(fit)

Description

Sample from a posterior using independence Metropolis-Hastings with a multivariate-normal proposal centred at the Laplace mode and scaled by the inverse Hessian. For posteriors that are well- approximated by a Gaussian near the mode this is dramatically cheaper than HMC: each iteration is one log-posterior evaluation plus one accept/reject step. Embarrassingly parallel across chains.

Use cases:

• Quick draws when Laplace is "almost right" but you want asymptotically-correct samples.

• Cross-validation, simulation studies, or other repeated-fit workflows where HMC startup cost dominates.

• Sanity check on the Laplace approximation: low IMH acceptance means the posterior is far from Gaussian and Laplace is biased.

Usage

imh_laplace(
  log_posterior,
  mode,
  hessian,
  n_iter = 2000L,
  warmup = n_iter%/%2L,
  scale = 1,
  init = NULL,
  thin = 1L,
  seed = NULL,
  verbose = FALSE
)

Arguments

Value

A list with class tulpa_fit carrying:

• draws: (n_iter - warmup) / thin x d matrix of draws.

• means: posterior means.

• accept_prob: per-iteration acceptance indicators.

• mean_accept: overall post-warmup acceptance rate.

• log_prob: per-draw log posterior values.

• inference_mode: "exact".

• inference_tier: 1L.

• backend: "imh_laplace".

Tier

Tier 1 (Exact). The MH accept/reject step makes the chain asymptotically correct under the standard MH conditions. Tier status does not depend on Laplace's quality – only its quality affects efficiency.

See Also

tulpa_laplace() for the mode + Hessian, bridge_sampling() for marginal-likelihood estimation on the resulting draws.

Examples

# Toy: Bernoulli logistic with one covariate.
set.seed(1)
n <- 100
x <- rnorm(n)
eta <- 0.3 + 1.2 * x
y <- rbinom(n, 1, plogis(eta))
X <- cbind(1, x)

lap <- tulpa_laplace(y, n_trials = rep(1L, n), X = X,
                     family = "binomial")

log_post <- function(beta) {
  eta <- as.numeric(X %*% beta)
  sum(y * eta - log1p(exp(eta))) +
    sum(dnorm(beta, 0, 10, log = TRUE))
}

fit <- imh_laplace(log_post, mode = lap$mode[1:2],
                   hessian = lap$H_beta, n_iter = 1000)
fit$mean_accept
colMeans(fit$draws)

Description

Displays information about the inference modes available in tulpa, including their tiers, guarantees, and appropriate use cases.

Usage

inference_mode_info()

Value

NULL, invisibly. Called for the side effect of printing the available inference modes, their tiers, guarantees, and use cases to the console.

Examples

inference_mode_info()

Description

Is an outer-grid setting marked as a default?

Usage

is_auto_grid(x)

Arguments

Value

TRUE when x carries the auto_grid() marker. This is the PROVENANCE question – whose choice the nodes are – and is TRUE whether or not the mark also asked for them to be integrated as written; that is auto_grid_place().

See Also

auto_grid(), auto_grid_place()

Examples

is_auto_grid(auto_grid(c(0.5, 1, 2)))
is_auto_grid(c(0.5, 1, 2))

Description

Wraps a user-defined latent block (currently a tgmrf() object) so the formula parser can recognise and route it to the inference layer. The call is structural – it is never executed at fit time. The parser evaluates the inner expression in the formula's environment, removes the latent(...) term from the fixed-effects formula, and attaches the resulting object to parsed$latent_blocks.

Usage

latent(block)

Arguments

Value

Returns block invisibly. Outside a formula context the call is a no-op pass-through.

Examples

n <- 8
blk <- tgmrf(
  Q = function(theta) Matrix::Diagonal(n, exp(-2 * theta[1])),
  prior = function(theta) stats::dnorm(theta[1], 0, 1, log = TRUE),
  init = c(log_sigma = 0))
parsed <- tulpa_parse_formula(y ~ x + latent(blk))
parsed$n_latent_blocks

Description

Define latent factors for capturing unmeasured shared structure between all processes. Factors are observation-level random effects that enter both linear predictors when shared = TRUE (default).

Usage

latent_factor(
  n_factors = 1L,
  prior = NULL,
  shared = TRUE,
  constraint = c("sum_to_zero", "first_zero"),
  scale = TRUE
)

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

Arguments

Details

Why Use Latent Factors?

When modeling multiple processes, they often share unmeasured confounders. For example, in relative abundance data, both the focal species count and total count might be affected by:

• Observer skill (unmeasured)

• Local microhabitat conditions (unmeasured)

• Weather on sampling day (unmeasured)

Without accounting for these shared drivers, estimates can be biased. Latent factors capture this shared structure without requiring the confounders to be measured.

Mathematical Model

For observation i with K latent factors, on each of two model arms (processes 1 and 2):

$\eta^{(1)}_i = X^{(1)}_i \beta^{(1)} + \sum_{k=1}^{K} f_{ik} \sigma_k + \ldots$

$\eta^{(2)}_i = X^{(2)}_i \beta^{(2)} + \sum_{k=1}^{K} f_{ik} \sigma_k + \ldots$

where:

• $f_{ik} \sim N(0, 1)$ are standardized factor scores

• $\sigma_k$ are factor standard deviations with PC prior

• Identifiability: $\sum_i f_{ik} = 0$ for each k

Because factors enter both linear predictors identically (when shared), they cancel in derived quantities (e.g., ratios, differences):

$\eta^{(1)}_i - \eta^{(2)}_i$

This means factors capture shared effects that would otherwise bias derived quantities.

Choosing n_factors

• Start with n_factors = 1 for simple unmeasured confounding

• Use n_factors = 2-3 if you suspect multiple independent confounders

• More than 3 factors is rarely needed and risks overfitting

• The PC prior provides regularization, shrinking unneeded factors toward zero

Relationship to Random Effects

Latent factors differ from random effects in several ways:

• Random effects are grouped (e.g., site-level), factors are observation-level

• Random effects require grouping structure, factors don't

• Factors capture residual correlation not explained by observed predictors

You can use both together: random effects for known grouping, factors for residual unmeasured confounding.

Value

A tulpa_latent object consumed by ratio / multi-arm model packages built on tulpa (e.g. tulpaRatio); not read by the single-response tulpa() front door (use latent(tgmrf(...)) for a single-response latent Gaussian block).

See Also

prior_pc() for prior specification; latent() / tgmrf() for a single-response in-formula latent Gaussian block

Examples

# Basic latent factor (single shared factor)
latent_factor()

# Two latent factors
latent_factor(n_factors = 2)

# Custom prior (more regularization)
latent_factor(n_factors = 1, prior = prior_pc(U = 0.5, alpha = 0.01))

# Numerator-only factor (not shared)
latent_factor(n_factors = 1, shared = FALSE)

# The spec is consumed by a ratio / multi-arm model package (tulpaRatio owns
# the two-arm `species | total ~ ...` formula and the negbin/negbin ratio
# family); it is that package's fitter, not the single-response tulpa() front
# door, that reads the `latent_factor()` object.

Description

Extract latent factor posteriors from fit

Usage

latent_factors(fit, ...)

## Default S3 method:
latent_factors(fit, summary = TRUE, probs = c(0.025, 0.5, 0.975), ...)

Arguments

Value

If summary = TRUE, a data frame with columns for observation, factor, and summary statistics. If summary = FALSE, a matrix of posterior draws.

Examples

# A fit sampled with latent factors comes from a model package (tulpaRatio
# reads the latent_factor() spec). The accessor reads two things off it:
# draw columns named factor[observation,factor] and the latent layout. Here
# both are simulated for three observations on two factors.
set.seed(1)
cols <- sprintf("factor[%d,%d]", rep(1:3, 2), rep(1:2, each = 3))
draws <- matrix(rnorm(200 * 6, mean = rep(c(-1, 0, 1), 2), sd = 0.3),
                200, 6, byrow = TRUE, dimnames = list(NULL, cols))
fit <- structure(
  list(draws = draws,
       .internal = list(latent_info = list(type = "latent", n_factors = 2L,
                                           n_obs = 3L))),
  class = "tulpa_fit")

latent_factors(fit)
dim(latent_factors(fit, summary = FALSE))

Description

What this returns depends on what the fit computed, and the returned object names it in a quantity attribute:

Usage

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

Arguments

Details

"log_posterior_mean"

a sampler fit: the mean over the draws of the log joint posterior density, prior included, in that backend's own parameterization: the unconstrained coordinates with their Jacobians for the ModelData samplers (hmc, ess, sghmc, sgld, mclmc, smc, vi), and the natural scale each quantity is sampled on for the RE-covariance and Polya-Gamma Gibbs samplers. Values are therefore comparable across fits of the same backend only.

"log_evidence"

a deterministic fit that estimated no hyperparameter from the data: the log marginal probability of the data under the model as specified. A hyperparameter the fit integrated counts, and so does one the caller supplied (a phi, a sigma_re, an outer grid laid at one value), which is part of the model rather than an estimate. A Laplace fit of a model with nothing to integrate reports its log marginal likelihood here, which already is that quantity.

"log_marginal_likelihood"

a deterministic fit that estimated some hyperparameters by maximising over the same data (empirical Bayes, or estimate_phi = TRUE): the log marginal likelihood at those estimates. The conditioned_on attribute names them.

"log_likelihood"

a fit that maximised over every parameter it reports, with the random effects integrated out (agq_fit()): the maximised log-likelihood, with df the number of maximised parameters.

On a nested-Laplace fit the value is the log evidence of its outer grid, ⁠log sum_k exp(log_marginal_k + log_cell_k)⁠, where log_marginal already carries each axis's hyperprior density on its integration coordinate and log_cell_k is the cell's absolute volume there. Every default outer axis carries a proper prior (a PC prior on each scale and range, a uniform on each bounded axis), so the value is the evidence under that prior and does not move with where the nodes were laid or how many there are, provided the grid covers the posterior. A fit carrying an axis with no proper prior (a flat random-effect covariance prior, a dispersion whose family has no PC prior yet, an axis on a block that does not carry its coordinates) reports NA with declined = "improper_hyperprior" and names the axes in a declined_axes attribute. A central-composite (CCD) design reproduces moments and carries no cell volume, so a fit integrated on one reports NA with a declined attribute.

A sampler fit whose producer kept no per-draw log posterior reports NA with quantity = "log_posterior_mean" and declined = "no_log_posterior_recorded". The Polya-Gamma Gibbs routes that recentre a proper field every sweep (the NNGP and multiscale-GP fields, and the negative-binomial kernel's iid random-effect block) leave no written density invariant and record none; a fit carrying neither draws nor a log marginal declines with "no_goodness_quantity_recorded". A value of NA always carries a declined attribute.

Values with a different quantity or conditioned_on are not on one scale, and compare_models(criterion = "loglik") refuses such a set. Only a "log_likelihood" is what AIC and BIC penalise, so AIC.tulpa_fit() and BIC.tulpa_fit() refuse every other quantity; compare those fits by their evidence, or by compare_models(criterion = "waic") / "loo", which score the pointwise predictive density.

Value

A logLik object with quantity and conditioned_on attributes, and a declined attribute naming the reason when no value could be read.

Description

Sample from a posterior using MALA: a Langevin proposal driven by the gradient of the log posterior, corrected by a Metropolis- Hastings accept/reject step. Each iteration is one log-posterior + one gradient evaluation. The natural stepping stone between random-walk MH and HMC: cheaper per iteration than HMC (no leapfrog integration), better-mixing than RWMH because the drift term moves toward higher density.

Step size epsilon is adapted via dual averaging during warmup to target an acceptance rate of 0.574 (Roberts & Rosenthal 1998 optimal for high-dimensional Gaussians).

Usage

mala(
  log_posterior,
  grad_log_posterior,
  init,
  n_iter = 2000L,
  warmup = n_iter%/%2L,
  epsilon = 0.1,
  target_accept = 0.574,
  mass_diag = NULL,
  mass_matrix = NULL,
  thin = 1L,
  seed = NULL,
  verbose = FALSE
)

Arguments

Value

A list with class tulpa_fit carrying:

• draws: post-warmup draws.

• means, n_params, n_samples, log_prob.

• accept_prob: per-iteration accepted indicators (post-warmup).

• mean_accept: post-warmup acceptance rate.

• epsilon: final adapted step size.

• inference_mode: "exact".

• inference_tier: 1L.

• backend: "mala".

Tier

Tier 1 (Exact). The MH step makes the chain asymptotically correct. Mixing depends on whether the gradient gives useful local geometry – poor for posteriors with very different scales across dimensions (use a preconditioner or HMC instead).

References

Roberts, G. O., & Tweedie, R. L. (1996). Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, 2(4), 341-363.

Roberts, G. O., & Rosenthal, J. S. (1998). Optimal scaling of discrete approximations to Langevin diffusions. JRSS B, 60(1), 255-268.

Examples

log_post <- function(t) -0.5 * sum((t - c(1, 2))^2)
grad <- function(t) -(t - c(1, 2))
fit <- mala(log_post, grad, init = c(0, 0), n_iter = 1000)
colMeans(fit$draws)  # near c(1, 2)
fit$mean_accept      # should adapt toward 0.574

Description

Returns a fit's posterior draws only when they form a genuine MCMC chain (⁠$draws_kind == "chain"⁠, or an untagged legacy fit); for an i.i.d. / approximation fit (nested Laplace, VI, SMC, ...) it returns NULL, because chain diagnostics do not apply. This is the accessor diagnostics() gates on. For the provenance-agnostic posterior sample used by summaries, see posterior_sample().

Usage

mcmc_draws(fit)

Arguments

Value

The chain draws matrix/array, or NULL for a non-chain fit.

See Also

posterior_sample(), diagnostics()

Examples

set.seed(1)
d <- data.frame(x = rnorm(60))
d$y <- rpois(60, exp(0.5 + 0.3 * d$x))
fit <- tulpa(y ~ x, data = d, family = "poisson", mode = "hmc",
             control = list(n_iter = 400L, warmup = 200L, n_chains = 2L,
                            seed = 1L))
dim(mcmc_draws(fit))

Description

Combine fitted values from several models using native model weights computed from the pointwise PSIS-LOO (or WAIC) elpd via tulpa_psis() – no loo package dependency. "loo" / "waic" give stacking weights (the simplex-optimal predictive combination); "pbma" / "pbma+" give pseudo-BMA(+) weights. Every model must carry an ⁠[n_draws x n_obs]⁠ pointwise log-likelihood (fit$draws$log_lik) over the same observations.

Usage

model_average(
  ...,
  weights = c("loo", "waic", "pbma", "pbma+"),
  fitted_fn = fitted
)

Arguments

Value

A list with averaged (the weighted fitted vector), weights (the named model weights), and comparison (the compare_models() table).

References

Yao, Vehtari, Simpson & Gelman (2018). Using stacking to average Bayesian predictive distributions. Bayesian Analysis 13(3):917-1007.

See Also

compare_models().

Examples

set.seed(1)
df <- data.frame(x = rnorm(120))
df$y <- rpois(120, exp(0.4 + 0.5 * df$x))
f1 <- tulpa(y ~ x, data = df, family = "poisson", mode = "hmc",
            control = list(n_iter = 500L, warmup = 250L, seed = 1L))
f2 <- tulpa(y ~ 1, data = df, family = "poisson", mode = "hmc",
            control = list(n_iter = 500L, warmup = 250L, seed = 1L))
ma <- model_average(full = f1, null = f2, weights = "waic")
ma$weights

Description

Tests whether residuals exhibit spatial structure after model fitting. Supports inverse-distance and k-nearest-neighbour weight matrices. No external dependencies (uses normal approximation for inference).

Usage

moran_i(
  object,
  coords,
  weights = c("inverse", "knn"),
  k = 10L,
  resid_type = "pearson",
  alternative = c("two.sided", "greater", "less")
)

Arguments

Value

An htest object with Moran's I, expected I, and p-value

Examples

set.seed(1)
coords <- cbind(runif(50), runif(50))
resid  <- rnorm(50)
moran_i(resid, coords)

Description

Counts divergent transitions recorded by an HMC/NUTS fit, reading whichever field the backend populated (the top-level ⁠$divergent⁠ flag vector every sampler in this package writes, or ⁠$diagnostics$divergent⁠, ⁠$diagnostics$divergent_idx⁠, ⁠$diagnostics$n_divergent⁠, ⁠$n_divergent⁠). It reads the same record plot_divergences() locates the divergent rows from, so the two always agree on one fit.

Usage

n_divergent(fit)

Arguments

Value

Integer count of divergent transitions (0 if none are recorded).

See Also

plot_divergences(), check_diagnostics()

Examples

set.seed(1)
d <- data.frame(x = rnorm(60))
d$y <- rpois(60, exp(0.5 + 0.3 * d$x))
fit <- tulpa(y ~ x, data = d, family = "poisson", mode = "hmc",
             control = list(n_iter = 400L, warmup = 200L, n_chains = 2L,
                            seed = 1L))
n_divergent(fit)

Description

Number of observations in a tulpa fit

Usage

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

Arguments

Value

Integer observation count.

Description

Translate a vector of original cell identifiers into the 1-based node indices of a tulpa_adjacency graph, by key (not by row order). Use it to add the node-index column the model's spatial grouping bar needs to the observation data, which typically has many rows per cell and a different row order than the graph.

Usage

node_index(graph, ids)

Arguments

Value

An integer vector the same length as ids, giving each one's node index in graph (NA for identifiers absent from the graph).

See Also

adjacency().

Examples

grid <- expand.grid(x = 1:3, y = 1:3)
grid$cell <- paste0("c", seq_len(nrow(grid)))
g <- adjacency(grid, id = "cell")

obs <- data.frame(cell = c("c5", "c1", "c9"))
obs$cell_idx <- node_index(g, obs$cell)
obs

Description

Single-path Pathfinder (Zhang, Carpenter, Gelman, Vehtari 2022): run L-BFGS toward the posterior mode, fit a Gaussian at the optimum using the inverse-Hessian estimate, and report draws plus the ELBO. Cheap, derivative-only, embarrassingly parallel – meant as an HMC warm-start, an initialiser for imh_laplace(), or a quick sanity check on the Laplace approximation.

This implementation is single-path only. The full multi-path Pathfinder (K parallel L-BFGS runs + mixture proposal + Pareto- smoothed importance reweighting) is a follow-on. The single-path version is what most users want as a Laplace-equivalent diagnostic.

Usage

pathfinder(
  log_posterior,
  init,
  grad_log_posterior = NULL,
  n_draws = 1000L,
  max_iter = 100L,
  tol = 1e-06,
  seed = NULL,
  verbose = FALSE
)

Arguments

Value

A list with class tulpa_fit carrying:

• draws: ⁠n_draws x d⁠ matrix of Gaussian draws at the mode.

• means: posterior means (= mode for a Gaussian fit).

• mode: the L-BFGS mode.

• cov: the proposal covariance (solve(-hessian)).

• elbo: Monte-Carlo estimate of ⁠E_q[log p - log q]⁠.

• n_iter: L-BFGS iterations.

• converged: logical.

• inference_mode: "structured".

• inference_tier: 2L.

• backend: "pathfinder".

Tier

Tier 2 (Structured). The output is a Gaussian approximation, not samples from the exact posterior – same epistemic class as tulpa_laplace(). Pair with imh_laplace() for an exact-tier upgrade.

References

Zhang, L., Carpenter, B., Gelman, A., & Vehtari, A. (2022). Pathfinder: parallel quasi-Newton variational inference. Journal of Machine Learning Research, 23(306), 1-49.

See Also

imh_laplace() for an exact-tier MH using the Pathfinder Gaussian as proposal; bridge_sampling() for marginal-likelihood estimation on the resulting draws.

Examples

# Toy: 2-D conjugate normal.
y <- c(0.5, -0.7)
log_post <- function(t) {
  sum(dnorm(y, t, 1, log = TRUE)) +
    sum(dnorm(t, 0, sqrt(10), log = TRUE))
}
pf <- pathfinder(log_post, init = c(0, 0), n_draws = 2000)
pf$mode      # near c(0.45, -0.64)
pf$elbo

Description

For each observation, computes the quantile of the observed value within the posterior predictive distribution. If the model is correct, PIT residuals follow Uniform(0, 1).

Usage

pit_residuals(object, ...)

## Default S3 method:
pit_residuals(object, observed = NULL, nsim = 250L, seed = 123L, ...)

Arguments

Details

For integer-valued responses, a randomisation step avoids discrete artefacts: the residual is drawn uniformly between P(sim < obs) and P(sim <= obs).

Value

Numeric vector of length n_obs with values in ⁠[0, 1]⁠

Examples

# From a matrix of replicates (n_obs x nsim) and the observed response:
set.seed(1)
y <- rpois(30, 3)
sims <- matrix(rpois(30 * 100, 3), nrow = 30)
r <- pit_residuals(sims, observed = y)
summary(r)

# From a fitted model with a simulate() method:
set.seed(1)
d <- data.frame(x = rnorm(80))
d$y <- rpois(80, exp(0.3 + 0.5 * d$x))
fit <- tulpa(y ~ x, data = d, family = "poisson", mode = "laplace")
head(pit_residuals(fit, nsim = 50))

Description

Creates autocorrelation function (ACF) plots for selected parameters. High autocorrelation indicates slow mixing and low effective sample size.

Usage

plot_acf(fit, pars = NULL, lags = 25, n_pars = 6)

Arguments

Details

Ideal ACF plots show rapid decay to zero. Slow decay indicates high autocorrelation, which reduces effective sample size and may indicate poor mixing.

Value

A ggplot object (if ggplot2 available) or base R plot (invisible).

See Also

plot_ess(), diagnostics()

Examples

set.seed(123)
df <- data.frame(x = rnorm(60))
df$y <- rpois(60, exp(0.5 + 0.3 * df$x))
fit <- tulpa(y ~ x, data = df, family = "poisson", mode = "hmc",
             control = list(n_iter = 500L, warmup = 250L, n_chains = 2L,
                            seed = 1L))
plot_acf(fit)
plot_acf(fit, lags = 10)

Description

Creates visualizations to investigate divergent transitions. Parallel coordinates and scatter plots highlight where in parameter space divergences occur.

Usage

plot_divergences(fit, pars = NULL, type = c("parcoord", "scatter"))

Arguments

Details

Divergent transitions indicate regions of high posterior curvature that the sampler cannot efficiently explore. Common causes:

• Very narrow funnels (hierarchical models)

• Strong correlations

• Multi-modality

Value

A ggplot object or base R plot (invisible).

See Also

plot_pairs(), n_divergent()

Examples

set.seed(123)
df <- data.frame(x = rnorm(60))
df$y <- rpois(60, exp(0.5 + 0.3 * df$x))
fit <- tulpa(y ~ x, data = df, family = "poisson", mode = "hmc",
             control = list(n_iter = 500L, warmup = 250L, n_chains = 2L,
                            seed = 1L))
plot_divergences(fit)
plot_divergences(fit, type = "scatter")

Description

Creates overlaid histograms of marginal energy and energy transition, along with the E-BFMI statistic. Low E-BFMI indicates poor exploration of the posterior.

Usage

plot_energy(fit)

Arguments

Details

E-BFMI (Energy Bayesian Fraction of Missing Information) compares the distribution of energy levels to energy transitions. Values below 0.3 indicate the sampler may not be exploring the full posterior.

Value

A ggplot object (if ggplot2 available) or base R plot.

See Also

diagnostic_summary(), check_diagnostics()

Examples

set.seed(123)
df <- data.frame(x = rnorm(60))
df$y <- rpois(60, exp(0.5 + 0.3 * df$x))
fit <- tulpa(y ~ x, data = df, family = "poisson", mode = "hmc",
             control = list(n_iter = 500L, warmup = 250L, n_chains = 2L,
                            seed = 1L))
plot_energy(fit)

Description

Creates a visual display of effective sample size (ESS) for all parameters, expressed as a ratio of ESS to total samples. Low ESS indicates high autocorrelation.

Usage

plot_ess(fit, type = c("bulk", "tail"), threshold = 400, pars = NULL)

Arguments

Details

ESS/iter ratio interpretation:

• Green: ESS >= threshold

• Yellow: threshold/2 <= ESS < threshold

• Red: ESS < threshold/2

Value

A ggplot object (if ggplot2 available) or base R plot (invisible).

See Also

plot_rhat(), diagnostic_summary(), diagnostics()

Examples

set.seed(123)
df <- data.frame(x = rnorm(60))
df$y <- rpois(60, exp(0.5 + 0.3 * df$x))
fit <- tulpa(y ~ x, data = df, family = "poisson", mode = "hmc",
             control = list(n_iter = 500L, warmup = 250L, n_chains = 2L,
                            seed = 1L))
plot_ess(fit)
plot_ess(fit, type = "tail")

Description

Create publication-ready maps from a tulpa spatial fit: the fitted response surface or the prediction uncertainty (credible-interval width), on the response scale.

Usage

plot_map(
  x,
  coords = NULL,
  what = c("fitted", "uncertainty"),
  summary = c("median", "mean", "q2.5", "q97.5", "sd"),
  newdata = NULL,
  title = NULL,
  palette = "viridis",
  points = FALSE,
  point_color = "grey30",
  point_size = 0.5,
  na_color = "transparent",
  crs = NULL,
  legend_title = NULL,
  ...
)