PaperlayerView exact revision on arXiv

Source-linked AI summary

An intuitive Bayesian spatial model for disease mapping that accounts for scaling

Andrea Riebler, Sigrunn H. Sørbye, Daniel Simpson, Håvard Rue

arXiv:1601.01180v1stat.ME

TL;DR

Disease mapping models must separate structured and unstructured variation while assigning hyperpriors whose meanings do not depend on the spatial graph. The paper evaluates a scaled BYM2 parameterisation with PC priors and finds good learning and shrinkage, while noting interpretability limitations in competing models.

  • Problem

    Existing BYM formulations confound structured and unstructured effects, while unscaled precision parameters and hyperpriors lack graph-independent, interpretable meanings.

  • Method

    The paper evaluates a scaled BYM2 parameterisation using PC priors to control total variance and distribute it between structured and unstructured components.

  • Results

    The BYM2 model shows good shrinkage and learning in simulations, with hyperparameters that are directly interpretable and performance at least comparable to alternatives.

  • Takeaways & Limitations

    Scaling and PC priors provide epidemiologically intuitive hyperparameter specifications that remain meaningful across spatial applications.

  • Takeaways & Limitations

    Comparisons with other models require care because their priors and parameter interpretations are not equivalent, and most are not properly scaled.

Abstract

from arXiv · show

In recent years, disease mapping studies have become a routine application within geographical epidemiology and are typically analysed within a Bayesian hierarchical model formulation. A variety of model formulations for the latent level have been proposed but all come with inherent issues. In the classical BYM model, the spatially structured component cannot be seen independently from the unstructured component. This makes prior definitions for the hyperparameters of the two random effects challenging. There are alternative model formulations that address this confounding, however, the issue on how to choose interpretable hyperpriors is still unsolved. Here, we discuss a recently proposed parameterisation of the BYM model that leads to improved parameter control as the hyperparameters can be seen independently from each other. Furthermore, the need for a scaled spatial component is addressed, which facilitates assignment of interpretable hyperpriors and make these transferable between spatial applications with different graph structures. We provide implementation details for the new model formulation which preserve sparsity properties, and we investigate systematically the model performance and compare it to existing parameterisations. Through a simulation study, we show that the new model performs well, both showing good learning abilities and good shrinkage behaviour. In terms of model choice criteria, the proposed model performs at least equally well as existing parameterisations, but only the new formulation offers parameters that are interpretable and hyperpriors that have a clear meaning.

1. Introduction

Disease mapping commonly uses Bayesian hierarchical models to obtain more reliable regional risk estimates by borrowing strength across neighboring areas. The paper motivates a scaled BYM2 parameterisation with interpretable parameters and hyperpriors.

  • Motivation: Bayesian hierarchical models improve regional disease estimates by borrowing strength across neighboring regions.They model spatial random effects alongside available covariates to capture unobserved heterogeneity or spatial correlation.
  • Model challenges: Combining structured and unstructured spatial effects creates confounding over how variability is divided between components.Reparameterised models replace separate precisions with a common precision and mixing parameter.
  • Model challenges: Unscaled spatial components make precision depend on graph structure and remain confounded with the mixing parameter.Consequently, fixed priors are not transferable across applications with different graphs, and hyperprior choice remains difficult.
  • Proposed approach: The proposed BYM parameterisation maps precision to marginal standard deviation, making parameters interpretable and hyperpriors meaningful.It applies penalised-complexity priors and treats the model as an extension of simpler base models.
  • Implementation: The study implements models with integrated nested Laplace approximations for full Bayesian inference without Markov chain Monte Carlo.The implementation uses the R-package INLA.
  • Evaluation: A simulation study investigates the new model’s properties and compares it with commonly used spatial models.The paper also includes an application following the simulation comparisons.

2. Disease Mapping

Disease mapping estimates geographical disease burden across areas while accounting for unstable crude rates. Bayesian hierarchical models smooth these estimates by borrowing strength across neighboring regions and incorporating latent risk variation.

  • Purpose: Disease mapping studies the geographical distribution of disease burden across non-overlapping areas.The framework defines expected populations at risk and area-specific relative risks.
  • Data model: Rare, non-contagious disease counts are modeled using an assumption stated for the number of cases or deaths in each region.The supplied passage introduces this assumption but does not provide the full displayed equation.
  • Data model: The underlying true area-specific relative risk is denoted by θi.The index i ranges over the study’s areas.
  • Estimation challenge: Directly mapping standardised mortality ratios can mislead because extreme estimates often occur where expected counts are small.Sampling variability is larger in such areas, with Var(θ̂i) = θi/Ei.
  • Smoothing: Bayesian hierarchical models smooth raw estimates by borrowing strength across neighboring regions.Their latent random effects provide local smoothness while accounting for extra-Poisson variation or spatial correlation.
  • Model formulation: The general latent formulation decomposes log risk ηi = log(θi) into covariate effects and a random effect.The random effects account for latent or unmeasured risk factors through extra-Poisson variation or spatial correlation.

3. Modelling the spatial dependency structure

Spatial disease models represent regional effects through structured and unstructured components, but their parameterisation, identifiability, and scaling determine whether priors are interpretable across applications.

  • Neighbourhood structures define which regions are spatially related, commonly by shared borders, with δ_i denoting region i’s neighbours.
  • Besag model: The Besag model uses an intrinsic GMRF whose conditional mean follows neighbouring effects, but it cannot represent purely unstructured variability.Its precision matrix is singular, so sum-to-zero constraints are imposed on connected subgraphs to prevent confounding with the intercept.
  • BYM model: The BYM model adds unstructured variation to the structured Besag component, yet the two components are not independently identifiable and require difficult hyperprior choices.The unstructured precision has a marginal interpretation, whereas the structured precision controls conditional variability given neighbouring effects.
  • Leroux model: The Leroux model makes the balance between unstructured and structured variation explicit through a mixing parameter φ, with φ = 0 giving pure overdispersion and φ = 1 the Besag model.
  • Dean model: Dean’s reparameterisation replaces component precisions with a common precision and mixing parameter, distributing variance between structured and unstructured effects.
  • Need for scaling: Unscaled structured components make precision depend on graph structure, so identical priors can imply different marginal variability across applications.Graphs with the same number of regions can have different generalized variances; examples give approximately 0.40 for Arnsberg, 0.56 for Germany, and 4.78 for continental U.S. counties.
  • Need for scaling: A precision hyperprior controls local deviation from the constant-level null space, with excessive precision oversmoothing and insufficient precision risking overfitting.
  • Need for scaling: Scaling makes the generalized variance equal to one, so τ_b represents marginal precision independently of the underlying graph.Sparse-matrix implementation is supported through null-space constraints and the INLA inla.scale.model function.

3.3. A modified BYM model accounting for scaling

The modified BYM2 parameterisation uses a scaled structured component and separates marginal precision from the mixing parameter, making both interpretable across graph structures. An augmented formulation preserves sparse precision matrices for efficient computation.

  • 3.3. A modified BYM model accounting for scaling: The BYM2 model uses a scaled structured component u⋆ and maps precision to marginal variability, addressing BYM identifiability and scaling issues.The mixing parameter φ represents the proportion of marginal variance attributed to spatial structure.
  • 3.3. A modified BYM model accounting for scaling: 0 ≤φ ≤1 controls the compromise between pure overdispersion and spatially structured correlation, with φ = 0 and φ = 1 yielding the two limiting models.The model reduces to pure overdispersion at φ = 0 and to the Besag model at φ = 1.
  • 3.3. A modified BYM model accounting for scaling: Scaling makes the prior on τb transferable between graphs because τb and φ are no longer confounded.The Leroux model cannot be scaled in the same way because its scaling would depend on φ.
  • 3.4. Parameterisation preserving sparsity: The covariance formulation introduces a generalised inverse Q⋆−, whose inverse precision matrix is no longer sparse.This loss of sparsity would reduce computational efficiency in INLA and block-updating MCMC methods.
  • 3.4. Parameterisation preserving sparsity: An augmented parameterisation restores sparsity while retaining the correct marginal distribution, with the structured component represented directly by the second half of w.This makes the formulation suitable for sparse computational methods.

4. Hyperprior choice: Penalised-complexity priors

The paper uses penalised-complexity priors to assign interpretable hyperpriors to total variability and its spatially structured fraction. These priors shrink toward simpler base models while allowing prior statements to be transferred across applications with scaled graphs.

  • 4. Hyperprior choice: Penalised-complexity priors: PC priors distribute total variance between the modified BYM components through a prior on τb and a separate prior on φ.τb controls marginal variance, while φ allocates variance between unstructured v and structured u⋆ effects.
  • 4. Hyperprior choice: Penalised-complexity priors: PC priors prefer simpler models, measure complexity with KLD, penalise distance at a constant rate, and use a user-defined probability statement.The supplied passage introduces these four construction principles, although the fourth principle is truncated.
  • 4.1. Controlling the total variability: The first-level base model has constant relative risk across regions, representing no spatial variation and infinite precision.The prior is constructed from the information-theoretic distance between this base model and the flexible model.
  • 4. Hyperprior choice: Penalised-complexity priors: Figure 2 compares the PC prior for τb with two gamma priors and the PC prior for φ with a uniform prior on (0,1).The displayed PC-prior settings are U = 1, α = 0.01 for τb and U = 0.5, α = 2/3 for φ.
  • 4.1. Controlling the total variability: τb is interpreted as marginal precision, allowing prior statements about total residual relative risk; for example, Prob((1/√τb) > 1) = 0.01 corresponds approximately to residual risks below 2 with probability 0.99.The passage cautions that gamma priors cannot attain the zero-distance base model because they place no density mass there.
  • 4.2. Controlling the spatial variance contribution: The second-level base model sets φ = 0, meaning only unstructured noise remains; increasing φ blends in spatial dependency.The model component explaining most variance shifts from v toward u⋆ as φ increases.
  • 4.2. Controlling the spatial variance contribution: φ directly represents the fraction of total variance attributed to spatial dependency, enabling statements such as Prob(φ < 0.5) = 2/3.This choice places more prior mass below 0.5 and conservatively favors unstructured variability.
  • 4.2. Controlling the spatial variance contribution: The PC prior for φ is computed on the logit scale for a specific Besag matrix Q, with graph scaling implemented by default.INLA can produce the prior in tabulated form for use in other software packages.

5. Applications

Applications show that BYM2 can learn the relevant spatial structure and shrink toward appropriate base models, while achieving model-choice performance comparable to existing approaches. However, competing parameter estimates require careful interpretation because priors and scaling differ across models.

  • Simulation study: 200 simulations across nine scenarios assessed BYM2 learning and shrinkage under constant, unstructured, and structured risk surfaces with varying expected case counts.Expected cases were set to 15, 60, or 200, using the neighbourhood structure of 366 Sardinian districts.
  • Simulation results: Under constant risk, the estimated structured standard deviation was close to zero and decreased with higher prevalence, while the prior dominated the mixing parameter.These data provided no information about the mixing parameter φ.
  • Simulation results: For independent area-specific risk, BYM2 estimated mixing parameters near zero, shrinking variation toward the unstructured base model.The posterior mean parameter estimates were close to the true values.
  • Simulation results: For spatially structured risk, the model estimated the structured standard deviation well and produced mixing estimates closer to one at higher expected case counts.The lower estimate when Ei = 15 compared with Ei = 60 and Ei = 200 indicates learning of a more complex model when supported by the data.
  • Prior sensitivity: Changing the mixing-parameter prior had little effect for unstructured or structured risk, but constant-risk estimates moved toward 0.5 under a uniform prior.In the constant-risk case, the mixing-parameter estimate was dominated by the new prior.
  • Comparison to other spatial latent models: BYM2 performed at least as well as competing models by model-choice criteria, although parameter estimates are not directly comparable because priors and parameterisations differ.Average criteria values were similar across models, with DIC slightly favouring the two BYM2 models.

6. Discussion

The BYM2 model addresses scaling and interpretability problems in Bayesian disease mapping by separating marginal precision from the mixing parameter. Its practical advantages accompany model-choice performance at least as good as existing methods, while further prior-sensitivity and broader-application work remains.

  • Model properties: Scaling makes hyperpriors transferable across applications with different graph structures and gives the hyperparameters direct interpretations.Without scaling, precision depends on the underlying graph and is confounded with mixing.
  • Model properties: The BYM2 model uses marginal precision and a mixing parameter that separately control overall variability and its allocation between unstructured and structured components.This parameterisation addresses graph-dependent interpretation and confounding concerns.
  • Model performance: The BYM2 model shrinks toward constant or spatially unstructured risk surfaces and learns well when risk is spatially structured across simulated disease prevalences.The authors interpret these results as evidence that the model does not overfit.
  • Model performance: For constant risk, BYM2 is slightly favored over Leroux and Dean by DIC and logarithmic score; for unstructured risk, model-choice characteristics are almost indistinguishable.The paper prioritizes interpretable parameters and hyperpriors rather than optimizing model-choice criteria alone.
  • Implementation: The model can be fitted efficiently with INLA, whose implementation preserves the sparse spatial structure and is illustrated with R code.The implementation uses a graph file, a BYM2 latent-effect specification, and posterior hyperparameter refinement through inla.hyperpar.
  • Open issues: Further work is needed for detailed prior-sensitivity analysis and for distributing variance across spatial components and other linear-predictor parameters in covariate or space-time models.The paper also notes potential applications beyond disease mapping, including genetics.

Appendix A. Supplementary tables and figures

The supplementary material reports simulation summaries and model-comparison quantities across risk surfaces, disease prevalences, and competing models. The experiments use repeated simulations and track estimation and predictive criteria with lower values preferred for the latter criteria.

  • Simulation tables: Table 4 varies risk surface type and disease prevalence for BYM2 simulations, reporting means and standard deviations over 200 simulations per setting.The settings include constant, spatially varying unstructured, and spatially varying structured risks, with true parameter values provided.
  • Simulation figures: Figure 5 displays logarithmic-score boxplots for all six models across nine simulation settings.The figure summarizes score distributions rather than only mean performance.
  • Simulation figures: Figure 6 displays DIC boxplots for six models across nine settings, separating constant, iid, and spatially structured risks by columns and using E = {15, 60, 200} across rows.Each setting contains 200 simulations, and y-axis range widths are kept constant only within columns.
  • Application comparison: Table 5 compares RMSE, DIC, and logarithmic score for different models analyzing Sardinian IDDM counts under four priors for φ.The priors include three PC-prior specifications and a uniform prior on (0, 1).

Appendix B. INLA R-code to implement the BYM2 model

The appendix demonstrates an INLA workflow for fitting BYM2 to Sardinian disease-count data. It covers graph and data input, model specification, fitting, hyperparameter refinement, and posterior inspection.

  • Data preparation: The code loads the INLA package and reads Sardinian data containing counts, expected counts, and standardized mortality ratios.The dataset is read from sardinia.dat with columns y, E, and SMR.
  • Spatial structure: The neighbourhood graph encodes 366 regions, their neighbour counts, and neighbour indices, which INLA transforms into the spatial matrix Q.The graph is supplied through a file such as sardinia.graph.
  • Model specification: The formula object specifies an intercept and BYM2 latent effect, while the graph and model scaling are supplied in the INLA specification.The R code uses the f(.) function to define the regional random effect.
  • Model fitting: The inla call receives the formula, data, Poisson likelihood, and expected counts, and requests posterior estimates for ηi and θi.The workflow then calls inla.hyperpar for improved hyperparameter estimates.
  • Output inspection: Posterior marginals for fixed effects and hyperparameters can be accessed through result$summary.fixed and result$summary.hyperpar.Overall results can also be inspected with summary(result) or plot(result).

C.1. Precision parameter τb

The PC-prior construction for τb begins from a high-precision base model, measures complexity through a Kullback–Leibler distance, and transforms an exponentially penalized distance prior into a prior for τb.

  • Construction: An exponential prior is assigned to the distance d(τb) with decay rate λ, implementing constant-rate penalisation before a change-of-variables transformation.The transformation maps the prior on distance to a prior on τb.
  • Resulting prior: In the limit τ0 →∞, the resulting prior for τb is a type-2 Gumbel distribution.This is the limiting form obtained from the stated PC-prior derivation.

C.2. Mixing parameter φ

The PC prior for φ is derived by comparing base and flexible covariance models through their Kullback–Leibler divergence. Computational shortcuts exploit sparse precision structure and the singularity of Q⋆.

  • The base and flexible covariance matrices are defined as Σbase = I and Σflex(φ) = (1 −φ) I + φQ−.
  • The PC prior derivation uses the Kullback–Leibler divergence between the base and flexible covariance models.
  • The trace computation is quick when Q⋆ is sparse.
  • The determinant can be evaluated from the eigenvalues of Q⋆, using ˜γi = 1/γi for positive eigenvalues and ˜γi = 0 otherwise.
  • Because Q⋆ has rank deficiency 1, one eigenvalue is zero; eigenvalue calculations can be avoided in high-dimensional settings using the cited appendix method.
Loading 1601.01180v1…