tulpa_eb: Empirical-Bayes random-effect covariances

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

Empirical-Bayes random-effect covariances

Description

Estimate one or more random-effect covariances Sigma by maximizing the Laplace marginal likelihood over them (plus the hyperprior), then report the fixed effects conditional on the maximizer. This is the plug-in ("ML-II" / empirical-Bayes) counterpart of tulpa_re_cov_nested(), which integrates over Sigma instead of fixing it at the maximizer.

Usage

tulpa_eb(
  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,
  n_quad = 1L,
  marginal = FALSE,
  estimate_phi = FALSE,
  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; see tulpa_re_cov_nested() for the per-term fields.
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, never estimated: estimate_phi covers phi alone.
prior_sigma, eta	Hyperparameters of the PC + LKJ prior used when hyperprior = "proper" (see re_cov_pc_lkj_prior()); NULL (the defaults) is c(3, 0.01) and 2. Ignored when hyperprior = "flat" or log_prior_theta is supplied. When active, the prior is part of the maximized objective, so it regularizes the estimate: with few groups it is what keeps a block off the sigma = 0 boundary.
hyperprior	"proper" (default) or "flat". "flat" maximizes with log_prior_theta the zero function – an unpenalized restricted (REML-type) estimate, beta integrated out rather than profiled (see Details), which can reach the sigma = 0 boundary on small designs (see the "lower end of the search bracket" warning). "proper" builds the PC + LKJ prior from prior_sigma / eta, regularizing the estimate away from that boundary. Ignored when log_prior_theta is supplied. Must match hyperprior on the paired tulpa_re_cov_nested() call for the two to share theta_hat.
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)).
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.
n_quad	Quadrature order for the inner marginal. 1 (default) uses the joint-field Laplace inner solve. ⁠> 1⁠ refines it with n_quad-point adaptive Gauss-Hermite quadrature, which requires a single shared grouping factor; see tulpa_re_cov_nested().
marginal	Report fixed-effect intervals that carry the hyperparameter uncertainty, instead of the intervals conditional on theta_hat. The posterior for theta is taken as Gaussian around the maximizer with covariance solve(H_theta), the inner mode is linearized in theta, and the law of total variance adds ⁠J solve(H_theta) J'⁠ to the conditional covariance, where ⁠J = d mode / d theta⁠. Both H_theta and J come from one central-difference stencil over the outer objective, costing ⁠1 + 2k^2⁠ further inner solves for k hyperparameter coordinates (k is 1 for a scalar (1 | g) block and 3 for a correlated (1 + x | g) one). Default FALSE. Widens intervals; never narrows them. When the variance components themselves are the target, or the correction's two approximations look strained (a strongly skewed variance-component marginal), integrate with tulpa_re_cov_nested() instead.
estimate_phi	Estimate the family's dispersion alongside the random-effect covariances, instead of conditioning on phi. When TRUE the supplied phi is the starting value and fit$phi is the estimate, with fit$phi_estimated distinguishing the two cases. log(phi) joins the maximization as one further coordinate, carrying the exact derivative of the Laplace log-marginal with respect to it, so the cost is one more coordinate for BFGS and not a second optimization. The dispersion enters unpenalized – the hyperprior covers the covariance coordinates only – so this is the ML-II estimate of phi, not a MAP under an undeclared prior. Like the covariances it is read off the likelihood marginal over beta, so it is the restricted (REML-type) estimate: for a gaussian response, the REML residual variance. Available for every family carrying a dispersion, which is every front-door family except poisson, binomial and truncated_poisson – those have no free dispersion at all, so estimating one is a category error rather than a missing feature, and it is refused. Needs n_quad = 1. Alongside X_zi both mixture kinds are covered. A hurdle (a zero-truncated base) has zero branch log(pi), which carries no dispersion, so the base family's registered derivative is already the mixture's. Genuine zero inflation has zero branch ⁠log(pi + (1 - pi) P(Y = 0))⁠, which depends on phi through P(Y = 0) and couples it to both linear predictors; that branch is supplied for neg_binomial_2, the only untruncated mixture family here with a free dispersion. Other untruncated bases are refused rather than handed the base derivative under a model it does not describe.
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 maximization is over the same covariance coordinates – the mixture changes the inner solve, not the outer objective's parameters. 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 knobs: max_iter, tol, n_threads (inner-solve controls, see tulpa_laplace()), and outer_maxit (iteration budget for the maximization over Sigma, default 500; applies to the Nelder-Mead simplex used from two parameters up, since the one-parameter case is bracketed by Brent). Exhausting the budget warns. outer_reltol sets that maximization's convergence tolerance (default 1e-10 for the gradient-driven methods, 1e-8 for the simplex, which cannot resolve as finely); it is converted to L-BFGS-B's factr on the bounded path, so one request means the same thing whichever method runs. sigma_init supplies the starting random-effect SD – a scalar, or one per coefficient across all blocks – replacing the method-of-moments guess taken from a pilot fit at Sigma = I. Worth setting when the true scale is far from 1, where that pilot starts the search on a flat stretch, and when a run should be reproducible from its inputs rather than from a pilot fit. It is diagonal: it sets each coefficient's scale and leaves any correlation to be fitted. Two further knobs tune marginal = TRUE and are inert without it: marginal_step (the stencil step in theta space, default 1e-3) and marginal_richardson (default FALSE; evaluate the stencil at step and step / 2 and extrapolate, turning the O(step^2) truncation error into O(step^4) at twice the solves – worth it only when the inner solver's own noise sits well below the truncation error, i.e. a tight tol).

Details

Blocks, coordinates and the default hyperprior are exactly those of tulpa_re_cov_nested() – a correlated (1 + x | g) term is a full ⁠Sigma = L L'⁠ in log-Cholesky coordinates, an uncorrelated (1 + x || g) term is diagonal in log-SD coordinates, and a scalar (1 | g) term is the degenerate one-coefficient block. Both functions call the same outer objective and the same optimizer, so tulpa_eb()$theta_hat and tulpa_re_cov_nested()$theta_hat are the same estimate on the same data – which requires hyperprior to default the same way on both: "proper", the PC + LKJ prior (see re_cov_pc_lkj_prior()) at the anchor every scale axis of the engine defaults to, which at small G keeps a block off the sigma = 0 boundary. hyperprior = "flat" maximizes the marginal likelihood alone.

The marginal likelihood is marginal over the fixed effects too: the inner Laplace approximation integrates beta out together with the random effects (under beta_prior, or the inner solve's vague N(0, 100^2) default) rather than profiling it. hyperprior = "flat" is therefore a RESTRICTED, REML-type estimate, not maximum likelihood: for a gaussian response it is REML, and for a non-gaussian one it agrees with glmmTMB(..., REML = TRUE), not with the ML fit of lme4::glmer() or glmmTMB(..., REML = FALSE). On a binomial random intercept with 40 groups of 6 the three read RE SDs of 0.7931 (tulpa_eb), 0.7931 (glmmTMB REML) and 0.7594 (ML); the fixed effects at theta_hat match the REML fit's as well. The restriction is what removes ML's downward bias in the variance components when beta is estimated from the same data, and it carries over to n_quad > 1 and to estimate_phi = TRUE. Validating against glmer() should expect exactly that gap.

The reported fixed-effect covariance is the conditional one at theta_hat (solve(H_beta)). It does not include the hyperparameter uncertainty that tulpa_re_cov_nested() integrates over, so EB intervals are narrower – increasingly so as the number of groups falls. Use the nested integrator when the variance components themselves, or calibrated fixed-effect intervals, are the target; use EB when the point estimate is, or as a fast starting fit.

Value

A tulpa_fit with:

• mode, H_beta: the fixed-effect (and, on the n_quad = 1 path, random- effect) mode and the fixed-effect precision at theta_hat, driving coef/confint/vcov/summary.

• map: the Sigma / sigma / rho summary at theta_hat (a single list for one block, a named list of them for several).

• Sigma: the estimated covariance (a matrix for one block, a named list of matrices for several).

• theta_hat, log_marginal, layout, n_blocks, n_coefs.

• converged: whether the inner Newton solve at theta_hat converged.

• outer_convergence: optim's code for the maximization over Sigma (0 on success). A non-zero value also warns.

• With marginal = TRUE and a correction that formed: cov_marginal (the widened fixed-effect covariance, which vcov / summary / confint then report), cov_conditional (the solve(H_beta) they would otherwise have reported, kept so the two are comparable on one fit), H_theta and theta_cov (the outer Hessian at theta_hat and its inverse, named by theta_names), and the marginal_step / marginal_richardson actually used. All absent when the correction was not requested or could not be formed, so is.null(fit$cov_marginal) tests whether the reported intervals are marginal. A requested correction that fails warns and leaves the conditional covariance in place.

References

Casella (1985). An introduction to empirical Bayes data analysis. The American Statistician 39(2):83-87. Rue, Martino & Chopin (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. JRSS-B 71(2):319-392.

See Also

tulpa_re_cov_nested() to integrate over Sigma rather than fix it; tulpa_laplace() for the inner solve.

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)
y <- rpois(n, exp(0.3 + 0.5 * x + b[grp]))
re_term <- list(idx = grp, n_groups = G, n_coefs = 1L)
fit <- tulpa_eb(y, NULL, cbind(1, x), re_term, family = "poisson")
fit$map$sigma          # empirical-Bayes RE standard deviation

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