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:

[Download]

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.