Expanding the ICAR model with smoothing component to a spatiotemporal framework

Dear all,

I hope this message finds you well. Sorry to disturb you folks with the frequent posting of questions. I’m making the transition from R-INLA to RStan and getting to grips with the coding.

[Problem]: I have implemented a Spatial ICAR model (with the ICAR component only [see incredibly helpful tutorials [HERE] (and insightful paper [HERE]) by author(s) Morris et al, 2019] which is a lot easier than the version with the structured and unstructured errors) to cross-sectional data on road accidents in the UK to estimate the area-specific relative risks and probability exceedance (i.e., RR > 1). I have shared a full reproducible example of the Stan & R-code and dataset (.csv & .shp) to churn the outputs you see below:

Image result [RR, Significant areas, Exceedance probabilities (RR>1)]

Stan code:

functions {
	real icar_normal_lpdf(vector phi, int N, array[] int node1, array[] int node2) {
		return -0.5 * dot_self(phi[node1] - phi[node2]);
data {
	int<lower=0> N;
	int<lower=0> N_edges;
	array[N_edges] int<lower=1, upper=N> node1;
	array[N_edges] int<lower=1, upper=N> node2;
	array[N] int<lower=0> Y;                                         // dependent variable i.e., number of road accidents
	vector<lower=0>[N] X;                                            // independent variable i.e., deprivation score
	vector<lower=0>[N] E;                                            // estimated number of expected cases of road accidents
transformed data {
	vector[N] log_offset = log(E);                                  // use the expected cases as an offset and add to the regression model
parameters {
	real alpha;                                                     // define the intercept (overall risk in population)
	real beta;                                                      // define the coefficient for the deprivation score variable
	real<lower=0> sigma;                                            // define the overall standard deviation producted with spatial effect smoothing term phi
	vector[N] phi;                                                  // spatial effect smoothing term or spatial ICAR component of the model 
model {
	phi ~ icar_normal(N, node1, node2);
	Y ~ poisson_log(log_offset + alpha + beta*X + phi*sigma);       // likelihood function i.e., spatial ICAR model
	alpha ~ normal(0.0, 1.0);                                       // prior for intercept   (weak/un- informative prior)
	beta ~ normal(0.0, 1.0);                                        // prior for coefficient (weak/un- informative prior)
	sigma ~ normal(0.0, 1.0);                                       // prior for SD          (weak/un- informative prior)
	sum(phi) ~ normal(0, 0.001*N);
generated quantities {
	vector[N] eta = alpha + beta*X + phi*sigma;                     // exponentiate only the alpha + beta*X + phi*sigma to get the estimated relative risk of areas 
	vector[N] mu = exp(eta);                                        // derive areas-specific relative risk ratios (RRs)
// end of script: Spatial ICAR model with 'ICAR' component only for spatial smoothing of risks

Here is access to the R-codes for compiling the Stan script and map production, and datasets:

I have longitudinal dataset of the road accidents from 2015 to 2020 and would like to implement a spatiotemporal version of the above model. Performing this in R-INLA is quite easy but I would like to use Stan. Is there a way of expanding the shared code into a spatiotemporal version? Any guidance or help would be greatly appreciated.

PS: This work is basically for an MSc computer tutorial. Last year, I taught the same thing with these materials in R-INLA. I am migrating all the content from R-INLA to RStan.

That’s what this website is for! :)

Can you describe the spatiotemporal model that you want to fit? Seperable spatiotemporal models, where the spatial component and the temporal components are represented with separate terms that decompose additively, are relatively easy extensions to your existing code (i.e. add a temporal term to your existing model). Complexity increases from models that look like y ~ spatial_term + temporal_term, through y ~ spatial_term + temporal_term[spatial_unit] and y ~ spatial_term[temporal_unit] + temporal_term, and up through y ~ nonseparable_spatiotemporal_term.

1 Like

Hello Jacob, thank you very much for the reply. Yes, the spatiotemporal model that I trying to implement has this structure (it’s Poisson) with a parametric trend, it’s a null model that I want to use where u_i is the spatial effect and \tau is the overall time trend and \delta_i is the area-specific time trend; and the interaction with the t causes the whole thing to vary through time and geographic areas.

log(\lambda_{it}) = \alpha+u_i + (\tau+\delta_i)\cdot t

(Sorry, my maths here is incredibly wishy-washy…)

The above model is a separable. In R-INLA, I last used this model code on the shared road accidents data to compute the area-time-specific relative risks:

model.formula  <- casualties ~ f(OBJECTID, model="bym", graph=g) + f(OBJECTID2, idyear, model="iid") + idyear
risk.estimates <- inla(model.formula, family = "poisson", data = analysis_data_temp, E = expected, control.predictor = list(compute = TRUE))

Honestly, not sure how to go about coding this in Stan. I will have a go at this but please I’m also open to any examples and suggestions you may have - which will be much appreciated.

I don’t know INLA syntax well, but others here might. I might be able to help further if you’re able to describe what the temporal model is. The first question is whether idyear gets treated as continuous, categorical, or something else (ordinal?).

1 Like

Thank you again for the reply.

Yes to respond to your first question, the idyear was indeed treated as ordinal numbers from 1 to 6 (where 1 represents 2015,… and 6 is 2020).