PaperlayerView exact revision on arXiv

Source-linked AI summary

Non-stationary Gaussian models with physical barriers

Haakon Bakka, Jarno Vanhatalo, Janine Illian, Daniel Simpson, Håvard Rue

arXiv:1608.03787v3stat.APstat.ME

TL;DR

Stationary spatial fields can smooth across physical barriers, and small changes in boundary polygons can sharply alter shortest-distance methods. The paper proposes the Barrier model, which reduces correlation across land; it has the same theoretical computational cost as the stationary model and changes the Archipelago application’s spatial estimates.

  • Problem

    Stationary fields do not account for physical barriers, and shortest-distance approaches can be sensitive to how the boundary polygon is drawn.

    Source: §1.1 · p. 2, §1.3 · p. 2

  • Method

    The Barrier SGF uses a Matérn field with a near-zero range on land to reduce correlation across the physical barrier.

    Source: §2.3 · p. 8

  • Results

    The Barrier SGF has the same theoretical computational cost as the stationary model, and the sparse Archipelago data yield substantially different posterior fields from the stationary SGF.

    Source: §6 · pp. 20–21

  • Takeaways & Limitations

    The authors recommend the Barrier model for applications where spatial estimation maps are important.

    Source: §6 · p. 21

    The land-to-normal range fraction is arbitrary, and the model cannot handle infinitely thin physical barriers.

    Source: §6 · p. 21

Abstract

from arXiv · show

The classical tools in spatial statistics are stationary models, like the Mat\'ern field. However, in some applications there are boundaries, holes, or physical barriers in the study area, e.g. a coastline, and stationary models will inappropriately smooth over these features, requiring the use of a non-stationary model. We propose a new model, the Barrier model, which is different from the established methods as it is not based on the shortest distance around the physical barrier, nor on boundary conditions. The Barrier model is based on viewing the Mat\'ern correlation, not as a correlation function on the shortest distance between two points, but as a collection of paths through a Simultaneous Autoregressive (SAR) model. We then manipulate these local dependencies to cut off paths that are crossing the physical barriers. To make the new SAR well behaved, we formulate it as a stochastic partial differential equation (SPDE) that can be discretised to represent the Gaussian field, with a sparse precision matrix that is automatically positive definite. The main advantage with the Barrier model is that the computational cost is the same as for the stationary model. The model is easy to use, and can deal with both sparse data and very complex barriers, as shown in an application in the Finnish Archipelago Sea. Additionally, the Barrier model is better at reconstructing the modified Horseshoe test function than the standard models used in R-INLA.

1 Introduction

Stationary spatial fields can smooth across physical barriers, while shortest-water-distance and boundary-condition approaches have limitations. The introduction motivates criteria for a practical, robust alternative and outlines the paper’s development and evaluation of the Barrier model.

1.1 Background

  • When physical barriers or holes matter, stationary fields can impose unrealistic dependence because their structure does not account for those features.The Finnish Archipelago example has many islands and peninsulas; a stationary field would smooth across them, creating an unrealistic dependence structure for aquatic animals.Source: §1.1 · p. 2, §1.2 · p. 2

1.2 Motivating example

  • The coastline problem includes physical barriers beyond coastlines, such as roads, power lines, mountains, and areas with different land use.The introduction assumes water is the normal area and land the barrier, reversing those roles when modelling data on land.Source: §1.2 · p. 2

1.3 Expanding the problem formulation

  • Different polygon approximations can alter shortest-water distances sharply, so models should not produce dramatically different results from slight boundary changes.The East–West channel disappears in one coastline approximation, changing the shortest distances between locations.Source: §1.3 · p. 2, Figure 2 · p. 3

1.4 Literature review

  • Shortest-distance methods can be cumbersome and sensitive to polygon choice, while boundary-condition methods can impose unrealistic assumptions or require known boundary values.Neumann conditions may strongly affect results; Dirichlet conditions can yield narrow intervals near boundaries and depend on boundary values being correct.Source: §1.4 · p. 4

1.5 Requirements for a new solution

  • The proposed solution is expected to be robust to polygon choice, near-stationary in computational cost, practical to use, flexible, and free of unrealistic new assumptions.The stated requirements also call for numerical stability on complex coastlines and suitability for sparse or rich data without extra tuning or assessment steps.Source: §1.5 · pp. 5–6

1.6 Outline of the paper

  • The paper defines and constructs the Barrier model, compares it on a modified Horseshoe reconstruction, and applies it to Finnish Archipelago fish data.The supplementary material provides additional model details and tutorials.Source: §1.6 · p. 6

2 The Barrier SGF

The Barrier SGF builds on the Matérn field but changes how dependence behaves across land. It interprets dependence through SAR paths, weakens paths crossing barriers, and uses an SPDE formulation to define the field.

2.1 Background

  • The Barrier SGF is based on the two-dimensional Matérn model and reduces to it in the absence of land; the smoothness parameter is fixed at ν = 1.The paper fixes ν because it can be unidentifiable in some cases and says this value is convenient and gives reasonably smooth fields.Source: §2.1 · pp. 6–7

2.2 Motivating the new approach

  • Interpreting the Matérn field as a SAR model makes dependence a collection of paths, which the approach weakens when paths cross land.The paths are not replaced by a shortest distance around the barrier; crossing paths are weakened to almost zero rather than completely removed.Source: §2.2 · p. 7

2.3 The Barrier SGF

  • On land, the Barrier SGF uses a Matérn field with the same σ and a range close to zero; the land range is fixed as a fraction of the normal-area range.The land and normal areas form the study area, and the paper gives rb = r/10 as an example to avoid adding a parameter requiring tuning, inference, or assessment.Source: §2.3 · p. 8
  • The field is chosen to be continuous across the boundary, with derivatives defined on each side of that boundary.The paper describes these as natural matching conditions and specifies the derivatives as orthogonal to the boundary, just inside the normal and land areas, respectively.Source: §2.3 · p. 8

3 Finite dimensional representation of u(s)

The field is represented with piecewise-linear finite elements on a mesh, yielding Gaussian coefficients with a sparse precision matrix. Barrier-specific local dependencies, mesh handling, and spatially varying uncertainty shape the resulting model.

3 Finite dimensional representation of u(s)

  • Piecewise-linear finite elements approximate the field as a sum of basis functions, with Gaussian coefficients whose precision matrix is Q.The linear finite-element approach avoids extra interpolation after computing the solution.Source: §3 · p. 9
  • Barrier and stationary SGFs differ in their Dq matrices, which encode local dependencies across the normal and barrier areas.The example matrices distinguish structural zeros where triangles do not connect nodes from numerical zeros arising from derivative inner products.Source: §3 · p. 10
  • The Barrier SGF is robust to boundary polygon selection, and the example mesh extends the outer boundary away from the region of interest.The blue rectangles mark land; the second plot zooms in on the first.Source: §3 · p. 10, Figure 4 · p. 10
  • The Barrier SGF has the same computational cost, O(n^3/2), as the stationary SGF because its precision matrix has the same sparsity structure.The repeated Cholesky factorization is the demanding operation; n is the number of finite-element mesh nodes.Source: §3 · pp. 10–11
  • The chosen marginal standard deviations vary spatially and are higher in inlets than under a stationary SGF.The authors choose this option for its interpretation, viewing inlet uncertainty as reasonable a priori.Source: §3 · p. 11

4 The modified Horseshoe reconstruction problem

The modified Horseshoe reconstruction compares Barrier, stationary, and Neumann spatial fields, using a test surface that does not satisfy typical boundary conditions. Simulations favor the Barrier model overall, while the Neumann model’s performance depends on sample size.

4 The modified Horseshoe reconstruction problem

  • The test uses the modified Horseshoe surface because the original satisfies the Neumann condition, giving Neumann-compatible approaches an unfair advantage.The three hierarchical Bayesian models use Gaussian observation noise and are fitted with R-INLA.Source: §4 · p. 12
  • On the modified Horseshoe surface, MB best matches the true reconstruction, while MS smooths over the gap and MN has height curves orthogonal to the boundary.Figure 8 shows posterior-mean reconstructions of the SGF.Source: §4 · p. 12, Figure 8 · p. 13
  • For σϵ = 0.1, the Barrier SGF’s RMSE is less than half the stationary SGF’s across the simulations, and MB is reported as best overall.The comparison uses 1000 runs per scenario; the authors report statistical and practical significance.Source: §4 · p. 14
  • MN performs better than MS with sparse data but significantly worse with richer data, because its zero-uncertainty Neumann boundary prior cannot be updated by data.The authors conclude that MN can perform worse than a stationary SGF.Source: §4 · p. 14

5 Modelling fish larvae in the Finnish Archipelago

The Archipelago application tests the Barrier SGF on sparse fish-larvae observations and a complex coastline. Its prior and posterior behavior differs from the stationary model around land, while cross-validation generally does not distinguish the spatial models.

Application context

  • The application addresses fish larvae in a complex coastline setting where sparse observations make inference challenging and habitat maps relevant to management.The data contain 198 observations per species.Source: §5 · p. 14

Prior behavior

  • The Barrier SGF smooths around land and is hindered by islands, whereas the stationary SGF smooths across land.The prior comparison used range 18 km in water and 3.6 km on land for the Barrier SGF.Source: §5 · p. 14, Figure 10 · p. 16

Model setup

  • The smelt analysis compares six model specifications per species, varying covariate inclusion and using stationary, Barrier, or no spatial field.The application models smelt, perch, and pikeperch larvae, with 198 observations for each species.Source: §5 · pp. 14–15

Posterior spatial field

  • For smelt larvae without covariates, the Barrier and stationary posterior means differ in how they smooth across land, including around a north-eastern peninsula.The Barrier SGF has higher uncertainty in south-eastern inlets, where the coastline separates them from nearby observations.Source: §5 · p. 15

Fixed effects

  • Changes in smelt fixed-effect estimates and intervals between stationary and Barrier models support recommending the Barrier SGF for future applications.The authors note that the differences may not justify a very resource-intensive non-stationary field, but describe the Barrier model as easy to use.Source: §5 · p. 15

Model comparison

  • LOOCV comparisons were mostly inconclusive: models with spatial fields were approximately equally good, while intercept-plus-iid models performed poorly.Scores used negative log predictive density and bootstrapped mean score differences; each model was run 198 times per dataset.Source: §5 · pp. 15, 17

6 Discussion

The discussion emphasizes the Barrier model’s computational efficiency and practical use for physical barriers, while qualifying its scope and the evidence for comparing methods. It also identifies parameter and geometry limitations and directions for further work.

Computational efficiency

  • The Barrier SGF has the stationary model’s theoretical computational cost because its precision matrix has the same sparsity structure.The current R implementation runs at different speeds from the stationary model implemented in C; the precision matrices support fast R-INLA fits at 100,000 observations and space-time Kronecker models.Source: §6 · p. 20

Comparison scope

  • The modified Horseshoe comparison is not a general basis for ranking coastline methods; the authors call for varied test functions and boundary conditions.A single example’s winner may depend on which method’s prior or low-penalty space most resembles the test function’s boundary.Source: §6 · p. 20

Spatial estimation maps

  • For spatial maps in the sparsely observed Archipelago, the stationary field smooths over a peninsula while the Barrier field does not, motivating the authors’ recommendation for map-focused applications.The authors describe the posterior spatial field as changing substantially between the two models in this setting.Source: Figure 12 · p. 18, §6 · p. 21

When to use the model

  • For coefficient-focused applications where the spatial field is a nuisance, coastline-aware and simple models may yield almost negligible differences.The authors nevertheless describe fitting effort for the Barrier model as almost the same as for the stationary model.Source: §6 · p. 21

Limitations

  • The range fraction rb/r is arbitrary, and infinitely thin barriers are unsupported; the authors suggest thickening such barriers to at least one mesh triangle and using a very small range fraction.Sensitivity analyses found little effect of the chosen fraction on parameter estimation or predictive performance; the authors advise against estimating it.Source: §6 · p. 21

Future research

  • Future work includes non-stationary models with k parameters, time-varying physical barriers, and extending boundary-condition discussions to SAR/CAR models..Source: §6 · p. 21

7 Supplementary materials

The supplementary materials are referenced for further details; the paper also points to code examples and mesh-construction guidance.

Implementation resources

  • The paper directs readers to supplementary material for additional details and lists online code examples and mesh-construction guidance for the R-INLA Barrier implementation.The implementation is named inla.barrier.pcmatern; fitting requires specifying the physical-barrier area, for example as a polygon.Source: §7 · p. 21

A Proof of theorem 1

The proof establishes operator properties needed to apply prior results and obtain existence and uniqueness for the SPDE solution. It also states solution regularity and notes that the same results can study the FEM approximation.

Operator properties

  • The proof uses that L is linear, self-adjoint, positive definite, and has a compact inverse to meet the conditions required by Bolin et al. (2017).The compact inverse follows from decay of the eigenvalues.Source: A Proof of theorem 1 · p. 24

Existence and uniqueness

  • Bolin et al. (2017) yields a unique continuous extension of L between the stated Sobolev spaces, which supports existence and uniqueness of the SPDE solution.The extension is an isometric isomorphism from ˙H^s to ˙H^{s−2} for any s.Source: A Proof of theorem 1 · p. 24

Solution regularity

FEM approximation

B The fish larvae dataset

The Archipelago Sea case study uses fish-larvae counts from 198 sampling locations, environmental covariates, and spatial coordinates. Sampling year was excluded, while RiverFlow was recoded as a factor.

B The fish larvae dataset

  • The data comprise larval counts and sampling effort at 198 Archipelago Sea locations, with six environmental variables available at 50 m resolution and spatial coordinates.The subset was collected in 2007 and 2011.Source: B The fish larvae dataset · p. 24
  • Sampling year was excluded because only two years were available and exploratory analysis found no differences between them.The authors note that a temporal component would be natural with more years.Source: B The fish larvae dataset · pp. 24–25
  • RiverFlow was recoded as a binary factor because 193 of 198 sites had only two distinct values, limiting the basis for a meaningful linear relationship.Zero represents values below average and one values above average.Source: B The fish larvae dataset · p. 25

C A small code example

The code example shows how users switch from a stationary spatial field to the Barrier SGF in the model formula. The barrier model also requires indices for mesh triangles covering land.

C A small code example

  • The example presents the stationary model formula as the reference code before showing the Barrier SGF substitution.Source: C A small code example · p. 25
  • For the Barrier SGF, replace the stationary model in the formula with barrier.model.The example otherwise retains the same formula structure and iid effect.Source: C A small code example · p. 25
  • The Barrier SGF implementation uses barrier.triangles, a list of indices for mesh triangles covering land.Source: C A small code example · p. 25

0 Supplementary material

The supplementary material specifies priors and predictor components for the horseshoe reconstruction and Archipelago analysis. Its hyperparameter priors use shrinkage toward a simple base model, with choices informed by model interpretation and study scale.

0.1 The prior model for the horseshoe reconstruction

  • The horseshoe reconstruction uses the same prior rationale as the Archipelago application, with λϵ = λσu = 1.5 and prior median spatial range 1.The simulated surface is the test surface to reconstruct, with locations sampled at random.Source: §0.1 · p. 26

0.2.1 The predictor

  • The predictor combines covariate effects, a spatial random effect, and an iid random effect; the spatial effect is compared as stationary, Barrier, or absent.The iid effect models extra variation, while the spatial field models spatial variability and unaccounted spatial structure.Source: §0.2.1 · p. 26

0.2.2 The priors for the hyper-parameters (r, σu, σϵ)

  • The priors are framed as shrinkage priors toward a simple base model, with the degree of shrinkage allowed to vary across applications to reflect measurement scale and competing model components.The base model is σϵ = 0, σu = 0, and r = ∞.Source: §0.2.2 · p. 27
  • The authors prefer fixed effects over random effects for explaining log-intensity variation, so σu and σϵ receive shrinkage priors.Smaller λ permits greater use of random effects and may weaken conclusions about fixed effects.Source: §0.2.2 · pp. 27–28
  • The iid-effect prior sets λϵ = 1.5 so its 1% and 99% quantiles are approximately -3 and 3 logarithmic units.The authors set these bounds so the iid effect could explain a datapoint on its own, though with low prior probability.Source: §0.2.2 · p. 28
  • The λ choice is somewhat arbitrary, but estimates remain reasonably close across λ values from 0.9 to 2; the Barrier penalty also varies over space.The authors report that this spatial variation matters less than variation among reasonable λ choices.Source: §0.2.2 · p. 28
  • The range prior favors longer spatial scales because the spatial effect is intended to capture non-local variation, leaving very short-scale variation to the iid effect.The 1/r parameterization addresses uneven behavioral differences across values of r.Source: §0.2.2 · p. 28
  • For the Archipelago case, the prior median range is 30 km; ranges around 4 km are unresolved by the mesh unless it is refined.This follows from a 60 km typical study-area length.Source: §0.2.2 · p. 29

0.3 Pike and pikeperch in the Archipelago

The perch and pikeperch analyses show species-specific spatial patterns and differences between stationary and Barrier models. For perch, the results resemble smelt and suggest similar covariate-confounding patterns; for pikeperch, model differences depend on the region and abundance across barriers.

0.3 Pike and pikeperch in the Archipelago

  • The supplementary section presents additional pike and pikeperch results beyond the selected smelt results in the main paper.Source: §0.3 · p. 29

0.3.1 Analysis of the perch data

  • Perch larvae show a distribution similar to smelt, with an additional south-east hotspot.Source: §0.3.1 · p. 29, Figure 2 · p. 31
  • Perch fixed-effect changes from the no-field to stationary and Barrier models resemble smelt’s, suggesting confounding with the spatial effect or shared missing covariates.The estimates’ covariate impacts are smaller than smelt’s, while significance levels are similar.Source: §0.3.1 · p. 29, Figure 1 · p. 30

0.3.2 Analysis of the pikeperch data

  • Pikeperch fixed effects differ from smelt for RiverFlow and ShoreDens, whose significance is lower, and show larger differences between stationary and Barrier models.Source: §0.3.2 · p. 29, Figure 4 · p. 33
  • For pikeperch, the stationary and Barrier spatial estimates are similar in the northern inlets because abundance is high on both sides of the land barrier.The figure caption reports only minor differences in spatial estimates without covariates.Source: §0.3.2 · p. 29, Figure 5 · p. 34

0.4 Model comparison for the archipelago

The archipelago comparisons assess stationary, Barrier, and IID models, then test Barrier predictions across mesh resolutions and barrier ranges. Barrier models outperform stationary models in simulation, while the real-data comparisons do not provide strong evidence that one model outperforms another.

0.4 Model comparison for the archipelago

  • The comparison ran 3,564 exact LOOCV fits across stationary, Barrier, and IID models, each with or without covariates.Source: §0.4 · p. 36
  • Figure 8 reports individual LOOCV scores with bootstrapped medians and 95% equal-tailed intervals; Figures 9–11 show pairwise-difference intervals.Source: §0.4 · p. 36, Figure 8 · p. 38, Figure 9 · p. 39
  • M5 is a statistically significant bad model; for pikeperch, spatial models outperform non-spatial models, but no result is strong evidence of one model outperforming another.Source: §0.4 · p. 36, Figure 11 · p. 40

0.5 Simulation study in the archipelago

  • The simulation compared 12 mesh and barrier-range variants using 2-fold cross-validation and RMSE on held-out predictions.The authors note that RMSE is not a proper scoring rule in the 5,000-dimensional joint prediction space.Source: §0.5 · p. 36
  • All Barrier models significantly outperform all stationary models, with similar predictive accuracy across barrier ranges and meshes.Coarser meshes perform consistently less well, but the differences are small relative to uncertainty.Source: §0.5 · p. 36, Table 1 · p. 38
Loading 1608.03787v3…