Nested Laplace approximation for latent Gaussian models — print.tulpa_nested

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

Supported priors:

Spatial (areal): "icar" (1D grid on tau), "bym2" (2D on (sigma, rho)), "car_proper" (2D on (tau, rho); rho lives in the eigenvalue interval (1/lambda_min, 1/lambda_max) of D^{-1} W).
Spatial (continuous): "nngp" (2D on (sigma2, phi_gp)), "hsgp" (2D on (sigma2, lengthscale)).
Temporal: "rw1", "rw2" (1D grid on tau), "ar1" (2D on (tau, rho))
SPDE continuous spatial: see cpp_nested_laplace_spde() (separate entry, rebuilds Q via SPDE Q-builder).

Usage

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

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

Arguments

x

A nested-Laplace fit.

...

Ignored.

y

Response vector.

n_trials

Trial sizes (binomial). Pass 1L-vector otherwise.

X

Fixed-effects design matrix.

prior

A list describing the latent prior block. Required field type one of {"icar", "bym2", "car_proper", "rw1", "rw2", "ar1"}. Type-specific fields:

icar: spatial_idx, n_spatial_units, adj_row_ptr, adj_col_idx, n_neighbors (CSR adjacency, 0-based); optional tau_grid.

Any block accepts an optional svc_weight, one weight per observation, which makes it a varying coefficient: observation i contributes svc_weight[i] * f_i instead of f_i, where f_i is the block's field at that row (z[spatial_idx[i]] for an areal block, (A u)_i for an SPDE one). Absent, the field enters unweighted.

bym2: same adjacency; scale_factor; optional sigma_grid, rho_grid.
car_proper: same adjacency; optional tau_grid, rho_grid, rho_bounds = c(lower, upper) (defaults to (0, 1)).
rw1/rw2: temporal_idx (1-based), n_times; optional tau_grid, cyclic (default FALSE).
ar1: temporal_idx, n_times; optional tau_grid, rho_grid (each rho_grid value strictly inside (-1, 1)).

An icar, bym2, car_proper or nngp block may carry projector, an N x n_spatial_units (N x n_spatial for nngp) matrix (dense or a Matrix sparse matrix), in place of spatial_idx: observation i then receives (A z)_i, so restricted spatial regression is A = S P_perp with S the observation-to-unit incidence. The block is fitted by the multi-block driver, alone or in a list of blocks. An intrinsic field's level is placed from A 1: removed with no fold where it is zero, folded into the intercept where it is the intercept column, otherwise left to the augmentation. The nngp and car_proper priors are proper and carry their level themselves.

Grids of a two-axis type (bym2, car_proper, ar1) are PAIRED when both are supplied: the i-th entries of tau_grid / sigma_grid and rho_grid form the i-th outer-grid cell, so the two must have equal length, and a tensor grid is written as expand.grid() columns. A single supplied axis is kept as an axis (its distinct values) and crossed with the default of the other.

A default grid axis is a starting axis, not a hard ceiling: for icar (tau_grid) and bym2 (sigma_grid) a posterior mode that rails a boundary node (pareto_k_regime = "collapsed_edge") triggers one mode-Hessian recenter-and-refit, reported through outer_grid_placement / outer_grid_recenter_declined – see tulpa_nested_laplace_joint()'s return docs. An axis the caller pinned is never moved; mark a grid your own code defaulted with auto_grid() to keep the recenter live on it.

spec

Optional tulpa_temporal or tulpa_spatial spec object (output of temporal_rw1(), temporal_rw2(), temporal_ar1(), spatial_car(), spatial_bym2(), etc.). When supplied alongside data, the prior list is built automatically via prior_from_spec() – pass either prior or spec, not both.

data

Data frame used to validate spec and resolve time/group/site indices. Required when spec is supplied.

re_idx

Optional 1-based RE group index per obs (defaults to no RE).

n_re_groups

RE group count (default 0).

sigma_re

RE standard deviation (default 1).

family

"binomial", "poisson", "neg_binomial_2", etc.

phi

Dispersion passed to the family, held fixed. One convention at every door: for gaussian / lognormal this is the residual VARIANCE (the SD is sqrt(phi)), for neg_binomial_2 the size, gamma the shape, beta the precision, t the scale; binomial and poisson ignore it. The compiled kernels parameterize the two variance families by the residual SD and are handed sqrt(phi) at the boundary.

Defaulted, it conditions at 1 and says so with a warning, since for a family that reads a dispersion that is a modelling choice rather than a neutral value. fit$phi_estimated records whether the value on the fit was estimated or conditioned on.

offset

Optional per-observation offset added to the linear predictor (length length(y)), as an offset() term in a model formula. Every single- and multi-block kernel carries it, and fitted_eta includes it.

likelihood

Optional model-supplied likelihood, replacing the built-in family. Pass an external pointer to a tulpa::NestedLikelihood (built in a model package's own C++ from a LikelihoodSpec); the inner Laplace solve then reads the per-observation score, Fisher weight, and log-likelihood from that spec instead of family, so family/phi are ignored. Used by model packages to fit a custom response without adding a family to tulpa – for example tulpaObs threads its marginalized single-season occupancy likelihood (a scaled Bernoulli, with the latent occupancy state integrated out) through this. Multi-block prior only. Default NULL (use family).

hyperprior

The prior an outer hyperparameter axis carries when the call states no density for it, "proper" (default) or "flat". "proper" folds a normalised density on the axis's integration coordinate into log_marginal: the PC prior P(sigma > 3) = 0.01 on a standard deviation, variance or precision axis, the PC range prior at 0.2 times the coordinates' bounding-box diagonal on a range or lengthscale axis, a uniform on a bounded axis, and the PC + LKJ prior over a free-covariance block. An axis with no sourced density is named in log_hyperprior_declined. "flat" folds no density of the engine's own, so such an axis is integrated under its cell measure alone, is named in log_hyperprior_declined as "flat_hyperprior", and the fit's log_evidence declines with "improper_hyperprior". A density the call states – a prior_* argument, a block's rho_prior, prior_range or prior_sigma, the copy coefficient's slab, a tgmrf() block's own prior – applies under either choice. The same two choices, under the same names, are offered by tulpa(), tulpa_eb(), tulpa_re_cov_nested() and fit_spde().

control

Optional list of perf/numerical tuning knobs (statistical arguments stay top-level), following the control convention of tulpa(). Recognised elements (defaults in parentheses):

max_iter (50L), tol (1e-6) – inner Newton iteration budget and tolerance.
n_threads (1L) – inner-loop OpenMP threads.
x_init (NULL) – warm-start for the first grid point's inner solve.
keep_grid_hessians (TRUE) – when TRUE, retain per-grid-point fixed-effects marginal Hessian $H_\beta$ and mode $\hat{\beta}$ on the return list as $grid_hessians (list of dense $p\times p$ matrices) and $grid_modes (list of length-$p$ vectors). Used downstream by simplified-Laplace (SLA) callers to assemble skew-aware marginals – see the cumulant pooling in rubins_pool().
diagnose_k (TRUE), k_samples (500L) – compute the outer Pareto-$\hat{k}$ accuracy diagnostic ($pareto_k) by importance sampling the hyperparameter posterior against the Gaussian proposal fitted to the grid, drawing k_samples extra inner-marginal evaluations. k_samples is a precision knob: the GPD tail size is held at the fraction the default budget implies, so a larger budget sharpens the same k-hat instead of moving it to a deeper quantile of the weight distribution (gcol33/tulpa#631). k_tail_points overrides that tail size directly (capped at 20% of the draws). Computed for a single-block, single positive-scale-axis grid; left NA (with the grid's quadrature ESS as the fallback diagnostic) for multi-block, multi-axis, or bounded-parameter grids. See tulpa_psis().
diagnose_skew (TRUE), skew_idx (NULL) – compute the inner-Laplace skewness diagnostic ($inner_skew, gamma_3, Rue Martino & Chopin 2009 Sec 3.2.3) at the fitted MAP grid cell: one extra Newton solve, scoring the p fixed-effects latent indices by default (pass skew_idx, 1-based latent indices, to score additional ones, e.g. specific spatial units – the full latent field is not scored by default since it costs one linear solve per index). This is the complementary layer to diagnose_k: the outer diagnostic scores the hyperparameter-grid integration around a FIXED inner Laplace, this scores whether that inner Gaussian approximation is itself a good fit to the latent-field conditional posterior. See diagnostics() for the combined whole-fit verdict.
within_cell ("log_quadratic") – the WITHIN-CELL construction the reported per-axis hyperparameter intervals are read with. The outer grid's weights say how much mass each cell holds; they do not say how it is spread inside the cell, and a quantile needs both. "log_quadratic" joins each row's node log densities by overlapping quadratics, exact for a Gaussian at any spacing, and continues the end quadratic past the outer node instead of stopping at the box edge; it runs on the rows that resolve their own conditional (h / sd <= 2) and reads every other row as boxes, declining to "box_uniform" where no row qualifies. "box_uniform" puts the cumulative FULL mass at each cell EDGE and interpolates between edges; "chord" puts the cumulative MID-mass at each cell coordinate and interpolates between coordinates. The default moved from "box_uniform" at gcol33/tulpa#932: a placed axis sits at 1.25 posterior SDs per cell, where the box read is narrow by construction, and against dense references of a two-arm ICAR fixture the log-quadratic read's summed error is 0.038 against the box read's 0.110 on the same grids. outer_grid_h_over_sd is how wide a cell is on each axis, and theta_within_cell is what each axis was actually read with. Only a "density" support admits the box and log-quadratic reads – a CCD design, a locally refined grid and a posterior sample are not cell partitions that tile – and an axis they decline on reports the read that ran with a reason rather than erroring. Nothing else moves: point estimates, moments, draws and weights are untouched.
skew_correct (TRUE) – consume the inner-Laplace expansion instead of only grading it: report Cornish-Fisher marginal quantiles at each coefficient's own gamma_3, about the centre gamma_1 + gamma_3 / 2 that Rue, Martino & Chopin (2009) eq. (22) implies, from summary() / confint() wherever the combined inner band says the leading-order expansion is in its regime, and the Gaussian quantiles everywhere else. It is post-processing on the reported quantiles: draws, modes and weights are untouched, so a fit run with it off is bit for bit the fit it was before. A coefficient whose location term could not be formed declines rather than reading it as zero. The band that bounded the relocation itself (centre_unreliable) is off (Inf): scored over seven fixtures with an exact reference, every finite cutoff declines the coefficients the correction helps most, because a large centre carrying a small gamma_3 is uniformly weak correlation rather than an expansion out of its regime. $skew_correction records the per-coefficient gamma_3, gamma_1 and the centre they form, the band, the inner importance k-hat, the combined band, the eligibility and the reason behind it; the skew_applied attribute on summary() / confint() records what was actually used at the requested level. RMC fit a skew normal here instead; the series correction is the same-order alternative, and unlike a skew normal its skewness does not saturate inside the band it is applied on.

MEASURED. Against exact quadrature quantiles of rare-event binomial-logit posteriors it cuts total absolute endpoint error 69.2%, improving both endpoints in every case. Scored over the WHOLE marginal – paired CRPS against the exact posterior in a 400-replicate prior-predictive experiment – it reads -0.01643 against the uncorrected Laplace at t = -1.89, essentially all of the -0.01662 the exact posterior itself achieves, and its PIT re-enters the simultaneous SBC band. Applied about the Laplace mode instead of about gamma_1 + gamma_3 / 2 the same reshaping scored +0.00775 at t = +3.54, a net loss; the location term is what supplies that centre.

IT IS ON BY DEFAULT, so summary() / confint() on a nested-Laplace fit report the corrected quantiles wherever the combined inner band admits the coefficient; skew_correct = FALSE restores the uncorrected report exactly. Scored against the mixture read a correction-off fit gives, the flip is t = -1.895 on the rare-event intercept and -3.765 / -3.201 on the small-group Bernoulli design. Across twelve model classes read off one solve per seed, pooled 95% coverage moves 0.9510 -> 0.9542 at a standard error of 0.0070, with every class inside the acceptance the shipped gates use. A fit the correction cannot help – a coupled one, whose location term is unreachable – reports what it reported before, to the bit.
subspace_debias (FALSE) – correct only the latent directions the inner-layer diagnostics flagged, by exact Metropolis, and leave the rest at their Gaussian conditional. TRUE takes every default; a list overrides band (the inner-reliability floor a coordinate is selected at, default "ok"), idx (pin the corrected set explicitly, skipping the selector), closure / closure_max (grow the set by strongly coupled precision-graph neighbours: a coordinate strongly coupled to a member of the corrected set and left out of it is carried linearly, which is the error the correction removes, so requesting the closure also retains the modal cell's joint precision for the selector to read), the sampler budget n_iter / warmup / thin, and n_draws. The selector reads the per-index bands diagnose_skew already attached, so it costs no extra solve; the correction itself re-runs the settled grid once with the sampler on, and the fit then reports $draws – the per-cell Metropolis sample for the selected coordinates, the rest from the Gaussian conditional given them – instead of the Gaussian-mixture moments. An EMPTY selection leaves the fit bit-for-bit identical to the plain path. $subspace_debias records what was selected, the bands it was read from, and the per-cell acceptance rate. Requesting it turns keep_grid_hessians on, since the recombination reads exactly those per-cell pieces.
cila (FALSE) – corrected integrated Laplace, the second inner-layer debias (after Lai, Margossian and Sheldon, arXiv:2605.20345). Where subspace_debias selects coordinates and runs exact Metropolis on them, this selects nothing: at every outer cell it draws n_points points from the whole inner Gaussian, weights each by the exact joint density it came from, and reports the weighted particles. TRUE takes the defaults; a list overrides n_points (1024L), variant ("qmc", a Sobol net; "is" for iid draws, "rqmc" for the net under n_shift random shifts), n_shift (8L), n_draws and seed. Below 512 points a cell's particle set is too coarse to be a marginal at all and the request is refused. The corrected per-cell masses become the fit's own weights / log_marginal and weights_source reports "cila"; the pre-correction pair is kept as $cila$laplace. A cell whose inner solve factorized sparsely draws through the CHOLMOD factor's own triangular and permutation solves; an LDL' factor has no square root to draw with and is declined with "sparse_factor_not_ll".
auto_recenter (TRUE) – outer-grid placement policy. TRUE re-centres the movable default axes on the posterior mode and refits when the grid either RAILS (an axis's own marginal is maximal at one of its own endpoints) or does not RESOLVE its posterior (an axis's node spacing exceeds 2 posterior SDs in its own coordinate). Both tests read the weights the fit already stored, so a grid that brackets and resolves its mode costs nothing beyond them; when the pass does fire it is a second full grid solve plus a finite-difference mode/Hessian stencil. A recentred axis is mode +/- 2.5 sd over 5 nodes, a cell width of 1.25 posterior SDs by construction, against a census median of 3.9 on the fixed spans.

Measured over 200 fixed-truth seeds on each of six configurations (icar chain / icar lattice / rw1 / bym2 / iid / nngp), the default moves mean |coverage - nominal| from 0.043 to 0.030 at the 95% level, 0.171 to 0.084 at 80% and 0.243 to 0.129 at 50% against the rail-only policy, at 0.63 times the 95% interval width and 0.76 times the median bias, for 1.71 times the wall clock.

Three other values. "rail" is the rail test alone, which is what TRUE meant before the sizing measurement settled the default. FALSE integrates over the grid exactly as given, whatever it is, and records outer_grid_recenter_declined = "auto_recenter_disabled". "always" re-centres every movable default axis whatever the fit did, at 2.04 times the wall clock; it agrees with the default seed for seed on five of the six measured configurations and differs on the one whose default axes already resolve their posterior, where it takes 50% coverage to 0.135 against a nominal 0.5.

Which families carry movable axes is .NL_REGISTRY_AXIS_FIELD: icar, rw1, rw2, iid, bym2, nngp, hsgp and spde. car_proper, ar1 and hsgp_mo each carry a correlation axis with no guessable coordinate, and mcar / miid / tgmrf hold their axes in one matrix field; a fit of any of them records which through outer_grid_recenter_declined rather than passing in silence. The per-axis policy names are the standalone registry path only – tulpa_nested_laplace_joint() and fit_st_nested() refuse them rather than accept them and ignore them.
max_grid_cells (NULL, no ceiling) – a ceiling on a multi-block outer grid's cell count, one inner Newton solve per cell; a grid past it is refused before any solve, naming the axes that produced the count. Unset, a grid of any size is solved, and one past 2048 cells is announced with a warning first.
checkpoint (list(path =, resume =)) – grid-cell checkpoint / resume. Each solved outer cell is appended to path, keyed by its hyperparameter coordinate; resume = TRUE loads the finished cells and solves only the rest, and a fingerprint mismatch (different data, settings or grid) errors rather than resuming onto a stale result. See ?tulpa_nested_laplace "Checkpoint / resume".
progress, progress.every, progress.file, progress.throttle – the outer-grid progress reporter: whether to print, how often (in cells), where to, and the minimum seconds between lines.
prune (TRUE), prune_tol (1e-3), screen_iters (.NL_SCREEN$iters, 2 – the depth gcol33/tulpa#640 measured the screen down to) – cheap-pass screening of the outer grid. The driver first sweeps the lattice running a screen_iters-step inner Newton per cell, each warm-started from the previous screened cell's quasi-mode, computes a screening Laplace log-marginal, softmax-normalises it, and skips the full inner Newton on every cell whose screened weight is below prune_tol. The neighbour warm-start keeps each cheap mode near its cell's own mode, so the cheap ranking tracks the full-solve ranking even where the latent mode moves across the grid. Pruned cells get log_marginal = -Inf, n_iter = 0 and inherit the pilot mode; the pilot cell is never pruned, and at least prune_min_keep cells (5) are solved in full whatever the tolerance says – the highest-ranked dropped cells are restored up to that floor, because the outer grid placement pass reads a finite-difference curvature off the cells that were SOLVED and a kept set of one carries none. After the full pass the posterior mass the dropped cells may carry is bounded by lifting each by the worst screening error measured on a kept cell (prune_dropped_mass_bound); past 1% the driver solves the dropped cells that contribute most to it until it is under 1% (prune_n_repaired), and a safety gate falls back to the full grid (with a warning) if the bound still fails. The bound assumes a dropped cell is screened no worse than the worst kept one. On by default (gcol33/tulpa#925); FALSE solves every cell. prune_tol must be a finite numeric in [0, 1) and screen_iters a single integer >= 1; both are validated whatever prune says. prune / prune_tol reach both the single-block and the multi-block dispatch; the screening DEPTH is declared by the single-block kernels only, so a multi-block prior refuses a pinned screen_iters rather than accepting it and ignoring it.
prune_log_gap (unset) – the screening cut stated in nats instead of as a normalised weight: keep every cell within this many nats of the best screened cell. prune_tol is a weight, so what it cuts at is the gap -(log(prune_tol) + log(Z)), and the whole range a caller would type (1e-3 to 1e-12) is 7 to 28 nats of it – on an outer surface whose log-marginal spans thousands of nats every setting returns the same kept set. Stating the gap directly keeps the knob's resolution wherever the surface is steep. It reaches the kernel as prune_tol = exp(-gap), so the realised cut is the requested gap less log(Z), at most log(n_grid) nats narrower; the fit reports the realised value. Setting both this and prune_tol is an error – they state one cut in two units.
fitted_var (TRUE) – fill the per-row predictive variance fitted_eta_var, the companion of fitted_eta a caller marginalises the per-row linear predictor with. It costs one back-solve per distinct loading vector per cell, which on a design with few repeated rows is the dominant cost of a cell, so a caller reading only fitted_eta or the coefficient summaries sets FALSE and the fit then carries no fitted_eta_var. Every single-block field reports it, the NNGP, HSGP, SPDE and spatiotemporal fields included, each read off its cell's own precision at the mode; a multi-block prior refuses FALSE rather than accepting it and ignoring it.
var_of_means_consistency (TRUE) – the consistency pass of tulpa_nested_laplace_joint(): a field SD (tau of an icar, rw1 or rw2 block, sigma of an iid, bym2 or spde block) whose marginal has collapsed onto too few nodes to carry a spread has levels added in every row of the other axes that holds the posterior until it resolves, and the new cells are solved and carried like the rest (var_of_means_consistency_info records the pass). Other axes keep the nodes they were laid on.
adaptive_grid (FALSE), adaptive_grid_edge_thresh (0.02), adaptive_grid_max_passes (1L) – the opt-in boundary / interior pass, on the same axes.
axis_refine (NULL) – per-axis override of how far those axes may be moved, named by the grid's axes ("tau"; "b1.tau" on a multi-block prior): "none", "densify" (inside the declared span) or "extend". Nodes the caller wrote down densify by default and an engine-placed axis extends, as in tulpa_nested_laplace_joint().

Value

A list with:

theta_grid: matrix or vector of grid hyperparameter values.
log_marginal: log p(y, mode | theta_k) at each grid point.
weights: integration weights normalising to sum 1.
theta_mean, theta_sd: posterior moments per hyperparameter.
n_iter: inner Newton iterations per grid point.
modes: matrix [n_grid x n_x] of inner modes, when stored.
pareto_k, pareto_k_is_ess: outer Pareto-$\hat{k}$ and its importance-sampling ESS (NA when not computed for the grid; see control$diagnose_k).
pareto_k_proposal_source, pareto_k_first_pass: which proposal family the reported k-hat came from (the outer k-hat is scored against several candidates and the best is kept), and the k-hat of the FIRST pass – the proposal exactly as this backend placed it, before any candidate refined or replaced it. A large gap between the two says the nodes are badly scaled around the hyperparameter posterior even though the verdict is fine; no gap says the placement was already right.
inner_skew, inner_skew_idx, inner_skew_dropped: the inner-Laplace skewness diagnostic (gamma_3) at each scored latent index and its 1-based index, plus a count of (index, observation) contributions dropped for a non-finite third derivative (see control$diagnose_skew).
timing: named numeric of wall-clock seconds (total, setup, grid, postproc, diagnostics); the grid phase is the inner Laplace pass that scales with grid size. Surfaced one-line in print.
prune_cheap_log_marginal, prune_mask, prune_n_pruned, prune_tol: present when the grid was screened and the safety gate did not trip – the cheap-screen log-marginal per cell (the truncated solve's value plus half its Newton decrement, the second-order estimate of the converged value the cells are ranked by; the decrement itself is prune_screen_decrement), the logical mask of pruned cells, the pruned-cell count, and the threshold actually applied.
prune_log_gap_cut, prune_cheap_lm_spread, prune_min_keep, prune_n_floor_restored: the same screen read on the scale it operates on – how many nats below the best screened cell the tolerance cut at, how many nats the screened surface spans, the kept-cell floor applied, and how many cells that floor put back. A cut that is a sliver of the spread is a tolerance with no resolution on this grid.
prune_dropped_mass_bound, prune_n_repaired, prune_gate_mass: the posterior mass the dropped cells may carry, with each lifted by the worst screening error measured on a kept cell; how many dropped cells the driver solved in full to bring that under the gate; and the gate.
prune_fallback_triggered, prune_fallback_reason: present only when the gate did trip; the rest of the fit is then the full-grid result and the reason string records which check fired.
prior: echoed input.

Details

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

References

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

Examples

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