Running a hierachical spatio-temporal geostaistical model in STAN

Modeling
,
1

Dear STAN groups

I am a STAN beginner, so I would like to know whether STAN can handle a hierachical spatio temporal model.

y\left ( s,t \right )=z\left ( x,t \right )\beta+\xi\left ( x,t \right )+\epsilon
\xi\left ( x,t \right )=a\xi\left ( x,t-1 \right )+\omega\left ( x,t \right )
\omega follows a multivariate gaussian distribution with zero mean and a matern defined covriance matrix.

It is clear that this model can be done in INLA (grouped AR1 model), but I intend to run this model through a sampling approach rather than an approximation.

Could STAN handle this model?
More clearly, I do not know how to deal with the second stage of this model, how I multiply the separatable covariance matrices.

Best regards

2

(Just FYI, LaTeX now works on the discourse (!!!), so you can wrap equations in $s to get it to render maths

In the case where the noise is Gaussian, the INLA method is exact up to the integration over the hyperparameters, so the solution you get there will probably be as good (if not better) than the one you get out of Stan.

Stan is a programming language, so you can just program up the model in the same way you wrote it:

 model {
      for (t in 2:T)
        xi[n] ~ multi_normal_c( a * xi[t-1], K);
}

where K is the covariance matrix of the GP. You’ll probably have ot do some tricks (like the noncentering described in the GP chapter of the manual) to make it run well.

1 Like
3

OK， I see. thank you Daniel

4

I am running the same model, but the chains are not mixed. Is there any suggestion on the code?

5

Hello, could you share your code here? Maybe I can help you find whether the problem comes from data or coding.

7

I am trying to run a similar model in stan but no luck. I am not sure if multi_normal_c function is still available on stan. Is it a typo or what? Also, is xi[n] a typo?

8

I came up with this code but it is quite slow. Can someone help me to check if I am doing something wrong which makes it slow?

data {
  int<lower=1> N;  // number of the entire data
  int<lower=1> Nt;  // number of time point
  int<lower=1> Nloc; // number of locations or dimension of space
  int y[N];    // observation
  real x1[Nloc]; // location x
  real x2[Nloc];  // location y
  vector[N] d;  // covariate
}
  
  
// transform data
transformed data {
  row_vector[2] x_tot[Nloc];
  for (n in 1:Nloc) {
        x_tot[n, 1] = x1[n];
        x_tot[n, 2] = x2[n];
  }
}

.
parameters {
  real<lower=0> length_scale; // scale parameter
  real<lower=0> alpha; // variance 
  real beta; // intercept
  real beta1;  // regresion coeffcient
  real<lower=0, upper=1> rho;  // temporal correlation parameter
  matrix[Nloc, Nt] S_x_t; // spatio-temporal process
}

transformed parameters {
    matrix[Nloc, Nloc] K;  // the covariance matrix
    vector[Nloc * Nt] flattened_S_x_t; // the flattened process
    K = cov_exp_quad(x_tot, alpha, length_scale);
    for (n in 1:Nloc) K[n, n] = K[n, n] + 1e-12;  
    flattened_S_x_t = to_vector(S_x_t);
}


model {
  length_scale ~ gamma(2, 20);
  alpha ~ normal(0, 1);
  beta ~ std_normal();
  beta1 ~ std_normal();
  rho ~ beta(0.5, 0.5);
  S_x_t[,1] ~ multi_normal(rep_vector(0, Nloc), K);
  for (t in 2:Nt)  S_x_t[,t] ~ multi_normal(rho*S_x_t[,t-1], K);
  y ~ poisson_log(beta + d*beta1 +  flattened_S_x_t);
}

EDIT: @maxbiostat has edited this post for legibilty.