fbesag (Partitioned Besag)

Overview

The fbesag model splits a Besag (ICAR) graph into P partitions via a user-supplied id vector, then assigns each partition its own log-precision hyperparameter \(\theta_p = \log \tau_p\). This lets the smoothing strength vary across the graph, which is more flexible than a single-scale Besag.

Within each partition, the conditional structure is the standard Besag CAR; across partitions a partition-level Gaussian effect with prior standard deviation sd_gamma ties the precisions together. A PC prior controls deviation of the partition log-precisions from a common reference value via param = (p1, p2), parameterised so that \(P(|\theta| > p_1) = p_2\) (the same scale used by pyinla's other PC priors).

How it differs from Besag

• Hyperparameter count: Besag has one (the global log-precision). fbesag has P (one per partition).
• Use case: Besag assumes the whole graph smooths at the same scale. fbesag allows e.g. urban vs rural partitions to smooth differently.
• Accessed via: Besag is a native pyinla model ('model': 'besag'). fbesag is a cgeneric wrapper: it ships as a precompiled C library inside the pyinla.gmrfs subpackage and you call pyinla.gmrfs.fbesag(…) to build a spec dict.

Runnable example

Below is a complete, self-contained example: an 8-node 4×2 grid split into two partitions (top row vs bottom row). Copy-paste into a Python prompt with pyinla installed and it runs end to end.

import numpy as np
import pandas as pd
from pyinla import pyinla
from pyinla.gmrfs import fbesag

# 4x2 grid (8 nodes):
#   1 - 2 - 3 - 4
#   |   |   |   |
#   5 - 6 - 7 - 8
n = 8
adj = np.zeros((n, n), dtype=int)
for i, j in [(1,2),(2,3),(3,4),
              (5,6),(6,7),(7,8),
              (1,5),(2,6),(3,7),(4,8)]:
    adj[i-1, j-1] = adj[j-1, i-1] = 1

# Pass the Besag precision Q = diag(degree) - A (intrinsic CAR structure).
Q = np.diag(adj.sum(axis=1)) - adj

# Partition id per node (top row -> 1, bottom row -> 2).
id_part = [1, 1, 1, 1, 2, 2, 2, 2]

# Synthetic Gaussian data, one observation per node.
rng = np.random.default_rng(42)
df = pd.DataFrame({
    'y':    rng.standard_normal(n),
    'area': np.arange(1, n + 1),
})

spec = fbesag(
    id        = 'area',
    graph     = Q,
    id_part   = id_part,
    sd_gamma  = 0.2,
    param     = {'p1': 1, 'p2': 1e-5},
)

result = pyinla(
    model={'response': 'y', 'fixed': ['1'], 'random': [spec]},
    family='gaussian', data=df,
)
print(result.summary_hyperpar)

Expected output (excerpt): a posterior with one Gaussian-noise precision and two Theta entries for area, one per partition.

Wrapper signature

Reference

If you use this model in published work, please cite:

Abdul-Fattah, E., Krainski, E., Van Niekerk, J., and Rue, Håvard (2024). Non-stationary Bayesian spatial model for disease mapping based on sub-regions. Statistical Methods in Medical Research, 33(6), 1093–1111. https://doi.org/10.1177/09622802241253518

@article{abdul2024non,
  title   = {Non-stationary Bayesian spatial model for disease mapping based on sub-regions},
  author  = {Abdul-Fattah, Esmail and Krainski, Elias and Van Niekerk, Janet and Rue, H{\aa}vard},
  journal = {Statistical Methods in Medical Research},
  volume  = {33},
  number  = {6},
  pages   = {1093--1111},
  year    = {2024},
  publisher = {SAGE Publications Sage UK: London, England}
}