If you have your parameter vector err\_spatial \sim CAR(\mu, \rho, \tau^2), then you’re assigning a multivariate normal probability to them; so it makes sense that those terms will be correlated, but the scatter plots look suspicious to me too.

The model is:

y \sim Gau(\mu, (I - \rho C)^{-1} M)

where C is a connectivity matrix. Its typically really sparse, meaning most entries are zero. But it gets inverted in the covariance matrix, and it becomes a dense matrix. So neighboring areas aren’t the only ones that will be correlated.

If there’s a covariate that has a similar spatial trend across the map as your outcome, then the covariate will be competing for influence with the spatial autocorrelation (SA) term (as intended) and I suppose that could induce additional correlation into the entire set of SA terms, as they expand/contract in response to the influence of the covariate.

If you’re using a model similar to what’s in that case study, where you’re modeling a rare disease or some other small rate in a log-linear Poisson model, then I might wonder why those err\_spatial terms are taking on such large values. It looks like the visible patterns in the scatter plots are dominated by the extreme values (long tails), and we can’t see what’s happening in the bulk of the distribution (where b1\_base and the spatial terms are both near zero). I don’t know if you have any other issues or concerns with the model, but some of the unusual aspects of these scatter plots appear to be a result of that long tail on b1\_base.