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)]
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.
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.
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?).
Hi，could you please let me know if you have learned how to build spatiotemporal models in Stan? If so, could you share your insights? I’ve been trying to work on spatiotemporal models in Stan recently, but I’ve encountered some difficulties.
Hi @Claire_Francis, I have some work on this topic that I’m getting ready to share - I can probably have at least something preliminary to post within the next couple days. I’ll share a link here when its up
Could you share what kind of difficulties? Is there a density definition somewhere you’re trying to fit, or is the issue formulating the spatio-temporal model? I was going to recommend a paper by Leonhard Held, but I can’t find the paper and the only ones I can find are paywalled. The simplest thing to do is just additive spatial and temporal effects. The main issue is that you need to reduce the degrees of freedom for identifiability, either by pinning a value or a sum-to-zero constraint on coefficients. This gets harder when you interact spatial and temporal effects, especially if you also include spatial and temporal effects and maybe heterogeneous effects.
We have case studies on how to implement the ICAR and BYM2 models:
Hi @Bob_Carpenter,
Thank you so much for sharing those materials! They’ve been a great help in my learning of spatial modeling in Stan. As a beginner in Stan, I’m now eager to work on spatio-temporal models. I’m wondering how to organize the data and how to add simple time effects in the spatial model. Do you happen to have any examples of spatio-temporal modeling that could guide me through this process? Your expertise and guidance would be extremely valuable. Thanks again!
this BYM2 variant adds a temporal trend and a spatial trend and ordinary random effects. the amount of spatial, temporal, and random variance is accounted for by parameter \rho, which is given a dirichlet distribution.
@Claire_Francis Here’s a tutorial on space-time modeling using Stan. I’m not sure how helpful some of the descriptive text is but I’m trying to communicate with those who haven’t seen it all before.
(Sorry the web page is now loading really slowly, due to mathjax, I’ll try to fix it)
Hi @Claire_Francis, my apology for the late response! Unfortunately, I was unable to build such model in Stan. I stuck with R-INLA and JAGS instead as I’m still learning Stan. But thank you so much for posting a message and resurrecting this thread!
Hi @mitzimorris, thank you very much for sharing this paper in the thread. This is much appreciated - I also look forward to learning and gaining insights from it!
Hi @cmcd, many thanks for sharing this online resource in the thread as well! It looks fantastic - I also look forward to learning this and gaining insight on how to pull this off in Stan. Thanks!
I’m adding a link to this other old post by @mark_payne – he was asking about space-time models that involve a series of (proper) CAR models, one per time point, combined with a temporal auto-regression. If I’m following his description correctly, that’s the model I implement in this tutorial (a CAR-AR model)
