Spatio-temporal INLA modelling with disconnected regions: regional vs global spatial variance and model specification

[Pauline Pouchin]

Dear R-INLA community,

I am working on modeling disease counts using the R-INLA package. My goal is to understand the spatio-temporal variation of disease risk in a regionally structured area.

My spatial units are IRIS (small areas), which are grouped into several disconnected regions. Because of this, I model separate spatial structures for each region using BYM2 models defined on the adjacency graph of each region. This approach allows me to assess the spatial variance of risk within each region.

For example, my base model looks like this:

O ~ 1 + f(region_id_f, model = "iid", constr = TRUE) + f(struct_region1, model="bym2", graph=region1.adj, scale.model=TRUE) + ...

I then extend this model with temporal structures (iid, RW1, RW2) and space-time interactions, and I compare models using WAIC.

My questions are:

1. Model selection and regional effect:

   Should I include the regional effect (f(region_id_f, model="iid") for example)  in all models when comparing them using WAIC/DIC, or should I exclude it during model selection and add it afterwards?

2. Fixed vs. random region effect:

   I am interested in whether some regions have higher risk than others.

   • Should I model the region effect as fixed or random?
   • If random, should I impose a sum-to-zero constraint (constr=TRUE) for identifiability?

3. Global spatial variance:

   Since my regions are disconnected, I do not include a spatial effect across the entire study area in the initial model. If I added a BYM2 effect on the full area, it would overlap with the region-specific spatial effects, because IRIS are already nested within regions.

   • How can I then obtain the global spatial variance of risk across the whole study area?

Thank you very much for your advice!

Best regards,