Latent Models

Latent models provide a structured way to define commonly used Gaussian Markov Random Fields (GMRFs) by specifying their key components: hyperparameters, precision matrices, mean vectors, and constraints. This abstraction enables easy construction and combination of standard models used in spatial statistics, time series analysis, and Bayesian modeling.

Overview

The LatentModel interface standardizes GMRF construction through a simple pattern:

1. Define the model structure (e.g., temporal, spatial, independence)

2. Specify hyperparameters and their validation

3. Construct GMRFs with automatic constraint handling

using GaussianMarkovRandomFields

# Define a temporal AR1 model
ar1 = AR1Model(100)
hyperparameters(ar1)  # (τ = Real, ρ = Real)

# Construct GMRF with specific parameters
gmrf = ar1(τ=2.0, ρ=0.8)  # Returns GMRF or ConstrainedGMRF automatically

LatentModel Interface

GaussianMarkovRandomFields.LatentModel Type

LatentModel

Abstract type for latent variable models that can be used to construct GMRFs.

A LatentModel provides a structured way to define commonly used GMRFs such as AR1, RW1, and other temporal/spatial models by specifying:

1. The hyperparameters of the model

2. How to construct the precision matrix from hyperparameters

3. How to construct the mean vector from hyperparameters

4. Any linear constraints that should be applied

Interface

Each concrete subtype must implement:

• length(model): Return the size/dimension of the latent process

• hyperparameters(model): Return a NamedTuple describing the hyperparameters

• precision_matrix(model; kwargs...): Construct precision matrix from hyperparameter values

• mean(model; kwargs...): Construct mean vector from hyperparameter values

• constraints(model; kwargs...): Return constraint information or nothing

• model_name(model): Return a Symbol representing the preferred name for this model type

• (model)(; kwargs...): Instantiate a concrete GMRF from hyperparameter values

Usage

# Define a model
model = SomeLatentModel(n=100)

# Get hyperparameter specification
params = hyperparameters(model)  # e.g. (τ=Real, ρ=Real)

# Instantiate GMRF with specific parameter values
gmrf = model(τ=2.0, ρ=0.8)  # Returns GMRF or ConstrainedGMRF

source

GaussianMarkovRandomFields.hyperparameters Method

hyperparameters(model::LatentModel)

Return a NamedTuple describing the hyperparameters and their types for the model.

Returns

A NamedTuple where keys are parameter names and values are their expected types.

source

GaussianMarkovRandomFields.precision_matrix Method

precision_matrix(model::LatentModel; kwargs...)

Construct the precision matrix for the model given hyperparameter values.

Arguments

• model: The LatentModel instance

• kwargs...: Hyperparameter values as keyword arguments

Returns

A precision matrix (AbstractMatrix or LinearMap) for use in GMRF construction.

source

Statistics.mean Function

mean(model::LatentModel; kwargs...)

Construct the mean vector for the model given hyperparameter values.

Arguments

• model: The LatentModel instance

• kwargs...: Hyperparameter values as keyword arguments

Returns

A mean vector (AbstractVector) for use in GMRF construction.

source

GaussianMarkovRandomFields.constraints Function

constraints(model::LatentModel; kwargs...)

Return constraint information for the model given hyperparameter values.

Arguments

• model: The LatentModel instance

• kwargs...: Hyperparameter values as keyword arguments

Returns

Either nothing if no constraints, or a tuple (A, e) where A is the constraint matrix and e is the constraint vector such that Ax = e.

source

GaussianMarkovRandomFields.model_name Function

model_name(model::LatentModel)

Return a Symbol representing the preferred name for this model type.

This name is used for parameter prefixing in CombinedModel to avoid conflicts. For example, if two models both have a τ parameter, they become τ_ar1, τ_besag, etc.

Returns

A Symbol that will be used as the suffix in parameter names.

source

Available Models

Temporal Models

GaussianMarkovRandomFields.AR1Model Type

AR1Model(n::Int; alg=LDLtFactorization(), constraint=nothing)

A first-order autoregressive (AR1) latent model for constructing AR1 GMRFs.

The AR1 model represents a temporal process where each observation depends on the previous observation with some correlation ρ and precision τ.

Mathematical Description

For n observations, the AR1 process has:

• Zero mean: μ = 0

• Precision matrix Q with tridiagonal structure:

  • Q[1,1] = τ

  • Q[i,i] = (1 + ρ²)τ for i = 2,...,n-1

  • Q[n,n] = τ

  • Q[i,i+1] = Q[i+1,i] = -ρτ for i = 1,...,n-1

Hyperparameters

• τ: Precision parameter (τ > 0)

• ρ: Correlation parameter (|ρ| < 1)

Fields

• n::Int: Length of the AR1 process

• alg::Alg: LinearSolve algorithm for solving linear systems

• constraint::C: Optional constraint, either nothing or (A, e) tuple

Example

model = AR1Model(100)
gmrf = model(τ=2.0, ρ=0.8)  # Construct unconstrained AR1 GMRF

# With sum-to-zero constraint
model = AR1Model(100, constraint=:sumtozero)
gmrf = model(τ=2.0, ρ=0.8)  # Returns ConstrainedGMRF

# With custom constraint
A = [1.0 1.0 zeros(98)...]
e = [0.0]
model = AR1Model(100, constraint=(A, e))
gmrf = model(τ=2.0, ρ=0.8)

source

GaussianMarkovRandomFields.RW1Model Type

RW1Model(n::Int; regularization=1e-5, alg=LDLtFactorization(), additional_constraints=nothing)

A first-order random walk (RW1) latent model for constructing intrinsic GMRFs.

The RW1 model represents a non-stationary temporal process where each observation is the previous observation plus Gaussian noise. This creates a smooth trend model that's popular for temporal smoothing and time-varying effects.

Mathematical Description

The RW1 process defines increments: x[i+1] - x[i] ~ N(0, τ⁻¹) for i = 1,...,n-1. This leads to a singular precision matrix with the tridiagonal structure:

• Q[1,1] = 1, Q[n,n] = 1

• Q[i,i] = 2 for i = 2,...,n-1

• Q[i,i+1] = Q[i+1,i] = -1 for i = 1,...,n-1

Since this matrix is singular (rank n-1), we handle it as an intrinsic GMRF by:

1. Scaling by τ first, then adding small regularization (1e-5) to diagonal for numerical stability

2. Always adding sum-to-zero constraint: sum(x) = 0 (required for identifiability)

Hyperparameters

• τ: Precision parameter (τ > 0)

Fields

• n::Int: Length of the RW1 process

• regularization::Float64: Small value added to diagonal after scaling (default 1e-5)

• alg::Alg: LinearSolve algorithm for solving linear systems

• additional_constraints::C: Optional additional constraints beyond the required sum-to-zero

Example

model = RW1Model(100)
gmrf = model(τ=1.0)  # Returns ConstrainedGMRF with sum-to-zero constraint

# Add additional constraint (e.g., fix first element to 0)
A_add = reshape([1.0; zeros(99)], 1, 100)
e_add = [0.0]
model = RW1Model(100, additional_constraints=(A_add, e_add))
gmrf = model(τ=1.0)  # Returns ConstrainedGMRF with both sum-to-zero AND additional constraint

source

Spatial Models

GaussianMarkovRandomFields.MaternModel Type

MaternModel{F<:FEMDiscretization,S<:Integer,Alg}(discretization::F, smoothness::S, alg::Alg)

A Matérn latent model for constructing spatial GMRFs from discretized Matérn SPDEs.

The MaternModel provides a structured way to define Matérn Gaussian Markov Random Fields by discretizing the Matérn SPDE using finite element methods. The model stores the discretization and smoothness parameter, while the range parameter is provided at GMRF construction time.

Mathematical Description

The Matérn SPDE is given by: (κ² - Δ)^(α/2) u(x) = W(x), where α = ν + d/2

This leads to a Matérn covariance function with range and smoothness parameters.

Two Construction Modes

1. Direct construction: Pass a pre-built FEMDiscretization

2. Automatic construction: Pass points, automatically creates convex hull mesh

Hyperparameters

• range: Range parameter (range > 0) - controls spatial correlation distance

Fields

• discretization::F: The finite element discretization

• smoothness::S: The smoothness parameter (Integer, controls differentiability)

• alg::Alg: LinearSolve algorithm for solving linear systems

Examples

# Direct construction
disc = FEMDiscretization(grid, interpolation, quadrature)
model = MaternModel(disc; smoothness = 2)
gmrf = model(range=2.0)

# Automatic construction from points
points = [0.0 0.0; 1.0 0.0; 0.5 1.0]  # N×2 matrix
model = MaternModel(points; smoothness = 1, element_order = 1)
gmrf = model(range=2.0)

# With custom algorithm
model = MaternModel(disc; smoothness = 2, alg = LDLtFactorization())
gmrf = model(range=2.0)

source

GaussianMarkovRandomFields.BesagModel Type

BesagModel(adjacency::AbstractMatrix; regularization::Float64 = 1e-5, normalize_var::Bool = true, singleton_policy::Symbol = :gaussian, alg=CHOLMODFactorization())

A Besag model for spatial latent effects on graphs using intrinsic Conditional Autoregressive (CAR) structure.

The Besag model represents spatial dependence where each node's precision depends on its graph neighbors. This creates a graph Laplacian precision structure that's widely used for spatial smoothing on irregular lattices.

Mathematical Description

For a graph with adjacency matrix W, the precision matrix follows:

• Q[i,j] = -τ if nodes i and j are neighbors (W[i,j] = 1)

• Q[i,i] = τ * degree[i] where degree[i] = sum(W[i,:])

• All other entries are 0

Since this matrix is singular (rank n-1), we handle it as an intrinsic GMRF by:

1. Scaling by τ first, then adding small regularization (1e-5) to diagonal for numerical stability

2. Adding sum-to-zero constraint: sum(x) = 0

Hyperparameters

• τ: Precision parameter (τ > 0)

Fields

• adjacency::M: Adjacency matrix W (preserves input structure - sparse, SymTridiagonal, etc.)

• regularization::Float64: Small value added to diagonal after scaling (default 1e-5)

• alg::Alg: LinearSolve algorithm for solving linear systems

Example

# 4-node cycle graph - can use sparse, SymTridiagonal, or Matrix
W = sparse(Bool[0 1 0 1; 1 0 1 0; 0 1 0 1; 1 0 1 0])
model = BesagModel(W)
gmrf = model(τ=1.0)  # Returns ConstrainedGMRF with sum-to-zero constraint using CHOLMODFactorization

# Or specify custom algorithm
model = BesagModel(W, alg=LDLtFactorization())
gmrf = model(τ=1.0)

source

Independence Models

GaussianMarkovRandomFields.IIDModel Type

IIDModel(n::Int; alg=DiagonalFactorization(), constraint=nothing)

An independent and identically distributed (IID) latent model for constructing simple diagonal GMRFs.

The IID model represents independent Gaussian random variables with identical precision τ. This is the simplest possible latent model, equivalent to a scaled identity precision matrix.

Mathematical Description

Each element is independent: x[i] ~ N(0, τ⁻¹) for i = 1,...,n. The precision matrix is simply: Q = τ * I(n)

This model is useful for:

• Modeling independent effects or noise

• Baseline comparisons with structured models

• Teaching/demonstration purposes

Hyperparameters

• τ: Precision parameter (τ > 0)

Fields

• n::Int: Length of the IID process

• alg::Alg: LinearSolve algorithm for solving linear systems

• constraint::C: Optional constraint, either nothing or (A, e) tuple

Example

model = IIDModel(100)
gmrf = model(τ=2.0)  # Returns unconstrained GMRF

# With sum-to-zero constraint (common in INLA for separating global offset)
model = IIDModel(100, constraint=:sumtozero)
gmrf = model(τ=2.0)  # Returns ConstrainedGMRF with sum-to-zero constraint

source

GaussianMarkovRandomFields.FixedEffectsModel Type

FixedEffectsModel(n::Int; λ::Real = 1e-6, alg=DiagonalFactorization())

Weakly regularized fixed-effects latent component.

Encodes standard GLM-style fixed effects inside the latent vector with a small ridge precision λ * I(n).

• No hyperparameters (returns NamedTuple())

• No constraints

• Name: :fixed (used for parameter prefixing in CombinedModel)

Fields

• n::Int: Length of the fixed effects vector

• λ::Float64: Ridge regularization parameter

• alg::Alg: LinearSolve algorithm for solving linear systems

Example

model = FixedEffectsModel(10)
gmrf = model()  # Returns GMRF with precision λ * I(10) using DiagonalFactorization

# Or specify custom algorithm
model = FixedEffectsModel(10, alg=CHOLMODFactorization())
gmrf = model()

source

Model Composition

GaussianMarkovRandomFields.CombinedModel Type

CombinedModel(components::Vector{<:LatentModel}; alg=CHOLMODFactorization())

A combination of multiple LatentModel instances into a single block-structured GMRF.

This enables modeling with multiple latent components, such as the popular BYM model (Besag-York-Mollié) which combines spatial (Besag) and independent (IID) effects.

Mathematical Description

Given k component models with sizes n₁, n₂, ..., nₖ:

• Combined precision matrix: Q = blockdiag(Q₁, Q₂, ..., Qₖ)

• Combined mean vector: μ = vcat(μ₁, μ₂, ..., μₖ)

• Combined constraints: Block-diagonal constraint structure preserving individual constraints

Parameter Naming

To avoid conflicts when multiple models have the same hyperparameters (e.g., multiple τ), parameters are automatically prefixed with model names:

• Single occurrence: τ_besag, ρ_ar1

• Multiple occurrences: τ_besag, τ_besag_2, τ_besag_3

Fields

• components::Vector{<:LatentModel}: The individual latent models

• component_sizes::Vector{Int}: Cached sizes of each component

• total_size::Int: Total size of the combined model

• alg::Alg: LinearSolve algorithm for solving linear systems

Example - BYM Model

# BYM model: spatial Besag + independent IID effects
W = sparse(adjacency_matrix)  # Spatial adjacency
besag = BesagModel(W)         # Spatial component
iid = IIDModel(n)            # Independent component

# Vector constructor
bym = CombinedModel([besag, iid])
# Or variadic constructor (syntactic sugar)
bym = CombinedModel(besag, iid)

# With custom algorithm
bym = CombinedModel([besag, iid], alg=LDLtFactorization())

# Usage with automatically prefixed parameters
gmrf = bym(τ_besag=1.0, τ_iid=2.0)

source

Usage Examples

Basic Model Construction

# Temporal AR1 process
ar1 = AR1Model(100)
gmrf = ar1(τ=2.0, ρ=0.8)

# Spatial Matérn model from points
points = [0.0 0.0; 1.0 0.0; 0.5 1.0]  # N×2 matrix
matern = MaternModel(points; smoothness = 2)
gmrf = matern(range=1.5)

# Spatial Besag model
W = sparse_adjacency_matrix
besag = BesagModel(W)
gmrf = besag(τ=1.0)  # Returns ConstrainedGMRF with sum-to-zero constraint

Model Composition

# Classic BYM model: spatial + independent effects  
W = spatial_adjacency_matrix
bym = CombinedModel(BesagModel(W), IIDModel(n))

# Check combined hyperparameters
hyperparameters(bym)  # (τ_besag = Real, τ_iid = Real)

# Construct combined GMRF
gmrf = bym(τ_besag=1.0, τ_iid=2.0)

# Continuous spatial field + independent effects
points = generate_observation_points()
spatial_matern = MaternModel(points; smoothness = 1) 
independent_effects = IIDModel(length(points))
combined = CombinedModel(spatial_matern, independent_effects)
gmrf = combined(range=2.0, τ_iid=0.1)

Smart Parameter Naming

CombinedModel automatically handles parameter naming conflicts:

# Single occurrence: no suffix
CombinedModel(AR1Model(10), BesagModel(W))  # τ_ar1, ρ_ar1, τ_besag

# Multiple occurrences: numbered suffixes  
CombinedModel(IIDModel(5), IIDModel(10), IIDModel(15))  # τ_iid, τ_iid_2, τ_iid_3

Advanced Compositions

# Spatiotemporal model with three components
W = spatial_adjacency(n_regions)
model = CombinedModel(
    BesagModel(W),           # Spatial structure
    RW1Model(n_time),        # Temporal trend  
    IIDModel(n_regions * n_time)  # Independent effects
)

gmrf = model(τ_besag=1.0, τ_rw1=2.0, τ_iid=0.1)

Key Features

Automatic Constraint Handling

Models automatically determine GMRF type based on constraints:

ar1 = AR1Model(10)     # No constraints → GMRF
rw1 = RW1Model(10)     # Sum-to-zero constraint → ConstrainedGMRF
mixed = CombinedModel(ar1, rw1)  # Inherits constraints → ConstrainedGMRF

Custom Constraints

Latent models support adding custom linear equality constraints Ax = e. This is particularly useful for INLA-style modeling where you want to learn a global offset separately from the latent field.

Models with Optional Constraints

For models without built-in constraints (AR1, IID, FixedEffects, Matern), use the constraint parameter:

# Sum-to-zero constraint (shorthand)
ar1 = AR1Model(100, constraint=:sumtozero)
gmrf = ar1(τ=1.0, ρ=0.8)  # Returns ConstrainedGMRF

# Custom constraint: fix first two elements to sum to 1
A = [1.0 1.0 0.0 0.0 0.0]  # 1×n constraint matrix
e = [1.0]                   # Constraint value
iid = IIDModel(5, constraint=(A, e))
gmrf = iid(τ=2.0)          # x[1] + x[2] = 1.0

Intrinsic Models with Additional Constraints

Intrinsic models (RW1, Besag) already have built-in constraints. Use additional_constraints to add extra constraints without removing the built-in ones:

# RW1 already has sum-to-zero; add constraint to fix first element
A_extra = [1.0 0.0 0.0 0.0 0.0]  # Fix x[1] = 0
e_extra = [0.0]
rw1 = RW1Model(5, additional_constraints=(A_extra, e_extra))
gmrf = rw1(τ=1.0)  # Both sum(x) = 0 AND x[1] = 0

# Besag with additional constraint
besag = BesagModel(W, additional_constraints=(A_extra, e_extra))
gmrf = besag(τ=1.0)  # Component-wise sum-to-zero + extra constraint

Don't Remove Required Constraints

Intrinsic models like RW1 and Besag require their built-in constraints for identifiability. Attempting to use constraint=:sumtozero on RW1Model will throw an error since it already has this constraint. Always use additional_constraints for these models.

Parameter Validation

All models include built-in parameter validation:

ar1 = AR1Model(100)
ar1(τ=-1.0, ρ=0.8)   # ArgumentError: τ must be positive
ar1(τ=1.0, ρ=1.5)    # ArgumentError: |ρ| must be < 1

Efficient Sparse Operations

• Block-diagonal precision matrices are always sparse

• Input matrix structures (sparse, SymTridiagonal) are preserved

• Minimal memory allocations through size caching

See Also

• GMRF and ConstrainedGMRF for the underlying GMRF types

• Hard Constraints for constraint handling details

• Observation Models for linking GMRFs to data

Formula Terms

Prefer writing models with formulas? The latent components discussed here can be constructed via formula terms and assembled automatically into a combined model and design matrix. See the Formula Interface reference for details:

• Formula Interface — terms (IID, RandomWalk, AR1, Besag) and build_formula_components.

• For a worked example combining Besag + IID + fixed effects under a Poisson likelihood with offset, see the tutorial: Advanced GMRF modelling for disease mapping