Skip to contents

Nested-Laplace integration over random-effect covariances

Source: R/nested_laplace_re_cov.R
tulpa_re_cov_nested.Rd

For one or more random-effects terms (e.g. (1 + x | g), (1 + x || g), or several terms together), integrate the Laplace marginal likelihood over the random-effect covariances Sigma instead of fixing them at point estimates. Reports weighted posterior summaries (mean, SD, median, 2.5\ Sigma and its derived scale (sigma_i) and correlation (rho_ij) parameters, marginalizing the joint posterior over a Sigma-grid.

This corrects the plug-in-MAP ("summary") bias: the mode of a skewed variance-component marginal is biased low relative to its median, so the headline summary should be the marginalized median, not the mode.

Usage

tulpa_re_cov_nested(
  y,
  n_trials = NULL,
  X,
  re_terms,
  family = "binomial",
  phi = 1,
  phi2 = NULL,
  prior_sigma = NULL,
  eta = NULL,
  hyperprior = c("proper", "flat"),
  log_prior_theta = NULL,
  beta_prior = NULL,
  offset = NULL,
  weights = NULL,
  n_quad = 1L,
  X_zi = NULL,
  zi_prior_sd = 2.5,
  control = list()
)

Arguments

y, n_trials, X, family, phi

Passed to tulpa_laplace() for the inner solve. n_trials = NULL defaults to 1 (binary / single-trial).

re_terms

Either a single random-effect term or a list of them. Each term is a list with idx (1-based group index per observation), n_groups, n_coefs (c), Z (the n_obs x c RE design, e.g. cbind(1, x) for (1 + x | g); only required when c > 1), and correlated (TRUE for a full Sigma, FALSE for a diagonal one; defaults to TRUE). An optional label / group_var names the block in the output. Any L / cov / sigma field is ignored – Sigma is what this function integrates over.

phi2

Optional second dispersion, threaded into every inner tulpa_laplace() solve: the Student-t degrees of freedom (family = "t", default 4 when NULL) or the Tweedie variance power (family = "tweedie", required – a defaulted power would be a statistical decision the caller never made). A phi2 supplied for any other family errors rather than being ignored. It is conditioned on: the integration is over the random-effect covariances, not over phi2.

prior_sigma, eta

Hyperparameters of the PC + LKJ prior used when hyperprior = "proper" (see re_cov_pc_lkj_prior()): prior_sigma = c(U, alpha) with P(sigma_i > U) = alpha (NULL, the default, is c(3, 0.01)) and LKJ shape eta (NULL is 2). Ignored when hyperprior = "flat" or log_prior_theta is supplied.

hyperprior

"proper" (default) or "flat". "proper" builds the PC + LKJ prior from prior_sigma / eta, the proper prior every other scale axis of the engine carries by default. "flat" integrates with log_prior_theta the zero function (flat in log(theta)); that prior is improper, so the fit then reports no evidence. Ignored when log_prior_theta is supplied.

log_prior_theta

Optional function(theta) returning a scalar log prior density on the full stacked parameter vector, overriding hyperprior entirely. Default NULL, which defers to hyperprior.

beta_prior

Optional Gaussian prior on the fixed effects, threaded into every inner tulpa_laplace() solve (list(mean, sd)). NULL (default) keeps the weak built-in prior.

offset

Optional observation-level offset on the linear predictor (length length(y)), e.g. log(exposure) for a rate model. Not supported with n_quad > 1, which errors rather than dropping it.

weights

Optional non-negative observation weights (length length(y)), each scaling its row's log-likelihood in every inner tulpa_laplace() solve, so the covariance is integrated under the weighted likelihood. Not supported with n_quad > 1, which errors rather than dropping them.

n_quad

Quadrature order for the inner marginal. 1 (default) uses the joint-field Laplace inner solve (tulpa_laplace()). > 1 refines the inner marginal with n_quad-point adaptive Gauss-Hermite quadrature (the tulpa_re_aghq() debias applied inside the Sigma integration), reducing the small-cluster variance attenuation for binary / low-count data. AGHQ requires a single shared grouping factor (the per-group integral must factorize); with crossed RE terms n_quad > 1 errors. When AGHQ is used the fixed effects are integrated, so the reported fixed-effect posterior is the marginal (ML-II) one rather than the joint-mode (PQL) estimate.

X_zi

Optional zero-inflation design matrix (length(y) rows), making the model a two-process mixture: each observation is a structural zero with probability plogis(X_zi beta_zi) and otherwise follows family. Paired with a zero-truncated family it is the hurdle model. The random effects enter the count predictor only, and the integration runs over the same covariance coordinates – the mixture changes the inner solve, not the parameters being integrated over. The ZI coefficients are reported alongside the count ones in coef() / vcov(), so the fixed block is ncol(X) + ncol(X_zi) wide. Needs n_quad = 1: the adaptive Gauss-Hermite inner marginal runs through a single-predictor oracle.

zi_prior_sd

Prior SD on beta_zi, keeping the logit identified where a level carries no zeros (the likelihood alone would send it to -Inf). Ignored when X_zi is NULL.

control

A named list of numerical / tuning knobs (statistical arguments stay in the signature above). Recognized entries:

  • integration: node layout, "ccd" (default, central-composite design, scales to larger total parameter count) or "grid" (full tensor product).

  • n_per_axis: points per parameter axis in the tensor grid (default 5); used only when integration = "grid".

  • span: half-width of the tensor grid in posterior standard deviations per whitened axis (default 3); grid only.

  • n_draws: posterior draws of the fixed effects synthesized from the node mixture (default 2000), exposed as draws for the generic tulpa_fit methods. The Sigma posterior is summarized directly from the integration nodes in posterior, independent of n_draws.

  • seed: optional integer seed for the fixed-effect draw synthesis.

  • diagnose_k: if TRUE (default), compute the outer Pareto k-hat accuracy diagnostic for the Gaussian proposal over the hyperparameters, returned as pareto_k. Several proposal candidates are scored and the best is kept: pareto_k_proposal_source names which one produced the reported number and pareto_k_first_pass is the k-hat of the proposal exactly as the mode-find placed it, before refinement. A large gap between the two says the placement is poor even where the verdict is fine – on a small-group binary fit the first pass runs 15 to 49 where the reported k-hat is 0.3 to 0.8.

  • k_samples: importance draws for the diagnose_k estimate (default 500). It is a precision knob: the GPD tail size is held at the fraction that default budget implies, so raising it supplies more tail ratios for the SAME estimand rather than moving the fit to a deeper quantile of the weight distribution (gcol33/tulpa#631). Each draw is one inner Laplace solve, so this sets the diagnostic's cost; the draws are solved in one compiled batch, each started from the mode at the proposal centre.

  • k_threads: thread width for those draws (NULL, the default, follows n_threads; "auto" takes the physical performance cores; an integer sets it, 1L forcing serial). Each draw is solved single-threaded and from the same start, so the k-hat does not depend on the width.

  • k_tail_points: expert override for that tail size, in upper-tail order statistics. Silently capped at 20% of the draws, beyond which body ratios enter the tail and bias the shape.

  • max_iter, tol, n_threads: inner-solve controls (see tulpa_laplace()).

  • outer_maxit: iteration budget for the mode-finding step that centres the integration grid (default 500). Applies to the Nelder-Mead simplex used from two parameters up; the one-parameter case is bracketed by Brent. Exhausting the budget warns, since the nodes are then centred on wherever the optimizer stopped.

  • checkpoint: node checkpoint/resume spec list(path = , resume = ). Each completed CCD / grid node (one inner Laplace solve) is cached to path; a resume = TRUE run loads the finished nodes and re-solves only the rest. resume = FALSE starts fresh. A file written for different data, layout, or grid is rejected (fingerprint mismatch). Default NULL (off).

  • subspace_debias: subspace debias, FALSE by default. 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), probe (the latent indices scored, default the fixed effects), closure (FALSE, TRUE, or a partial-correlation threshold: grow the set by the precision-graph neighbours it is strongly coupled to), closure_max, and the sampler budget n_iter / warmup / thin. When the selected set is non-empty, each integration node reports the selected fixed-effect coordinates from a Metropolis sample of the exact conditional along the Gaussian-conditional-mean surface, and the rest from the Gaussian conditional given them; an EMPTY set leaves the fit bit-for-bit identical to the plain path. What was selected is recorded in subspace_debias on the returned fit.

Value

A list with:

  • posterior: data frame with one row per parameter and columns mean, sd, median, ci_lo, ci_hi. Parameter names are sigma_i, rho_ij, Sigma_ij for a single block, prefixed by the block label (g.sigma_1, ...) when there are several blocks. Diagonal blocks report no rho.

  • map: the plug-in-mode summary at theta_hat (a single list(Sigma, sigma, rho) for one block, or a named list of them).

  • Sigma_mean: the weighted posterior mean of Sigma (a matrix for one block, or a named list of matrices).

  • beta, draws, means, param_names, process_info: the fixed-effect posterior from the node mixture (drives coef/confint/vcov/summary).

  • re_nodes, re_var_nodes: the per-group random-effect posterior at each integration node – conditional mean and marginal variance of every (block, group, coefficient), one row per node. ranef() reports the weights-mixture of them. NULL at n_quad > 1, whose inner marginal integrates each group out instead of conditioning on it; that fit carries ranef_unavailable (the reason) in their place.

  • re_debias_draws, re_debias_idx: present when the subspace debias selected a random-effect coordinate. The sampled draws of those coordinates on the node mixture the fixed-effect draws use, and their positions within the random-effect block. ranef() reports those rows empirically and the rest from the Gaussian mixture, recording which in its source column.

  • theta_hat, theta_grid, weights, log_marginal, n_grid, layout, n_blocks, n_coefs (vector of per-block c).

  • subspace_debias: present only when control$subspace_debias was set. idx are the corrected latent coordinates, bands the per-probed-index reliability table they were read from, closure_added what the coupling closure added, and accept the per-node Metropolis acceptance rate.

Details

Each term is one covariance block. A correlated block ((1 + x | g)) is a full Sigma = L L' parameterized by its lower Cholesky factor in log-Cholesky coordinates (the log-diagonal and the strictly-lower entries of L, c(c+1)/2 values for a c-coefficient block), which keeps Sigma positive definite for every coordinate. An uncorrelated block ((1 + x || g)) is a diagonal Sigma parameterized by its c log standard deviations. A scalar (1 | g) term is the degenerate c = 1 block. Several blocks stack their parameters into one integration vector; a single-term model is the length-1 case.

Integration nodes live in the whitened stacked-parameter space, centred at the joint marginal-likelihood mode and rotated/scaled by the Cholesky of the mode's posterior covariance (solve(Hessian)), so points track the posterior ridge. Two node layouts are available via integration:

  • "ccd" (default): a central-composite design (ccd_grid()) of 1 + 2k + 2^(k-q) points for the total k = sum_blocks parameters, with the corrected R-INLA design weights (ccd_weights()). Scales polynomially in k, where the tensor grid is exponential.

  • "grid": the full n_per_axis^k tensor product with uniform cell weights – denser and more robust to a non-Gaussian whitened posterior, but only tractable for small k.

Each node k contributes integration weight proportional to Delta_k * exp(log_marginal(Sigma_k) + log_prior_theta(theta_k)), following the INLA convention int ~ sum_k Delta_k pi(theta_k).

The two layouts also decide how the reported median and 2.5\ of each derived quantity are read off the nodes. A tensor grid's uniform cells discretize the posterior density, so the cumulative node weights are a CDF and the summary is the weighted quantile. A CCD is a moment rule: its nodes sit where they reproduce the integrand's first two moments and carry no probability mass of their own, so the summary is moment-matched instead – the first two weighted moments on each quantity's own coordinate (log for a scale or a variance, atanh for a correlation, the identity for a covariance) define a Gaussian there whose quantiles are mapped back. Scale intervals are therefore positive and asymmetric, and correlation intervals stay inside (-1, 1). The mean and sd columns are the weighted moments under either layout.

By default (hyperprior = "proper") log_prior_theta is the weakly-informative PC + LKJ hyperprior, built per block by re_cov_pc_lkj_prior() and summed over blocks (PC prior on each marginal SD via prior_sigma, LKJ prior on each correlated block's correlation matrix via eta), expressed in the same parameterization with the exact change-of-variables Jacobian. hyperprior = "flat" makes it the zero function, flat in log(theta). The likelihood it multiplies is the inner Laplace marginal, which integrates the fixed effects out together with the random effects rather than profiling them, so under "flat" the node centre theta_hat is a restricted (REML-type) estimate – REML for a gaussian response, agreeing with glmmTMB(..., REML = TRUE) rather than with the ML fit of lme4::glmer() otherwise (see tulpa_eb()). Supply a custom log_prior_theta function to override either default (then prior_sigma / eta / hyperprior are ignored); it must act on the full stacked parameter vector. tulpa_eb() shares this same objective and the same default, so tulpa_eb()$theta_hat and tulpa_re_cov_nested()$theta_hat stay the same estimate on the same data under either setting.

References

Rue, Martino & Chopin (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. JRSS-B 71(2):319-392. Lewandowski, Kurowicka & Joe (2009). Generating random correlation matrices based on vines and extended onion method. Journal of Multivariate Analysis 100(9):1989-2001.

See also

tulpa_laplace() for the inner solve; tulpa_nested_laplace() for the analogous outer integration over spatial / temporal prior hyperparameters.

Examples

# \donttest{
set.seed(1)
G <- 20L; per <- 12L; n <- G * per
grp <- rep(seq_len(G), each = per); x <- rnorm(n)
b <- cbind(rnorm(G, 0, 0.7), rnorm(G, 0, 0.5))     # random intercept + slope
eta <- -0.2 + 0.5 * x + b[grp, 1] + b[grp, 2] * x
y <- rbinom(n, 1L, plogis(eta))
re_term <- list(idx = grp, n_groups = G, n_coefs = 2L, Z = cbind(1, x),
                correlated = TRUE)
fit <- tulpa_re_cov_nested(y, rep(1L, n), cbind(1, x), re_term,
                           family = "binomial")
fit$Sigma_mean        # marginalized RE covariance
# }