Hello everyone

First of all I am not a statistician but a simple PhD candidate in geography who use statistics and try to understand it. So feel free to correct my statistical vocabulary!

I am actually learning how to fit a model accounting for spatial correlation in Stan and my model seems very inefficient in accounting for it. I am using the meuse dataset as my toy dataset.

I will go through the steps I have been doing for the analysis - so you can understand better and eventually correct:

First I fitted a simple linear model which didn’t account for spatial correlation. A simple

lm(meuse$log_zinc ~ 1 + meuse$dist)

I checked if there was any correlation within the residuals by fitting a semi-variogram. Here is the result:

Form this graph I deduce that the range might be somewhere around 1000.

Because I wanted to compare the results with a software I know better I fit the model in INLA and the SPDE. With INLA I get great results, the standard error are 2 times bigger and when plotting the semivariogram on the residuals I do not notice any spatial correlation:

Because I am trying to switch from INLA to Stan (more flexible and help me to understand better what I am really doing) I try to fit the same model in Stan. I found a great tutorial of how to approximate the Gaussian Random Field: https://mc-stan.org/users/documentation/case-studies/nngp.html and try to use the code provided to fit my own model.

Note that in the Stan code I changed the correlation structure from

```
iNNdistM[j, k] = exp(- phi * NN_distM[(i - 1), h]);
```

to

```
iNNdistM[j, k] = exp(- NN_distM[(i - 1), h] / phi);
```

Because I have no clue about how to set a prior for phi I set phi = Normal(300,200), don’t blame me. I tried a lot of different priors but nothing conclusive (I guess because I tried randomly)

After running the model I compute the semi variogram and obtain this:

As I said I am not an expert but this semi-variogram doesn’t look good at all. What worries me more is that the standard deviation of the coefficients are almost the same as the non-spatial model. I guess if the model accounted properly for the spatial correlation I should have standard deviation more or less like the INLA ones.

If you could provide some insights I would be very grateful!

Here is the stan code for the model. I also join the script for fitting the model in INLA and in Stan.

```
functions{
real nngp_w_lpdf(vector w, real sigmasq, real phi, matrix NN_dist,
matrix NN_distM, int[,] NN_ind, int N, int M){
vector[N] V;
vector[N] I_Aw = w;
int dim;
int h;
for (i in 2:N) {
matrix[ i < (M + 1)? (i - 1) : M, i < (M + 1)? (i - 1): M]
iNNdistM;
matrix[ i < (M + 1)? (i - 1) : M, i < (M + 1)? (i - 1): M]
iNNCholL;
vector[ i < (M + 1)? (i - 1) : M] iNNcorr;
vector[ i < (M + 1)? (i - 1) : M] v;
row_vector[i < (M + 1)? (i - 1) : M] v2;
dim = (i < (M + 1))? (i - 1) : M;
if(dim == 1){iNNdistM[1, 1] = 1;}
else{
h = 0;
for (j in 1:(dim - 1)){
for (k in (j + 1):dim){
h = h + 1;
iNNdistM[j, k] = exp(- NN_distM[(i - 1), h] / phi);
iNNdistM[k, j] = iNNdistM[j, k];
}
}
for(j in 1:dim){
iNNdistM[j, j] = 1;
}
}
iNNCholL = cholesky_decompose(iNNdistM);
iNNcorr = to_vector(exp(- phi * NN_dist[(i - 1), 1:dim]));
v = mdivide_left_tri_low(iNNCholL, iNNcorr);
V[i] = 1 - dot_self(v);
v2 = mdivide_right_tri_low(v', iNNCholL);
I_Aw[i] = I_Aw[i] - v2 * w[NN_ind[(i - 1), 1:dim]];
}
V[1] = 1;
return - 0.5 * ( 1 / sigmasq * dot_product(I_Aw, (I_Aw ./ V)) +
sum(log(V)) + N * log(sigmasq));
}
}
data {
int<lower=1> N;
int<lower=1> M;
int<lower=1> P;
vector[N] Y;
matrix[N, P] X; // design matrix
int NN_ind[N - 1, M];
matrix[N - 1, M] NN_dist;
matrix[N - 1, (M * (M - 1) / 2)] NN_distM; // Distance between the neighbours
}
parameters{
vector[P] beta;
real<lower = 0> sigma;
real<lower = 0> tau;
real<lower = 0> phi;
vector[N] w;
}
model{
beta ~ normal(0,5);
phi ~ normal(1000, 200);
sigma ~ inv_gamma(2, 1);
tau ~ inv_gamma(2, 1);
w ~ nngp_w(sigma, phi, NN_dist, NN_distM, NN_ind, N, M);
Y ~ normal(X * beta + w, tau);
}
generated quantities{
real y_rep[N];
for(n in 1:N){
y_rep[n] = normal_rng(X[n] * beta + w[n], tau);
}
}
```

instead of

model{

vector[N] mu = alpha+beta*x;

y~normal(mu,sigma);

}

meuse_SPDE.R (3.2 KB)

meuse_NNGP.R (1.9 KB)

Helper_functions.R (6.1 KB)