How to accelerate spatial models on STAN

Modeling
, ,
1

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

The model reads,
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?

Thanks in advance.
Charles

2

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
3

Hi,

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

Thanks

4

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

5

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

6

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.

7

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
8

Thanks, let me run this code.