Skip to content

meuse-spde.md

Latest commit

 

History

History
255 lines (204 loc) · 8.31 KB

meuse-spde.md

File metadata and controls

255 lines (204 loc) · 8.31 KB
Raw

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)
nothing # hide

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)
nothing # hide

Fitting

res = inla(model, y)
nothing # hide

Posterior summaries

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

fixed_effects(model, res)

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

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)

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)
nothing # hide

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)),
)

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)),
)

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),
)

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.