Updates to case study on spatial models for areal data - Poisson CAR/IAR BYM and BYM2

I’ve updated the case study on ICAR models for spatial smoothing:

instead of using a hard sum-to-zero constraint on the ICAR component, I’m using a soft-centering via a prior that (strongly) constrains the mean to be zero, i.e., the ICAR component is coded up as:

data {
  int<lower=0> N;
  int<lower=0> N_edges;
  int<lower=1, upper=N> node1[N_edges];  // node1[i] adjacent to node2[i]
  int<lower=1, upper=N> node2[N_edges];  // and node1[i] < node2[i]
parameters {
  vector[N] phi;
model {
  target += -0.5 * dot_self(phi[node1] - phi[node2]);
  // soft sum-to-zero constraint on phi)
  sum(phi) ~ normal(0, 0.01 * N);  // equivalent to mean(phi) ~ normal(0,0.01)

So much for what I presented at StanCon.

In this, as in everything else, I’m lagging the work done by @MigueBeneito, who rightly noted that the ICAR models are very slow to fit in Stan - I’ve added Miguel and Paqui’s names to the Acknowledgements section. In the previous thread on this topic - Case study on spatial models for areal data - Poisson CAR/IAR - comment 59 or thereabouts - Miguel proposed this solution. As did Andrew - during our initial discussions back in April one of the things he said was: “sum-to-zero is a terrible constraint - a soft constraint would be better”. At the time I didn’t understand how to do that - heck, at the time, I understood practically nothing - Miguel totally nailed the Stan implementation long before I did.

the soft-centering is about twice as fast - that’s the good news.

the bad news is that the BYM2 model is still a difficult model to fit and will be until we get more INLA-style matrix magic from @anon75146577 et al.

the other change to the case study is much slicker maps - which is a topic for its own case study or a blog post on spatial neighborhoods and how to get there from here.