Spatial — Meuse zinc SPDE

The Meuse-zinc dataset (Pebesma, sp and gstat) is the canonical geostatistics benchmark: top-soil zinc concentration at 155 sample points along the Meuse river, with distance to the river as the only covariate. The model is Gaussian-on-log-zinc with an SPDE–Matérn spatial random field on a triangulated mesh:

\[\begin{aligned} \log y_i &= \alpha + \beta\,\mathrm{dist}_i + u(s_i) + \varepsilon_i, \\ u &\sim \text{SPDE-Matérn}(\rho, \sigma; \alpha = 2), \\ \varepsilon_i &\sim \mathcal{N}(0, 1/\tau_\varepsilon), \end{aligned}\]

with PC priors on (ρ, σ) (Fuglstad et al. 2019) and a PC prior on τ_ε.

Loading the fixture

The fixture ships the input data, a fmesher-generated mesh (points, triangle vertices tv), the spatial projector A_field, and the R-INLA posterior summary. Pinning the mesh isolates the oracle from the M3 mesh-parity gate so the test exercises the LGM + SPDE + INLA path only.

using JLD2, SparseArrays, LinearAlgebra
using GMRFs, LatentGaussianModels, INLASPDE

const FIXTURE = joinpath(@__DIR__, "..", "..", "..",
    "packages", "INLASPDE.jl",
    "test", "oracle", "fixtures", "meuse_spde.jld2")

fxt = load(FIXTURE)["fixture"]

y       = Float64.(fxt["input"]["y"])
dist    = Float64.(fxt["input"]["dist"])
points  = fxt["mesh"]["loc"]::Matrix{Float64}
tv      = fxt["mesh"]["tv"]::Matrix{Int}
A_field = SparseMatrixCSC{Float64, Int}(fxt["A_field"])

n_obs = length(y)
n_v   = size(points, 1)

Building the model

Latent vector layout: x = [α; β; u], where u ∈ ℝ^{n_v} is the mesh-vertex SPDE field. The projector stacks the intercept column, the covariate column, and the fmesher-built barycentric A_field:

spde = SPDE2(points, tv; α = 2,
    pc = PCMatern(
        range_U = 0.5, range_α = 0.5,
        sigma_U = 1.0, sigma_α = 0.5,
    ))

c_int  = Intercept(prec = 1.0e-3)
c_dist = FixedEffects(1; prec = 1.0e-3)

A = hcat(ones(n_obs, 1), reshape(dist, n_obs, 1), A_field)

ℓ     = GaussianLikelihood(hyperprior = PCPrecision(1.0, 0.01))
model = LatentGaussianModel(ℓ, (c_int, c_dist, spde), A)

Fitting

res = inla(model, y)

Posterior summaries

spde_user_scale maps the internal log-scale SPDE hyperparameters back to interpretable (range, σ):

fixed_effects(model, res)

2-element Vector{@NamedTuple{name::String, mean::Float64, sd::Float64, lower::Float64, upper::Float64}}:
 (name = "Intercept[1]", mean = 6.616720890179365, sd = 0.1560711371815215, lower = 6.310827082030758, upper = 6.922614698327972)
 (name = "FixedEffects[2]", mean = -2.8234421079650867, sd = 0.3996178018777923, lower = -3.6066786078580746, upper = -2.0402056080720987)

θ̂ = res.θ̂
ρ̂, σ̂ = spde_user_scale(spde, θ̂[2:3])
τ̂_noise = exp(θ̂[1])
(range = ρ̂, sigma = σ̂, tau_noise = τ̂_noise)

(range = 0.5786891628419261, sigma = 0.45266851597075974, tau_noise = 15.389168886803379)

Comparing to R-INLA

The fixture's summary_fixed, summary_hyperpar, and mlik carry R-INLA's posterior. Per plans/testing-strategy.md, the Meuse oracle tolerances are 1% on fixed-effect means relative to their posterior SD, 25% on (ρ, σ), and 30% on τ_ε (looser bands reflect mode-vs-mean and integration-scheme differences). The full assertion suite lives in test/oracle/test_meuse_spde.jl.

sf_rows = fxt["summary_fixed"]["rownames"]
sf_mean = Float64.(fxt["summary_fixed"]["mean"])
α_R = sf_mean[findfirst(==("intercept"), sf_rows)]
β_R = sf_mean[findfirst(==("dist"), sf_rows)]

fe = fixed_effects(model, res)
(α_julia = fe[1].mean, α_R = α_R, β_julia = fe[2].mean, β_R = β_R)

(α_julia = 6.616720890179365, α_R = 6.6193211692206635, β_julia = -2.8234421079650867, β_R = -2.8116029235543802)

Mesh parity (M3)

The mesh used here was generated by R-INLA's fmesher and shipped with the fixture, so this vignette is decoupled from our own mesher. For a full Julia-only run, swap points, tv, and the projector for the output of inla_mesh_2d and MeshProjector. The fmesher_parity test suite tracks vertex counts and minimum-angle bounds against R-INLA's fmesher outputs on the unit square, L-shape, and Meuse hull.

Predicting on a raster grid

The fit above lives on the mesh; users typically want the spatial field interpolated onto a regular pixel grid for visualisation, overlay against rasterised covariates, or downstream geostatistical workflows. That is exactly what INLASPDERasters.jl provides via predict_raster(model, res, template).

The raster path requires an SPDE2 component constructed from a retained INLAMesh (per ADR-036), so we rebuild the model with a Julia-native mesh and MeshProjector rather than the fixture's (points, tv) + A_field scaffolding:

using INLASPDERasters
using Rasters: Raster, X, Y
using Random: Xoshiro

locs = fxt["input"]["locations"]::Matrix{Float64}

mesh_jl = inla_mesh_2d(locs;
    max_edge = (0.4, 0.8), cutoff = 0.05, offset = (0.2, 0.5))

spde_jl = SPDE2(mesh_jl; α = 2,
    pc = PCMatern(
        range_U = 0.5, range_α = 0.5,
        sigma_U = 1.0, sigma_α = 0.5,
    ))

P_obs = MeshProjector(mesh_jl, locs)
A_field_jl = SparseMatrixCSC{Float64, Int}(P_obs.A)
A_jl = hcat(ones(n_obs, 1), reshape(dist, n_obs, 1), A_field_jl)

model_jl = LatentGaussianModel(ℓ, (c_int, c_dist, spde_jl), A_jl)
res_jl = inla(model_jl, y)

The template raster fixes the output extent, resolution, and CRS; predict_raster returns a Raster matching it cell-for-cell, with NaN outside the mesh:

xs = 178.5:0.05:181.5
ys = 329.5:0.05:333.7
template = Raster(zeros(length(xs), length(ys)),
    (X(collect(xs)), Y(collect(ys))))

r_mean = predict_raster(model_jl, res_jl, template;
    component = SPDE2, quantity = :mean)
r_lower = predict_raster(model_jl, res_jl, template;
    component = SPDE2, quantity = :lower)
r_upper = predict_raster(model_jl, res_jl, template;
    component = SPDE2, quantity = :upper)

interior(r) = filter(!isnan, parent(r))
(
    grid_size = size(r_mean),
    n_cells_inside_mesh = count(!isnan, parent(r_mean)),
    mean_extrema  = extrema(interior(r_mean)),
    lower_extrema = extrema(interior(r_lower)),
    upper_extrema = extrema(interior(r_upper)),
)

(grid_size = (61, 85), n_cells_inside_mesh = 4395, mean_extrema = (-0.9464622147131828, 0.9251520703262933), lower_extrema = (-1.4255077830479734, 0.38206214048406323), upper_extrema = (-0.4674166463783921, 1.4682420001685235))

r_mean, r_lower, r_upper are projections of the per-vertex posterior summaries reported by random_effects(model_jl, res_jl) onto the cell centres of template. They share the template's extent, resolution, and CRS — ready for direct plotting or overlay against a rasterised covariate.

For tail probabilities — for example, "the probability that the SPDE residual exceeds 0.5 log-mg/kg" — the Gaussian-approximation overload cannot help, since a per-cell P(η > c) is not a linear functional of the vertex marginals. The sample-based overload threads draws from the joint posterior through the same projector and reduces in-cell:

rng = Xoshiro(20260506)
r_exc = predict_raster(rng, model_jl, res_jl, template;
    component = SPDE2,
    quantity = Exceedance(0.5),
    n_samples = 500)
(
    n_cells_inside_mesh = count(!isnan, parent(r_exc)),
    exceedance_extrema  = extrema(interior(r_exc)),
)

(n_cells_inside_mesh = 4395, exceedance_extrema = (0.0, 0.928))

Each cell of r_exc is the empirical fraction of posterior draws in which the SPDE field at that cell exceeds 0.5. The Exceedance wrapper distinguishes "value to threshold against" from "credible-level quantile to compute" — passing a Real in [0, 1] instead computes the per-cell quantile of the sampled field.

Coming from R-INLA: the @lgm macro

The same fit expresses cleanly through LGMFormula.jl's @lgm macro. Coordinate-indexed terms f((east, north), spde) invoke the MeshProjector for you and place the resulting block in the design matrix at the position implied by formula order — the explicit constructor above and the macro both produce the same LatentGaussianModel:

using LGMFormula

df = (
    logzinc = y,
    dist    = dist,
    east    = locs[:, 1],
    north   = locs[:, 2],
)

model_macro = @lgm logzinc ~ 1 + dist + f((east, north), spde_jl) data=df family=GaussianLikelihood(hyperprior=PCPrecision(1.0, 0.01))
(
    A_size            = size(model_macro.mapping.A),
    components        = nameof.(typeof.(model_macro.components)),
    matches_explicit  = size(model_macro.mapping.A) == size(model_jl.mapping.A),
)

(A_size = (155, 251), components = (:Intercept, :FixedEffects, :SPDE2), matches_explicit = true)

The tuple-coordinate form requires INLASPDE to be loaded so the LGMFormulaINLASPDEExt weakdep extension is active; that is satisfied here by the using INLASPDE at the top of the vignette. See the migration guide for the full range of forms @lgm accepts.