Package {tulpa}

Title:	Templated Unified Library for Posterior Approximation in Bayesian Hierarchical Models
Version:	0.7.5
Description:	A general-purpose engine for fitting Bayesian hierarchical models with spatial fields, temporal effects, spatially varying coefficients, and multiple inference backends. Scalable spatial structure includes Hilbert space approximate Gaussian processes (HSGP; Riutort-Mayol et al. 2023 <doi:10.1007/s11222-022-10167-2>), nearest-neighbor Gaussian processes (NNGP; Datta et al. 2016 <doi:10.1080/01621459.2015.1044091>), intrinsic conditional autoregressive models (ICAR; Besag, York, and Mollie 1991 <doi:10.1007/BF00116466>), the reparameterized Besag-York-Mollie model (BYM2; Riebler et al. 2016 <doi:10.1177/0962280216660421>), and stochastic partial differential equation fields (SPDE; Lindgren, Rue, and Lindstrom 2011 <doi:10.1111/j.1467-9868.2011.00777.x>). Temporal structure covers random walks, autoregressive processes, and Gaussian processes. Inference is tiered by correctness guarantee: exact Hamiltonian Monte Carlo with the No-U-Turn sampler, Laplace and nested Laplace approximations with hyperparameter integration (Rue, Martino, and Chopin 2009 <doi:10.1111/j.1467-9868.2008.00700.x>), and variational inference. Model-specific packages plug observation likelihoods into the engine through a templated C++ callback interface.
License:	MIT + file LICENSE
Encoding:	UTF-8
Language:	en-US
SystemRequirements:	C++17
Depends:	R (≥ 4.1.0)
Imports:	Rcpp (≥ 1.0.12), Matrix, tulpaMesh (≥ 0.1.3), generics, stats, tools, utils, methods, graphics, grDevices
Suggests:	testthat (≥ 3.0.0), bayesplot, coda, ggplot2, sf, terra, knitr, rmarkdown, posterior (≥ 1.5.0), loo (≥ 2.7.0), rstantools, lme4, nlme, lmtest, numDeriv, spdep, MASS, betareg, glmmTMB, pscl, tweedie, fmesher, patchwork, stars, statmod, withr
LinkingTo:	Rcpp, RcppEigen, Matrix
Config/testthat/edition:	3
VignetteBuilder:	knitr
URL:	https://github.com/gcol33/tulpa, https://gillescolling.com/tulpa/
BugReports:	https://github.com/gcol33/tulpa/issues
Config/roxygen2/version:	8.1.0
NeedsCompilation:	yes
Packaged:	2026-10-19 14:14:50 UTC; Gilles Colling
Author:	Gilles Colling [aut, cre, cph], Frances Y. Kuo [ctb, cph] (Sobol direction numbers in src/sobol_direction_numbers.h, BSD-3-clause), Stephen Joe [ctb, cph] (Sobol direction numbers in src/sobol_direction_numbers.h, BSD-3-clause), Arm Limited [cph] (exp and log routines in src/fastmath.h and src/fastmath_data.cpp, MIT), The Walnutpie Developers [cph] (WALNUTS transition in inst/include/tulpa/walnuts.h, ported from walnutpie, MIT)
Maintainer:	Gilles Colling <gilles.colling051@gmail.com>
Repository:	CRAN
Date/Publication:	2026-10-19 16:00:02 UTC

tulpa: Templated Unified Library for Posterior Approximation in Bayesian Hierarchical Models

Description

A general-purpose engine for fitting Bayesian hierarchical models with spatial fields, temporal effects, spatially varying coefficients, and multiple inference backends. Scalable spatial structure includes Hilbert space approximate Gaussian processes (HSGP; Riutort-Mayol et al. 2023 doi:10.1007/s11222-022-10167-2), nearest-neighbor Gaussian processes (NNGP; Datta et al. 2016 doi:10.1080/01621459.2015.1044091), intrinsic conditional autoregressive models (ICAR; Besag, York, and Mollie 1991 doi:10.1007/BF00116466), the reparameterized Besag-York-Mollie model (BYM2; Riebler et al. 2016 doi:10.1177/0962280216660421), and stochastic partial differential equation fields (SPDE; Lindgren, Rue, and Lindstrom 2011 doi:10.1111/j.1467-9868.2011.00777.x). Temporal structure covers random walks, autoregressive processes, and Gaussian processes. Inference is tiered by correctness guarantee: exact Hamiltonian Monte Carlo with the No-U-Turn sampler, Laplace and nested Laplace approximations with hyperparameter integration (Rue, Martino, and Chopin 2009 doi:10.1111/j.1467-9868.2008.00700.x), and variational inference. Model-specific packages plug observation likelihoods into the engine through a templated C++ callback interface.

Author(s)

Maintainer: Gilles Colling gilles.colling051@gmail.com (ORCID) [copyright holder]

Authors:

• Gilles Colling gilles.colling051@gmail.com (ORCID) [copyright holder]

Other contributors:

• Frances Y. Kuo (Sobol direction numbers in src/sobol_direction_numbers.h, BSD-3-clause) [contributor, copyright holder]

• Stephen Joe (Sobol direction numbers in src/sobol_direction_numbers.h, BSD-3-clause) [contributor, copyright holder]

• Arm Limited (exp and log routines in src/fastmath.h and src/fastmath_data.cpp, MIT) [copyright holder]

• The Walnutpie Developers (WALNUTS transition in inst/include/tulpa/walnuts.h, ported from walnutpie, MIT) [copyright holder]

See Also

Useful links:

• https://github.com/gcol33/tulpa

• https://gillescolling.com/tulpa/

• Report bugs at https://github.com/gcol33/tulpa/issues

Information criteria on a tulpa fit

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.

Random-effect variances and correlations

Description

Summarize the random-effect covariance of a fit: the standard deviation of each random-effect coefficient, the correlations between them within a term, and the covariance matrix itself.

Whether the covariance was estimated, sampled, or merely conditioned on is reported alongside it, because the three are different claims. A fit from mode = "laplace" conditions on sigma_re, so its values are the ones supplied; mode = "eb" estimates them; a sampler tier integrates them.

Usage

VarCorr(x, sigma = 1, ...)

## S3 method for class 'tulpa_fit'
VarCorr(x, sigma = 1, ...)

Arguments

Value

A data frame with one row per random-effect coefficient: term, coef, sd, and source (one of "estimated", "sampled", "conditioned"). A gaussian fit adds a final term = "Residual" row (coef NA) with the residual SD sqrt(phi), labelled "estimated" under estimate_phi = TRUE and "conditioned" otherwise; it has no entry in "cov". Correlated terms additionally carry the covariance matrices in the "cov" attribute, one per term, each with a "correlation" attribute. Returns an empty data frame when the fit has no random effects.

See Also

ranef() for the per-level deviations, tulpa_eb() to estimate the covariance rather than condition on it.

Examples

set.seed(1)
G <- 30L; per <- 10L; n <- G * per
grp <- rep(seq_len(G), each = per); x <- rnorm(n)
b <- rnorm(G, 0, 0.8)
d <- data.frame(y = rpois(n, exp(0.3 + 0.5 * x + b[grp])), x = x,
                g = factor(grp))
fit <- tulpa(y ~ x + (1 | g), data = d, family = "poisson", mode = "eb")
VarCorr(fit)

Construct a spatial adjacency graph for areal models

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

Adaptive Gauss-Hermite quadrature for one-RE GLMMs

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)

Posterior draws in the posterior package's format

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))
}

Mark an outer-grid setting as a default rather than a pin

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))

May the placement pass move a marked outer-grid axis?

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))

Bayesian R-squared

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)

Bridge sampling for marginal likelihood

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)

Observation-level accessors for categorical fits

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()

Central composite design (CCD) grid for nested-Laplace integration

Description

Produces a structured set of standardised hyperparameter points z \in \mathbb{R}^k for use as integration nodes in a nested Laplace approximation when k \ge 3. CCD scales much better than the full tensor expand.grid() used by the 1D and 2D backends: a CCD has 1 + 2k + 2^{k - q} points (1 centre, 2k axial, 2^{k - q} factorial), versus m^k for an m-per-axis tensor product.

Point layout (with centre+axial+factorial scaling f_0):

• 1 centre point at the origin;

• 2k axial points at \pm f_0 along each coordinate axis;

• 2^{k - q} factorial points at corners of the hypercube, scaled to lie on a sphere of radius f_0.

For k \le 6 the factorial portion is the full 2^k design. For k \ge 7 a half-fraction (q = 1) using the defining word x_1 \cdots x_k keeps the point count reasonable while preserving Resolution V.

Used by tulpa_nested_laplace() for higher-dimensional hyperparameter blocks. The standardised z-coordinates are mapped to physical hyperparameters \theta via ccd_to_theta().

Usage

ccd_grid(k, f_0 = sqrt(k))

Arguments

Value

A list with components:

• z: numeric matrix ⁠[n_points x k]⁠ of standardised hyperparameter coordinates.

• n_points: integer; total grid size.

• kind: character vector labelling each point as "center", "axial", or "factorial".

• f_0: the sphere radius used.

See Also

ccd_to_theta() to map z-coordinates to physical theta.

Examples

g <- ccd_grid(3)
g$n_points
table(g$kind)

Map standardised CCD coordinates to physical hyperparameters

Description

Converts standardised CCD z-coordinates (in \mathbb{R}^k) to physical hyperparameters \theta via the affine map \theta = \hat\theta + L \cdot z with optional log-scale transform per component. L is typically a Cholesky factor of the negative Hessian inverse evaluated at the (working) mode \hat\theta – i.e. it scales z to one posterior standard deviation per axis.

Usage

ccd_to_theta(z, theta_hat, L, log_scale = FALSE)

Arguments

Value

Numeric matrix ⁠[n_points x k]⁠ of physical theta-values.

Examples

g <- ccd_grid(2)
ccd_to_theta(g$z, theta_hat = c(0, 1), L = diag(c(0.5, 0.2)))

Corrected R-INLA CCD integration weights

Description

Per-point design weights for a ccd_grid() used as nested-Laplace integration nodes. Implements the corrected R-INLA central-composite-design weights: the formula in Rue, Martino and Chopin (2009) has a typo, corrected by the R-INLA team to give a positive central weight. With m hyperparameters, np design points and standardized scaling f0, every non-centre point gets

w = 1 / ((np - 1)(1 + e^{-m f0^2 / 2}(f0^2 - 1)))

and the centre point gets w_0 = 1 - (np - 1) w. Used as \Delta_k: the integration weight of node k is \Delta_k \exp(\log\,\mathrm{marginal}_k), renormalized. On a standardized (whitened) hyperparameter posterior these reproduce the Gaussian moments exactly.

Usage

ccd_weights(ccd)

Arguments

Value

Numeric vector of length ccd$n_points: the design weight for each point (the centre weight w0 for the central point, w otherwise).

See Also

ccd_grid(), ccd_to_theta().

Examples

ccd_weights(ccd_grid(3))

Validate a spatial adjacency matrix

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)

Quick convergence check

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)

Diagnostic panel plot

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)

Fixed-effect coefficients

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().

Compare models by information criteria

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")

Compute nearest neighbors for NNGP

Description

Compute the k nearest neighbors for each observation using Euclidean distance. Returns in a format suitable for the NNGP likelihood.

Usage

compute_nngp_neighbors(coords, k)

Arguments

Value

List with:

• nn_idx: N x k matrix of neighbor indices (0 for obs with fewer neighbors)

• nn_dist: N x k matrix of distances to neighbors

• nn_order: Ordering of observations for NNGP (by coordinate)

Credible intervals for the fixed effects

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().

DIC, CPO, WAIC and PSIS-LOO on a 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)

Comprehensive diagnostic summary

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)

Posterior diagnostics for a fitted model

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)

Durbin-Watson test for temporal autocorrelation

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)

Energy Bayesian Fraction of Missing Information (E-BFMI)

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

Supported R-level family names

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 accepted as well and not listed here.

Usage

family_names()

Value

A character vector of family names.

Examples

family_names()

Spatial field effects per unit

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]); for BYM2 each draw's two unit-scale components are combined into the same sigma * (sqrt(rho) * s * phi + sqrt(1 - rho) * theta) the sampler path reports. 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")

Find all bar terms in a formula's parse tree

Description

Recursively walks the formula AST and collects all | and || nodes found inside parentheses. These are the random effect specifications.

Usage

findbars(term)

Arguments

Value

A list of language objects, each a | or || call

Examples

findbars(y ~ x + (1 | g) + (1 + x | h))

Fit a spatial model using SPDE Laplace approximation

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
}

Fit an additive spatiotemporal GLM by nested Laplace

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")

Fitted values (population level)

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.

Fixed-effect coefficients (lme4-compatible)

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)

Geweke convergence test

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)

Model-level summary statistics (broom-compatible)

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)

Describe one outer-grid hyperparameter axis

Description

Builds a single axis spec for tulpa_hyper_grid(). Carries the candidate values, the optional log-prior, and the metadata (log-scale, bounds, refinable flag) that the generic refinement / consistency passes need.

Usage

hyper_axis_spec(
  name,
  grid,
  log_prior = NULL,
  log_scale = FALSE,
  bounds = NULL,
  refinable = FALSE,
  atom_mass = NULL,
  slab_bounds = NULL,
  log_prior_coord = c("integration", "natural"),
  extend = TRUE
)

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

Arguments

Value

An object of class tulpa_hyper_axis_spec (a validated list with the fields listed above).

See Also

tulpa_hyper_grid().

Examples

hyper_axis_spec("sigma", grid = c(0.5, 1, 2), log_scale = TRUE,
                bounds = c(0, Inf))
# A zero level on a log-scale axis is a point mass with a declared share:
hyper_axis_spec("alpha", grid = c(0, 0.5, 1), log_scale = TRUE,
                bounds = c(0, Inf), atom_mass = 0.5)

Independence Metropolis-Hastings with Laplace proposal

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)

Print inference mode information

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()

Is an outer-grid setting marked as a default?

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))

Mark an expression as a latent block in a tulpa formula

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

Create a latent factor specification

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.

Extract latent factor posteriors from fit

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))

The fit's log-scale goodness quantity

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.

Metropolis-adjusted Langevin algorithm (MALA)

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

MCMC chain draws from a fit

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))

Model-averaged predictions

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

Moran's I test for spatial autocorrelation in residuals

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)

Number of divergent transitions

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)

Remove all bar terms from a formula's parse tree

Description

Recursively rewrites the formula AST, removing any | or || nodes found inside parentheses. Returns the fixed-effects-only formula.

Usage

nobars(term)

Arguments

Value

A language object with all bar terms removed, or NULL if nothing remains

Examples

nobars(y ~ x + (1 | g) + (1 + x | h))

Number of observations in a tulpa fit

Description

Number of observations in a tulpa fit

Usage

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

Arguments

Value

Integer observation count.

Map cell identifiers to graph node indices

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

Pathfinder: variational warm-start via L-BFGS + ELBO scoring

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

PIT (probability integral transform) residuals

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))

Plot fixed-effect posteriors

Description

Plot fixed-effect posteriors

Usage

## S3 method for class 'tulpa_fit'
plot(x, type = c("density", "trace", "pairs", "smooth"), term = NULL, ...)

Arguments

Value

The input x, returned invisibly. Called for the side effect of producing base-graphics plots of the fixed-effect posteriors.

Plot method for tulpa_prior_predict

Description

Density overlay of prior predictive draws. Requires bayesplot; falls back to base graphics matplot of a subset of draws otherwise.

Usage

## S3 method for class 'tulpa_prior_predict'
plot(x, process = 1L, max_draws = 50L, ...)

Arguments

Value

A ggplot object (via bayesplot) when bayesplot is installed; otherwise NULL invisibly, after drawing a base-graphics overlay. Called for the density overlay of the prior predictive draws.

Plot method for spatiotemporal effects

Description

Plot method for spatiotemporal effects

Usage

## S3 method for class 'tulpa_st_summary'
plot(x, type = "heatmap", ...)

Arguments

Value

A ggplot object when ggplot2 is installed; otherwise NULL invisibly, after drawing a base-graphics plot. Called for the side effect of visualizing the spatiotemporal interaction effects.

Plot method for tulpa_svc_posterior

Description

Plot method for tulpa_svc_posterior

Usage

## S3 method for class 'tulpa_svc_posterior'
plot(x, term = 1, type = "mean", ...)

Arguments

Value

A ggplot object when ggplot2 is installed; otherwise NULL invisibly, after drawing a base-graphics map. Called for the side effect of mapping the selected spatially-varying coefficient.

Plot method for tulpa_temporal_posterior

Description

Plot method for tulpa_temporal_posterior

Usage

## S3 method for class 'tulpa_temporal_posterior'
plot(x, component = NULL, type = "ribbon", ...)

Arguments

Value

A ggplot object when ggplot2 is installed; otherwise NULL invisibly, after drawing a base-graphics plot. Called for the side effect of plotting the temporal-effect posterior.

Plot method for tulpa_tvc_posterior

Description

Plot method for tulpa_tvc_posterior

Usage

## S3 method for class 'tulpa_tvc_posterior'
plot(x, term = 1, type = "ribbon", ...)

Arguments

Value

A ggplot object when ggplot2 is installed; otherwise NULL invisibly, after drawing a base-graphics plot. Called for the side effect of plotting the selected temporally-varying coefficient.

Plot autocorrelation functions

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)

Diagnostic plotting functions for tulpa models

Description

Visual diagnostic tools for MCMC convergence assessment. All functions provide base R fallbacks when ggplot2/bayesplot are unavailable.

Creates a combined diagnostic figure with Rhat, ESS, trace plot, and energy/ACF panels. Requires the patchwork package for layout.

Usage

plot_diagnostics(fit, pars = NULL)

Arguments

Details

Creates a 2x2 grid:

• Top left: Rhat plot

• Top right: ESS plot

• Bottom left: Trace for worst parameter

• Bottom right: Energy (HMC) or ACF (other)

Value

A combined plot (ggplot + patchwork) or NULL if requirements not met.

See Also

diagnostic_summary(), plot_rhat(), plot_ess()

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_diagnostics(fit)

Plot divergent transitions

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")

Plot energy diagnostic (E-BFMI)

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)

Plot effective sample size diagnostic

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")

Plot spatial predictions as a map

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,
  ...
)

Arguments

Details

This function provides a streamlined workflow for visualizing spatial predictions from tulpa models. It handles:

• Extracting predictions from tulpa_fit objects

• Converting to appropriate spatial format (stars/sf)

• Creating publication-quality maps with sensible defaults

• Uncertainty visualization via credible interval width

For custom maps or more control, extract predictions using predict() and use ggplot2 directly with geom_stars() or geom_sf().

Value

A ggplot2 object that can be further customized.

Required packages

This function requires ggplot2. For raster-style maps, stars and sf are also needed. Install with:

install.packages(c("ggplot2", "stars", "sf"))

See Also

predict.tulpa_fit() for predictions at new locations

Examples

if (requireNamespace("ggplot2", quietly = TRUE)) {
  set.seed(123)
  n_sites <- 20
  df <- data.frame(
    y = rbinom(n_sites, 20, 0.4),
    elevation = rnorm(n_sites),
    site = factor(seq_len(n_sites)),
    lon = runif(n_sites),
    lat = runif(n_sites)
  )
  adj <- matrix(0, n_sites, n_sites)
  for (i in 1:(n_sites - 1)) adj[i, i + 1] <- adj[i + 1, i] <- 1
  fit <- tulpa(
    y ~ elevation + spatial(site),
    data = df,
    family = "binomial",
    n_trials = rep(20L, n_sites),
    spatial = spatial_car(adj, group_var = "site"),
    mode = "laplace"
  )
  # Areal (CAR/ICAR) fits carry no point coordinates, so pass `coords`:
  cc <- df[, c("lon", "lat")]
  plot_map(fit, coords = cc)                      # fitted response surface
  plot_map(fit, what = "uncertainty", coords = cc)
  plot_map_panel(fit, coords = cc)                # both side by side
}

Plot multiple maps in a grid

Description

Create a multi-panel figure with the fitted-value map and the prediction-uncertainty map side by side.

Usage

plot_map_panel(x, newdata = NULL, ncol = 2, ...)

Arguments

Value

A patchwork object (if patchwork is installed) or a list of ggplots.

Examples

if (requireNamespace("ggplot2", quietly = TRUE)) {
  set.seed(123)
  n_sites <- 20
  df <- data.frame(
    y = rbinom(n_sites, 20, 0.4),
    elevation = rnorm(n_sites),
    site = factor(seq_len(n_sites)),
    lon = runif(n_sites),
    lat = runif(n_sites)
  )
  adj <- matrix(0, n_sites, n_sites)
  for (i in 1:(n_sites - 1)) adj[i, i + 1] <- adj[i + 1, i] <- 1
  fit <- tulpa(
    y ~ elevation + spatial(site),
    data = df,
    family = "binomial",
    n_trials = rep(20L, n_sites),
    spatial = spatial_car(adj, group_var = "site"),
    mode = "laplace"
  )
  cc <- df[, c("lon", "lat")]
  plot_map_panel(fit, coords = cc)
  plot_map_panel(fit, coords = cc, ncol = 1)
}

Plot bivariate parameter posteriors (pairs plot)

Description

Creates a pairs plot showing bivariate relationships between parameters. Divergent transitions (if present) are highlighted to help identify problematic posterior regions.

Usage

plot_pairs(
  fit,
  pars = NULL,
  highlight_divergent = TRUE,
  n_pars = 5,
  alpha = 0.3
)

Arguments

Details

Pairs plots help identify:

• Strong correlations between parameters (potential non-identifiability)

• Multimodality

• Regions where divergences cluster (indicating problematic geometry)

Value

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

See Also

plot_divergences(), 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_pairs(fit)
plot_pairs(fit, n_pars = 2)

Plot Rhat convergence diagnostic

Description

Creates a visual display of Rhat values for all parameters, with color-coding to highlight convergence issues. Rhat values > 1.01 indicate potential convergence problems.

Usage

plot_rhat(fit, threshold = 1.01, pars = NULL)

Arguments

Details

Color coding:

• Green: Rhat < 1.01 (converged)

• Yellow: 1.01 <= Rhat < 1.05 (borderline)

• Red: Rhat >= 1.05 (not converged)

Value

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

See Also

plot_ess(), diagnostic_summary(), diagnostics()

Examples

# Diagnostic plots require a fitted tulpa model
# See tulpa() examples for fitting models


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_rhat(fit)

Fit a post-hoc linear model on estimated parameters

Description

Useful for exploring drivers of occupancy/detection/abundance variation after model fitting. Fits a weighted linear model and optionally generates bootstrap confidence intervals.

Usage

post_hoc_lm(
  formula,
  data,
  weights = NULL,
  n_boot = 1000L,
  probs = c(0.025, 0.975)
)

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

Arguments

Value

A list of class "post_hoc_lm" with:

summary

data.frame of coefficient estimates and CIs

lm_fit

the underlying lm object

boot_coefs

matrix of bootstrap coefficient samples (if n_boot > 0)

R2

R-squared from the fitted model

Examples

# Explore drivers of per-site estimates after fitting a model.
site <- data.frame(
  psi_hat = c(0.2, 0.5, 0.8, 0.4, 0.6, 0.3),
  se      = c(0.05, 0.04, 0.06, 0.05, 0.03, 0.05),
  trait   = c(1.0, 2.5, 3.8, 1.9, 3.1, 1.2)
)
fit <- post_hoc_lm(psi_hat ~ trait, data = site,
                   weights = 1 / site$se^2, n_boot = 200L)
fit

Posterior predictive replicates

Description

Draw replicated responses from the posterior predictive distribution: the in-sample linear predictor is drawn with every component the fit estimated – fixed effects, formula random effects, the offset, and any spatial or temporal field – and pushed through the family's sampling distribution. The same draws give the pointwise log-likelihood behind compare_models(criterion = "waic") / "loo".

Where the draws come from follows what the fit carries:

• A ModelData sampler fit (mode = "hmc" and its siblings) evaluates the engine's own linear predictor at each draw, so a field is drawn jointly with everything else.

• A nested-Laplace fit draws from its outer grid: a replicate picks a cell by its weight and draws each observation from that cell's Gaussian for the linear predictor (fitted_eta, fitted_eta_var). The cell's joint covariance across observations is not retained, so observations within a replicate are independent given the cell. A fit run with control$fitted_var = FALSE carries no fitted_eta_var, and its replicates hold the across-cell spread only.

• Any other fit carries its coefficients rather than its linear predictor. Fits with posterior draws use them (fixed and random effects jointly per draw). A Laplace or EB fit (mode = "laplace" / "eb") samples the fixed effects from N(coef(fit), vcov(fit)) and each draw's random effects from their Gaussian conditional given those fixed effects under the joint Laplace precision of ⁠[beta | b]⁠ at the mode, so the two carry their posterior correlation (on a conditional Laplace fit this is the joint Gaussian draw; an EB fit is conditional on the estimated variance components). A conditional-Laplace spatial fit holds its random effects and areal field at their posterior mode, which over-disperses the linear predictor wherever the intercept is aliased with them; an SPDE field enters at its posterior mean.

At newdata the prediction is population level (random effects at zero), matching predict.tulpa_fit().

A zero-inflated fit (ziformula) draws the structural-zero logit from the same posterior draw as the count predictor; each replicate observation is a structural zero with probability ⁠plogis(X_zi beta_zi)⁠ and otherwise a draw from the family (a zero-truncated family gives the hurdle model).

Usage

posterior_predict(object, ...)

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

Arguments

Value

A ⁠ndraws x n_obs⁠ numeric matrix of replicated responses.

See Also

pp_check(), which uses these replicates; simulate.tulpa_fit().

Examples

set.seed(1)
d <- data.frame(y = rpois(100, 4), x = rnorm(100))
fit <- tulpa(y ~ x, data = d, family = "poisson", mode = "laplace")
yrep <- posterior_predict(fit, ndraws = 100)
dim(yrep)  # 100 x 100

Posterior parameter sample from a fit

Description

Returns a fit's posterior draws for summary purposes (quantiles, derived quantities, density plots), regardless of how they were produced – an MCMC chain, a nested-Laplace node mixture, or a variational sample all answer here. For the chain-only view used by convergence diagnostics, see mcmc_draws().

Usage

posterior_sample(fit)

Arguments

Details

A fit that carries no draws returns NULL with a message naming its backend, the posterior representation it carries instead, and the accessor that turns that representation into draws where one exists – an absent draws matrix is a property of the backend, not a failure of the accessor, and saying which is the difference between a diagnosable answer and a bare NULL.

Value

The posterior draws matrix/array, or NULL if the fit carries none.

See Also

mcmc_draws(), tulpa_posterior_draws(), diagnostics()

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")
dim(posterior_sample(fit))

Posterior predictive check

Description

Generate posterior predictive checks for a fitted tulpa model. Compares observed data to replicated data from the posterior predictive distribution.

Usage

pp_check(object, ...)

## S3 method for class 'tulpa_fit'
pp_check(
  object,
  type = c("dens_overlay", "scatter", "intervals", "stat"),
  component = NULL,
  stat = mean,
  ndraws = 50,
  ...
)

Arguments

Value

A ggplot object

Examples

# pp_check is a generic; model packages (e.g. tulpaObs, tulpaRatio) provide
# the posterior-predictive method for their fits.

set.seed(123)
n <- 200L
df <- data.frame(y = rpois(n, 5), x = rnorm(n),
                 site = factor(rep(1:10, each = 20)))
fit <- tulpa(y ~ x + (1 | site), data = df, family = "poisson",
             mode = "hmc", control = list(n_iter = 500, warmup = 250))
# Density overlay (dispatches to the model package's pp_check method).
if (requireNamespace("bayesplot", quietly = TRUE)) {
  pp_check(fit, type = "dens_overlay")
}

Predict at new covariate values (population level)

Description

Prediction of the linear predictor at newdata, on the link or response scale. The fixed-effect part is ⁠X beta⁠ with credible bounds from the fixed-effect covariance (vcov()). For a fit carrying a continuous spatial field, the posterior-mean field is interpolated (kriged) to the newdata coordinates and added to the linear predictor by default, so predict() gives the conditional (location-specific) prediction. Three continuous field families are supported: an SPDE Matern field (spatial_spde()), projected through the mesh; a Hilbert-space GP field (spatial_gp(approx = "hsgp")), where the Laplacian basis is re-evaluated at the new coordinates (with the training centring / boundary); and a GP / NNGP field (spatial_gp()), interpolated by the NNGP conditional mean at each new location's nearest training locations. The HSGP and GP/NNGP fields are marginalised over the hyperparameter grid (not plugged in at the posterior mean). Ordinary random effects are held at zero (population level); add group effects from ranef() when needed. An areal (ICAR / BYM2 / CAR) or temporal (RW1 / RW2 / AR1) field is held at zero in the same way, at the training design too; the in-sample linear predictor with every fitted component is what posterior_predict() draws and what compare_models(criterion = "waic") / "loo" score.

Usage

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

Arguments

Value

If se.fit = FALSE, a numeric vector. If se.fit = TRUE, a data frame with fit, se.fit (link scale), lower, upper on the requested scale.

Summarize and print a tulpa fit

Description

print() shows a fit in a few lines: how it was produced (backend and inference tier), the call, the observation count and the fixed-effect point estimates. summary() returns the full report as a summary.tulpa_fit object, which prints as a sectioned display and whose parts are read by name; coef() on it returns the fixed-effect table.

Usage

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

## S3 method for class 'tulpa_fit'
summary(object, level = 0.95, ...)

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

## S3 method for class 'summary.tulpa_fit'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

Arguments

Details

The fixed-effect table reads the same posterior coef(), confint() and vcov() read. Sampler tiers report empirical quantiles; the Laplace tier reports the Gaussian approximation, and so does any fit whose reported posterior is a closed-form Gaussian (expectation propagation, the multinomial Laplace fit), whose draws are then samples from that Gaussian and are not what the summary reads.

On a nested-Laplace fit the estimate and standard error are the hyperparameter-grid-marginalized moments, and the bounds invert the Gaussian mixture ⁠sum_k w_k N(mu_kj, V_kjj)⁠ that grid defines, rather than reading ⁠mu +/- z sigma⁠ off the single Gaussian matching those moments. The table's interval_source attribute records which read produced them ("mixture_cdf", "gaussian_moment", "skew_map_cell", or "skew_map_cell/mixture_cdf" when the two are in play on different coefficients) and interval_declined says why, whenever the mixture read did not run. Its retained_mass attribute gives the share of the grid weight whose cells retained a fixed-effect block: 1 on a complete grid, and below 1 on one that dropped a positive-weight cell, whose report is then the posterior conditional on the cells that remain.

Where no per-cell block was retained the estimate falls back to the grid-weighted average of the per-cell modes, restricted to the cells whose inner solve reached a mode. A fit where none did reports NA with interval_declined = "not_converged", rather than the vector its Newton started from as an estimate.

With control$skew_correct = TRUE a coefficient whose inner-Laplace gamma_3 is in the band it is valid on reports Cornish-Fisher quantiles instead, and the table's skew_applied attribute names which coefficients took the correction. That correction is measured at the MAP cell, so it is reported on its own and is not composed with the mixture read; a coefficient it declines keeps the mixture read rather than falling back further.

The reliability headline is read from diagnostics(), so it answers the question that applies to how the draws were produced: Rhat and effective sample size for a chain, the approximation-reliability band (outer Pareto-k, grid quadrature, inner Laplace, combined verdict) for i.i.d. draws from a deterministic approximation, and nothing for a point summary.

Value

print.tulpa_fit() and print.summary.tulpa_fit() return their input invisibly.

summary.tulpa_fit() returns a list of class summary.tulpa_fit:

call, formula, family

as the fit records them, NULL where it records none.

inference

how the fit was produced: a list of backend, inference_tier, inference_mode and selection_reason (why mode = "auto" chose that backend).

n_obs, n_draws

the observation count and the number of posterior draws the fixed-effect summaries read (NULL without draws).

level

the interval level of coefficients.

coefficients

data frame, one row per fixed effect named by the term: estimate, std.error and the lower / upper credible bounds, labelled as confint.tulpa_fit() labels them. Carries the interval_source, interval_declined, retained_mass and skew_applied attributes described in Details. Also returned by coef() on the summary.

random

the random-effect standard deviations as VarCorr() reports them, with their estimated / sampled / conditioned label; NULL on a fit with no random-effect terms.

hyperparameters

data frame of the integrated hyperparameters, one row each: parameter, mean, sd, median and the 2.5% / 97.5% bounds ci_lo / ci_hi. On a nested-Laplace fit these are the outer-grid axes; on the random-effect covariance integrators the derived standard deviations and correlations. NULL where the fit integrated none.

sigma

the residual / dispersion scale a fit reports as one number, else NULL.

smooth_terms

the fit's s(...) term metadata, NULL without smoothers; smooth_effects() returns the fitted curves.

reliability

a list: kind ("chain", "iid" or "point", the draws provenance diagnostics() routes on); on a chain rhat_max, ess_bulk_min and ess_tail_min over the parameters; on i.i.d. draws band, the one-row headline diagnostics() attaches as its summary attribute (outer pareto_k and its band, grid quadrature ESS, inner-Laplace scores, reliability), and moment_source, which read the per-parameter moments came from; verdict, the whole-fit reading ("converged" / "not converged" on a chain that records its convergence check, the combined reliability verdict on i.i.d. draws, NA otherwise); convergence, the fit's own convergence record; and the backend-specific quality numbers mean_accept (Metropolis acceptance rate), n_divergent (divergent transitions) and elbo (the variational objective), each NULL where the backend has none.

axis_fields_dropped

the grid axes the fit's own resolved path could not read, one row per dropped field (block, type, field, path, integrates, reason); NULL whenever every supplied axis was used. diagnostic_summary() reads the same record in sentences.

prior

a list whose beta is the Gaussian fixed-effect prior the fit ran under, as list(mean, sd): the engine default prior_normal(0, 2.5) whenever the caller supplied none, and NULL on the paths that express no Gaussian prior on the fixed effects.

coef.summary.tulpa_fit() returns the coefficients data frame.

See Also

coef.tulpa_fit(), confint.tulpa_fit(), VarCorr(), diagnostics(), tulpa_outer_grid()

Examples

set.seed(1)
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 = "laplace")
fit
s <- summary(fit)
s
coef(s)
s$inference$backend

Nested Laplace approximation for latent Gaussian models

Description

Generic outer-grid nested Laplace driver. Builds a grid over the hyperparameters of a single latent prior block (spatial or temporal), runs an inner Laplace at each grid point with warm-starting, and integrates over the grid to give proper hyperparameter marginals.

Supported priors:

• Spatial (areal): "icar" (1D grid on tau), "bym2" (2D on (sigma, rho)), "car_proper" (2D on (tau, rho); rho lives in the eigenvalue interval (1/lambda_min, 1/lambda_max) of ⁠D^{-1} W⁠).

• Spatial (continuous): "nngp" (2D on (sigma2, phi_gp)), "hsgp" (2D on (sigma2, lengthscale)).

• Temporal: "rw1", "rw2" (1D grid on tau), "ar1" (2D on (tau, rho))

• SPDE continuous spatial: see cpp_nested_laplace_spde() (separate entry, rebuilds Q via SPDE Q-builder).

Usage

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

tulpa_nested_laplace(
  y,
  n_trials,
  X,
  prior = NULL,
  spec = NULL,
  data = NULL,
  re_idx = NULL,
  n_re_groups = 0L,
  sigma_re = 1,
  family = "binomial",
  phi = 1,
  offset = NULL,
  likelihood = NULL,
  hyperprior = c("proper", "flat"),
  control = list()
)

Arguments

Details

print() shows the integrated hyperparameters, the outer-grid size, the outer Pareto-\hat{k} diagnostic when present and the wall-clock timing line ("fit in 5h 25m (grid 2h 09m)") when fit$timing is attached, then the fixed effects.

Value

A list with:

• theta_grid: matrix or vector of grid hyperparameter values.

• log_marginal: log p(y, mode | theta_k) at each grid point.

• weights: integration weights normalising to sum 1.

• theta_mean, theta_sd: posterior moments per hyperparameter.

• n_iter: inner Newton iterations per grid point.

• modes: matrix ⁠[n_grid x n_x]⁠ of inner modes, when stored.

• pareto_k, pareto_k_is_ess: outer Pareto-\hat{k} and its importance-sampling ESS (NA when not computed for the grid; see control$diagnose_k).

• pareto_k_proposal_source, pareto_k_first_pass: which proposal family the reported k-hat came from (the outer k-hat is scored against several candidates and the best is kept), and the k-hat of the FIRST pass – the proposal exactly as this backend placed it, before any candidate refined or replaced it. A large gap between the two says the nodes are badly scaled around the hyperparameter posterior even though the verdict is fine; no gap says the placement was already right.

• inner_skew, inner_skew_idx, inner_skew_dropped: the inner-Laplace skewness diagnostic (gamma_3) at each scored latent index and its 1-based index, plus a count of (index, observation) contributions dropped for a non-finite third derivative (see control$diagnose_skew).

• timing: named numeric of wall-clock seconds (total, setup, grid, postproc, diagnostics); the grid phase is the inner Laplace pass that scales with grid size. Surfaced one-line in print.

• prune_cheap_log_marginal, prune_mask, prune_n_pruned, prune_tol: present when the grid was screened and the safety gate did not trip – the cheap-screen log-marginal per cell (the truncated solve's value plus half its Newton decrement, the second-order estimate of the converged value the cells are ranked by; the decrement itself is prune_screen_decrement), the logical mask of pruned cells, the pruned-cell count, and the threshold actually applied.

• prune_log_gap_cut, prune_cheap_lm_spread, prune_min_keep, prune_n_floor_restored: the same screen read on the scale it operates on – how many nats below the best screened cell the tolerance cut at, how many nats the screened surface spans, the kept-cell floor applied, and how many cells that floor put back. A cut that is a sliver of the spread is a tolerance with no resolution on this grid.

• prune_dropped_mass_bound, prune_n_repaired, prune_gate_mass: the posterior mass the dropped cells may carry, with each lifted by the worst screening error measured on a kept cell; how many dropped cells the driver solved in full to bring that under the gate; and the gate.

• prune_fallback_triggered, prune_fallback_reason: present only when the gate did trip; the rest of the fit is then the full-grid result and the reason string records which check fired.

• prior: echoed input.

References

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

Examples

set.seed(1)
S <- 30L                                   # spatial units arranged in a chain
nb <- lapply(seq_len(S), function(s) setdiff(c(s - 1L, s + 1L), c(0L, S + 1L)))
nn <- lengths(nb)
field <- as.numeric(scale(cumsum(rnorm(S, 0, 0.4))))   # smooth spatial field
idx <- rep(seq_len(S), each = 6L); n <- length(idx); x <- rnorm(n)
y <- rbinom(n, 1L, plogis(-0.2 + 0.6 * x + field[idx]))
prior <- list(type = "icar", n_spatial_units = S, spatial_idx = idx,
              adj_row_ptr = c(0L, cumsum(nn)), adj_col_idx = unlist(nb) - 1L,
              n_neighbors = nn, tau_grid = c(0.5, 1, 2, 4, 8))
fit <- tulpa_nested_laplace(y, rep(1L, n), cbind(1, x), prior = prior,
                            family = "binomial")
fit$theta_mean        # marginalized ICAR precision

Joint multi-likelihood nested Laplace approximation

Description

Outer-grid nested Laplace driver for joint models – multiple response arms sharing one latent prior block, parameterized as a per-arm field amplitude (sigma) on a unit-precision latent.

Supported priors:

• "bym2" – outer grid over ⁠(sigma, rho [, alpha])⁠. Latent: phi (n_s) | theta (n_s) with unit-precision ICAR + iid components.

• "icar" – outer grid over (sigma [, alpha]). Latent: phi (n_s) with unit-precision ICAR.

• "car_proper" – outer grid over ⁠(sigma, rho_car [, alpha])⁠. Latent: phi (n_s) with Q = D - rho_car * W.

Other backends (NNGP, HSGP, RW1/2, AR1) follow the same interface and land under Phase 3.

Usage

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

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

tulpa_nested_laplace_joint(
  responses,
  prior,
  copy = NULL,
  phi_grid = NULL,
  prior_sigma = NULL,
  prior_alpha = NULL,
  prior_phi = NULL,
  hyperprior = c("proper", "flat"),
  cell_coupling = "separable",
  control = list()
)

Arguments

prior_alpha = list(
  list(block = 1, prior = list("pc.prec", c(4.0, 0.01))),
  list(block = 2, prior = list("half_normal", 2.0))
)

Value

A list of class c("tulpa_nested_laplace_joint", "tulpa_nested_laplace", "list") with:

• theta_grid, theta_names – outer-grid hyperparameter values (includes the alpha axis when copy is set).

• log_marginal, weights – per-grid-point log-marginal and integration weights (sum to 1).

• theta_mean, theta_sd – posterior moments per hyperparameter, including alpha when copy is set.

• theta_median, theta_ci_lo, theta_ci_hi – weighted-quantile median and 2.5/97.5 empirical CI per hyperparameter axis (same names as theta_mean). Recommended summary for right-skewed scale-like axes (alpha at small n_pos, sigma/range/phi near a boundary), where the posterior mean is pulled by the right tail away from the bulk and ⁠mean +/- 1.96 sd⁠ mis-states the uncertainty.

• theta_sd_source – per axis, what the reported SD was read off: "weighted" (the spread of the weights, on a resolved axis), "stencil" (a parabola at the modal node) or "within_cell" (the modal node's mass spread over its cell, the read the interval of an unresolved axis comes from). theta_sd_stencil_declined says why a parabola was not formed.

• modes – ⁠[n_grid x n_x]⁠ matrix of inner modes.

• fitted_eta – ⁠[n_grid x N]⁠ matrix of the linear predictor at each cell's own mode, offset and latent field included. Present on a single-arm fit, which is the one that has a single linear predictor to report; a cell that was pruned or whose block preparation failed reads NA. fitted_eta_var is its per-cell within-cell variance, present under control$fitted_var (the default). Together these are the per-cell Gaussian posterior_predict() draws a replicate from.

• n_iter – inner Newton iterations per grid point.

• arm_layout – list with per-arm beta_start, re_start, spatial offset(s) and n_x for decoding modes.

• prior, responses, copy – echoed inputs.

• timing – named numeric of wall-clock seconds: total plus the setup (validation / encoding / grid construction), grid (inner Laplace solves, including adaptive-refinement and consistency passes), postproc (weight / moment / marginal assembly) and diagnostics (outer Pareto-\hat{k}) phases. The grid phase is the one that scales with grid size and core count. Surfaced one-line in print.

• pareto_k, pareto_k_is_ess, pareto_k_scope – outer Pareto-\hat{k} accuracy diagnostic and its importance-sampling ESS (both NA when control$diagnose_k = FALSE or the fit declines; see the diagnose_k control knob). pareto_k < 0.7 indicates the nested integration is reliable; ⁠>= 0.7⁠ that the (skewed / heavy- tailed) hyperparameter posterior is misfit by the Gaussian grid and the fit should escalate to an exact debias.

• pareto_k_declined – when pareto_k is NA, WHY: "not_requested" (control$diagnose_k = FALSE; nothing is wrong), "unguessable_axis: <axis>" (a support the engine will not guess, e.g. car_proper's rho_car – a PERMANENT limitation of that family, so read the quadrature ESS instead), "draws_too_few", "grid_too_small", "no_varying_axis", "degenerate_proposal" (the outer mode curvature came back non-finite – a signal about the fit), "not_applicable", or "internal_inconsistency" (an engine bug worth reporting). NA on a fit whose \hat{k} WAS computed.

• pareto_k_proposal_source – how the outer importance proposal the \hat{k} scores was built: "mode_hessian" from the Laplace curvature at the hyperparameter mode (the CCD design's, or a finite-difference Hessian when a sharp posterior collapses the grid), "grid_moment" from the grid-weighted covariance, "moment_matched" when the moment-matching refinement improved on either, "grid_mixture" when the faithful mixture-over-grid-cells proposal won, or "skew_normal" when the skew-normal rescue did. NA when the diagnostic is off or declines. The mode-Hessian source keeps the \hat{k} meaningful when the grid concentrates on ~1 cell.

• pareto_k_regime – whether the outer grid actually integrated hyperparameter uncertainty: "spread", or "collapsed_interior" / "collapsed_edge" when the quadrature effective sample size falls below two cells so the integration has degenerated to a point evaluation at the modal hyperparameter. An interior collapse is benign (the grid bracketed the mode; the fit is empirical Bayes at that mode and only the integrated hyperparameter uncertainty is missing); an edge collapse means the mode sits at an extreme node, so the grid may be too narrow. Attached from the stored weights, so it is present even with control$diagnose_k = FALSE.

• pareto_k_grid_edge_axes, pareto_k_grid_edge_sides – for an edge collapse, which axes the dominant cell sits against and on which side ("lower" / "upper"); widen those axes and refit to confirm the mode is bracketed. Axes the grid pins to one value are excluded (pinned, not at a boundary).

• axis_span – per outer axis, the span the fit worked over and the two terms that separate it from the nodes it was given: nodes (the range of the axis's initial continuum nodes), declared (their support, i.e. those nodes plus the half node step the outermost cells own – for k equally spaced nodes that is k / (k - 1) times the node range on the axis's integration coordinate, so 2x at two nodes and 1.125x at nine), integrated (the region the final grid's cell measure integrates, the same interval axis_support reports), refine (the mode refinement ran under: "none" / "densify" / "extend") and n_nodes (initial and final continuum node counts).

• outer_grid_placement – "fixed" (the default sigma_grid axis was used as-is) or "auto_recentered" when a railed sigma axis triggered the mode-Hessian recenter-and-refit (see the prior argument above). outer_grid_recenter_attempts (integer) and outer_grid_prior_added (logical: whether the light default PC(U=3, alpha=0.01) sigma prior was applied) are set alongside it, plus outer_grid_prior_declined = "prior_pinned" when a second attempt ran but the caller's pinned prior_sigma held that prior back.

• outer_grid_recenter_declined – on a "fixed" placement, why the recenter did not run: "grid_not_collapsed" (the grid already brackets the mode – the common case, no refit needed), "axis_pinned" (the caller pinned sigma_grid; mark it with auto_grid() if it is a default rather than a choice), "default_axis_pinned" (a wrapper package declared the nodes with auto_grid(place = FALSE) and asked for them as written – nothing you pinned), or "no_usable_curvature" (the outer mode-find the recenter lays the axis from did not reach a mode with a usable curvature, e.g. a car_proper grid whose rho_car axis has unguessable support). Absent when the fit WAS recentred. A joint fit's field SD and its per-arm dispersion are placed by different passes over the same grid, and this slot reduces over them rather than holding the last to speak: an axis held because it was declared does not stand as the fit's answer while a pass with a movable axis has one.

• outer_grid_axis_declined – the same question PER AXIS, as a named character vector. The slot above holds one reason for the whole fit and is written only while the fit is unplaced, so on a fit where one axis moved and another did not it says auto_recentered and nothing about the axis that stayed – which is the axis grid_coarsest_axis then names. Written by the field-SD passes for sigma (⁠b<k>.sigma⁠ on a copy block) and by the per-arm dispersion pass for ⁠phi_<arm>⁠, and carried across a later pass placing a different axis; absent on a fit no pass declined an axis of. An axis a pass DID move that still rails reads "attempts_exhausted" (or the reason its last attempt could not be laid), and is the axis the railed-mode warning names.

• outer_grid_pilot – present only when control$recenter_pilot ran: the pilot's resolution (n_pilot), its cell count (cells), the axes it thinned (axes) and those it could not (axes_kept), and what the placement was DECIDED from – the pilot's own regime, ess_grid and edge_axes. The reported grid is never the grid the trigger was read on, so a two-pass fit is legible from the object rather than only from the progress log. outer_grid_pilot_declined = "pilot_placement_declined" records a pilot whose detection placed nothing and which was therefore followed by the full fit and the full fit's own detection.

• outer_grid_recenter_sd_clamp, outer_grid_recenter_sd_raw, outer_grid_recenter_sd_used – on a recentred fit, one entry per moved axis: which mode-SD bound the placement hit ("none" / "floor" / "ceiling"), the SD the finite-difference stencil MEASURED, and the SD the axis was actually laid from. A clamped axis was laid from a substituted spread rather than a measured one, and the reported interval is read off that span; the two used to be indistinguishable on the fit.

• pareto_k_outer_skew – per-axis skewness of the hyperparameter marginal in the proposal's whitened coordinate, present only when a \hat{k} above the good band triggered the skew-normal rescue pass. It is the EXPLANATION for an inflated \hat{k}: a symmetric Gaussian proposal against a right-skewed variance-component marginal has a heavy importance-ratio tail whatever the integration's quality. NULL when the Gaussian proposal already fit.

• pareto_k_se_boot, pareto_k_ci_low, pareto_k_ci_high, pareto_k_se_formula, pareto_k_tail_points, pareto_k_tail_points_requested, pareto_k_band_confident – the outer Pareto-\hat{k} uncertainty, present whenever the diagnostic ran. pareto_k_se_boot is the bootstrap SE of the k-hat and pareto_k_ci_low / pareto_k_ci_high its 2.5\ quantiles – the estimator's sampling spread GIVEN the proposal, not a posterior credible interval; pareto_k_se_formula is the closed-form GPD-shape MLE asymptotic SE cross-check. pareto_k_tail_points is the GPD tail size used and pareto_k_tail_points_requested the k_tail_points request (NA when automatic). pareto_k_band_confident is TRUE iff the bootstrap CI lies within one reliability band, NA when the bootstrap was off or could not fit. diagnose_draws and diagnose_cost_ratio (the diagnostic's draw budget and its wall-clock cost relative to the fit) are attached at the top level.

• pareto_k_by_arm, pareto_k_by_arm_is_ess, pareto_k_by_arm_scope – present only with control$diagnose_k = "by_arm". Named (by arm) outer Pareto-\hat{k} restricted to each arm's hyperparameter axes, the other arms held at their posterior mean, so a tail-heavy joint k can be localised to one arm. A per-arm entry is NA when that arm carries no varying axis. Each arm carries its own bootstrap uncertainty: pareto_k_by_arm_se_boot, pareto_k_by_arm_ci_low / pareto_k_by_arm_ci_high, pareto_k_by_arm_se_formula, pareto_k_by_arm_tail_points and pareto_k_by_arm_band_confident.

• k_quality_requested, k_quality_reached, k_quality_best, k_quality_miss, k_quality_reason, k_quality_rounds, k_quality_k_trace – the reliability verdict for the control$k_quality intent. k_quality_requested echoes the intent; k_quality_best is the band THIS fit reached ("good" / "ok" / "unreliable", or "uncertain" when the bootstrap CI crosses a boundary) – the band of the returned fit, not the best band seen while escalating, which k_quality_k_trace carries; k_quality_reached is TRUE/FALSE for an "ok" / "good" target (NA for "report" / "none"), never silently downgraded; k_quality_miss classifies a miss as "resolution" (the \hat{k} confidently outside the band – a grid deficiency) or "precision" (the bootstrap CI straddling a boundary – the estimator's variance), NA when there is no miss to classify, and is what the escalation reads once refinement is exhausted; k_quality_reason records why it stopped; k_quality_rounds is the number of escalation re-fits performed (0 when the first fit sufficed or escalation was off); and k_quality_k_trace is the outer \hat{k} round by round, starting at the first fit's, so a chase that did not descend monotonically is visible. The refinement rungs rebuild the proposal from the grid they just changed, so a rising trace is a re-founded proposal on a strictly larger grid rather than a worse fit, which is why the LAST round is returned rather than the lowest-\hat{k} one.

• adaptive_grid_info – when adaptive_grid = TRUE, a list with triggered_axes (character) and n_points_added (integer) describing the refinement passes. NULL otherwise.

• local_ccd_info – when local_ccd engaged (or k_refine = "ccd"), a list with n_cells_refined, n_nodes_added, the refined cells, the n_cells_before / n_cells_after grid sizes, n_design_nodes (the cells that carry an in-cell design weight rather than their own mass), design_mass (the share of the integration weight sitting on them) and cell_share (the share each refined cell held on the BASE grid, before any node was placed). The two shares differ: the replacement nodes sit nearer the peak than the cell's own coordinate did, so refining raises the share the refined region holds, and reading both separates how concentrated the base grid already was from how much the refinement concentrated it. Three per-cell readings come with it, on three orthogonal axes, none of them gating anything. Shape: misfit (the standardized cubic magnitude of the part of the cell the design cannot represent, the one reading skew_max gates on, with misfit_declined / cells_declined for the cells put back as their own mass atoms). Centring: offset (the norm of the whitened gradient, i.e. how far the cell's own peak sat from the cell's coordinate in units of its marginal spread) and mode_gain (the same displacement in the cell's own curvature units, ⁠0.5 * g' (-H)^-1 g⁠ in nats, read off the same fit's quadratic coefficients and comparable across cells whose curvature differs; NA where -H is not positive definite, a cell with no interior peak to be displaced from). Mass: log_mass_ratio (the cell's refined mass over the coarse atom it replaced, log_mass_refined - log_mass_coarse, both carried on the log scale with the cell's outer design weight in them) and max_node_weight (the share the single largest node takes of its own cell's refined mass). Each carries a ⁠_declined⁠ twin for the cells the gate put back, whose nodes were evaluated before the score was read. misfit certifies that the design can represent the cell's shape and nothing more: a cell can be exactly quadratic and still be read at a point that is not representative of it (offset / mode_gain), and refining it can still move how much mass it competes for against its unrefined neighbours (log_mass_ratio, the embedded-rule local error indicator of adaptive cubature). NULL when local CCD was off or declined (single-block, ⁠< 4⁠ axes, an active phi_grid, or no peaked interior cell).

• integration_requested, integration_declined – what integration asked for, and why the requested CCD or adaptive lattice did not run. ⁠$integration⁠ names the integrator that RAN, and .nl_node_support() keys the interval construction off it, so a caller who asked for a moment rule and received a density grid reads the reason here rather than inferring it. integration_declined is NA_character_ when nothing was declined, and otherwise one of "axis_count" (fewer transformable latent axes than the requested mode's threshold), "unguessable_axis" (a CAR_proper rho_car or a non-BYM2 rho), "degenerate_axis" (a single-valued axis), "modefind_ridge" / "modefind_boundary" / "modefind_degenerate" / "modefind_failed", "hessian_singular", "hessian_not_pd", "copy_atom_components" (more copy atoms than the design splits into, see copy_atom_mass), "copy_atom_mass" (a declared atom mass outside ⁠[0, 1)⁠) or "placement_budget" (placing the design would cost more inner solves than integrating the tensor grid it replaces; see the ccd_budget control knob). Under integration = "grid_adaptive" it is one of "adaptive_small_grid" (fewer outer cells than adaptive_grid_min_cells), "adaptive_degenerate_lattice", "adaptive_seed_failed" (no finite marginal in the coarse seed), "adaptive_diffuse_seed" / "adaptive_flood_cap" / "adaptive_kept_cap" (the posterior mass spreads over more than adaptive_grid_max_frac of the dense grid, read at the seed, during the flood, or on the closed kept region), or "adaptive_error" (the flood raised an error, reported as a warning carrying its message); the fit then integrates the dense tensor. Multi-block fits only.

• ccd_modefind_evals, ccd_modefind_budget, ccd_modefind_rounds, ccd_modefind_seconds – what PLACING the CCD cost: inner Newton solves spent by the mode-find (seed, step calibration and every round of every mixture component), the ceiling they were spent against, the rounds taken and the wall-clock seconds. Present on any multi-block fit that attempted a placement, including one that declined – a decline reports what it had already paid – and absent when no placement was attempted, which is not the same as a placement that cost nothing. Read alongside ⁠$integration⁠ to judge whether the CCD earned its design on this model rather than arguing it from the round caps.

• weight_kind – one entry per outer-grid cell, "mass" or "design". A tensor cell holds the mass of its own cell and a CCD node holds a design weight, so a fit integrated by one rule reports one value throughout; a locally refined grid carries both, and this says which per cell instead of leaving it to be read off integration.

• theta_interval_read, theta_interval_design_mass – what the reported per-axis theta_median / theta_ci_lo / theta_ci_hi were read off, and how much of the integration weight sits on nodes whose cumulative sum is not a CDF. "density" is a weighted quantile over cells that discretize the posterior (design_mass 0), "moment_rule" an interval from the moments a central-composite design delivers (design_mass 1), and "mixed" the locally CCD-refined grid, which carries both kinds at once. A mixed support still reports the weighted quantile, which measured best against a converged reference in both design_mass regimes; the share is the regime variable, since on that part of the support the quantile is bounded by the refined cells' own grid neighbourhoods rather than by the posterior. Every nested path stamps the pair, not only the multi-block driver.

• within_cell_requested, theta_within_cell, theta_within_cell_declined – the WITHIN-CELL construction the same intervals were read with. The kind above says what the integrator left; this says how each cell's mass was spread inside its own box when the grid was read back. "log_quadratic" is the default and interpolates each row's node log densities; "box_uniform" (control$within_cell) puts the cumulative full mass at each cell edge; "chord" puts the cumulative mid-mass at each cell coordinate, the same masses over the same boxes with the knots moved half a cell. The construction is recorded per axis, and an axis whose cell partition could not be built falls back to "chord" on its own with the reason in theta_within_cell_declined ("support_<kind>", "single_node", "boxes_do_not_tile", "no_usable_node"). A "log_quadratic" request falls back to "box_uniform" first, with "unresolved" (no row resolves its own conditional), "not_a_grid_axis" or "no_node_density" as the reason.

• theta_cell_edge_coord, theta_cell_edge_declined – per axis, the COORDINATE that axis's outer cell edges were mirrored in, and why the support the axis declares did not produce them. A declared support is authoritative: its mirrored edge is used when finite and inside it, and otherwise the axis reports its extreme grid coordinate – which the containment test has already placed inside the support – rather than falling through to a coordinate guessed from the node values, which is what put a lower bound of 0 on a positive axis and a bound above 1 on a unit one. An axis that declares no support at all keeps the guess by design, but its mirrored edge is checked for being a representable double that brackets the coordinates, and falls back to the extreme grid coordinate with "mirrored_edge_not_representable" when it is not. theta_cell_edge_declined is NA where nothing was declined.

• outer_grid_cell_width, outer_grid_axis_sd, outer_grid_h_over_sd – per axis, the cell width and the posterior SD in that axis's own coordinate, and their ratio. A within-cell reconstruction resolves an interval endpoint to within one cell, so how much of the reported width and of its realized coverage is a property of where the grid fell rather than of the posterior is governed by this ratio; below .nl_diag("grid_resolved") the cells are narrower than the posterior they discretize and the two constructions converge to each other. NA on an axis whose own marginal is maximal at an endpoint, which is the placement question outer_grid_railed_axes reports.

• outer_grid_resolution_declined – per axis, why that axis carries no ratio, from the closed vocabulary too_few_nodes, mode_at_edge, coord_not_finite, spacing_not_finite, stencil_degenerate, curvature_not_negative; NA on an axis that scored. The reasons are not interchangeable: mode_at_edge says the grid does not contain that axis's own posterior mode, which is a stronger statement about the fit than any ratio, and the whole-grid resolved verdict is withheld while one is present rather than being read off the axes that did score.

• prune_cheap_log_marginal, prune_mask, prune_n_pruned, prune_tol, prune_screen_iters – present only when prune = TRUE and the safety gate did not fall back. Cheap-pass log-marginals at every cell – the truncated solve's value plus half its Newton decrement, the second-order estimate of the converged value the screen ranks by, so a cell whose inner mode had far to move from its warm start is not ranked on where the truncation left it (gcol33/tulpa#919); that decrement is prune_screen_decrement – a logical mask of pruned cells, the pruned-cell count, the threshold actually applied, and the screening depth the driver actually ran at. Pruned cells have log_marginal = -Inf so they get zero weight under .nl_normalise_weights_safe.

• prune_log_gap_cut, prune_cheap_lm_spread, prune_min_keep, prune_n_floor_restored – the screen read on the scale it operates on: how many nats below the best screened cell the tolerance cut at, how many nats the screened surface spans, the kept-cell floor applied, and how many cells that floor put back. A cut that is a sliver of the spread is a tolerance with no resolution on this grid.

• prune_dropped_mass_bound, prune_n_repaired, prune_gate_mass – the posterior mass the dropped cells may carry, with each lifted by the worst screening error measured on a kept cell; how many dropped cells the driver solved in full to bring that under the gate; and the gate.

• outer_mode_carried – TRUE on a placement refit that read the outer mode it was laid from rather than finding it again, which it does unless the placed grid rails.

• prune_fallback_triggered, prune_fallback_reason – present only when a fallback to the full grid fired. The returned fit is the full-grid (unpruned) result; the reason string records whether the cheap-pass safety gate tripped or the screened grid left the placement pass no curvature.

• n_threads_outer_requested, n_threads_outer_realised – the outer-grid width asked for and the width the grid ran at. The driver clamps the request to the team the OpenMP environment hands out, so a fit launched under OMP_NUM_THREADS=1 with n_threads_outer = 10 runs serial; the pair is what a timing has to be recorded against.

• converged, n_iter, score_max – per cell, whether the inner Newton reached a mode within control$max_iter, the iterations it took, and the largest absolute component of the joint penalized score at the mode it reported (the residual a cell that stopped short stopped at). A cell that stopped without a mode is solved again from the mode of its nearest converged cell (nonconverged_resolve counts the cells that first stalled and those the re-solve recovered); one that still has no mode stays in the integration at the value it stopped at, is listed in nonconverged_cells with the posterior mass such cells carry in nonconverged_mass, and is named by diagnostics(). Past the screen's mass gate the fit warns.

(sigma, alpha) parameterization

Each arm's linear predictor reads

\eta_{arm} = X_{arm} \beta_{arm} + \sigma_{arm} \cdot z_s,

where z_s is a unit-precision latent (ICAR(tau=1) for ICAR/BYM2, or the BYM2 mix; CAR_proper uses the structure of D - \rho_{car} W). All arms share the donor amplitude \sigma from the outer-grid sigma_grid axis. The copy arm's amplitude is \sigma_{arm} = \alpha \cdot \sigma, where \alpha is a direct outer-grid axis taken from copy$alpha_grid. The Cartesian product is over ⁠(sigma, [rho/rho_car,] alpha[, phi])⁠. Direct \alpha as an outer-grid axis (rather than a post-hoc ratio \alpha = \sigma_{pos} / \sigma_{occ}) avoids plug-in bias on the weakly-identified ratio at small n_pos and lets a regularizing hyperprior land on \alpha directly.

References

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

tulpa_nested_laplace() for the single-arm engine.

Examples

set.seed(1)
S <- 25L                                   # shared spatial units (chain graph)
nb <- lapply(seq_len(S), function(s) setdiff(c(s - 1L, s + 1L), c(0L, S + 1L)))
nn <- lengths(nb)
field <- as.numeric(scale(cumsum(rnorm(S, 0, 0.4))))
mk_arm <- function(m, fam) {
  si <- sample(S, m, replace = TRUE); x <- rnorm(m)
  lin <- 0.2 + 0.5 * x + 0.8 * field[si]
  y <- if (fam == "binomial") rbinom(m, 1L, plogis(lin)) else lin + rnorm(m, 0, 0.5)
  list(y = as.numeric(y), n_trials = rep(1L, m), X = cbind(1, x),
       spatial_idx = si, family = fam, phi = if (fam == "gaussian") 0.5 else 1)
}
prior <- list(type = "icar", n_spatial_units = S,
              adj_row_ptr = c(0L, cumsum(nn)), adj_col_idx = unlist(nb) - 1L,
              n_neighbors = nn, sigma_grid = c(0.3, 0.7, 1.4))
fit <- tulpa_nested_laplace_joint(
  responses = list(occ = mk_arm(200L, "binomial"), pos = mk_arm(200L, "gaussian")),
  prior = prior)
fit$theta_mean        # shared field amplitude, integrated across both arms

Extract temporal effects from a fitted model

Description

Extract posterior distributions of temporal effects from a fitted tulpa model with temporal specification.

Usage

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

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

temporal(
  object,
  component = "all",
  summary = FALSE,
  probs = c(0.025, 0.5, 0.975),
  ...
)

## S3 method for class 'tulpa_fit'
temporal(
  object,
  component = "all",
  summary = FALSE,
  probs = c(0.025, 0.5, 0.975),
  ...
)

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

## S3 method for class 'tulpa_temporal_posterior'
summary(object, probs = c(0.025, 0.5, 0.975), ...)

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

Arguments

Details

temporal() is overloaded. Given a fitted model it is the accessor described here. Given a one-sided formula (or a named ⁠formula =⁠ / ⁠structure =⁠ argument) it is instead the inline varying-coefficient field constructor used in a tulpa() model formula, the temporal mirror of spatial(): temporal(formula = ~ 1 + x || time, structure = "rw1") declares a smooth temporal level (the intercept column) plus a temporally varying slope on each covariate column. structure is one of "rw1" (default), "rw2", or "ar1"; only the double bar || (independent fields) is supported.

On a sampler fit the posterior is the field's own draws. On a nested-Laplace fit it is the Gaussian mixture the outer grid defines: each hyperparameter cell contributes its conditional mode and marginal variance, weighted by the cell's posterior weight, and summary() reports that mixture's exact mean, SD and quantiles. The draws of such a fit are sampled from the same mixture, each from its cell's joint Gaussian, so they carry the correlation between time points as well as each point's marginal, and a quantity derived from them (a trend, a difference between two times) has the right spread. A field so large that its per-cell covariance was not retained (more than 4 million values over the grid) is drawn with independent time points within each cell instead. The time column of the summary holds the time values the field was fitted on.

Value

A tulpa_temporal_posterior object

summary() on a tulpa_temporal_posterior returns a tulpa_temporal_summary data frame with one row per time point (per component for multi-scale fits), holding the posterior mean, SD and requested quantiles of the temporal effect.

See Also

temporal_multiscale(), temporal_rw1()

Examples

set.seed(131)
df <- data.frame(year = 1:40, x = rnorm(40))
df$count <- rpois(40, exp(1 + 0.2 * df$x))

fit <- tulpa(
  count ~ x,
  data = df,
  family = "poisson",
  temporal = temporal_multiscale("year", trend = "rw2", seasonal = 12),
  mode = "exact",
  control = list(n_iter = 200L, n_warmup = 100L, seed = 1L)
)

# Extract all temporal effects
temp_post <- temporal(fit)
summary(temp_post)

Beta prior

Description

Specify a beta prior for a parameter bounded in (0, 1).

Usage

prior_beta(alpha = 1, beta = 1)

Arguments

Details

• alpha = beta = 1: Uniform(0, 1)

• alpha = beta = 2: Symmetric, peaked at 0.5

• alpha > beta: Skewed toward 1

• alpha < beta: Skewed toward 0

Value

A tulpa_prior object

Examples

prior_beta(1, 1)   # Uniform
prior_beta(2, 2)   # Symmetric, centered at 0.5
prior_beta(5, 2)   # Skewed toward 1 (for high autocorrelation)

Exponential prior

Description

Specify an exponential prior for a positive parameter.

Usage

prior_exponential(rate = 1)

Arguments

Value

A tulpa_prior object

Examples

prior_exponential(1)

Build a prior list for tulpa_nested_laplace() from a tulpa spec object

Description

Validates a tulpa_temporal or tulpa_spatial specification against data, then converts it to the prior list shape consumed by tulpa_nested_laplace(). Mainly an internal helper for callers that already have a fitted spec; users typically pass spec + data directly to tulpa_nested_laplace() instead.

Supported spec types:

• tulpa_temporal with ⁠type in {"rw1", "rw2", "ar1"}⁠

• tulpa_spatial, areal ⁠type in {"car", "icar", "car_proper", "bym2"}⁠. An areal spec returns the graph fields only – type, spatial_idx, n_spatial_units, adj_row_ptr, adj_col_idx, n_neighbors – plus scale_factor for bym2 and the eigenvalue-derived rho_bounds for proper CAR. No grid field is returned for any areal type: the outer axes are built later by the family's registry defaults() closure, and proper CAR's ⁠(tau, rho)⁠ grid is built there from those bounds.

• tulpa_gp / tulpa_hsgp, continuous ⁠type in {"gp", "nngp", "hsgp"}⁠: validated against data and routed to the nngp / hsgp nested kernel through the shared converter.

SPDE is the one continuous field not built here – it carries its own (range, sigma) FEM integrator; call fit_spde() with the spatial_spde() spec.

Usage

prior_from_spec(spec, data)

Arguments

Value

A prior list ready for tulpa_nested_laplace().

Examples

d <- data.frame(year = rep(1:10, each = 3))
pr <- prior_from_spec(temporal_ar1("year"), d)
str(pr)

Gamma prior

Description

Specify a gamma prior for a positive parameter.

Usage

prior_gamma(shape = 2, rate = 0.1)

Arguments