OK spatial statisticians, need your insights here.
Fixed all indexing problems with the BYM2_islands.stan program, and have put together an Python notebook (running CmdStanPy) to fit the models and data, along with a bunch of R helper functions and scripts to do all the indexing for the spatial data, plus use INLA’s secret sauce for computing the scaling factors.
I am unable to get the same estimates as INLA on the Scotland data, per the Freni-Sterrantino et al paper - [1705.04854] A note on intrinsic Conditional Autoregressive models for disconnected graphs. - the model does fit the data quickly, without any diagnostic errors; the estimates are close, but as they say in Scotland, a wee bit off. This is a dataset with one connected component and 3 islands. I remain mystified the recommendations in section 4 of Freni-Sterrantino - I tried writing it up in my own words as part of the writeup, and that let to a bunch of experiments but no real insights.
The scaling factor for the singleton comes in as data and is computed in file bym2_helpers.R.
The ICAR prior is computed in the Stan program, according to Dan Simpson, it is doing the right thing.
This should work no problem for areal data where you’ve got two different components - it’s singletons that are problematic.