Package {tulpa}

Title:	Template Unified Latent Process Architecture for Bayesian Hierarchical Models
Version:	0.2.0
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, lifecycle, stats, tools, utils, methods, graphics, grDevices
Suggests:	testthat (≥ 3.0.0), bayesplot, 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.0.0
NeedsCompilation:	yes
Packaged:	2026-08-27 23:30:17 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)
Maintainer:	Gilles Colling <gilles.colling051@gmail.com>
Repository:	CRAN
Date/Publication:	2026-09-09 15:30:02 UTC

tulpa: Template Unified Latent Process Architecture for 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]

See Also

Useful links:

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

• https://gillescolling.com/tulpa/

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

Reject non-finite fitting inputs.

Description

The design is built with na.action = na.pass, so a missing predictor or response survives into X / y. tulpa() does not drop incomplete cases, so an NA/NaN/Inf would propagate into the C++ kernels as a NaN estimate. Fail loudly with the offending row instead.

Usage

.assert_finite_model_inputs(X, y)

Resolve a supplied fixed-effect prior to its Gaussian ⁠(mean, sd)⁠ fields

Description

The one place that decides whether a beta_prior argument can be expressed as a Gaussian prior on the fixed effects at all. .normalize_beta_prior() (per-coefficient) and .beta_prior_ridge_sd() (scalar ridge) both resolve through it and differ only in the shape they recycle the fields to, so a prior one accepts is a prior the other accepts.

Usage

.beta_prior_fields(beta_prior)

Arguments

Details

A prior with no sd is an input the fitters cannot express: substituting the default for it would replace the user's modelling statement with a different one and leave no trace on the posterior.

Value

list(mean, sd), unrecycled.

Resolve a beta_prior = list(mean, sd) to a single ridge SD.

Description

For the fitters whose fixed-effect prior is a mean-zero Gaussian ridge with a scalar SD (EP, multinomial, ordinal): validates the unified beta_prior object, enforces a mean of 0, and returns the positive scalar SD. Keeps the shared prior interface (beta_prior = list(mean, sd)) while rejecting the options those fitters do not implement (a non-zero mean, a per-coefficient SD). Validation of the supplied object itself is .beta_prior_fields(), shared with .normalize_beta_prior().

Usage

.beta_prior_ridge_sd(beta_prior, default_sd = .tulpa_prior_sd("ridge"))

Arguments

Assemble NUTS initial positions and an inverse-mass diagonal

Description

Assemble NUTS initial positions and an inverse-mass diagonal

Usage

.build_warm_start(
  fit,
  layout,
  re_terms,
  n_chains,
  sigma_re = NULL,
  jitter = 1,
  metric = FALSE
)

Arguments

Value

⁠list(init = <n_chains x D matrix>, inv_metric_diag = <D vector or NULL>)⁠.

Resolve a family spelling to its canonical registry name.

Description

A no-op for a canonical name, for a ⁠<family>_<link>⁠ code, and for anything unrecognized, which reaches .family_or_stop() and errors there against the canonical list.

Usage

.canonical_family(family)

Arguments

Proper-CAR precision Q = tau * (D - rho W) at fixed (tau, rho).

Description

D is the diagonal degree matrix (neighbour counts) and W the adjacency, so Q has the same nonzero structure as the ICAR precision with rho scaling the off-diagonals. Matches tulpa::add_car_proper_prior / the nested CAR_proper precision builder.

Usage

.car_proper_precision_Q(spatial, tau, rho)

Families whose dispersion can be estimated

Description

Families whose dispersion can be estimated

Usage

.dispersion_families()

Dispersion derivatives for a family, or NULL

Description

Dispersion derivatives for a family, or NULL

Usage

.family_dphi(family)

Arguments

Value

A list with dloglik, dscore and dweight, or NULL when the family has no free dispersion or none has been derived. NULL is a refusal to estimate, never a signal to fall back to a fixed value silently.

Second-order dispersion derivatives for a family, or NULL

Description

Second-order dispersion derivatives for a family, or NULL

Usage

.family_dphi2(family)

Arguments

Value

A list with dloglik2, dscore2, dweight2 and dweight_deta, or NULL when the family's phi Hessian has not been derived. NULL is a refusal, not a fixed-value fallback: it hands the phi Hessian to the differencing stencil rather than reporting an inexact closed form.

Whether a family's phi Hessian additionally needs d(W_obs)/dphi

Description

TRUE whenever the Newton working weight is not the observed curvature, which is exactly when the border's Hinv and the mode motion's Hinv_mode are different matrices. For a family whose weight is also free of eta the two inverses multiply a zero channel and the answer would come out the same either way, but the entry is required regardless rather than resting on that second coincidence.

Usage

.family_dphi2_needs_obs(family)

Observed curvature (-d^2 log-lik / d eta^2 at the realized y), elementwise.

Description

Falls back to the y-free weight for families whose curvature carries no y – those where the response enters the log-likelihood linearly in eta (Poisson, binomial, and the zero-truncated Poisson), for which the observed and expected curvature coincide identically rather than only in expectation.

Usage

.family_obs_weight(eta, y, family, n_trials = NULL, phi = 1, phi2 = NULL)

Look up a family's operation set, with a clear error for unknown families.

Description

Look up a family's operation set, with a clear error for unknown families.

Usage

.family_ops(family)

Build a length-n_obs -> length-n_units indicator design matrix Z.

Description

The latent field lives on units – unique coordinates for an NNGP field, adjacency units for an areal one – and each observation belongs to exactly one, so Z[i, j] = 1 iff observation i sits on unit j.

Usage

.field_design_Z(obs_to_unit, n_units, n_obs)

Arguments

Details

The map is passed in rather than read off the spec, because the two specs that reach here name it differently: a GP spec carries obs_to_loc, an areal spec carries spatial_idx. Reading one name served only one caller; for the other the field was absent, the map silently became seq_len(n_obs), and the column index ran past n_units as soon as a location carried more than one observation.

Attach the standard tulpa_fit contract to a fitter's result.

Description

The dispatch layer and the fitters that return generic-accessor-facing objects route their return value through this helper (some special-purpose fitters – the logpost samplers, tulpa_ep, the categorical drivers – still stamp their class by hand), so a directly-called fitter and a tulpa()-dispatched one yield the same enriched object: the tulpa_fit class (so the generic S3 methods – coef / summary / vcov / confint / tidy / glance / ranef – dispatch), the fixed-effect layout (n_fixed / fixed_names / param_names, which the .fit_fixed_table summary path reads), and an explicit posterior-provenance tag (draws_kind) the chain-vs-iid diagnostic gate (.tulpa_draws_kind()) reads to decide whether Rhat/ESS apply. Each field is filled only when the fitter did not already set it, so a fitter that knows better wins and the helper is idempotent under tulpa_dispatch().

Usage

.finalize_fit(
  fit,
  backend = NULL,
  draws_kind = NULL,
  n_fixed = NULL,
  fixed_names = NULL,
  param_names = NULL,
  extra_class = NULL
)

Arguments

Details

draws_kind precedence is: a value already on the fit, then the explicit draws_kind argument, then the registry emits property for backend. Backends that are not registry keys (the ⁠tgmrf_*⁠ fitters) must pass draws_kind explicitly, since their emits cannot be looked up.

Value

The enriched fit, classed c(extra_class, ..., "tulpa_fit").

Connected components of an adjacency graph, as a list of node-index vectors.

Description

Iterative depth-first search, so it carries no recursion-depth risk on a large map. Returns what graph_partition (inst/include/tulpa/graph_components.h) returns – the actual component MEMBERSHIP, not just a count – so a genuine disconnected map (a mainland plus islands) pins each component's constant over that component's real nodes rather than an equal-size contiguous split. Takes the dense adjacency rather than a sparse one because symmetric sparse storage (dsCMatrix) keeps a single triangle, which would put an edge's two endpoints in different components.

Usage

.graph_components(adjacency)

Number of connected components of an adjacency graph.

Description

The count of .graph_components, single-sourcing the DFS.

Usage

.graph_n_components(adjacency)

Intrinsic ICAR field precision, augmented to full rank.

Description

L = D - W (degree minus adjacency) plus the sum-to-zero augmentation ⁠sum_c 1_c 1_c' / J_c⁠ that identifies each component's constant null direction (inst/include/tulpa/sum_to_zero.h). The components are the actual connected components of the graph (.graph_components), each pinned over its own node set of size J_c – matching for_each_icar_component (src/icar_kernel.h), so an unequal / non-contiguous disconnected map reconstructs the same precision the kernel penalizes with.

Usage

.icar_precision_Q(spatial)

Details

The conditional Laplace kernel passes tau_spatial = 1 on both the ICAR/CAR path (R/fit_laplace.R) and the BYM2 structured block (which carries sigma and rho in its d_fac instead), and the field enters the linear predictor with d_fac = 1, so this is exactly the field block of that fit's joint Hessian. test-marginal-se-areal.R pins it against the kernel's own log_prior_icar.

Marginal H_beta for a BYM2 field. The latent is the two-block ⁠[phi (structured), theta (unstructured)]⁠ convolution; each enters the linear predictor through the field indicator scaled by its d_fac (sigma * sqrt(rho) * scale_factor for phi, sigma * sqrt(1 - rho) for theta), with phi carrying the augmented ICAR precision .icar_precision_Q and theta ~ N(0, I). The conditional Laplace kernel hardcodes sigma = 1, rho = 0.5.

Description

Marginal H_beta for a BYM2 field. The latent is the two-block ⁠[phi (structured), theta (unstructured)]⁠ convolution; each enters the linear predictor through the field indicator scaled by its d_fac (sigma * sqrt(rho) * scale_factor for phi, sigma * sqrt(1 - rho) for theta), with phi carrying the augmented ICAR precision .icar_precision_Q and theta ~ N(0, I). The conditional Laplace kernel hardcodes sigma = 1, rho = 0.5.

Usage

.marginal_H_beta_bym2(
  mode,
  X,
  spatial,
  family,
  phi,
  n_trials,
  weights = NULL,
  offset = NULL,
  sigma_spatial = 1,
  rho = 0.5,
  re_idx = NULL,
  n_re_groups = 0L,
  sigma_re = 1
)

Negate an additive term (numeric literal -> its negative; anything else -> a unary-minus call), for flattening binary - in a bar LHS.

Description

Negate an additive term (numeric literal -> its negative; anything else -> a unary-minus call), for flattening binary - in a bar LHS.

Usage

.negate_term(t)

Build NNGP precision Lambda = (I - A)' D^-1 (I - A) at given hyperparameters.

Description

Mirrors the algorithm in src/gpu_nngp_laplace.h::batch_nngp_scatter plus apply_nngp_full_prior_dense – for each Vecchia row i (in NNGP order):

• resolve conditioning set N(i) via nn_idx / nn_order

• form the n_nb x n_nb cov matrix C and the n_nb cov vector c

• solve C alpha = c

• v_i = sigma^2 - c' alpha

• beta_i is 1 at i and -alpha_k at each neighbor -> Lambda += beta_i beta_i' / v_i

Usage

.nngp_precision_Q(spatial, sigma2_gp, phi_gp)

Details

Identifiers follow the C++ side: nn_idx[i, k] is a 1-based NNGP-order index, nn_order[j] maps NNGP-order j (1-based) to obs idx (1-based).

Normalize an optional fixed-effect Gaussian prior

Description

Validates beta_prior and recycles scalar mean / sd to length p. Returns NULL (use the built-in weak prior) or list(mean, sd) with both vectors of length p. Shared by tulpa_laplace() and the EM driver so the validation rules live in one place.

Usage

.normalize_beta_prior(beta_prior, p)

Arguments

Normalize a zero-inflation prior to a scalar SD

Description

The compiled zero-inflated kernels carry one mean-zero Gaussian prior over the whole beta_zi block (ModelData::zi_prior_sd on the sampler paths, the appended BetaPrior tail on the Laplace path), so the prior is a single SD rather than the per-coefficient list(mean, sd) that beta_prior takes. The list form is kept so the two priors read alike at the front door and so a per-coefficient or non-zero-mean ZI prior can be added without a signature change.

Usage

.normalize_zi_prior(zi_prior)

Arguments

Value

A scalar prior SD; the engine default when zi_prior is NULL.

Interpret an RE term's covariance specification

Description

Maps one re_list element to the two representations the Laplace path needs: pack, the value the C++ kernel consumes (a length-n_coefs marginal-SD vector for a diagonal / uncorrelated term, or a packed lower-triangular Cholesky of length n_coefs (n_coefs + 1) / 2 for a correlated one), and Q, the ⁠n_coefs x n_coefs⁠ RE precision Sigma^{-1} used to build the marginal fixed-effect SE. A correlated term is signalled by r$L (a lower- triangular Cholesky factor, ⁠Sigma = L L'⁠) or r$cov (the covariance matrix); when present these take precedence over r$sigma.

Usage

.re_cov_spec(r)

Arguments

Value

list(pack, Q, diagonal).

n_obs x n_re_groups indicator design for an iid RE block.

Description

re_idx is the 1-based per-observation group index. Observations whose group falls outside 1:n_re_groups contribute no RE column, matching the g >= 0 && g < n_re_groups guard in the C++ kernels.

Usage

.re_design(re_idx, n_re_groups, n_obs)

Resolve a warm_start request into sampler init / inverse-mass vectors

Description

Accepts the front door's warm_start: "eb" or "laplace" to run that fit first, or an already-fitted object from either. Returns NULL when no warm start was asked for, so the caller passes nothing and the kernel keeps its own defaults.

Usage

.resolve_warm_start(warm_start, args, re_terms, sigma_re, beta_prior, n_chains)

Arguments

Schur complement of the latent block out of the joint Hessian.

Description

Given the fixed-effect design X, the combined latent design D (the iid RE indicator block stacked with the spatial-field design), the latent prior precision Q_latent, and the GLM weights W, returns the marginal fixed- effect precision ⁠H_beta = X'WX - (X'WD) (D'WD + Q_latent)^{-1} (X'WD)'⁠. Shared by the SPDE and NNGP marginal-SE paths.

Usage

.schur_H_beta(X, D, Q_latent, W)

Redirect a backend selection: swap the backend and restamp mode / tier / tier name from the registry, recording the reason shown on the fit.

Description

A redirect off an EXPLICIT request (anything but mode = "auto") is an override, not a resolution: the fit does what the model structure requires instead of what the caller asked for. That is recorded on the selection – sel$overridden for callers, and a clause appended to the reason the fit reports – so an overridden fit is distinguishable from one that was never asked for a mode at all. Only the FIRST override is recorded: with several redirects in a chain, the request the user actually made is the one worth naming, not the intermediate backend a previous redirect chose.

Usage

.sel_redirect(sel, backend, reason, notify = TRUE)

Details

notify says whether tulpa() should additionally WARN. It is for the case where the requested backend could have fitted the model and the redirect takes a capability away – a smoother sending an eb / agq request to the nested kernels loses the random-effect SD those two would have estimated. It is off where the requested mode is not expressible for the structure at all and the redirect is the documented resolution rather than a loss: a random slope has no scalar sigma_re for mode = "laplace" to condition on, and an SPDE field redirected from nested_laplace to spde is the same mode and tier reaching its own integrator. Those are recorded, not warned about, so a documented route does not warn on every fit.

Assemble the fractional rSPDE precision and obs map at a given (range, sigma)

Description

Shared by every fractional-nu fit path (Laplace single-point, nested, NUTS, marginal SEs). Maps ⁠(range, sigma)⁠ to the SPDE hyperparameters ⁠kappa = sqrt(8 nu) / range⁠, builds the latent precision ⁠Q = Pl' C^{-1} Pl⁠ and field shift Pr (field ⁠u = Pr x⁠) from the validated R oracle .spde_rational_assemble(), then normalizes the field to marginal variance sigma^2. The rational construction is correct in spectral shape but its overall scale carries kappa-dependent constants (the l_max normalization, the per-factor conditioning rescalings); without the variance normalization the implied field variance is not sigma^2 and varies with the range, which biases the nested ⁠(range, sigma)⁠ integration. The normalization rescales Pr (hence ⁠A_eff = A Pr⁠) by ⁠sigma / sqrt(mean marginal variance)⁠; Q is untouched, so logdet_Q remains the correct prior normalizer for the variance-normalized model.

Usage

.spde_assemble_at(spatial, range, sigma, order = 2L)

Arguments

Value

A list with Q, Pr, A_eff, Pl (all CSC, sized to the non-orphan submesh), keep (1-based indices of the retained mesh nodes), n_mesh_full (the full mesh size), kappa, var_scale, l_max, and logdet_Q.

Build the sparse FEM stiffness / projector from a spatial_spde spec

Description

Prefers the stored Matrix objects (spatial$G, spatial$A); falls back to rebuilding them from the pre-extracted CSC slots so a spec carrying only the slots still assembles.

Usage

.spde_fem_matrices(spatial)

Fractional rSPDE single-point Laplace fit at a fixed (range, sigma)

Description

The fractional counterpart of the integer cpp_laplace_fit_spde branch in laplace_spde_at(). Assembles ⁠(Q, A_eff)⁠ via .spde_assemble_at() and runs the precomputed C++ solve, whose latent is the auxiliary weights x; the returned mode mesh block is mapped back to the field ⁠u = Pr x⁠ so every downstream consumer reads field-space mesh effects exactly as on the integer path. The log-marginal is the well-conditioned B / determinant-lemma marginal (cpp_spde_fractional_logmarginal()) for the RE-free case (comparable across the ⁠(range, sigma)⁠ grid); a fit with an iid RE block keeps the precomputed precision-space marginal.

Usage

.spde_laplace_fractional_at(
  y,
  n_trials,
  X,
  spatial,
  family,
  phi,
  range,
  sigma,
  re_idx,
  n_re_groups,
  sigma_re,
  max_iter,
  tol,
  n_threads,
  offset,
  weights = NULL,
  order = 2L
)

Mean marginal variance of the rSPDE field u = Pr x, x ~ N(0, Q^-1)

Description

Estimates ⁠mean_i [Pr Q^{-1} Pr']_ii = tr(Pr Q^{-1} Pr') / n⁠ by Hutchinson probing: for z ~ N(0, I), ⁠a = Pr' z⁠, ⁠Q v = a⁠, then ⁠E[a' v] = tr(...)⁠. The solve is against Q through the SAME sparse Cholesky the precomputed C++ fit uses, so the normalization is consistent with the fit even when the wide rational spectrum makes Q ill-conditioned: a shared solver makes the implied field covariance identical between the normalization and the likelihood, which is what the nested ⁠(range, sigma)⁠ integration needs (an inconsistent solver breaks the cross-grid marginal). The probe matrix is fixed across calls for a deterministic, grid-smooth normalization.

Usage

.spde_mean_marginal_var(Q, Pr, C0, n_probe = .SPDE_VARNORM_NPROBE)

Numerically stable Laplace log-marginal for a fractional rSPDE at (range, sigma)

Description

Delegates the well-conditioned B / matrix-determinant-lemma marginal to C++ (cpp_spde_fractional_logmarginal()). The precision-space ⁠0.5(log|Q| - log|H|)⁠ is corrupted in a range-dependent way by the rational precision's wide spectrum (cond(Q) ~ 1e13+), so the marginal is formed through the obs-space ⁠B = (A_eff Pl^{-1}) C (A_eff Pl^{-1})' + X X'/tau_beta⁠, built through the operator factor Pl (cond = sqrt cond(Q)), never an explicit Q inverse. Gaussian is the exact conjugate marginal; non-gaussian uses the det-lemma at the precomputed Laplace mode. phi is the Gaussian residual SD (variance phi^2), consistent with the integer path and the engine family convention.

Usage

.spde_nested_logmarginal_at(
  spatial,
  range,
  sigma,
  y,
  X,
  family,
  phi,
  n_trials,
  order,
  max_iter,
  tol,
  n_threads,
  offset,
  tau_beta = 1e-04
)

Value

A list with log_marginal, n_iter, converged, and the assembly.

Is an SPDE smoothness nu fractional (rational path) or integer (exact path)?

Description

Is an SPDE smoothness nu fractional (rational path) or integer (exact path)?

Usage

.spde_nu_is_fractional(nu)

Arguments

Value

TRUE when alpha is non-integer (the rational SPDE path).

Backend names belonging to a tier (derived from the registry).

Description

Backend names belonging to a tier (derived from the registry).

Usage

.tier_backends(tier_key)

Reject families the compiled kernels do not implement.

Description

Reject families the compiled kernels do not implement.

Usage

.validate_family_compiled(family)

Validate the response y for a count family.

Description

Errors when family is a count family and y carries negative or non-integer values, which the integer-casting kernels would silently floor into a biased likelihood. A no-op for continuous families.

Usage

.validate_family_counts(family, y, zi = FALSE)

Arguments

Validate the dispersion parameter phi for a family.

Description

Errors when family carries a dispersion / precision parameter and phi is not a positive finite scalar. A no-op for families without dispersion (binomial, poisson) and for unknown / model-package families. Shared by the tulpa() front door and tulpa_laplace() so the rule lives in one place .

Usage

.validate_family_phi(family, phi)

Validate the response against a family's support.

Description

The one support rule set, shared by the tulpa() front door and by the engine fitters (tulpa_laplace(), tulpa_laplace_beta(), tulpa_nuts_beta()). Holding it in a single function is what keeps a boundary response rejected with the same message wherever it enters: a per-fitter copy is only as good as the set of fitters that carry it, and tulpa() reaches beta-with-RE through a sampler that carried none.

Usage

.validate_family_support(family, y, n_trials = NULL, zi = FALSE)

Arguments

Details

Counts delegate to .validate_family_counts(); the continuous families and the binomial denominator add their rules here.

Reject a ZI request against a family with no atom at zero.

Description

Reject a ZI request against a family with no atom at zero.

Usage

.validate_family_zi(family)

Reject a ZI request the compiled kernels cannot fit.

Description

Reject a ZI request the compiled kernels cannot fit.

Usage

.validate_family_zi_compiled(family)

Validate a GLM design bundle (y, X, n_trials) at a fitter's front door.

Description

Shared by the flagship drivers (tulpa_nested_laplace(), tulpa_gibbs(), the re_cov fitters) so nrow(X) == length(y) and length(n_trials) == length(y) are enforced in one place – otherwise a mismatched n_trials (as.integer(NULL) -> integer(0)) reaches the C++ kernel silently.

Usage

.validate_glm_design(y, X, n_trials, where)

Arguments

Value

List list(N, n_trials) with n_trials coerced to a length-N integer vector.

Validate a 0-based/1-based random-effect index against its group count.

Description

re_idx addresses one grouping factor: 0 marks "no random effect" and ⁠1..n_re_groups⁠ a group. An out-of-range id would index out of bounds in the C++ kernel, so bound it in R.

Usage

.validate_re_idx(re_idx, n_re_groups, N, where)

Arguments

Value

re_idx coerced to an integer vector.

Validate the SPDE smoothness parameter

Description

nu must be a single finite positive number. Integer nu (1, 2, 3, ...) gives an exact FEM construction; fractional nu uses the rational SPDE approximation (see rational_spde_coefficients()).

Usage

.validate_spde_nu(nu)

Arguments

Value

Invisibly TRUE; raises an error on an invalid nu.

Split a Laplace / EB mode vector into its semantic blocks

Description

Returns the count-side fixed effects, the zero-inflation coefficients, and the per-term random-effect deviations, so the caller places each into the sampler's corresponding block rather than relying on the two flat vectors agreeing end to end (they do not, under zero inflation).

Usage

.warm_start_blocks(fit, re_terms)

Arguments

Per-term random-effect SDs from a warm-start source fit

Description

An EB fit estimated them, so they are read off map. A Laplace fit conditioned on values the caller supplied, so those are passed in.

Usage

.warm_start_sigmas(fit, sigma_re, n_terms)

Arguments

Build the zero-inflation design matrix from ziformula.

Description

A one-sided formula giving the fixed effects of the structural-zero logit. ~ 1 is the constant-probability model. Random effects in the ZI predictor are not supported: the compiled kernels share the random-effect block into the count process only, so a bar here would be silently dropped.

Usage

.zi_design(ziformula, data, n_obs)

All registered backends across every tier (derived from the registry).

Description

All registered backends across every tier (derived from the registry).

Usage

ALL_BACKENDS

Backend -> supported families, derived from the registry.

Description

Backend -> supported families, derived from the registry.

Usage

BACKEND_FAMILY_SUPPORT

Backend registry – single source of truth for the inference backends.

Description

One entry per backend. Adding a backend is a single entry here; tier membership, family support, the input contract, and R-level reachability all derive from this list – the single source of truth for per-tier backends and family support.

Fields:

• emits – the kind of posterior representation the backend's draws hold, which fixes whether MCMC chain diagnostics apply (independent of tier: an exact SMC sampler emits "iid" particles, a Tier-3 VI fit also emits "iid", while a nested-Laplace Tier-2 fit emits "iid" too):

  • "chain" : autocorrelated MCMC output – Rhat and autocorrelation-ESS are meaningful (diagnostics() computes them).

  • "iid" : exchangeable draws from a deterministic approximation, or resampled particles – split-Rhat is vacuous and ESS = n by construction, so they say nothing about approximation bias.

  • "point" : no posterior sample (mode + covariance only).

• tier – tier key ("exact", "structured", "optimized").

• input – the input contract the backend's fitter consumes:

  • "design" : design-matrix bundle (y, n_trials, X, RE structure).

  • "logpost" : a log_posterior(theta) closure plus dimension/init.

  • "modeldata" : a tulpa ModelData / LikelihoodSpec (C-ABI NUTS path).

  • "nested" : a design-matrix bundle plus one or more latent prior blocks (latent(tgmrf(...))), integrated over the block hyperparameters by the nested-Laplace driver.

  • "spde" : a design bundle plus a self-contained SPDE spec (mesh + FEM matrices); fit_spde()'s own CCD / grid hyperparameter engine.

• fitter – name of the R function implementing the backend, or NULL when only a C++ kernel exists with no R entry point yet. Stored as a string (not the function object) so the registry is independent of source-load order; resolved lazily via resolve_backend_fitter(). NULL => not selectable from R; dispatch fails loudly.

• families – character vector of supported family identifiers, or NULL for an unrestricted backend.

• cabi – the registered C-ABI callable backing the backend (the symbol a model package reaches via LinkingTo: tulpa, and the one an R wrapper would call), or NULL.

• note – optional human-readable note.

Family identity (for families) is checked against family$name, family$distribution, and family$numerator$distribution (tulpaRatio-style ratio families nest a per-process distribution).

Usage

BACKEND_REGISTRY

Inference tiers with their properties (derived from TIER_META + BACKEND_REGISTRY). Preserves the ⁠INFERENCE_TIERS$<tier>$backends⁠ interface relied on by downstream code and tests.

Description

Inference tiers with their properties (derived from TIER_META + BACKEND_REGISTRY). Preserves the ⁠INFERENCE_TIERS$<tier>$backends⁠ interface relied on by downstream code and tests.

Usage

INFERENCE_TIERS

Tier metadata: the epistemic guarantee attached to each tier.

Description

Tier membership (which backends live in a tier) is not stored here – it is derived from BACKEND_REGISTRY. This table holds only the per-tier description/guarantee so the two never drift.

Usage

TIER_META

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

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

Convert adjacency matrix to CSR for tulpa C++

Description

Convert adjacency matrix to CSR for tulpa C++

Usage

adjacency_to_csr_tulpa(adj)

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

Apply RSR projection to spatial effect

Description

Project spatial effect into the space orthogonal to covariates. Called during posterior computation.

Usage

apply_rsr_projection(w, P_perp)

Arguments

Value

Projected spatial effect (length n)

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

Error if a selected backend has no R-level fitter.

Description

Enforces the registry honesty contract: a backend may ship a C++ kernel reachable from model packages (via LinkingTo: tulpa) yet have no R entry point. Such a backend must never be silently selectable from R – selecting it errors with a precise message naming the C-ABI symbol, rather than pretending to dispatch.

Usage

assert_backend_reachable(backend)

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.

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

• 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.

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)

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.

See Also

is_auto_grid(), 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-select mode (Tier 1 or Tier 2 only)

Description

Implements the "auto" mode selection. Chooses the most reliable method that is expected to finish for the given model.

Critical rule: Auto never selects Tier 3 (Optimized).

Usage

auto_select_mode(
  family,
  n_obs,
  has_spatial,
  has_temporal,
  has_latent,
  temporal = NULL,
  spatial_type = NULL,
  has_re = FALSE
)

Is a backend reachable from R (does it have an R-level fitter)?

Description

Is a backend reachable from R (does it have an R-level fitter)?

Usage

backend_is_reachable(backend)

Test whether a backend supports the given family.

Description

Test whether a backend supports the given family.

Usage

backend_supports_family(backend, family)

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)

Build a GLMM log-posterior (and gradient) from a model-data bundle.

Description

Build a GLMM log-posterior (and gradient) from a model-data bundle.

Usage

build_glmm_logpost(
  bundle,
  family,
  sigma_re = NULL,
  n_trials = NULL,
  phi = 1,
  beta_prior = .tulpa_default_beta_prior("glmm_logpost"),
  weights = NULL,
  phi2 = NULL
)

Arguments

Value

A list with:

• log_posterior(theta) – scalar log-posterior (up to a constant).

• grad_log_posterior(theta) – gradient vector.

• dim – length of theta.

• init – a zero starting vector of length dim.

• unpack(theta) – list(beta, u = list of n_groups x n_coefs matrices).

• param_names – character labels aligned with theta.

Check if data is gridded (can use stars)

Description

Check if data is gridded (can use stars)

Usage

can_use_stars(plot_data)

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.

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.

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

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, symmetric, zero on the diagonal, 0/1 valued, and free of isolated nodes. adjacency() runs the same checks on the graphs it constructs.

Usage

check_adjacency(adjacency, ids = NULL)

Arguments

Value

Invisibly, a tulpa_adjacency_check list with the per-check results (square, symmetric, zero_diag, binary, isolated-node indices, 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)

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

Collect additive terms from a + chain

Description

Recursively flattens a + b + c into list(a, b, c).

Usage

collect_additive_terms(expr)

Arguments

Value

A list of language objects (individual terms)

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.

Usage

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

Arguments

Value

A data frame. For "loglik": model, n_params, logLik. 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 BYM2 scaling factor

Description

Compute the scaling factor for BYM2 following Riebler et al. (2016). This makes the spatial fraction parameter interpretable.

Usage

compute_bym2_scale(adjacency)

Arguments

Value

Scaling factor (scalar)

Compute valid bounds for rho in proper CAR

Description

For proper CAR, rho must be in the range (1/lambda_min, 1/lambda_max) where lambda are the eigenvalues of D^(-1)W. In practice, we typically restrict to (0, 1) for interpretability (positive spatial autocorrelation).

Usage

compute_car_rho_bounds(adjacency)

Arguments

Value

Named vector with lower and upper bounds for rho

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)

Compute RSR projection matrix

Description

Compute the orthogonal projection matrix P_perp = I - P_X that projects the spatial effect into the space orthogonal to the covariates.

Usage

compute_rsr_projection(X)

Arguments

Value

Projection matrix (n x n)

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

Create the ggplot2 map

Description

Create the ggplot2 map

Usage

create_map(
  plot_data,
  obs_coords,
  title,
  palette,
  point_color,
  point_size,
  na_color,
  crs,
  legend_title,
  ...
)

Create point-based map

Description

Create point-based map

Usage

create_point_map(plot_data, palette, na_color, legend_title)

Create raster-style map using stars

Description

Create raster-style map using stars

Usage

create_raster_map(plot_data, crs, palette, na_color, legend_title)

DIC and CPO

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

dic(object, ...)

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

cpo(object, ...)

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

Arguments

Details

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.

Value

A tulpa_criteria object.

See Also

tulpa_criteria() for every criterion at once and for what the LOO unit means.

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)

Decompose the LHS of a bar term into intercept flag + slope language objects

Description

Walks the additive chain of the LHS and separates numeric intercept indicators (0, 1, -1) from slope term expressions. Returns language objects, not deparsed strings.

Usage

decompose_bar_lhs(lhs)

Arguments

Value

A list with has_intercept (logical) and slope_terms (list of language objects)

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)

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

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 laplace_diagnostics() 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 '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 that function 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.

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.

See Also

sbc() for calibration, laplace_diagnostics() for the approximation-reliability table in full, tulpa_draws_array(), 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)

Dispatch a spatial Polya-Gamma Gibbs fit to the correct sampler

Description

The Gibbs analogue of dispatch_laplace_spatial(): routes on spatial$type to the matching ⁠cpp_pg_<family>_gibbs_<structure>⁠ sampler, building the neighbour-list / coordinate inputs each one needs. The binomial Polya-Gamma augmentation backs the full areal (icar/bym2/rsr) + continuous (gp/nngp/ multiscale_gp) family; neg_binomial_2 is backed by the single areal ICAR negbin sampler (cpp_pg_negbin_gibbs_spatial), the only negbin spatial kernel.

Usage

dispatch_gibbs_spatial(
  y,
  n_trials,
  X,
  re_group,
  n_re_groups,
  spatial,
  family,
  iter,
  warmup,
  thin = 1L,
  prior_beta_sd = .tulpa_prior_sd("gibbs"),
  prior_sigma_re_scale = 2.5,
  verbose = FALSE,
  n_threads = 1L
)

Dispatch a temporal Polya-Gamma Gibbs fit

Description

The temporal analogue of dispatch_gibbs_spatial(): maps a validated temporal_multiscale() spec onto the multiscale temporal Polya-Gamma sampler (cpp_pg_binomial_gibbs_temporal), which composes an additive RW1 trend + cyclic-RW1 seasonal + AR1/IID short-term decomposition. Binomial only. The C++ kernel implements an RW1 trend (rw2 is rejected here rather than silently downgraded).

Usage

dispatch_gibbs_temporal(
  y,
  n_trials,
  X,
  re_group,
  n_re_groups,
  temporal,
  family,
  iter,
  warmup,
  thin = 1L,
  prior_beta_sd = .tulpa_prior_sd("gibbs"),
  prior_sigma_re_scale = 2.5,
  verbose = FALSE,
  n_threads = 1L
)

Dispatch spatial Laplace to the correct C++ backend

Description

weights is the per-observation likelihood weight. It reaches the same BuiltinFamilyResponse::weights channel the non-spatial route uses, which scales each row's log-density, score and Fisher curvature by the same w_i, so the mode these kernels return and the marginal precision ⁠.marginal_H_beta_*()⁠ builds at it describe one model.

Usage

dispatch_laplace_spatial(
  y,
  n_trials,
  X,
  re_idx,
  n_re_groups,
  sigma_re,
  spatial,
  family,
  phi,
  max_iter,
  tol,
  n_threads,
  offset = NULL,
  weights = NULL
)

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

Arguments

Value

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

Extract coordinates from various formats

Description

Extract coordinates from various formats

Usage

extract_coords(x)

Extract coordinates from tulpa_fit object

Description

Extract coordinates from tulpa_fit object

Usage

extract_coords_from_fit(fit)

Elementwise log-likelihood for a family.

Description

Elementwise log-likelihood for a family.

Usage

family_loglik(eta, y, family, n_trials = NULL, phi = 1, phi2 = NULL)

Inverse-link mean for a family. phi reaches entries whose response mean depends on the dispersion (lognormal); the others ignore it via ....

Description

Inverse-link mean for a family. phi reaches entries whose response mean depends on the dispersion (lognormal); the others ignore it via ....

Usage

family_mean(eta, family, phi = 1)

Supported R-level family names.

Description

Supported R-level family names.

Usage

family_names()

Response-scale mean of y given eta for a family (trial-scaled where relevant).

Description

Response-scale mean of y given eta for a family (trial-scaled where relevant).

Usage

family_response_mean(eta, family, n_trials = NULL, phi = 1)

One response draw per element of eta (posterior predictive), elementwise.

Description

One response draw per element of eta (posterior predictive), elementwise.

Usage

family_sample(eta, family, n_trials = NULL, phi = 1, phi2 = NULL)

Score (d log-likelihood / d eta), elementwise.

Description

Score (d log-likelihood / d eta), elementwise.

Usage

family_score_eta(eta, y, family, n_trials = NULL, phi = 1, phi2 = NULL)

Response variance Var(y | eta) for a family, elementwise.

Description

Response variance Var(y | eta) for a family, elementwise.

Usage

family_variance(eta, family, n_trials = NULL, phi = 1, phi2 = NULL)

Laplace/IRLS working weight (-d^2 log-lik / d eta^2), elementwise.

Description

Laplace/IRLS working weight (-d^2 log-lik / d eta^2), elementwise.

Usage

family_weight(eta, family, n_trials = NULL, phi = 1, phi2 = NULL)

Find all latent(...) calls in a formula's parse tree

Description

Recursively walks the formula AST and collects all latent(...) calls. The matched calls are returned unevaluated; tulpa_parse_formula() resolves them in the formula's environment.

Usage

find_latent_terms(term)

Arguments

Value

A list of language objects, each a latent(...) call

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

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

• 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 an areal spatial field (icar / bym2 / car_proper) and a temporal field (rw1 / rw2 / ar1), integrating the spatial precision, 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"),
  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,
  control = list()
)

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, and theta_grid over ⁠(tau_spatial, tau_temporal, rho)⁠. 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").

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)

Format mode selection message

Description

Format mode selection message

Usage

format_mode_selection(selection)

Arguments

Value

Character string for display

Render a slope spec as a display string

Description

Used by the print method. Derives the textual form on demand from the stored language objects so we never store deparsed slope text.

Usage

format_re_lhs(re)

Format tier information for display

Description

Format tier information for display

Usage

format_tier_info(tier_info, verbose = FALSE)

Arguments

Value

Character string

Gauss-Hermite quadrature (probabilist's, exp(-z^2/2) weight)

Description

Computes nodes and weights for ⁠\int f(z) (2\pi)^{-1/2} exp(-z^2/2) dz⁠ via Golub-Welsch on the Hermite Jacobi matrix. Sums of ⁠w_k f(z_k)⁠ approximate ⁠E_{Z ~ N(0,1)}[f(Z)]⁠.

Usage

gauss_hermite_prob(n)

Arguments

Value

List with nodes and weights (each length n).

Get tier for a backend

Description

Get tier for a backend

Usage

get_backend_tier(backend)

Arguments

Value

List with tier information

Map mode to valid backends

Description

Map mode to valid backends

Usage

get_mode_backends(mode)

Arguments

Value

Character vector of valid backends for this mode

Get appropriate color scale

Description

Get appropriate color scale

Usage

get_palette_scale(palette, na_color, legend_title, geom = "fill")

Geweke Convergence Test

Description

Performs Geweke's convergence diagnostic, comparing the mean of the first portion of a chain to the last portion. Useful for single-chain diagnostics.

Usage

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

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.

Value

A data frame with columns: parameter, z_score, p_value.

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)

Compute GLM working weights for Laplace Hessian

Description

Thin wrapper over the family-ops registry (family_weight()) so the weight formulas live in exactly one place (R/family_loglik.R).

Usage

glmm_weights(eta, family, n_trials = NULL, phi = 1, phi2 = NULL)

Map a spatial_gp covariance spec to the Laplace cov_type integer

Description

The Laplace NNGP kernel (laplace_core.cpp) supports three covariance functions: 0 = exponential, 1 = Matern(nu=1.5), 2 = Matern(nu=2.5). Anything else is rejected with a clear error rather than silently falling back to a different covariance.

Usage

gp_cov_type_for_laplace(spatial)

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
)

Arguments

Value

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

See Also

tulpa_hyper_grid().

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

Inference Mode System

Description

tulpa uses an explicit tier system for inference that encodes epistemic guarantees, not just runtime characteristics.

This design makes the difference between inference methods first-class and unavoidable, rather than hiding them as implementation details.

The Three Tiers

Tier 1 - Exact: Asymptotically correct posterior inference (up to Monte Carlo error). Credible intervals are interpretable as posterior uncertainty. This is the reference standard.

Tier 2 - Structured: Accurate inference conditional on explicit structural assumptions. Typically requires latent Gaussian structure, conditional independence, smooth posteriors. Very fast for the right model class, but can be wrong outside that class. Failure modes are predictable and explainable.

Tier 3 - Optimized: No general correctness guarantee beyond empirical usefulness. Point estimates usually good, but uncertainty often underestimated. Tails and correlations unreliable. Failure is usually silent. This is optimization, not sampling.

Auto Mode

mode = "auto" chooses between Tier 1 and Tier 2 only. It will never silently choose Tier 3 (Optimized).

The contract: "Use the most reliable method that is expected to finish for this model."

Auto decisions are deterministic, explainable, and overrideable.

Implementation Rules

1. Modes change semantics, not just runtime. Intervals from Optimized do not mean the same as Exact.

2. The mode must always be visible in output.

3. No silent upgrading or downgrading. If Exact fails, we error - we do not switch to Structured.

4. Backends slot into tiers. Adding a backend never introduces a new epistemic promise.

Initialize latent factor parameters

Description

Initialize latent factor parameters

Usage

initialize_latent_params(latent_info, seed = NULL)

Arguments

Value

Numeric vector of initial parameter values

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.

See Also

auto_grid()

Examples

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

Check if adjacency matrix is connected

Description

Check if the spatial graph defined by the adjacency matrix is fully connected. A disconnected graph can cause identifiability issues.

Usage

is_connected(adjacency)

Arguments

Value

Logical; TRUE if connected

Proper-CAR Laplace at given hyperparameters

Description

Single-point Laplace for a proper-CAR areal field at a fixed ⁠(tau, rho)⁠, the conditional counterpart of the nested CAR_proper integrator. Reuses the shared make_car_proper_latent_blocks factory + dense spec solver via cpp_laplace_fit_car_proper, so the mode + log-marginal equal the nested kernel at that one grid cell. tau / rho default to the spec's fields (or tau = 1, rho = midpoint of the eigenvalue-derived rho_bounds) and are recorded on the result.

Usage

laplace_car_proper_at(
  y,
  n_trials,
  X,
  spatial,
  family = "binomial",
  phi = 1,
  tau = NULL,
  rho = NULL,
  re_idx = NULL,
  n_re_groups = 0L,
  sigma_re = 1,
  max_iter = 100L,
  tol = 1e-06,
  n_threads = 1L,
  offset = NULL,
  weights = NULL
)

Approximation-reliability diagnostics for a deterministic nested-Laplace fit

Description

Use diagnostics(), which returns this table for any fit whose draws are an i.i.d. approximation sample. The sections below document that table; they remain the reference for its columns and attributes.

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 that diagnostics() reports for a chain fit does not apply. This is the accessor 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.

Usage

laplace_diagnostics(fit, pars = NULL)

Arguments

Value

A data frame with one row per parameter – parameter, mean, sd, ess_bulk, ess_tail, rhat – carrying attributes:

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.

ess_grid, n_grid, rel_ess_grid, max_weight

grid quadrature reliability.

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.

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 absent on a spread grid.

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.

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.

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.

References

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.

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.

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

diagnostics() (the front door, which returns this table for i.i.d. fits), tulpa_psis().

Examples

set.seed(1)
n <- 200L; x <- rnorm(n)
y <- rbinom(n, 1, plogis(-0.2 + 0.6 * x))
# `mode = "laplace"` returns a mode + covariance and carries no draws; a
# sampled deterministic backend is what this table describes.
fit <- tulpa(y ~ x, data.frame(y = y, x = x), family = "binomial",
             mode = "smc")
diagnostics(fit)

NNGP Laplace at given hyperparameters

Description

Single-point Laplace approximation for a Matern/exponential GP spatial field at fixed (sigma2_gp, phi_gp). Used by dispatch_laplace_spatial when spatial$type == "gp". The neighbor structure is read straight off the validated spec – call validate_gp(spatial, data) first if constructing manually.

Usage

laplace_gp_at(
  y,
  n_trials,
  X,
  spatial,
  family = "binomial",
  phi = 1,
  sigma2_gp = NULL,
  phi_gp = NULL,
  re_idx = NULL,
  n_re_groups = 0L,
  sigma_re = 1,
  max_iter = 100L,
  tol = 1e-06,
  n_threads = 1L,
  offset = NULL,
  weights = NULL
)

Arguments

Value

The raw cpp_laplace_fit_gp result list, augmented with sigma2_gp, phi_gp, and the spatial spec.

HSGP Laplace at given hyperparameters

Description

Single-point Laplace for a Hilbert-space GP field at a fixed ⁠(sigma2, lengthscale)⁠, the conditional counterpart of the nested HSGP integrator. Reuses the shared make_hsgp_block factory + dense spec solver (DENSE_BASIS scatter) via cpp_laplace_fit_hsgp, so the mode equals the nested kernel at that one grid cell. sigma2 / lengthscale default to the spec's fields (or 1) and are recorded on the result.

Usage

laplace_hsgp_at(
  y,
  n_trials,
  X,
  spatial,
  family = "binomial",
  phi = 1,
  sigma2 = NULL,
  lengthscale = NULL,
  re_idx = NULL,
  n_re_groups = 0L,
  sigma_re = 1,
  max_iter = 100L,
  tol = 1e-06,
  n_threads = 1L,
  offset = NULL,
  weights = NULL
)

SPDE Laplace at given hyperparameters

Description

Single-point Laplace approximation for an SPDE spatial field at a fixed (range, sigma). Used by both dispatch_laplace_spatial (single-point path) and fit_spde (single-point branch) so the call site stays a single source of truth.

Usage

laplace_spde_at(
  y,
  n_trials,
  X,
  spatial,
  family = "binomial",
  phi = 1,
  range = NULL,
  sigma = NULL,
  re_idx = NULL,
  n_re_groups = 0L,
  sigma_re = 1,
  max_iter = 100L,
  tol = 1e-06,
  n_threads = 1L,
  offset = NULL,
  weights = NULL
)

Arguments

Value

The raw cpp_laplace_fit_spde result list (mode, log_det_Q, log_marginal, n_iter, converged), augmented with range, sigma, and the spatial spec for downstream prediction.

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.

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
)

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

## Not run: 
# `fit` is a ratio / multi-arm model fitted with latent factors by a consumer
# package (tulpaRatio owns the two-arm formula and the negbin/negbin ratio
# family and reads the latent_factor() spec):

# Get summary
factors <- latent_factors(fit)
head(factors)

# Get full posterior draws
factor_draws <- latent_factors(fit, summary = FALSE)

## End(Not run)

Links available as a ⁠<family>_<link>⁠ suffix.

Description

Links available as a ⁠<family>_<link>⁠ suffix.

Usage

link_names()

Family codes the link layer accepts beyond the bare registry names.

Description

The cross product of the link-capable base families and the links, minus each family's own canonical form (which is already a registry name).

Usage

linked_family_names()

Log-likelihood at the posterior mean

Description

Log-likelihood at the posterior mean

Usage

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

Arguments

Value

A logLik object.

Numerically stable log(mean(exp(x)))

Description

Numerically stable log(mean(exp(x)))

Usage

logmeanexp(x)

Numerically stable log(exp(a) + exp(b)) elementwise (vectors / scalars)

Description

Numerically stable log(exp(a) + exp(b)) elementwise (vectors / scalars)

Usage

logsumexp_pair(a, b)

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,
  thin = 1L,
  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

Marginal fixed-effect Hessian for spatial-field Laplace fits

Description

The raw fixed-effect block of the joint Hessian gives the conditional precision on \beta \mid u^*, which under-states uncertainty; this returns the marginal precision instead.

Details

The correct marginal precision comes from a Schur complement on the joint Hessian at the mode:

H_\beta^{\mathrm{marg}} = X'WX - X'WZ (Z'WZ + Q_u)^{-1} Z'WX

where Z is the obs->latent map (SPDE: projection matrix A; NNGP: indicator from obs to unique-location field) and Q_u is the spatial precision at the fitted hyperparameters.

These helpers rebuild the spatial precision in R from the spec and fitted hyperparameters, then solve via sparse Cholesky. The shape matches what cpp_laplace_fit_spde / cpp_laplace_fit_gp use internally – Q construction here is a 1:1 port of spde_qbuilder.h (orphan ridge included).

MCMC convergence diagnostics

Description

Use diagnostics(), which reads a fit's draws provenance and returns the diagnostic that applies – chain mixing for MCMC draws, approximation reliability for deterministic fits. The name mcmc_diagnostics() described only one of the two branches it already routed between.

Usage

mcmc_diagnostics(
  fit,
  pars = NULL,
  measures = c("rhat", "ess_bulk", "ess_tail"),
  probs = c(0.05, 0.95)
)

Arguments

Value

The value of diagnostics() for fit.

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

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 (⁠$diagnostics$n_divergent⁠, ⁠$diagnostics$divergent_idx⁠, ⁠$diagnostics$divergent⁠, or the top-level ⁠$n_divergent⁠ / ⁠$divergent⁠).

Usage

n_divergent(fit)

Arguments

Value

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

See Also

plot_divergences(), check_diagnostics()

Remove all latent(...) calls from a formula's parse tree

Description

Mirror of nobars() for latent(...) calls.

Usage

no_latent_terms(term)

Arguments

Value

A language object with all latent(...) calls removed, or NULL if nothing remains.

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

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

Parse a single random effect bar term

Description

Takes a | or || language object and extracts the grouping variable(s), effect terms (intercept, slopes), and correlation structure. The bar operator drives the correlated flag: | -> TRUE (LKJ-Cholesky on the joint slope vector), || -> FALSE (diagonal covariance, one sigma per coefficient). This matches lme4 / glmmTMB and lets downstream packages branch on correlated to choose between independent-sigma and Cholesky parameterizations.

Usage

parse_bar_term(bar_term)

Arguments

Details

Nested grouping (1 | a/b) is expanded into one spec per level (a, then a:b). || is preserved as a single spec rather than split into multiple | bars, so the original user intent (one logical RE term with diagonal covariance) round-trips through the parsed object.

Value

A list of RE specs (one per grouping level). Each spec has:

• group_var: character, display label (colon-joined for nested)

• group_vars: character vector, each element a column name

• group_expr: language or NULL (set when grouping is a non-name expr)

• slope_terms: list of language objects (slope LHS terms)

• has_intercept: logical

• correlated: logical (TRUE for |, FALSE for ||)

• original: the original bar language object

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,
  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]⁠

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)

Base R ACF plot

Description

Base R ACF plot

Usage

plot_acf_base(draws, pars, lags)

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

Parallel coordinates plot for divergences

Description

Parallel coordinates plot for divergences

Usage

plot_divergences_parcoord(draws_df, pars)

Base R parallel coordinates

Description

Base R parallel coordinates

Usage

plot_divergences_parcoord_base(scaled, pars, divergent)

Scatter plot matrix for divergences

Description

Scatter plot matrix for divergences

Usage

plot_divergences_scatter(draws_df, pars)

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)

Base R energy plot

Description

Base R energy plot

Usage

plot_energy_base(energy, energy_diff, e_bfmi, status)

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

Base R ESS plot

Description

Base R ESS plot

Usage

plot_ess_base(diag, threshold, type)

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)

Base R pairs plot

Description

Base R pairs plot

Usage

plot_pairs_base(draws_df, pars, has_divergent, alpha)

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)

Base R Rhat plot

Description

Base R Rhat plot

Usage

plot_rhat_base(diag, threshold)

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

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 linear predictor is rebuilt per posterior draw (fixed effects, formula random effects, offset, and a posterior-mean SPDE field when present) and pushed through the family's sampling distribution.

Fits carrying posterior draws use them directly (fixed and random effects jointly per draw). The Laplace tier samples the fixed effects from the Gaussian approximation N(coef(fit), vcov(fit)) and holds the random effects at their posterior mode, so its replicates understate the RE posterior uncertainty. At newdata the prediction is population level (random effects at zero), matching predict.tulpa_fit().

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

Single component pp_check

Description

Single component pp_check

Usage

pp_check_single(y, yrep, type, stat, ndraws, title_suffix, ...)

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.

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.

Prepare latent factor data for HMC

Description

Prepare latent factor data for HMC

Usage

prepare_latent_for_hmc(latent, N)

Arguments

Value

List with latent factor data for C++

Prepare map data from data frame

Description

Prepare map data from data frame

Usage

prepare_map_data_from_df(x, coords, what, summary)

Prepare map data from tulpa_fit object

Description

Prepare map data from tulpa_fit object

Usage

prepare_map_data_from_fit(fit, what, summary, newdata, coords = NULL)

Print method for diagnostic summary

Description

Print method for diagnostic summary

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the diagnostic summary to the console.

Print method for Geweke test

Description

Print method for Geweke test

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the Geweke convergence diagnostic to the console.

Print method for tulpa_gp

Description

Print method for tulpa_gp

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the Gaussian-process spatial specification to the console.

Print method for tulpa_hsgp

Description

Print method for tulpa_hsgp

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the Hilbert-space GP spatial specification to the console.

Print method for tulpa_latent

Description

Print method for tulpa_latent

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the latent factor specification to the console.

Print method for tulpa_multiscale

Description

Print method for tulpa_multiscale

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the multi-scale spatial specification to the console.

Print method for nested-Laplace fits

Description

Compact one-screen summary of a tulpa_nested_laplace() or tulpa_nested_laplace_joint() fit: the hyperparameters integrated over, the outer-grid size, the outer Pareto-\hat{k} accuracy diagnostic when present, and the wall-clock timing line ("fit in 5h 25m (grid 2h 09m)") when fit$timing is attached. Inherited by the single-block joint and multi-block joint subclasses.

Usage

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

Arguments

Value

x, invisibly.

Print method for tulpa_prior

Description

Print method for tulpa_prior

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the prior specification to the console.

Print method for tulpa_prior_predict

Description

Print method for tulpa_prior_predict

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing a summary of the prior predictive draws to the console.

Print method for tulpa_priors

Description

Print method for tulpa_priors

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the full prior specification to the console.

Print method for tulpa_rsr

Description

Print method for tulpa_rsr

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the spatial specification and its restricted-spatial-regression modifier to the console.

Print method for tulpa_simulate

Description

Print method for tulpa_simulate

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing a summary of the simulated datasets to the console.

Print method for tulpa_spatial

Description

Print method for tulpa_spatial

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the spatial specification to the console.

Print method for tulpa_svc

Description

Print method for tulpa_svc

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the spatially-varying-coefficient specification to the console.

Print method for tulpa_svc_posterior

Description

Print method for tulpa_svc_posterior

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing a summary of the spatially-varying-coefficient posterior to the console.

Print method for tulpa_temporal

Description

Print method for tulpa_temporal

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the temporal specification to the console.

Print method for tulpa_temporal_gp

Description

Print method for tulpa_temporal_gp

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the Gaussian-process temporal specification to the console.

Print method for tulpa_temporal_multiscale

Description

Print method for tulpa_temporal_multiscale

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the multi-scale temporal specification to the console.

Print method for tulpa_temporal_posterior

Description

Print method for tulpa_temporal_posterior

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing a summary of the temporal-effect posterior to the console.

Print method for tulpa_tvc

Description

Print method for tulpa_tvc

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing the temporally-varying-coefficient specification to the console.

Print method for tulpa_tvc_posterior

Description

Print method for tulpa_tvc_posterior

Usage

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

Arguments

Value

The input x, returned invisibly. Called for the side effect of printing a summary of the temporally-varying-coefficient posterior to the console.

Format a single prior for printing

Description

Thin wrapper around format() for tulpa_prior objects. Each ⁠prior_*()⁠ constructor attaches a subclass (e.g. tulpa_prior_normal) so adding a new distribution is one new ⁠format.tulpa_prior_<dist>()⁠ method – no central if/else to extend.

Usage

print_prior(prior, indent = "")

Arguments

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

Gamma prior

Description

Specify a gamma prior for a positive parameter.

Usage

prior_gamma(shape = 2, rate = 0.1)

Arguments

Details

Mean = shape/rate, Variance = shape/rate^2.

Value

A tulpa_prior object

Examples

prior_gamma(2, 0.1)  # Mean = 20, weakly informative
prior_gamma(1, 1)    # Exponential(1)

Half-Cauchy prior

Description

Specify a half-Cauchy prior for a positive parameter. Has heavier tails than half-normal, often used for variance parameters.

Usage

prior_half_cauchy(scale = 2.5)

Arguments

Value

A tulpa_prior object

Examples

prior_half_cauchy(2.5)

Half-normal prior

Description

Specify a half-normal (truncated at 0) prior for a positive parameter.

Usage

prior_half_normal(sd = 1)

Arguments

Value

A tulpa_prior object

Examples

prior_half_normal(1)

Normal prior

Description

Specify a normal prior for a parameter.

Usage

prior_normal(mean = 0, sd = 2.5)

Arguments

Value

A tulpa_prior object

Examples

prior_normal(0, 2.5)
prior_normal(0, 1)

Penalized complexity (PC) prior

Description

Specify a PC prior for a positive parameter (typically a standard deviation). PC priors shrink toward simpler models by penalizing deviation from a base model.

Usage

prior_pc(U = 1, alpha = 0.01)

Arguments

Details

The PC prior is specified via: P(sigma > U) = alpha

• U is the upper bound you consider "large"

• alpha is the probability of exceeding it

This implies an exponential prior with rate = -log(alpha) / U.

Value

A tulpa_prior object

References

Simpson, D., Rue, H., Riebler, A., Martins, T. G., & Sorbye, S. H. (2017). Penalising model component complexity: A principled, practical approach to constructing priors. Statistical Science, 32(1), 1-28.

Examples

prior_pc(U = 1, alpha = 0.01)    # P(sigma > 1) = 0.01
prior_pc(U = 0.5, alpha = 0.05)  # Tighter, P(sigma > 0.5) = 0.05

Prior predictive simulation

Description

Draw datasets from the prior predictive distribution: parameters are sampled from their priors (no data conditioning) and pushed through the model's linear predictor and the family's simulator.

Useful for checking whether priors imply plausible data ranges before fitting.

Usage

prior_predict(
  formula,
  family,
  data,
  priors = NULL,
  n_draws = 100,
  seed = NULL,
  ...
)

Arguments

Value

A tulpa_prior_predict object: a list with

• y: list of length n_draws, each element the simulated response for that draw (matching whatever shape family$simulate_fn returns).

• theta: list of length n_draws of parameter draws (beta, sigma, RE coefficients u, family-specific extras).

• linpred: list of length n_draws, each a list of linear predictor vectors per process.

• family: the family used.

• n_draws, n_obs.

Examples

# Toy Gaussian family for illustration
fam <- tulpa_family(
  name = "gaussian",
  simulate_fn = function(eta, params, n_obs, ...) {
    rnorm(n_obs, eta[[1]], params$sigma_y)
  },
  extra_params = list(sigma_y = prior_half_normal(1))
)
df <- data.frame(y = rep(0, 20), x = rnorm(20))
pp <- prior_predict(y ~ x, fam, df, n_draws = 50, seed = 1)
length(pp$y)  # 50

Show default priors for a tulpa family

Description

Display the default prior specifications used for each model family. Useful for understanding what priors are applied before fitting and as a starting point for customization.

Usage

priors_default(family = NULL, spatial = FALSE, temporal = FALSE)

Arguments

Details

Default priors in tulpa follow these principles:

• Fixed effects (beta): Normal(0, 2.5) - weakly informative, allows coefficients roughly in [-5, 5] on the link scale.

• Random effect SD (sigma): PC prior with P(sigma > 1) = 0.01 - favors simpler models with smaller variance components.

• Overdispersion (phi): PC prior with P(phi > 10) = 0.01 on the NB2 size phi (larger phi is less overdispersion; phi -> Inf is the Poisson limit). The default keeps phi finite, allowing overdispersion while penalising extreme values.

• Temporal correlation (rho): Beta(2, 2) - symmetric prior centered at 0.5, appropriate for AR(1) correlation.

• Spatial mixing (rho_spatial): Beta(1, 1) = Uniform(0, 1) - no prior preference for structured vs. unstructured spatial variation.

Value

Invisibly returns a tulpa_priors object with the defaults. Primarily called for its side effect of printing.

See Also

tulpa_priors() for creating custom priors

Examples

# Defaults for all families (no family argument)
priors_default()

# Family-specific defaults take a tulpa_family object. Model packages
# (e.g. tulpaRatio) register rich families; a minimal one is enough here.
fam <- tulpa_family(
  name = "poisson_gamma",
  simulate_fn = function(eta, params, n_obs, ...) rpois(n_obs, exp(eta[[1]]))
)
priors_default(fam)

# Including spatial parameters
priors_default(fam, spatial = TRUE)

# Use as a starting point for customization
my_priors <- priors_default(fam)
my_priors$beta <- prior_normal(0, 1)  # Tighter prior on fixed effects

Random-effect summaries

Description

Random-effect summaries

Usage

ranef(object, ...)

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

Arguments

Details

What each backend reports for a group effect follows what it computes:

• Sampler tier and the RE-covariance Gibbs debias (tulpa_re_cov_gibbs(), reached by tulpa(..., control = list(re_cov = "gibbs"))) draw the random effects jointly with everything else, so estimate / sd / bounds are the empirical posterior summaries.

• The RE-covariance integrator (tulpa_re_cov_nested()) carries a Gaussian per-group posterior at each Sigma node; the reported summaries are the exact moments and quantiles of the weighted mixture of those, so they carry both the within-node curvature and the Sigma uncertainty. A group effect the subspace debias selected (control$subspace_debias) is moved by the Metropolis sampler at every node instead, and is reported from those draws.

• The Laplace tier reports the conditional mode with no spread (sd and the bounds are NA), which is the only per-group quantity it forms.

The source column says per row which of these produced it: "sampled" for a posterior draw summary, "mixture" for the node mixture, "mode" for a conditional mode. A fit whose backend never forms a per-group posterior at all (the adaptive Gauss-Hermite inner marginal integrates each group out by quadrature) errors with that reason rather than returning an empty table, which would be indistinguishable from a model with no random effects. A model that genuinely has none returns a zero-row data frame.

Value

Data frame with one row per random-effect coefficient: term (the group level, and the coefficient for a random slope), estimate, sd, the 2.5% / 97.5% bounds conf.low / conf.high, and source (which construction the row came from). sd and the bounds are NA on a backend that reports a point per group (see Details).

Examples

set.seed(1)
df <- data.frame(x = rnorm(100), g = factor(rep(1:10, 10)))
df$y <- rpois(100, exp(0.3 * df$x))
fit <- tulpa(y ~ x + (1 | g), data = df, family = "poisson")
ranef(fit)

Rational Approximation Coefficients for a Fractional SPDE

Description

Returns the coefficient descriptor for an SPDE Matern field of operator order alpha = nu + 1 (in 2D). For integer nu the construction is exact and no rational approximation is needed. For fractional nu it returns the rational-SPDE roots from the BRASIL best-rational approximation (see Details).

Usage

rational_spde_coefficients(nu, m = 4L, lambda_range = c(1e-04, 10000))

Arguments

Details

For integer nu (1, 2, 3, ...) the operator order alpha = nu + 1 is an integer and the precision is assembled directly from integer powers of the FEM operator ⁠L = kappa^2 C + G⁠ – an exact construction.

For fractional nu the field uses the operator-based rational SPDE approximation (Bolin & Kirchner 2020). With beta = alpha / 2 and m_beta = max(1, floor(beta)), the field-mode variance of the assembled precision ⁠Q = Pl' C^{-1} Pl⁠ (field ⁠u = Pr x⁠, x ~ N(0, Q^{-1})) tracks the Matern spectral density ⁠l^{-2 beta}⁠ when ⁠prod(1 - l rc) / (l^{m_beta - 1} prod(1 - l rb))⁠ approximates l^{-beta_rem} on the scaled spectrum, beta_rem = beta - (m_beta - 1). The roots come from the degree-⁠(m, m)⁠ best uniform (minimax) rational approximation of x^{-beta_rem}, computed by the BRASIL algorithm (Hofreither 2021); its numerator zeros map to rc = 1 / zero and denominator poles to rb = 1 / pole. The field assembly from these roots is .spde_rational_assemble(); it is validated against the Matern spectral density in test-spde-rational.R.