Scottish Lip Cancer Disease Mapping

Step 1: Prepare the Adjacency Graph

The adjacency structure defines which districts are neighbors. Three island districts (Orkney, Shetland, Western Isles) have no neighbors in the original graph and must be connected to their nearest mainland district to form a single connected component:

import json
from scipy.sparse import diags, csc_matrix
from scipy.sparse.linalg import inv as sparse_inv

# Load district centroids from GeoJSON
with open("scotland_districts.geojson") as f:
    gj = json.load(f)

from shapely.geometry import shape
coords = np.array([shape(feat['geometry']).centroid.coords[0] for feat in gj['features']])

def compute_icar_scaling(adj):
    """Geometric mean of ICAR marginal variances (scaling factor)."""
    n = adj.shape[0]
    D = diags(adj.sum(axis=1))
    Q = D - csc_matrix(adj)
    # Add small constant for invertibility (sum-to-zero constraint)
    Q_inv = sparse_inv(csc_matrix(Q + np.ones((n, n)) / n)).toarray()
    marginal_var = np.diag(Q_inv)
    return np.exp(np.mean(np.log(marginal_var)))

def parse_adjacency(filepath):
    """Parse INLA-format adjacency file into a symmetric matrix."""
    with open(filepath, 'r') as f:
        lines = f.readlines()

    n = int(lines[0].strip())
    adj = np.zeros((n, n), dtype=float)

    for line in lines[1:]:
        parts = [int(x) for x in line.split()]
        node = parts[0] - 1  # 0-based
        for nb in parts[2:]:
            adj[node, nb - 1] = 1.0
            adj[nb - 1, node] = 1.0

    return adj

adj_matrix = parse_adjacency("scotland.adj")

# Connect islands to nearest mainland neighbor by centroid distance
# Orkney    -> Caithness
# Shetland  -> Caithness
# Western Isles -> Skye-Lochalsh

num_nb = adj_matrix.sum(axis=1)
singleton_idx = np.where(num_nb == 0)[0]
connected_idx = np.where(num_nb > 0)[0]

for s in singleton_idx:
    dists = np.sqrt(((coords[s] - coords[connected_idx])**2).sum(axis=1))
    nearest = connected_idx[np.argmin(dists)]
    adj_matrix[s, nearest] = 1
    adj_matrix[nearest, s] = 1
    print(f"Connected: {data['county'].iloc[s]} -> {data['county'].iloc[nearest]}")

# Save connected adjacency matrix in INLA format
import os
os.makedirs("output", exist_ok=True)
with open("output/scotland_connected.adj", "w") as f:
    n = adj_matrix.shape[0]
    f.write(f"{n}\n")
    for i in range(n):
        nbs = np.where(adj_matrix[i] > 0)[0] + 1  # 1-based
        f.write(f"{i+1} {len(nbs)} {' '.join(map(str, nbs))}\n")

scaling_factor = compute_icar_scaling(adj_matrix)
print(f"ICAR scaling factor: {scaling_factor:.6f}")

Connected: Orkney -> Caithness
Connected: Shetland -> Caithness
Connected: Western Isles -> Skye-Lochalsh
All 56 nodes now in one connected component
ICAR scaling factor: 0.625549

ICAR Scaling Factor

The scaling factor $s$ is computed as the geometric mean of the marginal variances of the ICAR precision matrix. Dividing the ICAR component by $\sqrt{s}$ standardises its variance, making the two precision hyperparameters comparable.

Step 2: Fit the Model with pyINLA

The model is specified as a Python dictionary. The response field names the outcome variable, fixed lists the intercept and covariates, and random defines the spatial random effect. Inside the random effect we set model='bym' with scale.model=True to use the scaled BYM, and assign PC priors to both hyperparameters: theta1 controls the IID precision and theta2 controls the spatial (ICAR) precision. The control block requests model comparison criteria (DIC, WAIC) and sets vague Gaussian priors (precision 0.001, i.e. $\sigma = 1/\sqrt{0.001} \approx 31.62$) on the fixed effects. The expected counts E are passed separately as an offset.

from pyinla import pyinla

model = {
    'response': 'Y',
    'fixed': ['1', 'AFF'],
    'random': [{
        'id': 'idarea',
        'model': 'bym',
        'graph': 'output/scotland_connected.adj',
        'scale.model': True,
        'hyper': {
            'theta1': {'prior': 'pc.prec', 'param': [1, 0.01]},
            'theta2': {'prior': 'pc.prec', 'param': [1, 0.01]}
        }
    }]
}

result = pyinla(
    model=model,
    family='poisson',
    data=data,
    E=data['E'].values,
    control={
        'compute': {'dic': True, 'waic': True, 'config': True, 'return_marginals': True},
        'fixed': {'prec.intercept': 0.001, 'prec': 0.001},
    }
)

print("Fixed effects:")
print(result.summary_fixed[['mean', 'sd', '0.025quant', '0.975quant']])
print("\nHyperparameters:")
print(result.summary_hyperpar[['mean', 'sd', '0.025quant', '0.975quant']])

Fixed effects:
                 mean        sd  0.025quant  0.975quant
(Intercept) -0.264939  0.123755   -0.507560     -0.0208
AFF          4.217524  1.279915    1.662767      6.6975

Hyperparameters:
                                                mean           sd  0.025quant  0.975quant
Precision for idarea (iid component)      867.046703  6048.222795    8.009324  5845.17006
Precision for idarea (spatial component)    4.762687     1.855708    2.049507     9.23812

PyINLA fits the full scaled BYM model in approximately 0.63 seconds, including computation of marginals, DIC, and WAIC.