How does Markov random field (bs=mrf) in mgvc gam handle repeated measures on the spatial units?

Question

I am attempting a spatio-temporal model in mgcv gam.

I am using a factor smooth to define each of 27 areal units in a shapefile ("id") as subjects (essentially) which have undergone 23 repeated measurements in time ("year").

I am using the Markov random field to define the neighbor relationships of the area units.

I have combed through Wood (2017) and other resources and have not found an answer to my question, which is: how is the mrf handling the repeated measures? My model is:

model <- gam(response) ~ s(year) +
                         ti(year,id, bs="fs", m=1) +
                         s(id, bs = 'mrf', xt = list(nb = neighbors), k=5),
                         family=tw(), data=data)

Thank you.

Answer 1

The neighbors object defines, for the levels of id, which other locations are neighbours of the $i$th id. Your data object contains a an id for each observation. The MRF defines a spatial effect such that the effect of the MRF for any unit varies smoothly over the neighbours of that unit. This implies that the MRF is smooth over neighbouring units of those neighbours and so on.

How smooth the MRF is will depend on the rank of the MRF; if you estimate a low-rank MRF (use fewer knots than levels of id - 1) the field will be smooth, increasingly so the fewer the knots. If you fit a full-rank MRF then you will get a "rougher", with effectively one coefficient per unit.

How this plays out in the model, and you can see this if you use predict() with type = "terms", is that for a particular id the MRF smooth contributes the same value/effect to the fitted value for all of the 23 observations for that particular id.

You will get a similar effect if you put a fixed effect factor into the model; all the observations from id == 1 will get a contribution to the fitted value for belonging to that level of id (typically the intercept plus the coefficient for the $i$th level of id if observation is not from the reference level of the factor). Similarly, if you included a factor random effect; the fitted value for each observation is pulled up or down relative to the model intercept by some amount given by the "estimate" of the random effect for the unit that it belongs to.

The difference between the factor fixed effect or random effect and the MRF is that the MRF takes into account the spatial arrangement (or any other arrangement you can represent as a graph) of the levels of the factor.

Specifically for your model, the low-rank MRF effect is constant over the year observations. The factor smooth interaction given by the ti(year, id, bs = 'fs', m = 1) (this could really be s(year, id, bs = 'fs', m = 1)) implies a factor random effect (random intercept) plus a random smoother for year for each observation. The random intercept component of this smoother will account for non-spatial differences in the expected values due to the units.

If you wanted the MRF to vary over time, then you'd need to include it as a tensor product with the year term, for example:

te(year, id, bs = c("tp", "mrf"),
   xt = list(year = NULL, id = list(nb = neighbors)))

This doesn't include the random year smooths, which if you just wanted to estimate the smooths and not have the global smooth of year would be

te(year, id, bs = c("fs", "mrf"),
   xt = list(year = NULL, id = list(nb = neighbors)))

Doing the global smooth plus difference smooth of year and the MRF is trickier.