# How to accelerate spatial models on STAN

Hello, I’m working on spatial data by a GLMM.

y \sim binomial\left ( n,p \right )
p= logit(\eta)
\eta = \alpha + X\beta + \phi
\phi \sim MVN(0, \left [ \tau_{\mu}\left ( I-\lambda W \right )\left ( I-\lambda {W} \right )' \right ]^{-1} )

The stan code that I use

data{
int<lower=1> N;
int y[N];
int nTrials[N];
int<lower = 1> K;
matrix [N, K] X;
matrix<lower=0, upper=1>[N,N] I;
matrix<lower = 0 >[N, N] W;
}
transformed data{
vector[N] zeros;
zeros = rep_vector(0, N);
}
parameters {
real alpha;
vector[K] beta;
vector[N] phi;
real<lower = 0> tau;
real<lower = -1, upper = 1> lambda;
}
transformed parameters {
vector[N] eta;
eta = alpha + X * beta + phi;
}

model {
phi ~ multi_normal_prec(zeros, tau * crossprod(I - lambda * W));
alpha ~ normal(0,1000);
beta ~ normal(0,1000);
tau ~ gamma(2, 2);
lambda ~ uniform(-1,1);
y ~ binomial_logit(nTrials, eta);
}


It works, but it runs very slow on STAN. I wonder whether there are some ways (i.e. coding tricks) to accelerate this model on STAN?

Charles

hi,
is that the Leroux spatial model? yes, it will be very slow - that’s a problem with the Leroux model.
There’s a great paper by Riebler et al that goes over these models and motivates the BYM2 model proposed by Dan Simpson: https://arxiv.org/abs/1601.01180
Based on that, we’ve developed the BYM2 model in Stan: https://mc-stan.org/users/documentation/case-studies/icar_stan.html
You can also use BRMS, which has implementions for the CAR and ICAR models.

3 Likes

Hi,

It is not the usual CAR model. The variance-covariance matrix follows a SAR structure, rather than a CAR structure.

Thanks

interesting - not sure what to do about SAR models - @anon75146577 might have input?

maybe use BRMS? cf Spatial simultaneous autoregressive (SAR) structures in multilevel models

The cor_errorsar(W) in brm is not identical to the random intercept model (The SAR structure is imposed on the variance-covariance matrix.) mentioned here. You can see that in a cor_errorsar(W) model, the spatial structure is built on the error term, but there is not an error term in binomial regression.

With the caveat that I have not tried this, I would use the non-centred parameterization. Sadly, Stan does not have anything that speeds up the computations based on the Markov property (yet), so there’s a limit to how fast it can be.

transformed data{
vector[N] zeros;
zeros = rep_vector(0, N);
}
parameters {
real alpha;
vector[K] beta;
vector[N] z; // non-centred phi
real<lower = 0> tau_sd;  // Put this as the standard deviation NOT as precision
real<lower = -1, upper = 1> lambda;

}
transformed parameters {
vector[N] eta;
vector[N] phi;
// If I and W are symmetric then you don't need that transpose!
phi = (I - lambda * W)' \ z; // E(phi * phi') = ((I - lambda * W)' * (I - lambda * W))^{-1}
eta = alpha + X * beta + tau_sd * phi;
}

model {
z  ~ normal(0,1);
alpha ~ normal(0, 3);  // sensible priors on the logit scale
beta ~ normal(0, 3); // sensible prior on the logit scale
tau_sd ~ normal(0,3); // NB:  This is a different prior to the one in your model!
lambda ~ uniform(-1,1);
y ~ binomial_logit(nTrials, eta);
}


I also put tighter priors on everything because the linear predictor is on the logistic scale so nothing is going to be outside [-5,5]

1 Like

Thanks, let me run this code.