LatentGaussianModels.jl

LatentGaussianModels

Julia-native implementation of latent Gaussian models and their approximate-Bayesian inference (INLA-style). Consumes the sparse precision numerics from GMRFs.jl; owns the likelihood, hyperprior, component, and inference-strategy abstractions.

See CLAUDE.md and plans/plan.md for the design document.

DEFAULT_BYM2_PHI_ALPHA

Default right-tail probability α for the PC prior on BYM2's mixing parameter φ, given P(φ < U) = α with U = 0.5.

Set to 2/3 from Riebler et al. (2016). Current R-INLA may ship a different default — some releases use α = 0.5. This constant must be reconciled against the running R-INLA in scripts/verify-defaults/bym2_phi.R before v0.1 release; until then the divergence is documented in plans/defaults-parity.md.

AR1(n; precprior = PCPrecision(),
    ρprior = _NormalAR1ρ(0.0, 1.0),
    τ_init = 0.0,
    fix_τ::Bool = false)

Stationary AR(1) latent field of length n. Two hyperparameters by default:

• τ (marginal precision) on internal scale θ₁ = log(τ).
• ρ ∈ (-1, 1) (lag-1 correlation) on internal scale θ₂ = atanh(ρ) = ½ log((1 + ρ)/(1 - ρ)).

precprior is a scalar prior on θ₁ (default PC prior on τ). ρprior defaults to a Normal(0, σ=1) on θ₂ — the simple diffuse-on-Fisher-z choice used in many INLA tutorials. R-INLA's built-in f(., model="ar1") default differs: a Normal on logit(ρ) = 2 atanh(ρ) with precision 0.15 (i.e. σ ≈ 2.582 on logit(ρ), equivalently σ ≈ 1.291 on atanh(ρ)) — close to but not identical to our default. Pass ρprior = PCCor1(U, α) to opt into the textbook PC prior on ρ with reference at ρ = 1 (Sørbye-Rue 2017), or any other AbstractHyperPrior to override.

Parameterization follows Rue & Held (2005, Eq. 1.39) so that Var(x_i) = 1/τ at every index.

Fixing the precision

• τ_init: initial value of log τ on the internal scale. Used as the optimisation start, and as the held-fixed value when fix_τ = true.
• fix_τ: when true, drop the precision from the θ vector and hold it at τ = exp(τ_init). The component's hyperparameter vector is then just [atanh(ρ)] (length 1). Useful for separable space-time composition with KroneckerComponent, where R-INLA's f(field, group, control.group = list(model = "ar1")) matches the fix_τ = true, τ_init = 0 pattern (prec = 1 is implicit on the group dimension for identifiability against the spatial precision).

AbstractHyperPrior

Prior on a scalar hyperparameter. All priors are expressed on the internal unconstrained scale θ ∈ ℝ; the transformation from the user-facing parameter (precision τ > 0, correlation ρ ∈ (-1, 1), mixing φ ∈ [0, 1]) to θ is fixed by the component that owns the hyperparameter (see plans/defaults-parity.md).

Concrete subtypes implement:

• log_prior_density(prior, θ) — log π(θ) on the internal scale, including the Jacobian of the internal-to-user transformation.
• user_scale(prior, θ) — convert internal θ to the user-facing parameter for reporting.
• prior_name(prior) -> Symbol — short identifier for printing.

This type is for scalar priors. Multi-dimensional priors (e.g. the PC prior on (σ, range) for SPDE) live in INLASPDE.jl because they are inherently coupled.

AbstractIIDND <: AbstractLatentComponent

Umbrella type for the multivariate-IID family. Two concrete subtypes (see ADR-022):

• IIDND_Sep{N} — separable parameterisation: N marginal precisions
  • N(N-1)/2 correlations, each carrying its own scalar
  AbstractHyperPrior. This matches R-INLA's f(., model = "2diid"/"iid3d", ...) defaults.
• IIDND_Joint{N} — joint Λ parameterisation with a single AbstractJointHyperPrior (Wishart / InvWishart). Lands in PR-1c.

The latent vector has length n · N and stores the N blocks of length n consecutively: x = (x^{(1)}, x^{(2)}, …, x^{(N)}) where each x^{(j)} ∈ ℝ^n. Joint precision is Q = Λ ⊗ I_n for both subtypes.

AbstractInferenceResult

Return type for fit.

AbstractInferenceStrategy

Dispatch type for the fit entry point. Concrete strategies:

• Laplace — fit at fixed θ, return the Gaussian approximation x | θ, y ≈ N(x̂, (Q + A' D A)⁻¹) (D from the likelihood Hessian).
• EmpiricalBayes — plug-in estimate θ̂ = argmax π(θ | y) via the outer Laplace log-marginal, then Laplace at θ̂.
• INLA — the full thing (deferred).

AbstractIntegrationScheme

Strategy for numerically integrating functions of θ against the Gaussian approximation at the mode θ̂ with covariance Σ = H⁻¹ (where H is the Hessian of the negative log-posterior at θ̂).

A scheme produces a list of design points θ_k and weights w_k such that

∫ f(θ) π(θ | y) dθ  ≈  ∑_k w_k f(θ_k)

is a good quadrature rule when π(θ | y) is well approximated by a Gaussian in a log-posterior sense. Implementations return (points, weights) with points::Vector{Vector{Float64}} and weights::Vector{Float64} summing to 1.

AbstractLatentComponent

The Bayesian wrapper around a GMRF — the LGM equivalent of R-INLA's f(...) term. A concrete component composes a GMRFs.AbstractGMRF (or builds one lazily from θ) with an AbstractHyperPrior and any constraint metadata.

Required methods:

• length(c) — dimension of the component's latent field.
• initial_hyperparameters(c) -> Vector — internal-scale start values.
• nhyperparameters(c) -> Int — convenience: length(initial_hyperparameters(c)).
• precision_matrix(c, θ) — sparse precision matrix as a function of this component's slice of θ.
• log_hyperprior(c, θ) — log-prior density at θ (on internal scale).

Optional (with default implementations):

• prior_mean(c, θ) — defaults to zeros.
• constraints(c) — defaults to NoConstraint().
• graph(c) — the underlying GMRF graph if any.
• gmrf(c, θ) — the AbstractGMRF at θ, default built from precision_matrix(c, θ) via a Generic0GMRF wrapper.

AbstractLikelihood

The observation model y | η, θ. η is the linear predictor — the sum of latent-component contributions — and θ carries any likelihood-specific hyperparameters (Gaussian precision, NegBinomial size, Gamma shape). Likelihoods without hyperparameters ignore θ.

Required methods:

• log_density(ℓ, y, η, θ) -> Real — scalar log p(y | η, θ).
• ∇_η_log_density(ℓ, y, η, θ) -> AbstractVector{<:Real} — gradient w.r.t. η, same length as y.
• ∇²_η_log_density(ℓ, y, η, θ) -> AbstractVector{<:Real} — diagonal of the Hessian w.r.t. η, same length as y. The likelihood factorises across observations so the Hessian is diagonal.
• link(ℓ) -> AbstractLinkFunction.
• nhyperparameters(ℓ) -> Int — number of scalar hyperparameters attached to the likelihood (default 0).

Return-type conventions are load-bearing for the inner Newton hot path: ∇_η_log_density and ∇²_η_log_density return a plain AbstractVector{<:Real}, not Diagonal, not SparseMatrixCSC. The caller wraps with Diagonal(...) only where algebraically required.

AbstractLinkFunction

Link between the linear predictor η = x[1] + x[2] + … and the likelihood's mean parameter μ. Concrete subtypes implement:

• inverse_link(g, η) — μ as a function of η (scalar → scalar).
• ∂inverse_link(g, η) — derivative dμ/dη.
• ∂²inverse_link(g, η) — second derivative (needed for closed-form Gaussian-likelihood fit with a non-identity link).

Link functions are stateless singletons; the concrete structs carry no fields.

AbstractMarginalStrategy

Dispatch type for the per-coordinate posterior marginal density of x_i | y and, in the INLA integration stage, for the per-θ approximation of x_mean / x_var. Mirrors R-INLA's control.inla$strategy.

Concrete strategies (see ADR-026):

• Gaussian — Gaussian centred at the Newton mode, with constraint-corrected Laplace marginal variance. R-INLA's strategy = "gaussian".
• SimplifiedLaplace — Newton mode plus the Rue-Martino mean shift in the integration stage; Edgeworth first-order skewness correction in the per-coordinate density. Reduces to Gaussian when the likelihood third derivative ∇³_η log p is zero. R-INLA's strategy = "simplified.laplace".
• FullLaplace — per-x_i refitted Laplace via constraint injection. The reference quality strategy when the latent posterior is sharply non-Gaussian. R-INLA's strategy = "laplace". PR-3 wires it into posterior_marginal_x only; the integration-stage hook delegates to Gaussian until PR-4 ships the rank-1 update.

The mean-shift facet of SimplifiedLaplace (integration-stage) and its density-skew facet (per-coordinate marginals) are independent — they happen at different points in the pipeline and are dispatched on the same type but through different internal hooks. See ADR-016 and ADR-026.

AbstractObservationMapping

Abstract type for the linear map from the stacked latent vector x into the linear predictor η = mapping * x. The seam introduced by ADR-017 to support joint-likelihood / multi-response models without inheriting inla.stack's ergonomic baggage.

Concrete subtypes

• IdentityMapping — η = x. Areal models with A = I.
• LinearProjector — wraps an A::AbstractMatrix. The v0.1 default; constructed automatically when LatentGaussianModel is called with an AbstractMatrix.
• StackedMapping — block-row stack of per-likelihood mappings. Lands fully in Phase G PR2 (multi-likelihood); the struct ships now so the type hierarchy is complete.
• KroneckerMapping — separable A_space ⊗ A_time. Forward/adjoint multiplies via vec(A_time · X · A_space'); the dense Kronecker product is never materialized.

Required interface

Every concrete subtype implements:

apply!(η, mapping, x)            # η .= mapping * x, in place
apply_adjoint!(g, mapping, r)    # g .= mappingᵀ * r, in place
nrows(mapping) -> Int            # row count of the implicit matrix
ncols(mapping) -> Int            # column count = length(x)

Optional (default returns 1 for single-block mappings):

likelihood_for(mapping, i) -> Int   # which likelihood-block owns row i

Convenience

Defined here for any subtype satisfying the required interface:

• Base.:*(mapping, x) — allocating forward multiply.
• Base.size(mapping), Base.size(mapping, dim).
• as_matrix — concrete materialization for the ~3 inference sites that need a sparse/dense A (multi-RHS solves; A' D A Hessian quadratic form).

References

• ADR-017 in plans/decisions.md.

BYM(graph; hyperprior_iid = PCPrecision(), hyperprior_besag = PCPrecision(),
    scale_model = true)

Classical Besag–York–Mollié (1991) decomposition of a spatial random effect into an iid (unstructured) part and an intrinsic CAR (Besag) part — R-INLA's model = "bym". The latent vector has length 2n and is arranged as x = [v; b]:

• v — IID component with precision τ_v. Length n.
• b — intrinsic Besag component with precision τ_b. Length n, subject to a sum-to-zero constraint per connected component (Freni-Sterrantino et al. 2018).

The combined effect entering the linear predictor is v + b; the user projector A is responsible for that summation (typical pattern: A_block = [I_n | I_n]).

Internal hyperparameters are θ = [log τ_v, log τ_b] (two precisions, no mixing parameter — that is the BYM2 reparameterisation).

With scale_model = true (default) the Sørbye-Rue (2014) geometric- mean scaling is applied to the Besag block, matching R-INLA ≥ 17.06.

The joint precision is block-diagonal:

Q(τ_v, τ_b) = blockdiag(τ_v · I_n, τ_b · R)

where R = c · L is the (optionally scaled) Besag Laplacian.

BYM2(graph; hyperprior_prec = PCPrecision(), hyperprior_phi = nothing,
     U = 0.5, α = DEFAULT_BYM2_PHI_ALPHA)

Riebler et al. (2016) BYM2 reparameterisation of the spatial + unstructured random effect. The latent vector has length 2n and is arranged as x = [b; u]:

• b = (1/√τ)(√(1-φ) v + √φ u) — the combined mixed effect that enters the linear predictor. Length n.
• u — the scaled Besag spatial component with geometric-mean marginal variance 1 on its non-null subspace. Length n, subject to a sum-to-zero constraint per connected component (Freni-Sterrantino et al. 2018).

τ > 0 is the overall precision of b; φ ∈ [0, 1] is the mixing parameter (φ = 0 is pure IID, φ = 1 is pure spatial). Internal hyperparameters are θ = [log(τ), logit(φ)].

The joint precision of x = [b; u] is

Q(τ, φ) = [  (τ/(1-φ)) I_n           -(√(τφ)/(1-φ)) I_n        ;
             -(√(τφ)/(1-φ)) I_n      R_scaled + (φ/(1-φ)) I_n  ]

where R_scaled = c · L is the Sørbye-Rue–scaled Besag Laplacian. Derivation: change of variables from (v, u) (with independent quadratic forms v'v + u' R_scaled u) to (b, u) via v = (√τ b - √φ u)/√(1-φ).

If hyperprior_phi = nothing, a PCBYM2Phi(U, α, γ) is constructed internally using the eigenvalues γ of R_scaled on its non-null subspace. Defaults match Riebler et al. (2016); see plans/defaults-parity.md for the ⚠ unverified-default note on α.

Besag(graph; hyperprior = PCPrecision(), scale_model = true)

Intrinsic CAR (Besag) component on the supplied graph (AbstractGMRFGraph or adjacency matrix). One hyperparameter on the internal scale θ = log(τ).

With scale_model = true (default) the Sørbye-Rue (2014) geometric- mean scaling is applied, matching R-INLA ≥ 17.06. Carries one sum-to-zero constraint per connected component of graph (Freni-Sterrantino et al. 2018).

BetaBinomialLikelihood(n_trials; link = LogitLink(),
                       hyperprior = GaussianPrior(0.0, sqrt(2.0)))

y_i | n_i, η_i, ρ ∼ BetaBinomial(n_i, μ_i (1−ρ)/ρ, (1−μ_i)(1−ρ)/ρ) in R-INLA's mean-overdispersion parameterisation: mean E[y_i] = n_i μ_i with μ_i = g⁻¹(η_i) ∈ (0, 1), variance n_i μ_i (1−μ_i) (1 + (n_i−1) ρ). The overdispersion ρ ∈ (0, 1) is a single likelihood hyperparameter carried on the internal scale θ = logit(ρ). Equivalently, s := (1−ρ)/ρ = exp(−θ) is the Beta-prior "size" / effective sample, and the canonical α/β form is α = μ s, β = (1−μ) s.

Density (with α = μ s, β = (1−μ) s, s = exp(−θ)):

p(y | n, μ, ρ) = C(n, y) · B(y + α, n − y + β) / B(α, β)

Currently only the canonical LogitLink is supported. n_trials is a fixed vector of positive integers; pass per-observation trial counts at construction.

The default hyperprior matches R-INLA's family = "betabinomial" default (gaussian(0, prec = 0.5) on logit(ρ), i.e. mean 0 and σ = √2 on the internal scale). For PR-level oracle fits, override on both sides to keep the prior bit-for-bit identical.

BetaLikelihood(; link = LogitLink(), hyperprior = GammaPrecision(1.0, 0.01))

y_i | η_i, φ ∼ Beta(μ_i φ, (1 - μ_i) φ) in R-INLA's mean-dispersion parameterisation: mean μ_i = g⁻¹(η_i) ∈ (0, 1), variance μ_i (1 - μ_i) / (φ + 1). The dispersion φ > 0 is a single likelihood hyperparameter carried on the internal scale θ = log(φ).

Density (for y ∈ (0, 1)):

p(y | μ, φ) = Γ(φ) / [Γ(μ φ) Γ((1 - μ) φ)] ·
              y^(μ φ - 1) · (1 - y)^((1 - μ) φ - 1)

Currently only the canonical LogitLink is supported. Boundary values (y = 0 or y = 1) yield -Inf — these need a separate zero-/one- inflated family.

The default hyperprior matches R-INLA's family = "beta" default (loggamma(1, 0.01) on log(φ)). For PR-level oracle fits, override on both sides to keep the prior bit-for-bit identical.

BetaPrior(a, b)
BetaPrior(d::Distributions.Beta)

Generic Beta(a, b) prior on a user-scale bounded-ratio parameter p ∈ (0, 1). Internally expressed via the logit transform θ = logit(p) with Jacobian |dp/dθ| = p(1-p), so on the internal scale

log π_θ(θ) = a · log(p) + b · log(1 - p) - log B(a, b),
p          = σ(θ) = 1 / (1 + exp(-θ)).

This is a Distributions.jl-friendly entry point that complements LogitBeta: both implement the same prior with the same numerically stable softplus form, but BetaPrior accepts a Distributions.Beta(a, b) instance directly so authors can pass a shared distribution object across components. Functionally BetaPrior(a, b) ≡ LogitBeta(a, b); pick whichever reads cleanest at the call site.

Defaults a = b = 1 give the uniform-on-p prior (standard logistic on θ).

Example

using Distributions
hyper = BetaPrior(Beta(2.0, 5.0))                    # via Distributions
hyper = BetaPrior(2.0, 5.0)                          # equivalent
log_prior_density(hyper, 0.0)                        # at p = 1/2

BinomialLikelihood(n_trials; link = LogitLink())

y_i | η_i ∼ Binomial(n_trials[i], p_i) with p_i = g⁻¹(η_i). Canonical link is logit.

n_trials must be a vector of positive integers of length matching y at fit time.

CCD(; f0 = nothing)

Central composite design in the eigenbasis of H, after Rue, Martino & Chopin (2009, §6.5). For m = dim(θ) ≥ 2 the design consists of

• 1 center point,
• 2m axial points at ±f0 on each eigenaxis,
• 2^(m-p) fractional-factorial corner points with coordinates in {-1, +1} (we use p = 0 — full factorial — up to m = 5; above that callers should prefer Grid or GaussHermite).

Weights are chosen so that the rule integrates a quadratic + pure fourth-order term exactly against the standard normal. The default f0 = sqrt(m + 2) is the Rue-Martino-Chopin recommendation.

Falls back to Grid(n_per_dim = 7) for m ≤ 1.

Censoring

Per-observation censoring mode for survival likelihoods. One of NONE, RIGHT, LEFT, INTERVAL.

Constructable from Symbol for keyword-call ergonomics: Censoring(:none) == NONE. Storage is the enum (one byte per row), not the Symbol, so the inner Newton loop's per-row branch is a cheap integer compare.

See ADR-018 in plans/decisions.md for the full design rationale.

Copy(source_indices::UnitRange{Int};
     β_prior = GaussianPrior(1.0, 0.1),
     β_init  = 1.0,
     fixed   = false)

Specifier for one Copy contribution attached to a CopyTargetLikelihood. The receiving likelihood adds β * x[source_indices][i] to its η-row i (one-to-one row mapping between observations and source-component slots).

source_indices is the range of the source component inside the model's stacked latent vector x. For the common case where the source is the j-th declared component, that is m.latent_ranges[j] after model construction; users can compute it inline with cumsum(length, components_before) or the helper component_range once the model exists.

β_prior is on the unconstrained user scale (β = β-internal). Default is R-INLA's gaussian(mean=1, prec=100) — GaussianPrior(1.0, 0.1) in this package's parameterisation.

If fixed = true, β is held at β_init for the duration of inference (no hyperparameter slot is allocated). Useful for the closed-form regression test that pins β to 1.0 to recover the unscaled-share result, and for sensitivity studies that want to ablate β.

CopyTargetLikelihood(base::AbstractLikelihood, copies::Tuple{Vararg{Copy}})
CopyTargetLikelihood(base, copy::Copy)

Wrap base with one or more Copy contributions. Forwards every required AbstractLikelihood method to base after the inference loop folds the Copy contributions into η via add_copy_contributions!.

Each non-fixed Copy adds one β slot to the wrapper's hyperparameter vector, in the order they are supplied. The full layout is [θ_base..., β_1, β_2, ...], matching how the model constructor accounts for nhyperparameters(ℓ).

Example

# Two-block joint Gaussian + Poisson model where the Poisson arm
# *copies* the Gaussian arm's IID random effect with an estimated β.
ℓ_g = GaussianLikelihood()
ℓ_p_with_copy = CopyTargetLikelihood(
    PoissonLikelihood(; E = fill(1.0, n)),
    (Copy(2:(n + 1); β_prior = GaussianPrior(1.0, 0.5)),))
model = LatentGaussianModel(
    (ℓ_g, ℓ_p_with_copy),
    (Intercept(), IID(n)),
    StackedMapping(...))

CoxphAugmented(y, E, subject, interval, breakpoints, n_subjects, n_intervals)

Result of inla_coxph data augmentation. Contains the per-row augmented Poisson response, exposure offsets, subject and baseline-interval indices, the breakpoints used, and the original subject / interval counts.

Pass (y, E) to a PoissonLikelihood (E as the offset), and use coxph_design to build the joint design matrix that pairs the piecewise-constant baseline-hazard component with the original covariates.

EmpiricalBayes(; laplace = Laplace(), θ0 = nothing, optim_options = NamedTuple())

Empirical Bayes inference: plug-in estimate of the hyperparameters,

θ̂ = argmax_θ log π(θ | y)  ≈  argmax_θ [ log p(y | θ) + log π(θ) ],

using the Laplace approximation for log p(y | θ). The final result is the Laplace fit at θ̂ plus the mode and Hessian of the hyperparameter log-posterior (for reporting only — no integration is performed, by contrast with full INLA).

Optimisation is performed through the SciML Optimization.jl frontend; the default algorithm is OptimizationOptimJL.LBFGS(). Custom options may be supplied via optim_options (forwarded to Optimization.solve).

EmpiricalBayesResult <: AbstractInferenceResult

• θ̂::Vector{Float64} — mode of log π(θ | y) on the internal scale.
• laplace::LaplaceResult — Laplace fit at θ̂.
• log_marginal::Float64 — log p(y | θ̂).
• optim_result — raw Optimization.jl solution (for diagnostics).

ExponentialLikelihood(; link = LogLink(), censoring = nothing, time_hi = nothing)

y_i | η_i ∼ Exponential(rate = g⁻¹(η_i)). With the canonical LogLink, the rate is λ_i = exp(η_i), the mean is 1/λ_i = exp(-η_i), and y_i > 0 is required.

censoring is an optional AbstractVector{Censoring} of length length(y); when nothing (default), every observation is uncensored. For INTERVAL rows, time_hi[i] is the upper bound (with y[i] the lower bound); time_hi is otherwise unread and may be nothing.

No likelihood hyperparameters. θ is ignored.

Closed-form derivatives are provided for LogLink for all four censoring modes; see ADR-018 for the contract.

Example

# Right-censored survival data
ℓ = ExponentialLikelihood(censoring = [NONE, RIGHT, NONE, RIGHT])

# Interval-censored
ℓ = ExponentialLikelihood(
    censoring = [NONE, INTERVAL],
    time_hi   = [0.0, 5.0],
)

FixedEffects(p; prec = 1.0e-3)

p-dimensional block of fixed effects β with iid Normal prior β_j ~ N(0, prec⁻¹). No hyperparameters. prec = 1e-3 matches R-INLA's default control.fixed = list(prec = 0.001).

FullLaplace()

Per-x_i refitted Laplace marginal strategy via constraint injection. R-INLA's strategy = "laplace". Definition + helpers in inference/full_laplace.jl; this declaration lives in abstract.jl alongside Gaussian and SimplifiedLaplace so that method tables in inference/marginals.jl and inference/inla.jl can dispatch on ::FullLaplace at load time.

See AbstractMarginalStrategy and ADR-026.

GEVLikelihood(; link = IdentityLink(),
              precision_prior = GammaPrecision(1.0, 5.0e-5),
              shape_prior = GaussianPrior(0.0, 2.0),
              xi_scale = 0.1,
              weights = nothing)

Generalised Extreme Value (GEV) observation model in R-INLA's family = "gev" parameterisation. The CDF is

F(y; η, τ, ξ) = exp(− [1 + ξ √(τ s) (y − η)]^(−1/ξ))

defined on the half-line 1 + ξ √(τ s) (y − η) > 0, where η is the linear predictor (location), τ > 0 the precision (so the GEV scale is σ = 1/√(τ s)), s > 0 a fixed per-observation weight (weights[i], default 1), and ξ ∈ ℝ the shape parameter.

Two hyperparameters on the internal scale:

• θ[1] = log τ
• θ[2] = ξ / xi_scale — internal scaling of the shape parameter (default xi_scale = 0.1, matching R-INLA's gev.scale.xi).

Defaults match R-INLA's family = "gev":

• precision_prior = GammaPrecision(1, 5e-5) ↔ loggamma(1, 5e-5) on log τ.
• shape_prior = GaussianPrior(0, 2.0). R-INLA's documented default is gaussian(0, prec = 25), applied to the user-scale ξ, not to the internal θ[2]. Reparameterising (since ξ = xi_scale · θ[2]) gives a Gaussian on θ[2] with prec = 25 · xi_scale² = 0.25 for the default xi_scale = 0.1, i.e. σ = 2.0. If you change xi_scale, scale this σ by 0.1/xi_scale to keep the user-scale prior identical.
• xi_scale = 0.1.

R-INLA's gev family is marked "disabled" in current releases in favour of bgev. The two share the body parameterisation; bgev adds a smooth tail blend to guarantee finite support. Use this likelihood for body-of-distribution fits or replicating older R-INLA results.

Note: GEV is not globally log-concave. The inner Newton step may require damping when the initial η falls near the support boundary; ensure initial_η is well inside 1 + ξ √(τ s)(y − η) > 0 for all observations.

GammaLikelihood(; link = LogLink(), hyperprior = GammaPrecision(1.0, 5.0e-5))

y_i | η_i, φ ∼ Gamma(μ_i, φ) in R-INLA's mean-precision parameterisation: mean μ_i = g⁻¹(η_i), variance μ_i² / φ. The precision φ > 0 is carried as a single hyperparameter on the internal scale θ = log(φ).

Density (with shape a = φ, rate b = φ/μ):

p(y | μ, φ) = (φ/μ)^φ / Γ(φ) · y^(φ-1) · exp(-φ y / μ),  y > 0

The default hyperprior matches R-INLA's family = "gamma" default (loggamma(1, 0.00005) on log(φ)). Currently only the canonical LogLink is supported. The likelihood requires strictly positive observations.

GammaPrecision(shape, rate)

Classic Gamma prior on τ, τ ~ Gamma(shape, rate). Internal scale θ = log(τ). Included for backwards compatibility with older R-INLA workflows; PC priors are preferred.

GammaSurvLikelihood(; link = LogLink(), censoring = nothing,
                      time_hi = nothing,
                      hyperprior = GammaPrecision(1.0, 5.0e-5))

Gamma survival likelihood matching R-INLA's family = "gammasurv". Uses R-INLA's mean-precision parameterisation: with shape a = φ and rate b = φ/μ where μ = exp(η) under the canonical LogLink,

f(t)     = b^a t^(a-1) exp(-b t) / Γ(a)
S(t)     = Γ(a, b t) / Γ(a)             (regularised upper-incomplete)
F(t)     = γ(a, b t) / Γ(a)             (regularised lower-incomplete)

The single hyperparameter is the precision φ > 0, carried internally as θ_ℓ = [log φ]. The default hyperprior matches R-INLA's family = "gamma" default — loggamma(1, 5e-5) on log φ, encoded here as GammaPrecision(1.0, 5.0e-5).

censoring is an optional AbstractVector{Censoring} of length length(y); when nothing (default), every observation is uncensored. For INTERVAL rows, time_hi[i] is the upper bound (with y[i] the lower bound).

Closed-form derivatives are provided for LogLink for all four censoring modes via the unitless ratio g(x) = x · δ(x) / Γ(a, x) (the gamma analog of the inverse-Mills ratio), where δ(x) = x^(a-1) exp(-x) / Γ(a) and x = b t. The analogous h(x) = x · δ(x) / γ(a, x) covers the LEFT/INTERVAL branches. See ADR-018 for the contract.

GaussHermite(; n_per_dim = 5)

Tensor-product Gauss-Hermite quadrature in the eigenbasis of H. Exact for polynomial integrands of degree 2 n_per_dim - 1 against the Gaussian weight.

Gaussian()

Marginal strategy: per-θ Gaussian centred at the Newton mode, with constraint-corrected Laplace marginal variance. R-INLA's strategy = "gaussian". Default for both INLA (integration stage) and posterior_marginal_x (per-coordinate density).

See AbstractMarginalStrategy.

GaussianLikelihood(; link = IdentityLink(), hyperprior = PCPrecision(1.0, 0.01))

y_i | η_i, τ ∼ N(μ_i, τ⁻¹) with μ_i = g⁻¹(η_i) and observation precision τ carried as a single hyperparameter on the internal scale log(τ). The default hyperprior is the PC prior on standard deviation σ = τ^(-1/2) with default P(σ > 1) = 0.01 (matching R-INLA's family = 'gaussian' default).

With IdentityLink, η = μ, which is the common case and admits a closed-form posterior when paired with a Gaussian latent field.

GaussianPrior(μ, σ)

A Normal(μ, σ) prior on the internal scale θ (no transform). Matches R-INLA's gaussian(μ, prec) with σ = 1/√prec.

Used for hyperparameters that already live on an unconstrained scale — e.g. θ = logit(p) for the inflation probability of zero-inflated likelihoods (R-INLA default gaussian(0, prec=1)), or θ = log(α) for the type-2 inflation-intensity hyperparameter.

Defaults μ = 0, σ = 1 correspond to a unit-precision Gaussian on the internal scale — diffuse but proper.

Generic0(R; rankdef = 0, scale_model = false, hyperprior = PCPrecision(),
        constraint = nothing)

User-supplied structure-matrix component, R-INLA's model = "generic0". The latent vector has length n = size(R, 1) and precision Q = τ · R, with one hyperparameter on the internal scale θ = log(τ).

R must be symmetric and non-negative definite. rankdef is the dimension of the null space (rank deficiency); supply this when known because we cannot infer it cheaply for a general R. With scale_model = true the Sørbye-Rue (2014) geometric-mean scaling is applied to R once at construction (matching R-INLA's scale.model = TRUE).

constraint is an optional LinearConstraint attached to the component (e.g. sum-to-zero for an intrinsic R); if nothing, no constraint is applied. For intrinsic components the user is responsible for providing the appropriate constraint.

The default hyperprior is PCPrecision(). R-INLA's generic0 default is loggamma(1, 5e-5) — pass hyperprior = GammaPrecision(1.0, 5.0e-5) to match.

Generic1(R; rankdef = 0, hyperprior = PCPrecision(), constraint = nothing)

Like Generic0 but rescales R at construction so that its largest eigenvalue is 1. Matches the structure-matrix normalisation used by R-INLA's generic1 before multiplying by τ. Precision is Q = τ · R̃ with R̃ = R / λ_max(R) and τ = exp(θ[1]).

Unlike R-INLA's full generic1, this component does not introduce the mixing parameter β ∈ (0, 1) — only the eigenvalue rescaling. The β-mixing flavour is deferred; see packages/GMRFs.jl/plans/plan.md M2 note.

R must be symmetric and non-negative definite with a strictly positive largest eigenvalue. rankdef is the dimension of the null space; if non-zero the caller is responsible for supplying a matching constraint (e.g. sum-to-zero).

The default hyperprior is PCPrecision(). R-INLA's generic1 default is loggamma(1, 5e-5) — pass hyperprior = GammaPrecision(1.0, 5.0e-5) to match.

Generic2(C; rankdef = 0, scale_model = false,
        hyperprior_τv = PCPrecision(), hyperprior_τu = PCPrecision(),
        constraint = nothing)

R-INLA's model = "generic2" — hierarchical Gaussian model on a joint latent vector [u; v] of length 2n (where n = size(C, 1)):

v ~ N(0, (τ_v C)⁻¹),
u | v ~ N(v, (τ_u I_n)⁻¹).

The joint precision (eq. 1 of the R-INLA generic2 doc) is the 2n × 2n block matrix

Q = [[ τ_u I,         -τ_u I       ],
     [ -τ_u I,        τ_u I + τ_v C ]]

with two hyperparameters on the internal scale matching R-INLA's parameterisation (theta1, theta2):

θ[1] = log τ_v   (precision multiplying C),
θ[2] = log τ_u   (precision of the conditional noise u | v).

Used for shared/specific decompositions — the u block carries the linear-predictor-relevant signal, the v block carries the structured prior mediated by C. With C the ICAR Laplacian this is a joint-precision-parameterised cousin of BYM (which exposes u + v directly via summation rather than the joint precision).

Arguments

• C: symmetric, non-negative-definite "structure" matrix.
• rankdef: dimension of null(C). The rank of Q is then 2n − rankdef, and the caller is responsible for supplying a matching LinearConstraint when rankdef > 0.
• scale_model: if true, applies the Sørbye-Rue (2014) geometric-mean scaling to C once at construction (matches R-INLA's scale.model = TRUE).
• hyperprior_τv: scalar prior on θ[1] = log τ_v. Default PCPrecision(). R-INLA's default is loggamma(1, 5e-5) — pass GammaPrecision(1.0, 5.0e-5) to match.
• hyperprior_τu: scalar prior on θ[2] = log τ_u. Default PCPrecision(). R-INLA's default is loggamma(1, 1e-3) (looser shape than τ_v) — pass GammaPrecision(1.0, 1.0e-3) to match.
• constraint: optional LinearConstraint over the joint length-2n latent. For intrinsic C (e.g. ICAR Laplacian, rankdef = 1) the nullspace of Q is span([1_n; 1_n]) — a single sum-to-zero on either the u or the v block alone breaks the ambiguity.

Grid(; n_per_dim = 5, span = 3.0,
       stdev_corr_pos = nothing, stdev_corr_neg = nothing)

Tensor-product midpoint grid in the eigenbasis of H. n_per_dim points per dimension on [-span, span] (in units of standard deviations of the Gaussian approximation), weights derived from the standard-normal density.

Recommended only for dim(θ) ≤ 3; cost is n_per_dim^dim(θ).

Asymmetric skewness corrections

stdev_corr_pos / stdev_corr_neg are length-m = dim(θ) per-axis stretch factors applied to the eigen-basis grid coordinates z before mapping back to the native θ-parameterisation: along eigen-axis k, a point with z_k > 0 is placed at θ̂ + halfσ * (z_k * stdev_corr_pos[k]) e_k and a point with z_k < 0 at θ̂ + halfσ * (z_k * stdev_corr_neg[k]) e_k. Quadrature weights are unchanged from the standard-normal weights — the importance- sampling reweight π̂(θ_k|y) / q(θ_k) performed by the surrounding INLA pipeline absorbs any proposal-mismatch from the asymmetric placement, while the asymmetric placement itself concentrates points where π(θ|y) actually has mass on each side.

The reference values defining the stretch are σ⁺k / σ⁻k from Rue–Martino–Chopin (2009) §6.5: σ⁺_k is the positive-direction step along axis k at which log π̂(θ̂) - log π̂(θ̂ + σ⁺_k e_k) = ½, and similarly for σ⁻_k in the negative direction. For an exactly Gaussian posterior σ⁺k = σ⁻k = √λk (the corresponding eigenvalue of Σ), so the stretches reduce to 1. See [`computeskewnesscorrections](@ref) for the helper that estimates them, and theskewnesscorrectionkeyword on [INLA`](@ref) for the opt-in flag in the standard pipeline.

Pass nothing (the default) on either side to disable the correction on that side; mixed behaviour (only one side specified) is not allowed — use ones(m) to leave the other side at 1.0.

Group(components::AbstractVector{C}) where {C <: AbstractLatentComponent}
Group(factory, group_id::AbstractVector{<:Integer}; kwargs...)

Stack a list of inner components, all sharing a single hyperparameter block, with potentially different sizes per group. Generalises Replicate to the non-uniform case — Group([AR1(5), AR1(5), AR1(5)]) and Replicate(AR1(5), 3) produce equivalent priors, but Group([AR1(10), AR1(7), AR1(12)]) cannot be expressed as a Replicate.

Hyperparameter sharing

Like Replicate, the wrapper exposes the same nhyperparameters, initial_hyperparameters, and log_hyperprior as the first inner component, and the same θ slice flows verbatim into each group's precision_matrix(component, θ) call. All inner components must have the same number of hyperparameters; the constructor validates this. Beyond hyperparameter count, ensuring that the inner components share prior structure (i.e. they were all built with the same prior types and parameters) is the user's responsibility.

Public API

• Group(components): explicit per-group component vector. The Vector's element type is locked at construction.
• Group(factory, group_id; kwargs...): convenience form. group_id is an integer vector with labels 1:G (each present at least once); the constructor counts members per group and calls factory(s_g; kwargs...) for each group size s_g. Works for components whose primary constructor takes the size as a positional argument and prior options as keyword arguments — IID, AR1, RW1, RW2, Seasonal. Components needing more than a size (e.g. Besag, Generic0) must use the Vector form.

Mechanics

• length(g) = sum(length, components).
• precision_matrix(g, θ) = blockdiag(Q₁, Q₂, …, Q_G) where Q_g = precision_matrix(components[g], θ).
• prior_mean(g, θ) is the concatenation of per-group prior means.
• GMRFs.constraints(g) block-stacks each group's constraint at its column offset within the stacked latent; groups with NoConstraint contribute no rows.
• log_normalizing_constant(g, θ) = Σ_g log_normalizing_constant( components[g], θ).

Example

# Replicated AR1 with unequal panel lengths (subjects with different
# numbers of visits), all sharing (τ, ρ).
panels = [AR1(10), AR1(7), AR1(12)]
g      = Group(panels)
m      = LatentGaussianModel(GaussianLikelihood(), (g,), A)

IID(n; hyperprior = PCPrecision(),
    τ_init = 0.0,
    fix_τ::Bool = false)

IID Gaussian random effect of length n with precision τ. One hyperparameter on the internal scale θ = log(τ); the user-facing parameter is τ.

Composes GMRFs.IIDGMRF — the precision is τ I_n. No linear constraint.

Arguments

• n: number of independent slots.
• hyperprior: scalar prior on log τ. Default PCPrecision().
• τ_init: initial value of log τ on the internal scale. Used as the optimisation start, and as the held-fixed value when fix_τ = true.
• fix_τ: when true, drop the precision from the θ vector and hold it at τ = exp(τ_init). Useful for the multinomial-via-Poisson reformulation (ADR-024), where the per-row α-intercepts carry as a fixed-precision IID nuisance block (τ_init = -10, fix_τ = true, matching R-INLA's prec = list(initial = -10, fixed = TRUE)).

IIDND_Sep{N, PP, PC} <: AbstractIIDND

Separable parameterisation of an N-variate IID component. PP is the type of precpriors (an NTuple{N, AbstractHyperPrior}); PC is the type of corrpriors (an NTuple{N(N-1)/2, AbstractHyperPrior}).

Hyperparameter vector layout (length N + N(N-1)/2):

θ = (log τ_1, …, log τ_N,  atanh z_{2,1}, atanh z_{3,1}, atanh z_{3,2}, …,  atanh z_{N,N-1})

with the strictly-lower-triangular entries enumerated in row-major order ((i, j) with j < i). The z_{i,j} are the LKJ canonical partial correlations of the correlation matrix R (Lewandowski-Kurowicka-Joe 2009): for the first column they coincide with the direct correlations (z_{i,1} = ρ_{i,1}); for entries with i, j > 1 they are partial correlations of variables i and j controlling for variables 1, …, j − 1. The correlation matrix is reconstructed from the z's via stick-breaking, which guarantees R is positive-definite for any (z_{i,j}) ∈ (-1, 1)^{N(N-1)/2}.

PR-1a ships the N = 2 case (where z_{2,1} = ρ_{1,2} directly); PR-1b lifts the N = 3 case via the LKJ stick-breaking step.

INLA(; int_strategy = :auto, laplace = Laplace(),
      latent_strategy = Gaussian(),
      skewness_correction = false,
      θ0 = nothing, optim_options = NamedTuple())

Integrated Nested Laplace Approximation for a LatentGaussianModel.

The algorithm (Rue, Martino & Chopin 2009) proceeds in three stages:

1. Find the mode θ̂ = argmax_θ log π(θ | y) using the Laplace-approximated log marginal and Optimization.jl.
2. Compute the Hessian H of the negative log-posterior at θ̂ by finite differences. Let Σ = H⁻¹.
3. Integrate x | y, θ | y, and x_i | y against π(θ | y) using a design {θ_k, w_k} that is accurate for Gaussian integrands; reweight each point by π̂(θ_k | y) / q(θ_k) where q ~ N(θ̂, Σ).

int_strategy

• :auto (default) — Grid() for dim(θ) ≤ 2, CCD() otherwise.
• :grid / Grid() — tensor-product midpoint grid.
• :ccd / CCD() — central composite design per Rue-Martino-Chopin §6.5.
• :gauss_hermite / GaussHermite() — tensor-product Gauss-Hermite.
• Any subtype of AbstractIntegrationScheme can be passed directly.

latent_strategy

Controls the per-θ approximation of the latent posterior summary (x_mean, x_var) — this is R-INLA's control.inla$strategy, distinct from the θ-integration scheme above.

Accepts either an AbstractMarginalStrategy instance or the legacy symbols (:gaussian, :simplified_laplace). See [ADR-026].

• Gaussian() / :gaussian (default) — use the Newton mode x̂(θ) and the constraint-corrected Laplace marginal variance. Bit-for-bit unchanged from prior releases.
• SimplifiedLaplace() / :simplified_laplace — apply the Rue-Martino mean-shift correction Δx(θ) = ½ H⁻¹ Aᵀ (h³ ⊙ σ²_η) per integration point before accumulating x_mean / x_var. Reduces to Gaussian exactly when the likelihood third derivative ∇³_η log p is zero.

The mean-shift correction is independent of the density-shape skew correction in posterior_marginal_x(strategy = :simplified_laplace): that one expands the per-coordinate density p(x_i | y) around the unshifted Newton mode (Rue-Martino-Chopin 2009 §4.2). They can be toggled independently — see ADR-016.

R-INLA's full simplified.laplace also includes a variance correction; that piece is deferred to v0.3 (see ADR-006 amendment).

skewness_correction

Opt-in asymmetric per-axis stretches for the Grid quadrature design (Rue-Martino-Chopin 2009 §6.5). When true, after θ̂ and Σ are in hand the pipeline probes log π̂(θ|y) at θ̂ ± Σ^{1/2} e_k along each eigen-axis k and rebuilds the Grid scheme with per-axis stdev_corr_pos / stdev_corr_neg factors that align the design points with the actual log-posterior shape — useful when the hyperparameter posterior is heavy-tailed on one side (e.g. log-precision posteriors with sharp upper-curvature walls and slow decay below).

Quadrature weights stay as standard-normal weights — the asymmetric placement only relocates points; the IS reweight in step (3) absorbs the proposal mismatch.

The flag has no effect for non-Grid schemes (CCD / GaussHermite); for those a future PR can wire equivalent corrections. Default false matches R-INLA's control.inla$skew.corr.positive = 1.0, skew.corr.negative = 1.0 (off by default).

INLALogDensity(model::LatentGaussianModel, y; laplace = Laplace())

LogDensityProblems.jl conformance for the INLA hyperparameter posterior log π(θ | y) = log p(y | θ) + log π(θ).

For a fixed observation vector y, the bundle (model, y) defines a density over the internal-scale hyperparameter vector θ. Evaluating logdensity at θ runs a Laplace approximation of p(x | θ, y) and returns the marginalised log posterior up to an additive constant. The gradient is computed by central finite differences, which is adequate for NUTS / Pathfinder given that each logdensity call is already O(Laplace) in cost and dim(θ) is typically small.

LogDensityProblems interface

LogDensityProblems.dimension(ld)                 # n_hyperparameters(model)
LogDensityProblems.capabilities(INLALogDensity)  # LogDensityOrder{1}()
LogDensityProblems.logdensity(ld, θ)             # log p(y|θ) + log π(θ)
LogDensityProblems.logdensity_and_gradient(ld, θ)

Example

using LogDensityProblems
ld = INLALogDensity(model, y)
LogDensityProblems.dimension(ld)         # number of hyperparameters
ℓθ = LogDensityProblems.logdensity(ld, θ)

Downstream samplers (AdvancedHMC, Pathfinder, Turing via LGMTuring.jl) consume this interface unchanged.

INLAResult <: AbstractInferenceResult

• θ̂::Vector{Float64} — posterior mode of θ | y on the internal scale.
• Σθ::Matrix{Float64} — Gaussian covariance approximation at θ̂.
• θ_points::Vector{Vector{Float64}} — integration design points.
• θ_weights::Vector{Float64} — normalised posterior weights on the points.
• laplaces::Vector{LaplaceResult} — Laplace fit at each design point.
• x_mean::Vector{Float64} — posterior mean of x | y, integrated over θ.
• x_var::Vector{Float64} — posterior variance of x | y, integrated over θ.
• θ_mean::Vector{Float64} — posterior mean of θ | y (internal scale).
• log_marginal::Float64 — log p(y) estimate (marginal likelihood of the model).
• optim_result — raw Optimization.jl solution for θ̂.

IdentityMapping(n)

η = x; algebraically A = I_n. Used for areal models where the observation index equals the latent index. The dimension n is carried so shape assertions at construction are cheap.

Intercept(; prec = 1.0e-3, improper = true)

Scalar intercept. With improper = true (default, matching R-INLA's prec.intercept = 0), the prior is treated as flat (improper Lebesgue); prec is then only a numerical regulariser added to the joint precision to keep it invertible, and is dropped from the log normalising constant of the joint Gaussian. With improper = false the prior is the proper Normal α ~ N(0, prec⁻¹).

The latent field of this component has length 1.

KroneckerComponent(space::AbstractLatentComponent,
                   time::AbstractLatentComponent)

Separable space-time latent component composing two child AbstractLatentComponents into a Kronecker-structured GMRF with precision

Q = Q_space ⊗ Q_time.

R-INLA's f(field, model = spatial, group = time, control.group = list(model = "ar1")) builds the same separable form for SPDE-Matérn × AR(1); the Cameletti et al. (2013) PM₁₀ case study is the canonical oracle. Pairs with the KroneckerMapping projector so the flattening conventions agree end-to-end.

Hyperparameter layout

θ = [θ_space; θ_time] concatenated. The composer adds no hyperparameters of its own:

nhyperparameters(c) == nhyperparameters(c.space) + nhyperparameters(c.time)

log_hyperprior(c, θ) sums the two children's log-priors.

Latent flattening — matches KroneckerMapping

The latent vector x is vec(X) where X has shape (length(time) × length(space)). Equivalently: x[(s - 1) · n_t + t] is the value at space slot s, time slot t. This matches the KroneckerMapping(A_space, A_time) storage convention ((A_space ⊗ A_time) vec(X) = vec(A_time · X · A_space')), so a KroneckerComponent(space, time) paired with KroneckerMapping(A_space, A_time) projects correctly without any transposition.

Constraints

• Both children proper (NoConstraint): composer is proper.
• One child constrained, the other proper: the composer's constraint is the child's constraint lifted via Kronecker with the identity on the unconstrained axis.
• Both children constrained: rejected at construction. The joint-constraint case (rank, conditioning-by-kriging interaction with the Kronecker factor) is deferred — see ADR-029.

Example

using LatentGaussianModels: KroneckerComponent, AR1, IID

# AR1 in space × AR1 in time, sharing nothing.
space = AR1(50)
time  = AR1(20)
spt   = KroneckerComponent(space, time)
length(spt) == 50 * 20  # 1000

# Compose with a KroneckerMapping projector for observation-side
# Kronecker structure (e.g. SPDE2 ⊗ AR1 for Cameletti).

KroneckerMapping(A_space, A_time)

Separable space-time projector A = A_space ⊗ A_time. Forward and adjoint multiplies use the matrix identity (A_space ⊗ A_time) vec(X) = vec(A_time · X · A_space') so the dense Kronecker product is never materialized — critical when both factors are large (e.g. SPDE mesh × time series).

Storage layout

Treat x and η as column-major flattenings:

• reshape(x, ncols(A_time), ncols(A_space)) — column index ranges over space, row index over time.
• reshape(η, nrows(A_time), nrows(A_space)) — same convention.

Laplace(; maxiter = 50, tol = 1.0e-8)

Laplace approximation to p(x | θ, y) at a given θ. Finds the mode x̂ by Newton iteration and returns a Gaussian with precision H = Q(θ) - A' diag(∇²_η log p) A evaluated at x̂. The signs follow the convention that ∇²_η log p is negative semi-definite so that H is positive definite whenever Q is PD on the non-null subspace.

Hard linear constraints declared by components via GMRFs.constraints(c) are stacked into a model-level constraint C x = e and enforced via the Rue & Held (2005) §2.3 kriging correction — the rank-deficient prior is regularised along null(Q) with V V^T (V = C^T (C C^T)^{-1/2}), Newton is run on the PD regularised posterior, and the final x̂ is projected back to {x : C x = e}.

LaplaceResult

The output of a Laplace fit at fixed θ.

• mode::Vector{Float64} — x̂, on the constrained surface C x̂ = e.
• precision::SparseMatrixCSC{Float64,Int} — posterior precision at x̂ on the non-null subspace, regularised with the null-space bump V V^T so that it is PD and has the same quadratic form as the original H = Q + A'DA on null(C).
• factor::FactorCache — cached sparse Cholesky of precision.
• θ::Vector{Float64} — internal-scale hyperparameters used.
• log_joint::Float64 — log p(x̂, y | θ).
• log_marginal::Float64 — Laplace log marginal log p(y | θ).
• iterations::Int, converged::Bool.
• constraint::Union{Nothing, NamedTuple} — if constrained, the triple (C, U, W_fact) where U = H_reg^{-1} C^T and W_fact is the Cholesky of C U. Used downstream to subtract the kriging correction from conditional marginal variances.

LatentGaussianModel(likelihoods, components, mapping)
LatentGaussianModel(likelihood,  components, mapping_or_A)
LatentGaussianModel(likelihood,  component,  mapping_or_A)

Bayesian latent Gaussian model:

y | η, θ_ℓ ∼ p(y | η, θ_ℓ),     η = mapping * x,
x | θ     ∼ N(0, Q(θ)⁻¹),       Q(θ) = blockdiag(Q_1(θ_1), ...),
θ         ∼ π(θ).

Arguments

• likelihoods::Tuple{Vararg{AbstractLikelihood}} — one observation model per block. Each block applies to a contiguous slice of rows of η, with the slice ranges supplied by the mapping. For a single-likelihood model the convenience signature LatentGaussianModel(ℓ, components, mapping_or_A) accepts a scalar ℓ::AbstractLikelihood and promotes it internally to (ℓ,).
• components::Tuple{Vararg{AbstractLatentComponent}} — latent components, concatenated into a single x of length sum(length, components).
• mapping::AbstractObservationMapping — projector from the stacked latent vector x to the linear predictor η. The convenience signature LatentGaussianModel(ℓ, components, A::AbstractMatrix) wraps the matrix in LinearProjector automatically (the v0.1 default).

For multi-likelihood models the mapping must be a StackedMapping whose number of blocks equals length(likelihoods); the k-th block of mapping defines the row range owned by likelihoods[k].

Internally we store the latent slices (index ranges of each component in x), the per-component hyperparameter slices in the full θ, and the per-likelihood hyperparameter slices in the full θ. The full θ layout is [θ_ℓ_1, ..., θ_ℓ_K, θ_c_1, ..., θ_c_M] — likelihood hyperparameters first (in block order), component hyperparameters after.

Leroux(graph; hyperprior_tau = PCPrecision(),
       hyperprior_rho = LogitBeta(1.0, 1.0))

Leroux (1999) CAR component. The precision is

Q(τ, ρ) = τ · ((1 - ρ) · I_n + ρ · R)

where R = D - W is the combinatorial graph Laplacian. The mixing parameter ρ ∈ (0, 1) interpolates between pure IID (ρ = 0) and a pure ICAR (ρ = 1, rank-deficient). For ρ strictly inside (0, 1) the precision is positive definite and no sum-to-zero constraint is needed, distinguishing Leroux from Besag / BYM.

Two hyperparameters on the internal scale: θ = [log(τ), logit(ρ)]. Default priors: PC precision on τ (matching PCPrecision()), Beta(1, 1) on ρ (uniform), the latter expressed on the logit scale via LogitBeta.

R-INLA equivalent: model = "besagproper2" with the convex Leroux combination, or model = "leroux" in the proprietary inla build.

LinearProjector(A::AbstractMatrix)

Wraps an explicit row-to-latent matrix A. The v0.1 default — every existing LatentGaussianModel(...) call with a matrix A becomes LinearProjector(A) internally via the compatibility constructor.

Forward / adjoint multiplies dispatch directly onto LinearAlgebra.mul! on the wrapped matrix; no additional allocation.

LogNormalPrecision(μ, σ)

LogNormal prior on τ with parameters μ, σ on the log scale. Internal scale is θ = log(τ), so the prior is simply N(μ, σ²) on θ.

LogitBeta(a, b)

Beta(a, b) prior on a parameter ρ ∈ (0, 1) expressed on the internal logit scale θ = logit(ρ). Includes the Jacobian |dρ/dθ| = ρ(1 - ρ), so on the internal scale

log π(θ) = a · log(ρ) + b · log(1 - ρ) - log B(a, b)

where B(a, b) = Γ(a) Γ(b) / Γ(a + b). With a = b = 1 this is the uniform-on-ρ prior (i.e. a standard logistic on θ), matching R-INLA's default logit-Beta(1, 1) prior on Leroux's mixing parameter.

Numerically stable via the -softplus identities log σ(θ) = -log(1 + e^{-θ}), log(1 - σ(θ)) = -log(1 + e^{θ}).

LognormalSurvLikelihood(; link = IdentityLink(), censoring = nothing,
                          time_hi = nothing,
                          hyperprior = PCPrecision(1.0, 0.01))

Lognormal survival likelihood, matching R-INLA's family = "lognormalsurv". Under the canonical IdentityLink,

log T_i  ~ N(η_i, σ²),    σ² = 1/τ
f(t)     = (1 / (t σ √(2π))) exp(-(log t − η)² / (2σ²))
S(t)     = Φ̄((log t − η) / σ) = Φ((η − log t) / σ)
F(t)     = Φ((log t − η) / σ)

η_i is the mean of log T_i — equivalently, exp(η_i) is the median survival time. The single hyperparameter is the precision τ, carried internally as θ_ℓ = [log τ]. The default hyperprior is the PC prior on standard deviation σ = τ^(-1/2) with P(σ > 1) = 0.01, matching R-INLA-modern Gaussian conventions.

censoring is an optional AbstractVector{Censoring} of length length(y); when nothing (default), every observation is uncensored. For INTERVAL rows, time_hi[i] is the upper bound (with y[i] the lower bound).

Closed-form derivatives are provided for all four Censoring modes via the inverse-Mills ratio h(u) = φ(u)/Φ(u) (and its first/second derivatives). Stable evaluation in both tails uses Distributions.logpdf / logcdf / logccdf.

Example

ℓ = LognormalSurvLikelihood(censoring = [NONE, RIGHT, NONE])

MEB(values; scale = ones(length(values)),
    τ_u_prior = GammaPrecision(1.0, 1.0e-4),
    τ_u_init = log(1000.0))

R-INLA's model = "meb" — Berkson measurement-error component (ADR-023). The latent block x is one slot per supplied values entry with prior

x ~ N(values, (τ_u · diag(scale))⁻¹).

The Berkson tie x_i = w_i + u_i, u_i ~ N(0, (τ_u s_i)⁻¹) has been marginalised so that the latent prior carries the observed (proxy) w as a θ-constant prior mean.

Arguments

• values: length-n vector of per-slot prior means (the observed proxy w, deduplicated by the caller if appropriate).
• scale: length-n per-slot diagonal scaling. Matches R-INLA's scale = ... argument. Default is all-ones.
• τ_u_prior: scalar prior on θ[1] = log τ_u. Default GammaPrecision(1.0, 1.0e-4) (matches R-INLA's loggamma(1, 1e-4)).
• τ_u_init: initial value for θ[1]. Default log(1000.0), matching R-INLA's initial = log(1000).

Hyperparameters

The β scaling that multiplies x before it lands in η lives on the receiving likelihood as a Copy (per ADR-021/ADR-023), not on the component. The R-INLA β default is GaussianPrior(1.0, 0.001) on the user scale; users attach it as

c = MEB(w)
m = LatentGaussianModel(...)        # `c` placed in the component tuple
range_c = component_range(m, c_idx) # 1-indexed position in m.components
β_copy = Copy(range_c; β_prior = GaussianPrior(1.0, 0.001), β_init = 1.0)
target = CopyTargetLikelihood(receiving_likelihood, β_copy)

Note on gmrf(c, θ)

The gmrf(c::MEB, θ) accessor returns a GMRFs.Generic0GMRF carrying only the precision; the non-zero prior mean is exposed separately via prior_mean(c, θ). Inference inside the LGM package reads prior_mean(c, θ) through joint_prior_mean and is therefore correct regardless of how gmrf represents the mean. Direct rand(rng, gmrf(c, θ)) calls draw from the centred prior — add prior_mean(c, θ) if you want the shifted draw.

MEC(values; scale = ones(length(values)),
    τ_u_prior = GammaPrecision(1.0, 1.0e-4),
    μ_x_prior = GaussianPrior(0.0, 1.0e-4),
    τ_x_prior = GammaPrecision(1.0, 1.0e4),
    τ_u_init  = log(10000.0),
    μ_x_init  = 0.0,
    τ_x_init  = -log(10000.0),
    fix_τ_u::Bool = true,
    fix_μ_x::Bool = true,
    fix_τ_x::Bool = true)

R-INLA's model = "mec" — Classical measurement-error component (ADR-023). The latent block x is one slot per supplied values entry with prior x ~ N(μ_x · 1, (τ_x I)⁻¹) and an observed proxy w | x ~ N(x, (τ_u D)⁻¹) with D = diag(scale). Gaussian conjugacy absorbs the Berkson tie into the prior:

Q̂(θ) = τ_x I + τ_u D
μ̂(θ) = Q̂⁻¹ · (τ_x μ_x 1 + τ_u D · values)

so the LGM-level latent prior is x ~ N(μ̂(θ), Q̂(θ)⁻¹). Unlike MEB's θ-constant mean, MEC's prior_mean(c, θ) depends on θ through (τ_u, μ_x, τ_x).

Arguments

• values: length-n vector of observed proxy w (deduplicated by the caller if appropriate).
• scale: length-n per-slot diagonal scaling. Default all-ones.
• τ_u_prior, μ_x_prior, τ_x_prior: scalar priors. Defaults match R-INLA's mec.tex.
• τ_u_init, μ_x_init, τ_x_init: initial values. Defaults match R-INLA.
• fix_τ_u, fix_μ_x, fix_τ_x: toggle whether each slot is estimated. R-INLA's default is fix_* = true for all three — the model degrades to plain regression unless the user opts in.

Hyperparameters

Per ADR-023, the component carries up to three internal slots in canonical order (log τ_u, μ_x, log τ_x). Fixed slots are excluded from the θ vector and held at their *_init values. The β scaling that multiplies x before it lands in η lives on the receiving likelihood as a Copy (per ADR-021/ADR-023), not on the component. The R-INLA β default is GaussianPrior(1.0, 0.001) on the user scale; users attach it as

c = MEC(w)
m = LatentGaussianModel(...)        # `c` placed in the component tuple
range_c = component_range(m, c_idx) # 1-indexed position in m.components
β_copy = Copy(range_c; β_prior = GaussianPrior(1.0, 0.001), β_init = 1.0)
target = CopyTargetLikelihood(receiving_likelihood, β_copy)

Note on gmrf(c, θ)

gmrf(c::MEC, θ) returns a GMRFs.Generic0GMRF carrying only the precision Q̂(θ); the non-zero, θ-dependent prior mean is exposed separately via prior_mean(c, θ). LGM inference reads prior_mean(c, θ) through joint_prior_mean and is correct regardless of how gmrf represents the mean.

NegativeBinomialLikelihood(; link = LogLink(), E = nothing,
                             hyperprior = GammaPrecision(1.0, 0.1))

y_i | η_i, n ∼ NegBinomial(μ_i, n) with μ_i = E_i · g⁻¹(η_i) and size (overdispersion) parameter n > 0. The mean is μ and the variance is μ + μ² / n; as n → ∞ this recovers Poisson.

E is an optional offset / exposure vector (e.g. expected counts in disease mapping); defaults to 1. The size parameter is carried as a single hyperparameter on the internal scale θ = log(n). The default hyperprior matches R-INLA's family = "nbinomial" default (loggamma(1, 0.1) on log(n)).

Currently only the canonical LogLink is supported.

PCAlphaW(λ = 5.0)

Penalised-complexity prior on the Weibull shape parameter α > 0 (Sørbye-Rue 2017; Simpson et al. 2017). The reference model is the exponential α = 1, and the prior is exponential on the Kullback-Leibler distance

d(α) = √(2 · KL(Weibull(1, α) ‖ Exp(1)))
     = √(2 · [Γ(1 + 1/α) + log α - γ + γ/α - 1])

with rate λ. The corresponding density on α > 0 is

π(α) = ½ λ exp(-λ d(α)) |dd/dα|,

with the factor of ½ because α = 1 is interior to the parameter space and the unsigned distance d corresponds to two values of α either side of 1. On the internal scale θ = log α (where the Weibull likelihood carries α), the Jacobian |dα/dθ| = α cancels with the 1/α in |dd/dlogα| = α |dd/dα|, so

log π_θ(θ) = log(½ λ) - λ d(α) + log|dd/dlogα|.

Mirrors R-INLA's inla.pc.dalphaw exactly (cross-checked at α ∈ {0.5, 1, 2, 5} to ≈1e-12 nat).

The default λ = 5 matches R-INLA's pc.alphaw default — a moderately informative prior on shrinkage toward the exponential.

Example

hyper = PCAlphaW(5.0)
log_prior_density(hyper, 0.0)   # log π_θ(0) at α = 1 (the reference)

PCBYM2Phi{T}(U, α, γ)

PC prior on the BYM2 mixing parameter φ ∈ [0, 1], specified by P(φ < U) = α. The base model is φ = 0 (pure IID random effects); the alternative is φ > 0 (increasing weight on the scaled spatial component). See Riebler, Sørbye, Simpson & Rue (2016).

γ is the vector of non-zero eigenvalues of the scaled Besag structure matrix R_scaled = c · L, where L is the combinatorial graph Laplacian and c is the Sørbye-Rue (2014) scaling constant. The length of γ equals n - r, where n is the graph size and r is the number of connected components.

The prior is obtained from the exponential-on-distance PC construction, with

d(φ) = sqrt(2 · KLD(φ))
KLD(φ) = (1/2) Σ_k [ φ(1/γ_k - 1) - log(1 + φ(1/γ_k - 1)) ]

on the non-null subspace. Since d is bounded on φ ∈ [0, 1] (i.e. d(1) < ∞), the exponential density on d is truncated to [0, d(1)] and renormalised by Z = 1 - exp(-λ · d(1)). The rate λ > 0 is solved numerically so that

P(φ < U) = (1 - exp(-λ · d(U))) / (1 - exp(-λ · d(1))) = α.

Internal scale is θ = logit(φ); the Jacobian φ(1-φ) is included in log_prior_density.

Defaults U = 0.5, α = DEFAULT_BYM2_PHI_ALPHA match Riebler et al. (2016). See plans/defaults-parity.md for the ⚠ unverified-default note.

PCCor0(U = 0.5, α = 0.5)

Penalised-complexity prior on a correlation ρ ∈ (-1, 1) with reference model at ρ = 0 (independence). Mirrors R-INLA's pc.cor0 — used for the bivariate-IID correlation slot in IIDND_Sep{2} (and recursively for each pairwise correlation in IIDND_Sep{N} for N ≥ 3).

The Kullback-Leibler distance from a bivariate normal with correlation ρ and unit marginals to its ρ = 0 counterpart is

KLD(ρ) = -½ log(1 - ρ²),

so the natural PC distance is

d(ρ) = √(2 · KLD(ρ)) = √(-log(1 - ρ²)).

The PC prior is exponential on d(ρ) with rate λ; because d is symmetric under ρ ↔ -ρ, the density on ρ ∈ (-1, 1) carries a factor of ½:

π_ρ(ρ) = ½ λ exp(-λ d(ρ)) · |dd/dρ|,
|dd/dρ| = |ρ| / ((1 - ρ²) d(ρ))    (for ρ ≠ 0).

The user-facing parameters are (U, α) with P(|ρ| > U) = α, i.e. α = exp(-λ · d(U)), hence λ = -log(α) / √(-log(1 - U²)).

Internal scale is θ = atanh(ρ), so ρ = tanh(θ) and |dρ/dθ| = 1 - ρ². The Jacobian cancels the explicit 1 - ρ² in |dd/dρ|, leaving

log π_θ(θ) = log(λ) - log(2) - λ d(ρ)
             - ½ log( -log(1 - ρ²) / ρ² )

where the bracketed ratio is well-defined at ρ = 0 (limit equals 1) — handled by a Taylor short-circuit when ρ² is below a small threshold.

Defaults U = 0.5, α = 0.5 match R-INLA's pc.cor0(0.5, 0.5) default.

PCCor1(U = 0.7, α = 0.7)

Penalised-complexity prior on a correlation ρ ∈ (-1, 1) with reference model at ρ = 1 (perfect positive correlation). Mirrors R-INLA's pc.cor1 — the textbook PC prior on the lag-1 correlation of an AR(1) process, where ρ = 1 is the random-walk limit.

Density on the user scale (R-INLA inla.pc.dcor1):

π_ρ(ρ) = λ exp(-λ √(1-ρ)) / (1 - exp(-√2 λ)) · 1/(2√(1-ρ)),  ρ ∈ (-1, 1).

The "complexity distance" is μ(ρ) = √(1-ρ), monotonically decreasing from √2 (at ρ = -1) to 0 (at ρ = 1); the density is exponential in μ truncated to μ ∈ (0, √2). Calibration is via (U, α) with P(ρ > U) = α, i.e. λ is the unique positive root of

(1 - exp(-λ √(1-U))) / (1 - exp(-√2 λ)) = α,    α > √((1-U)/2).

The lower bound on α enforces λ > 0 — at α = √((1-U)/2) the rate collapses to zero and the prior becomes uniform on μ.

Internal scale is θ = atanh(ρ) (matching AR1's lag-1 correlation parameterisation), so ρ = tanh(θ) and |dρ/dθ| = 1 - ρ². Combining the density and the Jacobian:

log π_θ(θ) = log(λ) - λ √(1-ρ) - log(1 - exp(-√2 λ))
             - log 2 + ½ log(1-ρ) + log(1+ρ).

Computations of log(1-ρ), log(1+ρ), and √(1-ρ) route through softplus of ±2θ so that 1 - tanh(θ) underflowing to 0 for θ ≳ 19 does not break the prior at saturation.

Defaults

U = 0.7, α = 0.7 matches R-INLA's textbook example for pc.cor1 (also the AR(1) ρ-prior default suggested in INLA tutorials when the user opts out of the built-in Normal on the logit scale).

Notes

• Not a default for IID2D / IID3D / 2diid: the bivariate-IID correlation slot uses PCCor0 (reference at ρ = 0) — see ADR-022.
• R-INLA's actual f(., model="ar1") default is a Normal prior on 2 atanh(ρ) with precision 0.15, not pc.cor1. Pass ρprior = PCCor1(...) to AR1 to opt into this prior.

PCGevtail(λ = 7.0, interval = (0.0, 1.0); xi_scale = 0.1)

Penalised-complexity prior on the GEV tail (shape) parameter ξ ∈ interval (Opitz et al. 2018; R-INLA's pc.gevtail). The reference value is ξ = 0 (Gumbel), and the prior takes the linearised PC distance d(ξ) = ξ exact for ξ ≥ 0 close to the reference. The user-scale density is the truncated exponential

π_ξ(ξ) = λ exp(-λ ξ) / Z,   ξ ∈ [low, high],
Z      = exp(-λ low) - exp(-λ high),

matching inla.pc.dgevtail(xi, lambda, interval) in R-INLA.

The internal scale matches GEVLikelihood: θ = ξ / xi_scale, with |dξ/dθ| = xi_scale. On the internal scale the log-density is

log π_θ(θ) = log λ - λ · xi_scale · θ - log Z + log(xi_scale),

returning -Inf outside [low / xi_scale, high / xi_scale].

Defaults λ = 7, interval = (0, 1), xi_scale = 0.1 match R-INLA's pc.gevtail defaults composed with GEVLikelihood's default xi_scale = 0.1 (R-INLA's gev.scale.xi).

Example

hyper = PCGevtail(7.0, (0.0, 1.0); xi_scale = 0.1)
log_prior_density(hyper, 0.5)   # at θ = 0.5 ⇒ ξ = 0.05

PCPrecision(U, α)

The PC prior on precision (Simpson et al. 2017), specified via P(σ > U) = α where σ = 1/√τ. User-facing parameter is τ > 0; internal scale is θ = log(τ).

The marginal prior on σ is Exponential(λ) with λ = -log(α) / U. Under τ = σ⁻², the density on the internal scale θ = log(τ) is

π(θ) = π_σ(σ(θ)) · |dσ/dθ|
     = λ exp(-λ σ) · (σ / 2)       (using σ = exp(-θ/2))

where the σ/2 is |dσ/dθ| = ½ · exp(-θ/2) = σ/2.

Defaults U = 1, α = 0.01 match R-INLA's pc.prec.

POMLikelihood(n_classes; link = LogitLink(),
              dirichlet_concentration = 3.0)

Proportional-odds ordinal regression in R-INLA's family = "pom" parameterisation. With ordered response y_i ∈ {1, …, K} and linear predictor η_i (without an intercept — the cut points absorb it),

P(y_i ≤ k | η_i, α) = F(α_k − η_i),    k = 1, …, K − 1
P(y_i = k | η_i, α) = F(α_k − η_i) − F(α_{k−1} − η_i)

with α_0 = −∞, α_K = +∞, and F the standard logistic CDF. The K − 1 ordered cut points α_1 < α_2 < ⋯ < α_{K−1} carry as likelihood hyperparameters on the internal scale:

θ[1]   = α_1                    (real-valued, no transform)
θ[k]   = log(α_k − α_{k−1}),    k = 2, …, K − 1

so the increments α_k − α_{k−1} = exp(θ[k]) are strictly positive, guaranteeing the ordering for any θ ∈ ℝ^{K−1}. The struct stores n_classes = K; nhyperparameters(ℓ) = K − 1.

Prior — Dirichlet on the implied class probabilities

R-INLA's pom family hard-wires a single Dirichlet prior on the class probabilities implied by the cut points at η = 0,

π_k(α) = F(α_k) − F(α_{k−1}),    k = 1, …, K

with α_0 = −∞, α_K = +∞, so (π_1, …, π_K) ∈ Δ^{K−1}. The prior density is

p(π_1, …, π_K) = Γ(K γ) / Γ(γ)^K · ∏_k π_k^{γ−1}

for a single concentration γ > 0 (dirichlet_concentration, default 3 — matching R-INLA). Pushing this back to θ via the chain θ → α → π adds the Jacobian

log |det(dπ/dθ)| = Σ_{k=1}^{K−1} log f(α_k)  +  Σ_{k=2}^{K−1} θ_k

where f(x) = F(x) (1 − F(x)) is the logistic pdf.

The R-INLA pom doc records that R-INLA's internal Dirichlet prior is "correct only up to a multiplicative constant due to a missing correction in the log-Jacobian for the sum-to-zero constraint." This implementation includes the full Jacobian (i.e. is mathematically exact); the consequence is that the marginal log-likelihood reported by Julia and R-INLA on the same fit will differ by a fixed θ- independent additive constant, while every posterior moment of α, β, and θ matches.

Currently only LogitLink is supported (the cumulative link). The linear predictor η_i is fed through identity — there is no separate "link on η". Probit is the only alternative R-INLA supports for family = "pom" and could be added later via dispatch on link; the closed forms below are specific to logistic-CDF arithmetic.

Closed-form derivatives (logit case): with g_k = F(α_k − η_i) ∈ (0,1) and f_k = g_k (1 − g_k), f'_k = f_k (1 − 2 g_k) for k ∈ {1, …, K − 1}, and g_0 = f_0 = f'_0 = 0, g_K = 1, f_K = f'_K = 0,

log p_i              = log(g_{y_i} − g_{y_i − 1})
∂ log p_i / ∂η_i     = (f_{y_i − 1} − f_{y_i}) / p_i
∂² log p_i / ∂η_i²   = (f'_{y_i} − f'_{y_i − 1}) / p_i
                       − (f_{y_i − 1} − f_{y_i})² / p_i²

For the boundary classes the formulae collapse to:

y_i = 1:   ∂ log p / ∂η = −(1 − g_1),
           ∂² log p / ∂η² = −g_1 (1 − g_1)
y_i = K:   ∂ log p / ∂η = g_{K−1},
           ∂² log p / ∂η² = −g_{K−1} (1 − g_{K−1})

(∂³ inherits the abstract finite-difference fallback.)

The cumulative-logit likelihood is globally log-concave in η, so the inner Newton step is well-behaved without damping.

PoissonLikelihood(; link = LogLink(), E = nothing)

y_i | η_i ∼ Poisson(E_i · g⁻¹(η_i)). E_i is an optional offset / exposure vector (e.g. expected counts in disease mapping); defaults to 1. With the canonical LogLink, the mean is E · exp(η).

No likelihood hyperparameters. θ is ignored.

RW1(n; hyperprior = PCPrecision(), cyclic = false)

First-order random walk of length n. One hyperparameter on internal scale θ = log(τ). The component is intrinsic (rank deficient by 1) and carries a sum-to-zero constraint so that the intercept is identifiable.

NOTE: Sørbye-Rue (2014) scale.model = TRUE behaviour is not yet applied here — the unit-τ geometric-mean marginal variance is not forced to 1. Tracked in plans/defaults-parity.md.

RW2(n; hyperprior = PCPrecision(), cyclic = false)

Second-order random walk. Open version has rank deficiency 2 (linear trend); cyclic has rank 1. Defaults match RW1; the one-hyperparameter contract is identical.

Replicate(component::AbstractLatentComponent, n_replicates::Integer)

R-INLA's f(idx, model = M, replicate = id) — stack n_replicates independent copies of component, all sharing a single hyperparameter block. The latent vector is the concatenation [x⁽¹⁾; x⁽²⁾; …; x⁽ⁿ⁾] of length n_replicates · length(component), and the joint prior precision is blockdiag(Q, Q, …, Q) where Q = precision_matrix(component, θ).

Hyperparameter sharing — the load-bearing piece

Replicate exposes the same nhyperparameters, initial_hyperparameters, and log_hyperprior as the inner component; the same θ slice flows verbatim into each replicate's precision_matrix(component, θ) call. Without this, the replicate prior would silently degenerate to n_replicates independent priors — the test nhyperparameters(Replicate(c, n)) == nhyperparameters(c) is the public guarantee that hyperparameters are shared.

Prior mean, constraints, log-NC

• prior_mean(r, θ) repeats the inner mean n_replicates times.
• GMRFs.constraints(r) block-stacks the inner constraint (A_c, e_c): if A_c is (k × n), the replicate carries an (n_replicates · k) × (n_replicates · n) constraint with A_c on each diagonal block and the same e_c repeated.
• log_normalizing_constant(r, θ) = n_replicates · log_normalizing_constant(component, θ) — the blockdiag log-det factorises exactly across replicates.

Example

ar1   = AR1(20)                       # 20 time points per panel member
panel = Replicate(ar1, 50)            # 50 panel members share (τ, ρ)
m     = LatentGaussianModel(GaussianLikelihood(), (panel,), A)

Seasonal(n; period, hyperprior = PCPrecision())

Intrinsic seasonal-variation component of length n and period s = period. Thin LGM wrapper around GMRFs.SeasonalGMRF: the precision is Q = τ · B' B penalising every s-consecutive sum. The null space has dimension s - 1; the default constraint is a single 1 × n sum-to-zero row (matching R-INLA's model = "seasonal"), and the remaining s - 2 null directions are identified by the likelihood.

One hyperparameter on the internal scale θ = log(τ).

SimplifiedLaplace()

Marginal strategy: Rue-Martino mean shift in the integration stage, plus Edgeworth first-order skewness correction in the per-coordinate density. Reduces to Gaussian when the likelihood third derivative ∇³_η log p is zero. R-INLA's strategy = "simplified.laplace".

See AbstractMarginalStrategy and ADR-016.

SkewNormalLikelihood(; link = IdentityLink(),
                    precision_prior = GammaPrecision(1.0, 5.0e-5),
                    skew_prior = GaussianPrior(0.0, 1.0))

Skew-normal observation model in R-INLA's family = "sn" parameterisation. With z_i = (y_i − η_i) √τ we have z_i ∼ f where

f(z) = (2 / ω_α) · φ((z − ξ_α) / ω_α) · Φ(α (z − ξ_α) / ω_α)

is the standardised skew-normal density (mean 0, variance 1). The shape α, location ξ_α, and scale ω_α depend only on the standardised skewness γ ∈ (−0.988, 0.988) via

γ = (4 − π) / 2 · (δ √(2/π))³ / (1 − 2δ²/π)^{3/2},   δ = α / √(1 + α²)
ω_α = 1 / √(1 − 2δ²/π)
ξ_α = − ω_α · δ · √(2/π)

so the mean of y_i is η_i and the variance is 1/τ. Two likelihood hyperparameters on the internal scale:

θ[1] = log τ
θ[2] = logit-skew, mapped to γ via  γ = 0.988 · tanh(θ[2] / 2)

The user-facing skewness is γ; the upper magnitude 0.988 matches R-INLA's hard cap (slightly tighter than the theoretical SN limit).

Defaults match R-INLA's family = "sn":

• precision_prior = GammaPrecision(1, 5e-5) ↔ loggamma(1, 5e-5) on log τ
• skew_prior = GaussianPrior(0, 1) for θ[2] (override of R-INLA's PC prior pc.sn, which has no closed-form Julia equivalent yet); pin both sides explicitly when validating an oracle.

StackedMapping(blocks::Tuple, rows::Vector{UnitRange{Int}})

Block-row stack of per-likelihood mappings. Each block shares the same latent vector x (shared ncols); each block contributes a contiguous slice of observation rows defined by rows[k].

Invariants

• length(blocks) == length(rows).
• All blocks have the same ncols.
• rows[k] are disjoint, contiguous, and sorted; together they cover 1:nrows(mapping).

The contiguous-and-sorted invariant is documented as the v0.2 design point in ADR-017's "Open questions" — interleaved row indices are deferred unless a real user needs them.

StudentTLikelihood(; link = IdentityLink(),
                   precision_prior = GammaPrecision(1.0, 1.0e-4),
                   dof_prior = GaussianPrior(2.5, 1.0))

Scaled Student-t observation model:

y_i = η_i + σ ε_i,    ε_i ~ Student-t(ν),    σ = 1 / √τ

with mean η_i (so the canonical link is identity), precision τ > 0, and degrees of freedom ν > 2. Two likelihood hyperparameters, θ = (log τ, log(ν − 2)). The ν > 2 floor (rather than ν > 0) ensures finite variance — this matches R-INLA's family = "T".

Density:

p(y | η, τ, ν) = Γ((ν+1)/2) / [√(πν) Γ(ν/2)] · √τ
                 · [1 + τ (y − η)² / ν]^{−(ν+1)/2}

Defaults match R-INLA's family = "T" in spirit (informative-precision gamma + mildly-informative Gaussian on log(ν − 2)); for PR-level oracle fits, override on both sides to keep the prior bit-for-bit identical.

UserComponent(callable; n, θ0 = Float64[]) -> UserComponent

R-INLA rgeneric equivalent. Wraps a user callable θ ↦ NamedTuple into a fully-functional AbstractLatentComponent. The callable is invoked at every θ visited by the inference pipeline and must return a NamedTuple with at least one key:

• Q::AbstractMatrix — symmetric length(c) × length(c) precision matrix (sparse strongly preferred; dense input is converted).

Optional keys, all with defaults:

• log_prior::Real — log-prior density of θ on the internal scale, including any user→internal Jacobian (default 0).
• log_normc::Real — per-component log normalizing constant contributing to the marginal-likelihood formula via log_normalizing_constant (default 0). For a proper Gaussian model in R-INLA's "mlik up to a constant" convention this is ½ log|Q| − ½ n log(2π); for an intrinsic model the structural ½ log|R̃|_+ is dropped (see log_normalizing_constant for the GMRFLib convention).
• constraint::Union{NoConstraint, LinearConstraint} — θ-independent hard linear constraint (default NoConstraint()). Read once at construction time by invoking the callable at θ0; later evaluations may include the key but its value is ignored.

Constructor arguments:

• callable — the function θ → NamedTuple described above. Receives the component's slice of the global hyperparameter vector.
• n::Integer — dimension of the latent field for this component.
• θ0::AbstractVector{<:Real} — initial hyperparameter values on the internal unconstrained scale; defaults to Float64[] (zero hyperparameters).

The two extension paths are complementary (see ADR-025):

• UserComponent is the callable path — a one-line port of an R-INLA rgeneric model definition. Use it for prototyping and for the long tail of R-INLA components we have not yet ported natively (crw2, besag2, besagproper, clinear, z, ou, …).
• Subtyping AbstractLatentComponent directly is the power-user path. Pick this when you also need to override prior_mean(c, θ) (a θ-dependent shifted prior, ADR-023) or gmrf(c, θ) (a custom AbstractGMRF factorization). See docs/src/extending.md.

Notes on cgeneric (C-callable user models)

R-INLA exposes cgeneric for C-coded user models — primarily for performance. In Julia this layer is unnecessary: the callable here is JIT-compiled to native code by Julia, with no FFI overhead. Users who already have a C library implementing the precision matrix can call it directly via @ccall inside the Julia callable; there is no cgeneric wrapper to learn.

Examples

A minimal IID component (Q = τ I) implemented as a UserComponent:

using LatentGaussianModels, SparseArrays
using LinearAlgebra: I

n = 50
component = UserComponent(n=n, θ0=[0.0]) do θ
    τ = exp(θ[1])
    return (; Q = sparse(τ * I, n, n),
              log_prior = -θ[1],     # `Exponential(1)` on τ
              log_normc = -0.5 * n * log(2π) + 0.5 * n * θ[1])
end

WeakPrior()

A deliberately-improper "nearly-flat" prior on the internal scale. Returns 0.0 for any θ. Exists so that quick-and-dirty fits (and EB-with-point-mass-prior tests) have something to dispatch on.

WeibullCureLikelihood(; link = LogLink(), censoring = nothing,
                        time_hi = nothing,
                        hyperprior_alpha = GammaPrecision(1.0, 0.001),
                        hyperprior_p     = LogitBeta(1.0, 1.0))

Weibull mixture-cure survival likelihood, matching R-INLA's family = "weibullcure". A latent fraction p ∈ (0, 1) of the population is cured (will never experience the event); the remaining 1 − p follow a Weibull distribution in the proportional- hazards parameterisation. With shape α > 0, cure fraction p, and the canonical LogLink,

λ_i  = exp(η_i),    u_i = λ_i t^α
S(t) = p + (1 − p) exp(−u_i)
F(t) = (1 − p) (1 − exp(−u_i))
f(t) = (1 − p) α t^(α−1) λ_i exp(−u_i)

Two likelihood hyperparameters, carried internally as θ_ℓ = [log α, logit p]. Defaults: loggamma(1, 0.001) on log α matching weibullsurv, and a uniform prior on p (i.e. LogitBeta(1, 1) on the logit scale).

censoring is an optional AbstractVector{Censoring} of length length(y); when nothing (default), every observation is uncensored (each contributes f(t), not the cured-mass marginal — there is no information about cured-vs-event without censoring). For INTERVAL rows, time_hi[i] is the upper bound (with y[i] the lower bound); time_hi is otherwise unread and may be nothing.

η-derivatives for the NONE/LEFT/INTERVAL branches coincide with WeibullLikelihood because the (1 − p) factor in f/F/F_hi − F_lo is η-independent. Only the RIGHT branch differs: with v = exp(−u) and D = p + (1 − p) v, the chain rule gives ∂η log D = D'/D, ∂²η log D = D''/D − (D'/D)², and ∂³η log D = D'''/D − 3 D' D''/D² + 2 (D'/D)³, where

D'   = −(1 − p) u v
D''  =  (1 − p) u v (u − 1)
D''' =  (1 − p) u v (3u − 1 − u²)

These reduce to plain Weibull's RIGHT branch as p → 0. See ADR-018 for the contract.

WeibullLikelihood(; link = LogLink(), censoring = nothing, time_hi = nothing,
                   hyperprior = GammaPrecision(1.0, 0.001))

Weibull survival likelihood in the proportional-hazards (PH) parameterisation, matching R-INLA's family = "weibullsurv" (variant 0). With shape α > 0 and the canonical LogLink,

λ_i  = exp(η_i)
h(t) = α λ_i t^(α-1)
S(t) = exp(-λ_i t^α)
f(t) = α λ_i t^(α-1) exp(-λ_i t^α)

so the cumulative hazard at t is u(t) = λ_i t^α. The shape α is the single likelihood hyperparameter, carried internally as θ_ℓ = [log α]. The default hyperprior matches R-INLA's weibullsurv default — loggamma(1, 0.001) on log α, encoded here as GammaPrecision(1.0, 0.001). (PCAlphaW (Sørbye-Rue 2017) lands in a follow-up PR per ADR-018; until then, GammaPrecision is the placeholder.)

censoring is an optional AbstractVector{Censoring} of length length(y); when nothing (default), every observation is uncensored. For INTERVAL rows, time_hi[i] is the upper bound (with y[i] the lower bound); time_hi is otherwise unread and may be nothing.

Closed-form derivatives are provided for LogLink for all four censoring modes; the η-derivative shapes coincide with ExponentialLikelihood once u = λ t is replaced by u = λ t^α and time differences by (t_hi^α − t_lo^α). See ADR-018 for the contract.

Example

# Right-censored survival data with shape α
ℓ = WeibullLikelihood(censoring = [NONE, RIGHT, NONE, RIGHT])

# Interval-censored
ℓ = WeibullLikelihood(
    censoring = [NONE, INTERVAL],
    time_hi   = [0.0, 5.0],
)

ZeroInflatedBinomialLikelihood0(n_trials; link = LogitLink(),
                                hyperprior = GaussianPrior(0.0, 1.0))

R-INLA family = "zeroinflatedbinomial0" — hurdle parameterisation.

P(y = 0)        = π
P(y = k | k>0)  = (1 - π) · C(n,k) p^k (1-p)^{n-k} / (1 - (1-p)^n)

Hyperparameter θ = logit(π); default prior gaussian(0, 1) on the internal scale matches R-INLA. Currently only LogitLink is supported on the count component.

ZeroInflatedBinomialLikelihood1(n_trials; link = LogitLink(),
                                hyperprior = GaussianPrior(0.0, 1.0))

R-INLA family = "zeroinflatedbinomial1" — standard mixture parameterisation.

P(y = 0)  = π + (1 - π) (1 - p)^n
P(y = k)  = (1 - π) C(n,k) p^k (1-p)^{n-k}     (k ≥ 1)

Hyperparameter θ = logit(π); default prior gaussian(0, 1) matches R-INLA. Currently only LogitLink is supported.

ZeroInflatedBinomialLikelihood2(n_trials; link = LogitLink(),
                                hyperprior = GaussianPrior(0.0, 1.0))

R-INLA family = "zeroinflatedbinomial2" — intensity-modulated mixture with π_i = 1 - p_i^α, p_i = sigmoid(η_i). Larger Bernoulli success probabilities push more mass into the count component; smaller ones inflate the zero. Single hyperparameter α > 0 carried internally as θ = log α. Default prior gaussian(0, 1) matches R-INLA.

ZeroInflatedNegativeBinomialLikelihood0(; link = LogLink(), E = nothing,
                                        hyperprior_size = GammaPrecision(1.0, 0.1),
                                        hyperprior_zi = GaussianPrior(0.0, 1.0))

R-INLA family = "zeroinflatednbinomial0" — hurdle parameterisation of the negative-binomial.

P(y = 0)        = π
P(y = k | k>0)  = (1 - π) · NB(k | s, μ) / (1 - NB(0 | s, μ))

with μ = E · exp(η) and size s (overdispersion). Two hyperparameters: θ[1] = log s (size, identical to plain NB) and θ[2] = logit(π). Defaults match R-INLA: loggamma(1, 0.1) on log s, gaussian(0, 1) on logit π.

ZeroInflatedNegativeBinomialLikelihood1(; link = LogLink(), E = nothing,
                                        hyperprior_size = GammaPrecision(1.0, 0.1),
                                        hyperprior_zi = GaussianPrior(0.0, 1.0))

R-INLA family = "zeroinflatednbinomial1" — standard mixture parameterisation:

P(y = 0)  = π + (1 - π) · NB(0 | s, μ)
P(y = k)  = (1 - π) · NB(k | s, μ)        (k ≥ 1)

with NB(0 | s, μ) = (s/(s+μ))^s. Hyperparameters and defaults match type 0.

ZeroInflatedNegativeBinomialLikelihood2(; link = LogLink(), E = nothing,
                                        hyperprior_size = GammaPrecision(1.0, 0.1),
                                        hyperprior_zi = GaussianPrior(0.0, 1.0))

R-INLA family = "zeroinflatednbinomial2" — intensity-modulated parameterisation: π_i = 1 - (μ_i / (1 + μ_i))^α. Two hyperparameters: θ[1] = log s (size) and θ[2] = log α (intensity-modulation exponent). Defaults: loggamma(1, 0.1) on log s, gaussian(0, 1) on log α.

ZeroInflatedPoissonLikelihood0(; link = LogLink(), E = nothing,
                                 hyperprior = GaussianPrior(0.0, 1.0))

R-INLA family = "zeroinflatedpoisson0" — hurdle parameterisation.

P(y = 0)        = π
P(y = k | k>0)  = (1 - π) · μ^k e^{-μ} / k! / (1 - e^{-μ})

The point mass at zero π replaces the count component's zero entirely; the count component is the Poisson distribution truncated to y ≥ 1.

Hyperparameter θ = logit(π) (one scalar, identical to R-INLA's internal scale). The default hyperprior is the R-INLA default gaussian(mean = 0, prec = 1) on the internal scale (σ = 1).

E is an optional offset / exposure vector. With LogLink, μ = E · exp(η). Currently only LogLink is supported.

ZeroInflatedPoissonLikelihood1(; link = LogLink(), E = nothing,
                                 hyperprior = GaussianPrior(0.0, 1.0))

R-INLA family = "zeroinflatedpoisson1" — standard mixture parameterisation.

P(y = 0)  = π + (1 - π) · e^{-μ}
P(y = k)  = (1 - π) · μ^k e^{-μ} / k!     (k ≥ 1)

The Poisson zero is not removed; π is an additional mixing probability of an extra mass at zero. This is the most common zero-inflated parameterisation in disease-mapping and ecology applications.

Hyperparameter θ = logit(π); default prior gaussian(0, 1) on the internal scale matches R-INLA.

ZeroInflatedPoissonLikelihood2(; link = LogLink(), E = nothing,
                                 hyperprior = GaussianPrior(0.0, 1.0))

R-INLA family = "zeroinflatedpoisson2" — intensity-modulated mixture, where the zero-inflation probability is a function of the mean:

π_i = 1 - (μ_i / (1 + μ_i))^α,    μ_i = E_i · exp(η_i)
P(y = 0)  = π_i + (1 - π_i) · e^{-μ_i}
P(y = k)  = (1 - π_i) · μ_i^k e^{-μ_i} / k!     (k ≥ 1)

Smaller intensities push more mass into the zero point — useful when the excess-zero process scales with the latent risk. The single hyperparameter is α > 0, carried internally as θ = log α. Default prior gaussian(0, 1) on θ matches R-INLA.

Internal-scale Normal prior on the Fisher-transformed correlation. Kept local to this file because the Normal-on-θ case recurs only here.

GMRFs.constraints(c::AbstractLatentComponent) -> AbstractConstraint

Default hard linear constraint attached to the component. Intrinsic components override; proper ones inherit NoConstraint.

GMRFs.constraints(c::BYM2) -> LinearConstraint

Sum-to-zero constraint on the spatial block u (positions n+1:2n), with one row per connected component of the graph. The b block is left unconstrained.

GMRFs.constraints(c::BYM) -> LinearConstraint

Sum-to-zero constraint on the Besag block b (positions n+1:2n), with one row per connected component of the graph. The v block is left unconstrained.

IID2D(n; hyperprior_precs = (PCPrecision(), PCPrecision()),
          hyperprior_corr  = PCCor0())

Bivariate IID component — ergonomic alias for IIDND(n, 2; …). Mirrors R-INLA's f(idx, model = "2diid", …) defaults.

IID3D(n; hyperprior_precs = (PCPrecision(), PCPrecision(), PCPrecision()),
          hyperprior_corrs = (PCCor0(), PCCor0(), PCCor0()))

Trivariate IID component — ergonomic alias for IIDND(n, 3; …). The three correlation slots are LKJ canonical partial correlations on atanh scale (see IIDND_Sep's docstring): θ entries 4 and 5 are direct correlations ρ_{1,2} and ρ_{1,3} on atanh scale; entry 6 is the partial correlation of variables 2 and 3 controlling for variable 1, also on atanh scale.

The closest R-INLA analogue is f(idx, model = "iid3d", …), which defaults to a Wishart prior on the joint precision; the separable default here uses PCCor0 per pairwise CPC and so will not match R-INLA's defaults exactly. Per ADR-022, the Wishart alternative lands in PR-1c.

IIDND(n, N = 2; hyperprior_precs, hyperprior_corr, hyperprior_corrs)

Multivariate IID latent component of dimension n × N. See ADR-022 for the parameterisation choice; PR-1a/PR-1b ship N ∈ {2, 3}. ADR-022 caps the separable form at N ≤ 3 (the atanh-of-each-pairwise-corr parameterisation is not injective onto positive-definite correlation matrices for N ≥ 4).

For N = 2, pass hyperprior_corr (a single AbstractHyperPrior on atanh ρ). For N ≥ 3, pass hyperprior_corrs as an NTuple of length N(N-1)/2.

Defaults: each marginal precision uses PCPrecision(); each correlation uses PCCor0() (R-INLA's pc.cor0(0.5, 0.5) default).

_bym2_params(θ) -> (τ, φ)

Map internal θ = [log τ, logit φ] to user-facing (τ, φ).

Distance d(φ) = sqrt(2 · KLD(φ)).

Internal: KLD on the non-null subspace between the alternative BYM2 covariance (1-φ) I + φ R_scaled^{-1} and the base I, in terms of the eigenvalues γ_k > 0 of R_scaled.

KLD(φ) = (1/2) Σ_k [ φ(1/γ_k - 1) - log1p(φ(1/γ_k - 1)) ]

Uses log1p for numerical accuracy when u = φ(1/γ_k - 1) is small.

Internal: derivative d'(φ) / φ — pulling out the leading factor of φ so the ratio is well-defined at φ → 0. From

KLD'(φ) = (φ/2) Σ_k (1/γ_k - 1)^2 / (1 + φ(1/γ_k - 1))
d'(φ)   = KLD'(φ) / d(φ)

we get, at finite φ > 0,

d'(φ) = KLD'(φ)/d(φ).

At φ → 0, d(φ) ~ A φ with A^2 = (1/2) Σ_k (1/γ_k - 1)^2, and d'(φ) → A.

_coerce_censoring(c)

Internal helper: accepts nothing, AbstractVector{Censoring}, or AbstractVector{Symbol} (e.g. [:none, :right, :none]) and returns the canonical Vector{Censoring} storage form. Anything else throws.

_constrained_marginal_variances(H_reg, constraint_data) -> Vector{Float64}

Per-coordinate conditional variances under the hard constraint:

Var(x_i | y, θ, C x = e) = (H_reg^{-1})_{ii} - (U W^{-1} U^T)_{ii}

For constraint_data === nothing, returns diag(H_reg^{-1}) unchanged. H_reg must be PD; callers supply the regularised posterior precision produced by laplace_mode.

_is_bad_theta_failure(err) -> Bool

Classify an exception thrown from inside laplace_mode as a "bad-θ" failure (numerical pathology that should activate the smooth penalty) versus a genuine bug (which must propagate). Treating only specific numerical exceptions as bad-θ keeps real bugs unmasked: a typo in a user-defined likelihood that throws MethodError should surface, not silently land in the penalty region.

Currently PosDefException, SingularException, and DomainError are recognised — these cover Cholesky failure on a near-singular H, LAPACK-side singular factorisations, and out-of-domain log-densities or component domain violations (e.g. AR1GMRF rejecting ρ = 1.0 when LBFGS overshoots atanh(ρ)). Add new types here when a new numerical-failure mode is identified; do not widen to a bare Exception catch.

By convention, ArgumentError is reserved for static shape/programming bugs (n ≥ 2, mismatched matrix dimensions) and is NOT in the bad-θ class — those should never be reachable from a θ-step and indicate a real bug. Component constructors should raise DomainError for parametric domain violations.

_apply_kriging!(x, C, e, cache) -> x

Project x onto {x : C x = e} using the factored H_reg in cache:

x ← x - U (C U)^{-1} (C x - e),   U = H_reg^{-1} C^T.

Returns x. Also returns the (U, W_fact) pair for reuse by the posterior-variance correction.

_neg_log_posterior_θ(m, y, strategy) -> function

Return a closure (θ, _p) -> -log π(θ | y) evaluated by Laplace. The returned function signature matches Optimization.jl's OptimizationFunction.

When the inner Laplace step throws (typically a PosDefException from Cholesky on an ill-conditioned Q + GᵀWG at extreme θ) or returns a non-finite log_marginal, the closure returns a smooth large penalty 1e10 + 1e3 · ‖θ‖² rather than Inf. This is load-bearing for the outer LBFGS line search (LineSearches.HagerZhang asserts isfinite(ϕ_c)), which can otherwise crash mid-bracket when its trial step lands in an infeasible region — see ADR-022 / IIDND_Sep{2} where ρ → ±1 saturation triggers it. The penalty's 1e10 floor is much larger than any feasible -log π(θ | y) we have ever observed (≲ 1e6 on Phase F–I oracles), so the optimum is unaffected, and the smooth quadratic in θ gives a usable FD gradient pointing back to the origin.

Only exceptions classified as "bad-θ" by _is_bad_theta_failure trigger the penalty; everything else is rethrown so genuine bugs in the model definition surface unmasked.

_null_bump(C::AbstractMatrix) -> SparseMatrixCSC

Return C^T (C C^T)^{-1} C as a sparse matrix. This is V_C V_C^T for the orthonormalised V_C = C^T (C C^T)^{-1/2} and adds a rank-k bump along range(C^T). PD-ness of Q + V_C V_C^T + A' D A is the responsibility of the caller — if null(Q) ⊋ range(C^T), the residual rd - k null directions must be covered by A' D A (typically true for informative observations).

_resolve_marginal_strategy(s) -> AbstractMarginalStrategy

Backwards-compatibility shim accepting either an AbstractMarginalStrategy instance (returned as-is) or a symbol from the legacy whitelist (:gaussian, :simplified_laplace, :full_laplace). Mirrors _resolve_scheme(::Symbol, ::Int) for the integration-scheme keyword.

Throws ArgumentError for unknown symbols.

_sla_mean_shift(lp::LaplaceResult, model::LatentGaussianModel, y)
    -> Vector{Float64}

Simplified-Laplace mean-shift correction at the Laplace fit lp:

Δx = ½ H⁻¹ Aᵀ (h³ ⊙ σ²_η)

where:

• H is the constraint-regularised Laplace Hessian factored in lp.factor.
• h³_i = ∇³_η log p(y_i | η_i, θ_ℓ) evaluated at η̂ = A x̂.
• σ²_η_i = (A H⁻¹ Aᵀ)_ii is the conditional variance of the linear predictor under the Laplace at θ, with the Rue-Held kriging correction applied when the model carries hard linear constraints.

Returns a zero vector when h³ ≡ 0 (Gaussian-likelihood collapse) so that latent_strategy = SimplifiedLaplace() reduces exactly to the Gaussian() path on quadratic-in-η likelihoods.

Cost per call: one multi-RHS sparse triangular solve (H⁻¹ Aᵀ, n_x × n_obs) plus one vector solve. With FactorCache reuse this is roughly one inner Newton iteration per integration point.

Internal: solve (1 - exp(-λ·d_U)) / (1 - exp(-λ·d_max)) = α for λ > 0. This is strictly monotone in λ: as λ → 0⁺ the ratio approaches d_U/d_max; as λ → ∞ it approaches 1. So a solution with λ > 0 exists iff α > d_U/d_max; when α ≤ d_U/d_max we return the limit λ = 0, which corresponds to a uniform prior on d.

_θ_mode_and_hessian(m, y, strategy) -> (θ̂, H, optim_result)

Find the θ-mode via LBFGS and compute the Hessian of the negative log posterior at the mode by finite differences.

add_copy_contributions!(η_block, ℓ::AbstractLikelihood, x, θ_ℓ) -> η_block

Add this likelihood's Copy-component contributions to its slice of the linear predictor, in place, and return η_block.

Copy (R-INLA f(..., copy=...)) shares a latent component across linear predictors with an estimated scaling β: the receiving block gets η_target[i] += β * x_source[k(i)]. ADR-021 places β on the receiving likelihood (rather than on the projection mapping), so each likelihood that opts in stores its (source_indices, β_index) pairs and reads β from its own θ_ℓ slice.

The default no-op covers every likelihood without copies — Gaussian, Poisson, Binomial, NegBinomial, Gamma, the survival likelihoods on arms that don't share latents — and leaves η_block unchanged. The inner Newton loop calls this hook after every η = mapping * x evaluation, before computing likelihood derivatives.

Implementing likelihoods receive an η_block view spanning their observation rows, the full latent vector x (so they can index any source component), and their hyperparameter view θ_ℓ (where β lives).

apply!(η, mapping, x) -> η

In-place forward multiply: η .= mapping * x. Concrete subtypes must specialise.

apply_adjoint!(g, mapping, r) -> g

In-place adjoint multiply: g .= mappingᵀ * r. Concrete subtypes must specialise.

as_matrix(mapping) -> AbstractMatrix

Concrete materialization of the mapping as an AbstractMatrix. Used by inference sites that need to multi-RHS solve through the projector or build the Aᵀ D A Hessian quadratic form. Sparse where the underlying mapping is sparse.

component_range(model, i::Int) -> UnitRange{Int}

Index range of the i-th component in the stacked latent vector x.

compute_skewness_corrections(log_post, θ̂, Σ;
                             threshold = 0.05,
                             max_stretch = 5.0)
    -> (stdev_corr_pos, stdev_corr_neg)

Estimate the per-eigen-axis Rue–Martino–Chopin (2009) §6.5 stretches that align a Gaussian midpoint grid with the actual density-shape asymmetry of π(θ|y).

Along eigen-axis k of Σ, the stretch on the positive side is stdev_corr_pos[k] = √(½ / Δ⁺_k) where Δ⁺_k = log π(θ̂) - log π(θ̂ + halfσ_k e_k) and halfσ is Σ^{1/2} expressed in the eigenbasis. Symmetric correction: for an exactly Gaussian posterior the drop in log-density across one nominal std along any axis is exactly ½, so the stretch is 1.

log_post is the callable θ → log π̂(θ|y) (up to a θ-independent constant — the constant cancels out of the ratios). Pass any function of θ only (no extra arguments). Internally the helper queries log_post(θ̂) once and log_post(θ̂ ± halfσ e_k) for each k, so the total cost is 2m + 1 evaluations.

Stretches are floored at threshold and capped at max_stretch:

• Δ < threshold (very flat axis or numerical noise) → stretch left at 1.0,
• otherwise stretch is clamped to [1/max_stretch, max_stretch] to avoid runaway proposals when one side of the posterior is a stiff wall (Δ becomes large) or pathologically flat.

Example

f = θ -> -negative_log_post(θ)        # log π̂(θ|y) up to a constant
stdev_pos, stdev_neg = compute_skewness_corrections(f, θ̂, Σθ)
scheme = Grid(n_per_dim = 21, span = 4.0,
              stdev_corr_pos = stdev_pos,
              stdev_corr_neg = stdev_neg)

coxph_design(aug::CoxphAugmented, X::AbstractMatrix) -> SparseMatrixCSC

Assemble the joint design matrix for the augmented Cox-PH-as-Poisson fit: a sparse [B | Xrep] block where B is the one-hot baseline-interval indicator and Xrep repeats each subject's covariate row across that subject's augmented rows.

X must have aug.n_subjects rows.

coxph_design(aug::CoxphAugmented) -> SparseMatrixCSC

Baseline-only design matrix (no covariates).

cpo(rng, res, model, y; n_samples = 1000)
  -> @NamedTuple{CPO, log_CPO, log_pseudo_marginal}

Held-Schrödle-Rue CPO via the harmonic-mean identity

CPO_i⁻¹ = E[ 1 / p(y_i | η_i, θ_ℓ) ]

where the expectation is over the posterior (η_i, θ_ℓ) | y. The pseudo-marginal likelihood Σ log CPO_i is returned as a summary log-predictive score.

The harmonic-mean estimator is known to be noisy when a few observations have very small likelihood values; R-INLA flags these as "failures" via a failure indicator. We do not replicate that flag yet — outliers will produce large log_CPO_i variance across RNG seeds. Increasing n_samples stabilises the estimator.

dic(res, model, y) -> @NamedTuple{DIC, pD, D_bar, D_mode}

Spiegelhalter DIC, computed without sampling:

• D_bar = E_{x,θ | y}[-2 log p(y | η, θ_ℓ)] estimated by averaging the per-θ_k Laplace approximation with a second-order Taylor correction in the posterior variance of η_i.
• D_mode = -2 log p(y | η̂, θ̂_ℓ) evaluated at the posterior means.
• pD = D_bar - D_mode (effective number of parameters).
• DIC = D_bar + pD.

The second-order correction is

E[-2 log p(y_i|η_i, θ_ℓ)] ≈ -2 log p(y_i|η̂_i, θ_ℓ)
                           + (-∇²_η log p)_i · Var(η_i | θ_k)

which is exact for Gaussian identity-link and a good approximation for Poisson/Binomial near the mode. For more complex likelihoods the MC-based waic is a robust alternative.

empirical_bayes(m, y; kwargs...)

Alias for fit(m, y, EmpiricalBayes(; kwargs...)).

fit(m::LatentGaussianModel, y, strategy::EmpiricalBayes) -> EmpiricalBayesResult

Empirical-Bayes fit. Convenience alias: empirical_bayes(m, y; kwargs...).

fit(m, y) -> EmpiricalBayesResult

Convenience: default strategy is EmpiricalBayes().

fixed_effects(model, res; level = 0.95) -> Vector{@NamedTuple{...}}

Posterior summaries for each scalar "fixed-effect"-shaped component — those whose latent dimension is 1 (e.g. Intercept, scalar slopes). Each element has fields (name, mean, sd, lower, upper) with the (1-level)/2-quantile intervals computed from the Gaussian at (x_mean, x_var).

Components of length > 1 are surfaced through random_effects.

hyperparameters(model, res; level = 0.95) -> Vector{@NamedTuple{...}}

Gaussian-approximation summaries of the hyperparameters on the internal scale: one row per entry of θ, ordered as (likelihood, component_1, ..).

Each row has (name, mean, sd, lower, upper). For display on user scale, apply the relevant user_scale transforms (e.g. exp for log-precisions).

initial_hyperparameters(c::AbstractLatentComponent)

Initial values on the internal unconstrained scale.

initial_hyperparameters(ℓ::AbstractLikelihood) -> Vector

Initial values on the internal (unconstrained real-valued) scale. A likelihood with a Gaussian precision τ uses log(τ) internally; this function returns [0.0] (i.e. τ = 1 at init).

initial_hyperparameters(m::LatentGaussianModel) -> Vector{Float64}

Stack the initial internal-scale hyperparameters of all likelihoods and components, in the canonical [θ_ℓ; θ_c] order.

inla(m, y; kwargs...)

Convenience alias for fit(m, y, INLA(; kwargs...)).

inla_coxph(time, event;
           breakpoints = nothing, nbreakpoints = 15) -> CoxphAugmented

Cox proportional-hazards data augmentation, mirroring R-INLA's inla.coxph. Each subject is expanded into one row per baseline-hazard interval [τ_k, τ_{k+1}) they survive into; the augmented response becomes a binary indicator and the exposure offset becomes the time spent in that interval. The augmented data are then suitable for a Poisson regression with linear predictor

η_ik = γ_k + x_iᵀ β

where γ_k is the piecewise-constant baseline log-hazard for interval k (typically given an RW1 prior). The piecewise-exponential likelihood for time-to-event under right-censoring is recovered exactly by this Poisson augmentation (Holford 1980; Laird & Olivier 1981; Clayton 1991).

Arguments

• time::AbstractVector{<:Real} — observation times (>0).
• event::AbstractVector{<:Integer} — event indicators (1 = uncensored event, 0 = right-censored).

Keyword arguments

• breakpoints::AbstractVector{<:Real} — explicit interval boundaries. Must be sorted, start at 0, and end at or above maximum(time). When nothing (default), quantile-based breakpoints over the observed event times are constructed (see nbreakpoints).
• nbreakpoints::Integer = 15 — number of interior quantile-based breakpoints when breakpoints is not given.

Example

aug = inla_coxph(time, event)

ℓ           = PoissonLikelihood(E = aug.E)
c_baseline  = RW1(aug.n_intervals; hyperprior = PCPrecision(1.0, 0.01))
c_beta      = FixedEffects(size(X, 2))
A           = coxph_design(aug, X)

model = LatentGaussianModel(ℓ, (c_baseline, c_beta), A)
res   = inla(model, aug.y)

inla_summary([io], model, res; level = 0.95) -> Nothing

Print an R-INLA-shaped summary of a fitted INLAResult: fixed effects, random effects (one block per vector-valued component), hyperparameters on the internal scale, and the marginal log-likelihood. io defaults to stdout.

This mirrors the layout of summary.inla() without the parts (DIC / WAIC / CPO) that require diagnostics — those are computed on demand by dic, waic, cpo, pit.

integration_nodes(scheme, θ̂, Σ) -> (points, log_weights)

Return the design points in the native θ-parameterisation and their log-weights. The weights carry the Gaussian base measure already; caller multiplies by exp(Δ log π) to get weights under the true posterior, then renormalises. Returning log-weights avoids underflow at high dimension.

inverse_link(g::AbstractLinkFunction, η) -> μ

Apply the inverse link. Works scalar-in / scalar-out and broadcasts.

joint_apply_copy_contributions!(η, m::LatentGaussianModel, x, θ) -> η

Apply each likelihood's add_copy_contributions! hook to its own slice of η, in place, and return η. Called after every η = mapping * x evaluation in the inner Newton loop and posterior sampler so that Copy-scaled latent shares fold into the linear predictor before likelihood derivatives are computed.

The default add_copy_contributions! is a no-op, so single-likelihood models without copies pay only the per-block dispatch loop — concrete Tuple{Vararg{<:AbstractLikelihood}} likelihoods union-split at compile time, so the no-op specialisation lowers to a no-op.

joint_effective_jacobian(m::LatentGaussianModel, θ) -> AbstractMatrix

Effective design Jacobian dη/dx, including any Copy contributions on the receiving likelihoods.

For models without copies this returns as_matrix(m.mapping) directly (no allocation overhead). For models with copies, the Copy share η_target[i] += β * x[source_indices[i]] adds a β entry to J[block_rows[k][i], source_indices[i]], returned as A + B(θ).

Used by every inference site that needs dη/dx — the inner Newton Hessian/gradient, simplified-Laplace mean-shift, posterior-marginal skewness, and DIC.

joint_log_density(m::LatentGaussianModel, y, η, θ) -> Real

Sum of per-block likelihood log-densities, routed through m.likelihoods and m.block_rows. Reduces to a single log_density(m.likelihoods[1], y, η, θ_ℓ) call when there is exactly one block.

joint_pointwise_cdf(m::LatentGaussianModel, y, η, θ) -> AbstractVector

Length-n_obs per-observation predictive CDF F(y_i | η_i, θ_ℓ). Used by PIT diagnostics. Each block must implement pointwise_cdf on its likelihood.

joint_pointwise_log_density(m::LatentGaussianModel, y, η, θ) -> AbstractVector

Length-n_obs per-observation log p(y_i | η_i, θ_ℓ). Sums to joint_log_density.

joint_precision(m::LatentGaussianModel, θ) -> SparseMatrixCSC

Block-diagonal latent precision Q(θ) = blockdiag(Q_1, ..., Q_K) where Q_i = precision_matrix(components[i], θ[θ_ranges[i]]).

joint_prior_mean(m::LatentGaussianModel, θ) -> Vector

Stacked latent prior mean μ(θ) = [μ_1(θ_1); ...; μ_K(θ_K)] where μ_i = prior_mean(components[i], θ[θ_ranges[i]]). The default per component is zeros, so this returns a zero vector for every model that does not use measurement-error or other shifted-prior components.

ADR-023 promotes prior_mean from documented-optional to load-bearing in the inner Newton loop: the latent quadratic is ½ (x - μ)' Q(θ) (x - μ), not ½ x' Q x.

joint_∇_η_log_density(m::LatentGaussianModel, y, η, θ) -> AbstractVector

Length-n_obs gradient of Σ_k log p(y_k | η_k, θ_ℓ_k) w.r.t. η, assembled by routing each block through its ∇_η_log_density.

joint_∇²_η_log_density(m::LatentGaussianModel, y, η, θ) -> AbstractVector

Length-n_obs diagonal of the η-Hessian. Same routing pattern as joint_∇_η_log_density.

joint_∇³_η_log_density(m::LatentGaussianModel, y, η, θ) -> AbstractVector

Length-n_obs diagonal of the third-derivative tensor of the joint likelihood w.r.t. η. Used by simplified-Laplace mean-shift and posterior-marginal skewness.

laplace(m, y, θ; kwargs...)

Laplace fit at fixed θ. Returns a LaplaceResult.

laplace_mode(model, y, θ;
             strategy = Laplace(), x0 = nothing,
             extra_constraint = nothing) -> LaplaceResult

Find the mode x̂ of log p(x | θ, y) ∝ log p(y | A x, θ_ℓ) - ½ x' Q(θ) x by Newton iteration. Returns a LaplaceResult.

extra_constraint, if supplied, is a NamedTuple{(:rows, :rhs)} whose rows are stacked onto the model-level constraint produced by model_constraints. Used by FullLaplace to fix a single coordinate x_i = a for the per-x_i Laplace refit; any caller-supplied linear equality is acceptable as long as the augmented constraint matrix is full-rank.

laplace_mode_fixed_xi(model, y, θ, i, a;
                      strategy = Laplace(),
                      x0 = nothing) -> LaplaceResult

Constrained Laplace fit with the additional hard linear constraint x_i = a stacked onto the model-level constraints declared by model_constraints. Forwards to laplace_mode with extra_constraint = (rows = e_i^T, rhs = [a]).

x0 (optional) seeds the inner Newton iteration. Used by FullLaplace to warm-start adjacent grid-point refits from the previous fit's mode — for a smooth posterior the conditional mode at a + Δa differs from the mode at a by O(Δa), so 1-2 Newton steps suffice in place of a fresh ~5-10 step descent.

Used by FullLaplace inside posterior_marginal_x.

likelihood_for(mapping, i) -> Int

Block index of the likelihood applied to observation row i. Defaults to 1 for any single-block mapping; multi-block mappings (e.g. StackedMapping) override.

link(ℓ::AbstractLikelihood) -> AbstractLinkFunction

The link function attached to this likelihood instance.

log_density(ℓ::AbstractLikelihood, y, η, θ) -> Real

Scalar log-density log p(y | η, θ).

log_hyperprior(c::AbstractLatentComponent, θ)

Log-prior density on the internal scale; includes any Jacobians from the user-to-internal transformation. Zero-hyperparameter components return 0.

log_hyperprior(ℓ::AbstractLikelihood, θ) -> Real

Log-prior density for the likelihood's hyperparameters, evaluated on the internal scale. Zero when nhyperparameters(ℓ) == 0.

log_hyperprior(m::LatentGaussianModel, θ) -> Real

Sum of log-priors on the internal scale for all likelihoods and components.

log_marginal_likelihood(res::INLAResult) -> Float64

Estimate of log p(y) (model evidence) from the importance-sampling normalising constant.

log_normalizing_constant(c::AbstractLatentComponent, θ) -> Real

Per-component log normalizing constant of the prior, in the R-INLA convention. Used by the Laplace marginal-likelihood formula

mlik ≈ log p(y|x*) - ½ x*'Q x* + ½ (n_x - r) log(2π)
       - ½ log|H_C| + Σ_i log_normc_i(θ_i)

where log|H_C| = log|H| + log|C H⁻¹ C^T| - log|C C^T| and the sum runs over components. The convention follows R-INLA's GMRFLib: each proper component contributes the full Gaussian log-NC -½ d log(2π) + ½ log|Q_i|, while intrinsic components drop the structural log-determinant ½ log|R̃|_+ (it is absorbed into the global ½ (n_x - r) log(2π) and log|H_C| corrections).

Default: zero. Components must override to participate correctly in the marginal likelihood.

log_prior_density(prior::AbstractHyperPrior, θ) -> Real

Log-prior density at internal-scale value θ, including Jacobian.

log_prior_density(p::PCBYM2Phi, θ) -> Real

Log-density of the PC prior on φ = σ(θ) transformed to the internal scale θ = logit(φ). Comprises

• the PC density on φ: log π_φ(φ) = log λ - λ d(φ) + log|d'(φ)|
• the Jacobian |dφ/dθ| = φ(1-φ).

At θ → ±∞, the prior mass goes to zero: log π(θ) → -∞.

logsubexp(a, b)

log(exp(a) - exp(b)) for a > b, computed via log1p for digit retention near a ≈ b.

model_constraints(m::LatentGaussianModel) -> AbstractConstraint

Assemble the model-level hard constraint by stacking each component's GMRFs.constraints(c) block into the full (k_total × n_x) constraint matrix. Returns NoConstraint() if no component declares a constraint.

multinomial_design_matrix(helper, X; reference_class = helper.K)

Build the class-specific covariate block of the long-format design matrix from an (n_rows, p) covariate matrix X and the layout helper returned by multinomial_to_poisson.

Returns a sparse (n_long, (K - 1) * p) matrix A_β. The rows match the order of helper.row_id, helper.class_id; the reference class's block of p coefficients is dropped to identify the model (corner-point parameterisation, Agresti 2010 §8.5).

For long-format index idx with i = helper.row_id[idx], k = helper.class_id[idx]:

• if k != reference_class: columns ((k′ - 1) * p + 1) : (k′ * p) carry x_i^⊤, where k′ ∈ 1:(K-1) is the index of k after dropping the reference.
• if k == reference_class: the row is all zeros (the linear predictor for the reference class is just the per-row α_i).

The per-row nuisance intercept α_i is not included here; attach it as a separate IID block — see the multinomial_to_poisson docstring for the recipe.

multinomial_to_poisson(Y; class_names = 1:K)

Reshape a multinomial-counts matrix Y::AbstractMatrix{<:Integer} of shape (n_rows, K) into the long-format independent-Poisson layout per ADR-024. Returns a NamedTuple:

• y::Vector{Int} — length n_rows * K count vector with row-major layout y[(i-1)*K + k] = Y[i, k].
• row_id::Vector{Int} — same length; row_id[(i-1)*K + k] = i.
• class_id::Vector{Int} — same length; class_id[(i-1)*K + k] = k.
• n_rows::Int, K::Int, n_long::Int = n_rows * K.
• class_names::Vector — pass-through of the supplied class names.

The reformulation is the Multinomial-Poisson trick (Baker 1994; Chen 1985): the multinomial likelihood Y_i ~ Multinomial(N_i, π_i) is equivalent to Y_ik ~ Poisson(λ_ik) with λ_ik = exp(α_i + x_i^⊤ β_k), where α_i is a per-row nuisance intercept attached as IID(n_rows; τ_init = -10.0, fix_τ = true) (matching R-INLA's prec = list(initial = -10, fixed = TRUE)).

Use multinomial_design_matrix to build the class-specific covariate block; combine with the IID α block via hcat/StackedMapping to get the full LGM design matrix.

n_hyperparameters(m) -> Int

Number of hyperparameters (all likelihoods + all components).

n_latent(m) -> Int

Length of the stacked latent vector x.

n_likelihood_hyperparameters(m) -> Int

Total number of likelihood-attached hyperparameters across all blocks.

n_likelihoods(m) -> Int

Number of likelihood blocks. Equal to length(m.likelihoods).

n_observations(m) -> Int

Number of observation rows. Equal to nrows(m.mapping).

ncols(mapping) -> Int

Column count of the implicit matrix; equal to length(x).

nhyperparameters(c::AbstractLatentComponent) -> Int

nhyperparameters(ℓ::AbstractLikelihood) -> Int

Number of scalar hyperparameters attached to the likelihood. Defaults to 0; Gaussian/NegBinomial/Gamma override.

nrows(mapping) -> Int

Row count of the implicit matrix; equal to length(η).

pit(rng, res, model, y; n_samples = 1000) -> Vector{Float64}

Posterior-predictive probability integral transform

PIT_i = E[ F(y_i | η_i, θ_ℓ) ]

evaluated via Monte-Carlo. For a well-specified model PIT values are approximately uniform on [0, 1]; deviations flag miscalibration. For discrete likelihoods this is the naive (non-randomised) variant — the CDF is evaluated at the observed integer, so ties at the endpoints appear but summary statistics still diagnose gross miscalibration.

Requires the likelihood to implement pointwise_cdf.

pointwise_cdf(ℓ, y, η, θ) -> AbstractVector{<:Real}

Per-observation CDF F(y_i | η_i, θ) = P(Y_i ≤ y_i | η_i, θ), length length(y). Needed for PIT diagnostics. Not all likelihoods have a closed-form CDF — the default raises, and concrete likelihoods that support PIT implement this method.

pointwise_log_density(ℓ, y, η, θ) -> AbstractVector{<:Real}

Per-observation log-density log p(y_i | η_i, θ), length length(y). Likelihoods factorise across observations, so the sum of this vector equals log_density(ℓ, y, η, θ).

Needed by WAIC, CPO, DIC, and PIT diagnostics. Default implementation falls back on log_density evaluated on singletons — correct but inefficient; concrete likelihoods should override with a vectorised form.

posterior_marginal_x(res::INLAResult, i::Int;
                     strategy = Gaussian(),
                     model = nothing, y = nothing,
                     grid_size = 75, span = 5.0,
                     grid = nothing) -> @NamedTuple{x::Vector, pdf::Vector}

Posterior density of the i-th latent coordinate, evaluated on a grid.

Returns a named tuple (x, pdf). The density is the θ-mixture

p(x_i | y) ≈ Σ_k w_k · π_k(x_i)

where π_k depends on strategy (an AbstractMarginalStrategy instance, or one of the legacy symbols :gaussian / :simplified_laplace — see ADR-026):

• Gaussian() / :gaussian (default) — π_k = φ(x_i; mode_k[i], var_k[i]).
• SimplifiedLaplace() / :simplified_laplace — π_k = φ · (1 + γ_k / 6 · H₃(s)) with s = (x_i - mode_k[i]) / σ_k, H₃(s) = s³ - 3s, and posterior skewness γ_k = κ_3(x_i|θ_k) / σ_k³. The third cumulant is assembled from ∇³_η_log_density and the Laplace precision at θ_k; for a Gaussian likelihood this collapses to the Gaussian strategy. Requires model and y to be supplied.
• FullLaplace() / :full_laplace — at each grid point a in xs, refit the joint Newton mode under the augmented constraint e_i^T x = a (via laplace_mode_fixed_xi) and read the constrained Laplace log-marginal. Per-θ density is the renormalised ratio exp(log p̂(y | θ_k, x_i = a) - log p̂(y | θ_k)); the mixture is again renormalised on the grid. Requires model and y.

If grid is supplied it is used verbatim; otherwise a grid of grid_size equally spaced points spanning ±span · √posterior_var about the posterior mean is generated.

posterior_marginal_θ(res::INLAResult, j::Int;
                     grid_size = 75, span = 5.0,
                     grid = nothing) -> @NamedTuple{θ::Vector, pdf::Vector}

Gaussian marginal of the j-th hyperparameter on the internal scale, centred at res.θ̂[j] with standard deviation √res.Σθ[j,j]. This is the "gaussian" strategy in R-INLA. A numerically integrated density over the INLA design is a future extension.

posterior_predictive(rng, res, model, mapping; n_samples = 1000)
  -> @NamedTuple{x::Matrix{Float64}, θ::Matrix{Float64}, η::Matrix{Float64}}

Posterior predictive samples of the linear predictor η_new = A_new x at a new observation mapping. mapping may be either an AbstractObservationMapping (e.g. LinearProjector, IdentityMapping, StackedMapping) or an AbstractMatrix A_new — matrices are wrapped in LinearProjector automatically, mirroring LatentGaussianModel's convenience constructor.

The function returns the joint draws (x, θ) from posterior_sample plus η::Matrix{Float64} of size nrows(mapping) × n_samples. Each column is η_s = mapping * x_s.

η is the foundation for downstream predictive inference: applying the likelihood's inverse link gives μ posterior samples, and sampling y_new ∼ p(y | η, θ) per likelihood gives full posterior predictive y samples. The latter is left to per-likelihood sampling code (Phase K follow-up); this function ships the likelihood-agnostic part.

mapping must satisfy ncols(mapping) == n_latent(model).

Example

res = inla(model, y)
# Predict at new covariate rows X_new (rows match the original design):
draws = posterior_predictive(rng, res, model, X_new; n_samples = 500)
μ_new = inverse_link.(link(model.likelihood), draws.η)

posterior_predictive_y(rng, res, model; n_samples = 1000)
  -> @NamedTuple{x::Matrix{Float64}, θ::Matrix{Float64},
                 η::Matrix{Float64}, y_rep::Matrix{Float64}}

Posterior-predictive samples on the response scale. Extends posterior_predictive with a per-likelihood sample_y draw, returning the simulated replicate observations y_rep of size n_obs × n_samples alongside the underlying joint draws (x, θ, η).

For multi-likelihood models the response sampler dispatches per block: each likelihood m.likelihoods[k] draws a replicate of length length(m.block_rows[k]) from p(y | η_block, θ_ℓ_k) using its own metadata (Binomial n_trials, Poisson/NegBin offset E, …) read directly from the likelihood instance.

y_rep is always returned as Matrix{Float64} regardless of the likelihood family — counts arrive as integer-valued floats, continuous likelihoods as real-valued floats — to keep multi-likelihood blocks typed-uniformly.

Used by posterior-predictive checks (Gelman et al. 2014; bayesplot's pp_check). To overlay densities of y_rep against the observed response, draw n_samples ≈ 200–1000 and plot a subset of the columns of y_rep against y_obs.

Example

res = inla(model, y)
draws = posterior_predictive_y(rng, res, model; n_samples = 400)
# draws.y_rep is n_obs × 400 — each column a posterior-predictive
# replicate of the full observation vector.

Note

Likelihoods without sample_y defined (currently: censored survival families and zero-inflated likelihoods) raise ArgumentError. Closed-form samplers ship for GaussianLikelihood, PoissonLikelihood, BinomialLikelihood, NegativeBinomialLikelihood, and GammaLikelihood.

posterior_sample(rng, res, model; n_samples = 1000)
  -> @NamedTuple{x::Matrix{Float64}, θ::Matrix{Float64}}

Joint draws (x, θ) ∼ π̂(· | y) from the INLA posterior. The output matrices have one column per draw:

• x is n_x × n_samples — full latent vector (component blocks concatenated in joint order).
• θ is n_θ × n_samples — full hyperparameter vector (likelihood block first, then per-component blocks; same ordering as hyperparameters and res.θ_points).

The sampler picks θ_k from the discrete integration design with probabilities res.θ_weights, then draws x | θ_k ∼ N(mode_k, H_k⁻¹) using the cached Cholesky factor at θ_k. If the component stack declared hard linear constraints, the kriging correction is applied so each x sample satisfies C x = e to working precision.

Use posterior_samples_η instead when you want draws of η = A x and the likelihood-only hyperparameters θ_ℓ (the form that drives WAIC / CPO / PIT). posterior_sample is the right building block for joint-x summaries (e.g. random-effect contrast posteriors) and Stan/NUTS triangulation.

When n_hyperparameters(model) == 0 (the dim(θ)=0 fast path used by ADR-024's multinomial-via-Poisson), the returned θ matrix has zero rows.

posterior_samples_η(rng, res, model; n_samples = 1000)
  -> @NamedTuple{η::Matrix{Float64}, θℓ::Matrix{Float64}}

Joint draws of the linear predictor η = A x and the likelihood hyperparameters θ_ℓ from the INLA posterior. Returns an n_obs × n_samples matrix η and an n_ℓ × n_samples matrix θℓ (empty rows when the likelihood has no hyperparameters).

The sampler picks θ_k from the discrete integration design with probabilities res.θ_weights, then draws x | θ_k ∼ N(mode_k, H_k⁻¹) using the cached Cholesky factor at θ_k. If the component stack declared hard linear constraints, the kriging correction is applied to each draw so that C x = e to working precision.

pp_check(rng, res, model, y_obs; n_samples = 400)
  -> @NamedTuple{y::AbstractVector, y_rep::Matrix{Float64}}

Posterior-predictive-check data: the observed response y_obs paired with n_samples replicate response vectors y_rep drawn under the fitted model. Convenience wrapper over posterior_predictive_y that drops the (x, θ, η) joint draws and keeps just the data needed for pp-check overlays.

n_samples defaults to 400 — enough for stable density estimates, small enough that overlaying every column of y_rep is reasonable.

Plotting (Makie)

using GLMakie, LatentGaussianModels
ck = pp_check(rng, res, model, y_obs; n_samples = 400)

fig = Figure(); ax = Axis(fig[1, 1])
for j in axes(ck.y_rep, 2)[1:50]                         # 50 replicates
    density!(ax, ck.y_rep[:, j]; color = (:gray, 0.4))
end
density!(ax, ck.y; color = :black, linewidth = 3)
fig

The y-axis covers min(y_rep) to max(y_rep) automatically; a clear visual gap between the observed and replicated densities indicates mis-specification (location, scale, or shape).

precision_matrix(c::AbstractLatentComponent, θ)

Sparse symmetric precision matrix at θ (this component's slice).

prior_mean(c::AbstractLatentComponent, θ) -> AbstractVector

Prior mean vector. Default: zeros of length length(c).

prior_name(prior::AbstractHyperPrior) -> Symbol

psis_loo(rng, res, model, y; n_samples = 1000)
  -> @NamedTuple{elpd_loo, looic, pointwise_elpd_loo, pointwise_p_loo,
                 p_loo, pareto_k}

Pareto-smoothed importance sampling estimate of the leave-one-out expected log pointwise predictive density (Vehtari, Gelman & Gabry 2017; Vehtari, Simpson, Gelman, Yao & Gabry 2024).

Returns:

• elpd_loo — sum of pointwise LOO elpd (higher is better).
• looic = -2 · elpd_loo — information-criterion form.
• pointwise_elpd_loo — per-observation elpdloo, length `nobs`.
• pointwise_p_loo — per-observation effective parameter contribution.
• p_loo — sum of pointwise effective parameters.
• pareto_k — vector of fitted Pareto shape parameters k̂_i. Values above 0.7 indicate the LOO importance ratios for that observation have a heavy tail and the PSIS estimate is unreliable; refit-based LOO is the gold-standard remedy. Values above 0.5 (but below 0.7) are a soft warning.

PSIS-LOO supplements cpo: both target the per-observation predictive distribution but PSIS-LOO is more robust to outliers and provides the diagnostic k̂ for failure detection.

This function requires the PSIS.jl weakdep:

using PSIS, LatentGaussianModels
res = inla(model, y)
loo = psis_loo(rng, res, model, y; n_samples = 2000)

Without PSIS loaded a MethodError is raised.

References

• Vehtari, A., Gelman, A. & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27(5), 1413-1432.
• Vehtari, A., Simpson, D., Gelman, A., Yao, Y. & Gabry, J. (2024). Pareto smoothed importance sampling. JMLR, 25(72), 1-58.

random_effects(model, res; level = 0.95)
  -> Dict{String, @NamedTuple{mean, sd, lower, upper}}

Per-component posterior summaries for vector-valued components (length > 1). Each entry is a NamedTuple of vectors — one entry per latent coordinate within that component.

refine_hyperposterior(res, model, y;
                      n_grid = 15, span = 4.0,
                      skewness_correction = true,
                      laplace = Laplace(),
                      latent_strategy = Gaussian())
  -> INLAResult

Re-evaluate the INLA posterior on a denser hyperparameter grid given an existing fit. R-INLA equivalent: inla.hyperpar.

Reuses res.θ̂ and res.Σθ from the input fit (skipping the outer LBFGS / FD-Hessian stage) and re-runs only the integration stage of the INLA pipeline against a fresh Grid(n_grid, span). The denser quadrature improves IS estimates of log p(y), x_mean, x_var, and θ_mean, and gives smoother plotted hyperparameter densities. On a well-conditioned posterior the refined log_marginal should agree with the original within IS error.

Skewness correction is on by default here (versus inla()'s false default): the typical reason a user calls refine_hyperposterior is that the hyperparameter posterior is heavy-tailed enough that the default 5×5 design under-resolves the tails, and asymmetric per-axis stretches are exactly the right knob for that case.

Throws ArgumentError if n_hyperparameters(model) == 0 — the fixed-θ fast path has no hyperparameter design to refine.

sample_conditional(model, θ, y; rng = Random.default_rng(),
                   laplace = Laplace()) -> Vector{Float64}
sample_conditional(model, θ, y, n_samples; rng, laplace) -> Matrix{Float64}

Draw samples of the latent vector x from the Laplace approximation to p(x | θ, y). Paired with INLALogDensity for INLA-within-MCMC (ADR-009) — each NUTS / HMC step on θ is followed by a draw from the conditional surface at that θ.

Algorithm

1. Run laplace_mode at θ to get mode x̂ and posterior precision H.
2. Draw z ∼ N(0, I) and solve L'ᵀ x₀ = z where L is the sparse Cholesky factor of H (so Cov(x₀) = H⁻¹).
3. If the component stack declared hard linear constraints Cx = e, subtract the kriging correction so the draw satisfies the constraint to working precision.

Return value

• Single-sample call: Vector{Float64} of length n_latent(model).
• n_samples positional argument: Matrix{Float64} of size (n_latent(model), n_samples).

Example

rng = Random.Xoshiro(42)
θ = initial_hyperparameters(model)
x = sample_conditional(model, θ, y; rng)
X = sample_conditional(model, θ, y, 100; rng)

sample_y(rng, ℓ::AbstractLikelihood, η, θ) -> Vector{Float64}

Draw a posterior-predictive replicate y_rep ∼ p(y | η, θ) for each row in η. The returned vector has length length(η) and matches the likelihood's own observation rows — its block, in a multi-likelihood model. Per-observation metadata that lives on the likelihood (Binomial n_trials, Poisson/NegBin offset E) is read directly from ℓ, so the caller does not pass it explicitly.

Used by posterior_predictive_y for posterior-predictive checks. The default raises ArgumentError; concrete likelihoods that admit closed-form sampling (Gaussian, Poisson, Binomial, NegativeBinomial, Gamma) override.

user_scale(prior::AbstractHyperPrior, θ) -> Real

Map internal-scale θ to the user-facing parameter (e.g. exp(θ) for a log-precision prior).

validate_censoring(censoring, time_hi, y)

Public-boundary validation for survival-likelihood inputs. Asserts that censoring and time_hi (when non-nothing) match y in length, and that every INTERVAL row has time_hi[i] > y[i].

Called once per fit; assertions are not in the inner Newton loop.

waic(rng, res, model, y; n_samples = 1000)
  -> @NamedTuple{WAIC, lpd, pWAIC, elpd_WAIC}

Watanabe's widely-applicable information criterion (WAIC-2):

elpd_WAIC = Σ_i (lpd_i - pWAIC_i)
WAIC      = -2 · elpd_WAIC

with

lpd_i   = log E[ p(y_i | η_i, θ_ℓ) ]
pWAIC_i = Var[ log p(y_i | η_i, θ_ℓ) ]

Expectations are Monte-Carlo under the posterior samples.

∂inverse_link(g::AbstractLinkFunction, η)

Derivative dμ/dη.

∂²inverse_link(g::AbstractLinkFunction, η)

Second derivative d²μ/dη².

∇_η_log_density(ℓ::AbstractLikelihood, y, η, θ) -> AbstractVector

Gradient of log p(y | η, θ) w.r.t. η.

∇²_η_log_density(ℓ::AbstractLikelihood, y, η, θ) -> AbstractVector

Diagonal of the Hessian of log p(y | η, θ) w.r.t. η. Length equals length(y).

∇³_η_log_density(ℓ::AbstractLikelihood, y, η, θ) -> AbstractVector

Diagonal of the third derivative tensor of log p(y | η, θ) w.r.t. η. Length equals length(y). The likelihood factorises across observations, so only the diagonal entry ∂³ log p(y_i|η_i) / ∂η_i^3 is non-zero.

Needed for Rue-Martino-Chopin (2009) simplified Laplace / skew correction of posterior marginals π(x_i | y). Concrete likelihoods should override with closed forms; the default falls back on a central finite-difference of ∇²_η_log_density per coordinate.