1What the method is for

Handling the rival explanation — “is it not simply that neighbouring settlements resemble one another?” — inside the model.

Data from sites laid out across a geography almost always carry spatial structure. Adjacent settlements have people moving between them, belong to the same school district, and share a similar history. So the premise that observations from different settlements are mutually independent does not hold as it stands. A regression that assumes independence counts this overlapping information as if it were new, and underestimates the uncertainty on the coefficients.

There are two ways to deal with it. One is to remove the spatial structure before regressing; the other is to estimate the spatial structure explicitly as part of the model. ICAR is the latter. A spatial random effect \(\phi_i\) is placed at each site, and \(\phi\) is given a prior saying that neighbours take nearby values.

What this study wanted to ask was whether the association between population growth rate and vocabulary loss could be explained away by spatial proximity. Nearby settlements have similar population dynamics, so the point of the check is whether the coefficient survives once a spatial effect is included.

2Origin — a method that came from cancer maps

Epidemiology built it in order to map incidence rates by region.

The prototype is the model published in 1991 by Julian Besag (1945–2010, English statistician) with York and Mollié. Its original purpose was to draw maps of cancer incidence by region. Computing a rate naively for each region gives a small denominator wherever the population is small, so that a few extra cases by chance send the rate leaping upward. The map then depicts sampling variation rather than reality.

So the knowledge that “of course a region resembles its neighbours in incidence” is brought in as a prior, pulling each region’s estimate towards its neighbours. The smaller the population, the more strongly it is pulled, and the map stabilizes. Rather than removing spatial autocorrelation as a nuisance, this model uses it as information that helps the estimation.

From epidemiology to language research

The epidemiological picture, in which rural areas near a city show rates resembling the city’s, transfers directly to the picture in linguistic geography, in which neighbouring settlements show similar vocabulary change. It is no accident that this study uses rarefaction from ecology and spatial models from epidemiology: the shape of the problem — aggregated data whose sample size differs from site to site — is shared.

3Where it comes from mathematically

A spatial version of the Markov property, written as a random field on a graph.

The mathematical root of ICAR is the same as that of Markov chains (Markov 1906). The Markov property in a chain over time said that once the present is fixed, the past is irrelevant to the future. Moved into space, it becomes “once the values at adjacent sites are fixed, sites further away are irrelevant”. A random field with this property is called a Markov random field, and the Hammersley–Clifford theorem gives the conditions under which specifying the conditional distributions determines the joint distribution.

ICAR is therefore a model that, instead of writing the joint distribution directly, writes only the conditional distribution at each site. Specifying local relations without specifying the whole is its characteristic feature.

4Definitions and assumptions

Conditional autoregression (CAR), and ICAR as its limit.

CAR — defined by the conditional distribution at each site
\[ \phi_i \mid \phi_{-i} \;\sim\; \mathcal{N}\!\left(\alpha\,\frac{\displaystyle\sum_{j} w_{ij}\,\phi_j}{\displaystyle\sum_{j} w_{ij}},\;\; \frac{\sigma^{2}}{\displaystyle\sum_{j} w_{ij}}\right) \]

\(w_{ij}\) … an entry of the adjacency matrix; 1 if adjacent, 0 otherwise
\(\phi_{-i}\) … the spatial effects at every site other than \(i\)
The mean … the average over adjacent sites multiplied by \(\alpha\). A site’s own value is pulled towards the average of its neighbours
The variance … smaller for sites with more neighbours. The more information, the more firmly the value is pinned down

ICAR — the limit at \(\alpha = 1\)
\[ \phi_i \mid \phi_{-i} \;\sim\; \mathcal{N}\!\left(\bar{\phi}_{\partial i},\;\; \frac{\sigma^{2}}{n_i}\right), \qquad p(\boldsymbol{\phi}) \;\propto\; \exp\!\left[-\frac{1}{2\sigma^{2}}\sum_{i \sim j}\bigl(\phi_i - \phi_j\bigr)^{2}\right] \]

\(\bar{\phi}_{\partial i}\) … the mean of \(\phi\) over the sites adjacent to site \(i\). \(n_i\) … the number of neighbours of site \(i\)
\(i \sim j\) … a sum over adjacent pairs. The penalty takes the form of a sum of squared differences between neighbours, so that the further apart adjacent sites are, the lower the probability

ICAR is improper as a prior

Look closely at the expression on the right: adding the same constant to every element of \(\boldsymbol{\phi}\) leaves the differences \(\phi_i - \phi_j\) unchanged. In other words this density fixes the overall level not at all. Being flat in the constant direction, its integral diverges and it cannot be normalized as a probability distribution. That is what “improper” means.

There are two practical consequences. First, as it stands it is incomplete as a prior. Second, if the model has an intercept, the level of the intercept and the level of the spatial effect cannot be told apart, and identifiability is lost. The standard remedy is to impose a sum-to-zero constraint \(\sum_i \phi_i = 0\), which pushes the level onto the intercept and leaves the spatial effect carrying only the differences between sites. The car(type = "icar") in brms is the implementation with this constraint.

The adjacency matrix is decided by the analyst, not by the data

The \(w_{ij}\) appearing in the expressions above does not emerge automatically from anywhere. What counts as a neighbour is a design decision. Sites within a fixed distance? The \(k\) nearest? Regions sharing an administrative boundary? A Delaunay triangulation? The choice changes the distribution of neighbour counts, which changes the strength of the pull, which changes the spatial effect that is estimated. The data do not settle this choice.

5How it is used in this research

Neighbours were defined as sites within 3 km, and models with and without a spatial term were compared.

  1. Computing distances　Great-circle distances (haversine formula, Earth radius 6371 km) are computed from the latitude and longitude of each settlement.
  2. Building the adjacency matrix　A binary matrix \(W\) is formed, taking pairs at a distance greater than 0 and at most 3 km as adjacent (1) and everything else as 0. The 3 km threshold follows the value used in the earlier 2016 analysis.
  3. Checking the adjacency structure　The resulting \(W\) is 25×25, with a mean of 4.80 neighbours and no isolated sites. ICAR becomes unstable if isolated sites are present, so this check is essential.
  4. Entering it into the model　car(W, gr = settlement, type = "icar") is added to the hierarchical beta regression. The estimation framework is the same as on the Bayesian estimation page: four chains, 8,000 draws, adapt_delta = 0.97.
  5. Judging whether the spatial term is needed　The models with and without the spatial term are compared by approximate leave-one-out cross-validation, to see whether the data support the spatial term.
Comparisonelpd_looelpd_diffse_diffVerdict
Loss rate, no spatial term (base)33.80.00.0this one is selected
Loss rate, with ICAR (sp)33.1−0.70.3the spatial term does not improve fit

Figures are quoted from the run log of R/12_brms_loss.R. Adding the spatial term barely moved the coefficient on population growth rate: \(\beta = +0.228\), 95% credible interval [+0.04, +0.40], \(P(\beta>0) = 0.992\) (against \(+0.231\) without the spatial term).

The spatial variance collapsed to near zero

In the diversity model of the same study (51 neighbourhood associations) the spatial variance sdcar can be read directly. For the Kansai stratum it is 0.04, 95% credible interval [0.00, 0.16]; for the total trait count, 0.02, [0.00, 0.09]. In both the lower bound touches 0, and no spatial signal is detected. In this situation the posterior takes a funnel shape, in which each individual \(\phi_i\) shrinks as the spatial variance approaches 0, and the sampler struggles. In this model max_treedepth was in fact reached every time, the effective sample sizes were as low as 179 and 57, and \(\hat{R}\) was 1.02 and 1.07.

Conclusion: not in the main model, but used as a robustness check

On this basis the study does not include the spatial term in its main model. ICAR is used not to “correct for spatial structure” but to confirm that spatial structure does not account for the finding. Loading 25 spatial random effects onto a sample of 25 sites is excessive, and LOO does not support it either.

6Limits — what the method cannot say

Change the definition of adjacency and the conclusion can change

This is not a general remark; it is what actually happened within this group of studies. The 2016 analysis defined adjacency at distance thresholds of 1, 2 and 3 km and, using Moran’s I and Geary’s C, detected no significant spatial autocorrelation either in the loss and gain rates or in the regression residuals. The 2017–2018 analysis, however, redefined adjacency in several ways, including Delaunay triangulation and \(k\) nearest neighbours, and obtained indicators pointing to the presence of spatial autocorrelation.

The working record leaves this discrepancy flagged as a point requiring re-examination, and which reading is correct has not been settled. What can be stated here is the fact itself: that for the same data, changing the definition of adjacency produced different indicators as to whether spatial autocorrelation is present. The 3 km used on this page therefore has to be read as one design decision among others. A sensitivity analysis at thresholds of 1 km and 2 km remains an outstanding task in the working record.

It does not amount to proof that a spatial effect is absent

Even with LOO selecting the non-spatial model and sdcar pinned at 0, what that shows is that no spatial signal could be detected under this definition of adjacency and at this sample size. With a sample of 25 sites the power itself is low. The 2016 analysis also recorded the reservation that “Moran says no, but power is low in a small sample”, and the Bayesian CAR re-examination was introduced in the first place in order to answer that reservation.

In addition, the LOO computation for the spatial model produced one observation (4.0%) with a Pareto \(k\) above 0.7, and relative effective sample sizes as low as 0.1–0.3. The ground on which the comparison stands is not entirely stable either.

Spatial proximity is only one kind of proximity

What ICAR expresses is geographical adjacency alone. Dialect spreads along road networks, school catchments, marriage circles and the barriers of terrain; it is not settled by straight-line distance. A 3 km circle is only a proxy, letting distance stand in for those routes. This is the same situation as the point made on the Mantel test page (in Japanese), that spatial distance is a designed object rather than a measurement.

7Sources and data

The subject is a three-generation dialect vocabulary survey conducted in settlements running from Maibara in Shiga Prefecture to Sekigahara in Gifu Prefecture. The original survey was carried out in 2015 by Habiro Net, a local citizens’ group, and this research is a secondary analysis of its anonymized data. The construction of the adjacency matrix and the estimation of the spatial models are written in R/03_brms_car_template.R, R/09_brms_diversity.R and R/12_brms_loss.R. The paper is in preparation and was not published at the time of writing.

References on the method

  • Besag, J., York, J. & Mollié, A. (1991) Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics, 43(1), 1–20. doi:10.1007/BF00116466
  • Bürkner, P.-C. (2017) brms: An R package for Bayesian multilevel models using Stan. Journal of Statistical Software, 80(1), 1–28. doi:10.18637/jss.v080.i01
  • Carpenter, B., Gelman, A., Hoffman, M. D. et al. (2017) Stan: A probabilistic programming language. Journal of Statistical Software, 76(1), 1–32. doi:10.18637/jss.v076.i01
  • Ferrari, S. & Cribari-Neto, F. (2004) Beta regression for modelling rates and proportions. Journal of Applied Statistics, 31(7), 799–815.

Which perspective this method belongs to

ICAR writes knowledge about spatial proximity into the model in the form of a prior. The term penalizing differences between neighbours is itself the log density of the prior, and estimation proceeds from the posterior formed by combining the data with it. It belongs to perspective 3 (Bayesian) in Seven Perspectives in Equations. The Bayesian estimation page in the same perspective covers the estimation framework as a whole; this page covers how spatial structure is written within it.

This page sets out, for readers in the same field, a method Onohara uses in actual research. Every figure is quoted from the run logs and output files of the analysis scripts; no estimated or approximate values are given.
The whole picture of where each method sits is in Seven Perspectives in Equations, and which teaching session falls under which perspective is in Methods & Mathematics.