Modeling SPDE - pyINLA

Overview

This page walks you through fitting an SPDE model in pyINLA, from mesh creation to extracting results. For the theory behind SPDEs (Matérn covariance, FEM matrices, precision matrix Q), see Understanding SPDEs.

# Quick import of all SPDE functions
from fmesher import (
    fm_mesh_2d,       # Create triangular mesh
    fm_fem,           # Compute FEM matrices (C, G)
    spde2_pcmatern,   # Create SPDE model with PC priors
    spde_make_A,      # Build projector matrix
)

The SPDE Workflow

Here's the step-by-step process for fitting an SPDE model with pyINLA:

Step 1: Create Mesh

Build a triangular mesh covering your spatial domain using the fmesher package.

# Using the fmesher package (wraps R's fmesher)
from fmesher import fm_mesh_2d
import numpy as np

coords = np.column_stack([lon, lat])  # Your observation locations
mesh = fm_mesh_2d(
    loc=coords,
    max_edge=[0.5, 1.0],   # Max edge length [inner, outer]
    offset=[0.1, 0.3],      # Extension distance
    cutoff=0.05             # Min distance between points
)

# mesh.n = number of vertices
# mesh.n_triangle = number of triangles
# mesh.loc = vertex coordinates
# mesh.tv = triangle-vertex indices

Step 2: Create SPDE Model

Create the SPDE model object with PC priors on range and sigma. The FEM matrices (C, G) are computed internally using fm_fem.

# Using the fmesher package (pure Python)
from fmesher import spde2_pcmatern

spde = spde2_pcmatern(
    mesh=mesh,
    alpha=2,                      # Smoothness: alpha = nu + d/2
    prior_range=(0.3, 0.5),       # P(range < 0.3) = 0.5
    prior_sigma=(1.0, 0.01)       # P(sigma > 1.0) = 0.01
)

# The SPDE object contains:
# - spde.n_spde: number of mesh vertices
# - spde.fem: FEM matrices (c0, c1, g1, g2)
# - spde.precision(theta): compute Q for given [log_range, log_sigma]
# - spde.log_prior(theta): compute PC prior log-density

PC Priors: The Penalized Complexity prior is set via two conditions: prior_range=[r0, p] means P(range < r0) = p, and prior_sigma=[s0, p] means P(sigma > s0) = p. Choose r0 based on your domain size and s0 based on expected variability.

Step 3: Build Projector Matrix

Link mesh vertices to observation locations using barycentric interpolation.

# Using the fmesher package (pure Python)
from fmesher import spde_make_A

A = spde_make_A(mesh=mesh, loc=coords)
# A is a sparse matrix: n_obs x n_vertices
# Each row has at most 3 non-zero values (barycentric coordinates)
# Each row sums to 1

Step 4: Create Data Stack

Use inla.stack() to organize data, projector matrices, and effects together. This handles the different dimensions (n observations vs m mesh vertices).

# Using R's inla.stack via rpy2
from pyinla.fmesher_rpy2 import inla_stack

n_obs = len(response)
stk = inla_stack(
    data={'y': response},           # Response variable
    A=[A, 1],                       # Projector matrices: A for spatial, 1 for intercept
    effects=[
        {'i': range(1, spde.n + 1)},  # Spatial effect indices
        {'beta0': [1] * n_obs}        # Intercept (vector of 1s)
    ],
    tag='est'                       # Label for this stack
)

Step 5: Fit Model

Call INLA with the SPDE random effect. The formula removes the default intercept and adds it as a covariate so all terms use projector matrices.

# Using R-INLA via rpy2
from pyinla.fmesher_rpy2 import inla_fit, inla_stack_data, inla_stack_A

result = inla_fit(
    formula='y ~ 0 + beta0 + f(i, model=spde)',
    data=inla_stack_data(stk),
    control_predictor={'A': inla_stack_A(stk)}
)

# Access results
print(result.summary_fixed)        # Intercept posterior
print(result.summary_hyperpar)     # Precision, range, sigma
print(result.summary_random['i'])  # Spatial field at vertices

          Key Functions:
          inla.stack() - Combines data, projector matrices, and effects
inla.stack.data() - Extracts data from stack for inla()
inla.stack.A() - Extracts combined projector matrix
f(i, model=spde) - Specifies the SPDE random effect in the formula

        

What You Get

After fitting, you can extract these results from the pyINLA result object:

Result	pyINLA Access	Description
Fixed effects	`result.summary_fixed`	Posterior of \(\beta_0\) and other covariates (mean, sd, quantiles)
Spatial field	`result.summary_random['i']`	Posterior mean/sd at each mesh vertex
Hyperparameters	`result.summary_hyperpar`	DataFrame with range, sigma, nugget posteriors
Marginal likelihood	`result.mlik`	Log marginal likelihood for model comparison
Fitted values	`result.summary_linear_predictor`	Posterior of linear predictor at each observation

Example: Extracting Results

# Fixed effects (intercept, covariates)
print(result.summary_fixed)
#          mean       sd  0.025quant  0.5quant  0.975quant
# beta0   1.234    0.052       1.132     1.234       1.336

# Spatial random effects at mesh vertices
spatial_mean = result.summary_random['i']['mean']
spatial_sd = result.summary_random['i']['sd']

# Hyperparameters (range, sigma, nugget precision)
print(result.summary_hyperpar)
# Range and sigma are in practical units if PC priors were used

# Project spatial field to a prediction grid
A_pred = spde_make_A(mesh=mesh, loc=pred_coords)
field_pred = A_pred @ spatial_mean

# Compute precision matrix Q for given hyperparameters
theta = [np.log(range_est), np.log(sigma_est)]
Q = spde.precision(theta)  # Sparse precision matrix

Next Steps

Understanding SPDEs - Learn the theory behind SPDE models
Creating Meshes - Detailed guide on mesh construction
Priors - Setting PC priors for range and sigma