tulpa_laplace: Fit a model via Laplace approximation

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

Fit a model via Laplace approximation

Description

General-purpose Laplace approximation for latent Gaussian models. Finds the mode of the latent field (beta + random effects) and returns the Laplace-approximated marginal likelihood.

This is the public API for model packages (tulpaGlmm, tulpaObs, etc.) to call tulpa's Laplace engine. As the low-level engine entry point (not a front-door fitter), it keeps its numerical controls (max_iter, tol, n_threads, return_hessian) inline in the signature rather than in a control list, so callers assembling many solves pass them positionally.

Usage

tulpa_laplace(
  y,
  n_trials,
  X,
  re_list = list(),
  family = "binomial",
  phi = 1,
  phi2 = NULL,
  spatial = NULL,
  weights = NULL,
  offset = NULL,
  max_iter = 100L,
  tol = 1e-06,
  n_threads = 1L,
  return_hessian = TRUE,
  beta_prior = NULL,
  return_re_cov = FALSE,
  X_zi = NULL,
  zi_prior_sd = 2.5,
  return_joint_hessian = FALSE,
  compute_skew = FALSE,
  skew_idx = NULL,
  debias = NULL
)

Arguments

y	Response vector (integer for binomial/poisson/negbin, numeric for gaussian)
n_trials	Trial sizes (integer vector, used for binomial only)
X	Fixed-effects design matrix
re_list	List of RE specifications. Each element is a list with: idx: integer vector of group indices (1-based) n_groups: number of groups n_coefs: coefficients per group (1 = intercept-only, >1 = random slopes) sigma: per-coefficient RE standard deviation(s), a diagonal covariance (uncorrelated, lme4 (x || g)). A scalar is recycled to n_coefs. Z: slope design matrix (n_obs x n_coefs) when n_coefs > 1; NULL means intercept-only. L / cov: optional ⁠n_coefs x n_coefs⁠ covariance for a correlated term (lme4 (1 + x | g)) – supply either a lower-triangular Cholesky factor L (covariance = ⁠L L'⁠) or the covariance matrix cov. When present these take precedence over sigma; the off-diagonal enters both the joint Hessian (mode finding) and the marginal fixed-effect SE.
family	Character: "binomial", "poisson", "neg_binomial_2", "gaussian"
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.
phi2	Optional second dispersion: the Student-t degrees of freedom (family = "t"; default 4 when NULL). Non-spatial path only.
spatial	Optional spatial specification (tulpa_spatial object)
weights	Optional observation weights (numeric vector, length length(y)). Scales each observation's log-density, score and Fisher curvature by the same w_i, on the spatial route as well as the non-spatial one, so the mode and the marginal precision H_beta describe one weighted model. NULL (default) uses 1.
offset	Optional observation-level offset on the linear predictor (numeric vector, length length(y)). NULL (default) uses 0.
max_iter	Maximum Newton iterations (default 100)
tol	Convergence tolerance (default 1e-6)
n_threads	Number of threads (default 1)
return_hessian	Logical: return the fixed-effect Hessian block? (default TRUE)
beta_prior	Optional Gaussian prior on the fixed effects. NULL (default) keeps the weak built-in prior beta ~ N(0, 100^2). Otherwise a list with element sd (prior standard deviation, required) and optional mean (prior mean, default 0). Each may be a scalar (applied to every coefficient) or a length-ncol(X) vector. Adds sum((beta - mean)^2 / (2 * sd^2)) to the negative log-posterior, so the mode is the penalized (MAP) estimate. A coefficient's sd may be +Inf, which sets its precision to 0 (no penalty on that coefficient). Not supported on the spatial path.
return_re_cov	If TRUE, additionally return per-group marginal posterior covariance blocks Cov(u_g | y, Sigma) – one ⁠n_coefs x n_coefs⁠ matrix per (RE term, group), with the fixed effects and other groups marginalized out (each block is a diagonal block of the full inverse Hessian, not the inverse of a diagonal block). Used by the EM M-step for a full random-effect covariance. Non-spatial multi-RE path only.
X_zi	Optional zero-inflation design matrix (length(y) rows). When supplied the latent fixed-effect block becomes ⁠[beta_count | beta_zi]⁠ and the family's compiled zero-inflated kernel is used. Non-spatial path only.
zi_prior_sd	Prior SD on the zero-inflation coefficients, beta_zi ~ N(0, zi_prior_sd^2) (default 2.5, matching the samplers' ModelData::zi_prior_sd). It is what keeps the logit identified when a level contributes no zeros, where the likelihood alone drives beta_zi to -Inf. +Inf removes the penalty. Ignored when X_zi is NULL.
return_joint_hessian	If TRUE, additionally return H_joint: the full joint posterior precision of the latent field ⁠[beta | random effects]⁠ at the mode, as a symmetric sparse matrix. This is the matrix the Laplace approximation takes the determinant of, so it is what an exact derivative of the log-marginal has to differentiate through; H_beta is only its fixed-effect Schur complement. Costs one extra copy of the Hessian, so it is off by default. Non-spatial multi-RE path only. With return_hessian = TRUE it also returns H_latent, the same precision at the observed curvature (identical to H_joint except for the families whose Newton weight is the expected one); H_beta is its Schur complement, and posterior_predict() draws the random effects from it.
compute_skew	If TRUE, additionally return the inner-Laplace reliability material at skew_idx: inner_skew (gamma_3, Rue Martino & Chopin 2009's cubic term), and the importance curve inner_is_z / inner_is_log_joint the inner Pareto-k-hat is fitted from. Costs one linear solve plus a fixed batch of objective evaluations per probed index. Non-spatial path only.
skew_idx	1-based latent indices to probe (in the ⁠[beta | random effects]⁠ layout of mode). NULL with compute_skew = TRUE probes every latent index.
debias	Subspace debias: a list with idx (1-based latent indices to correct by Metropolis along the Gaussian-conditional-mean surface through the mode) and optional n_iter / warmup / thin. The result then carries debias_draws (⁠n_kept x length(idx)⁠, the sampled x_S - mode_S), debias_sigma_ss (the inner Laplace's own marginal covariance of x_S), debias_accept and debias_idx. NULL (default) or an empty idx leaves the solve bit-for-bit as it was and consumes no random number. Non-spatial path only.

Value

A list with:

• mode: full mode vector (beta, then RE values per term)

• log_marginal: Laplace-approximated log-marginal likelihood

• n_iter: number of Newton iterations

• converged: logical, whether the stopping rule was met

• score_max: largest absolute component of the joint penalized score at the returned mode – the residual the solve actually achieved, which is a different question from converged and the one anything differentiating through the mode depends on

• log_det_Q: log-determinant of the Hessian

• H_beta: fixed-effect block of the Hessian (if return_hessian = TRUE)

• H_joint, H_latent: the joint latent precision, at the kernel's and at the observed curvature (if return_joint_hessian = TRUE)

• cov_blocks: list of per-group posterior covariance matrices, one per (RE term, group) in term-major then group order (if return_re_cov = TRUE)

Examples

set.seed(1)
n <- 200L
X <- cbind(1, rnorm(n))
eta <- X %*% c(-0.3, 0.8)
y <- rbinom(n, 1, plogis(eta))
fit <- tulpa_laplace(y, rep(1L, n), X, family = "binomial")
fit$mode          # posterior mode of the fixed effects

tulpa documentation built on Oct. 19, 2026, 5:09 p.m.