That’s the ICAR prior - pairwise distance - and it’s not identifiable - you can add any constant you want to sai and it washes out of the pairwise difference. adding a sum-to-zero constraint fixes this problem. we recommend a soft sum-to-zero constraint.
You can code up the ICAR prior as a function in Stan as follows:
functions {
real icar_normal_lpdf(vector phi, int N, int[] node1, int[] node2) {
return -0.5 * dot_self(phi[node1] - phi[node2])
+ normal_lpdf(sum(phi) | 0, 0.001 * N);
}
}
When used in a sampling statement, you call it as follows:
Many thanks for your comments @mitsimorris. Yes, I added it to the codes before but forgot to write it here. soft sum-to-zero constraint on sai in codes is:
it looks like the BYM2 part is OK.
the complexity of your model and the fact that I haven’t yet implemented a spatio-temporal model
make it hard to see what’s wrong.
the standard debugging advice is to simplify - could you set up a basic spatio-temporal regression based on your code and share?
Thanks a lot for your reply @mitsimorris. I wrote the code more simply by distributing gamma. The result is the same as before! The codes:
data {
int<lower=0> N;
int <lower=0> K;
int<lower=1> N_edges;
int<lower=1, upper=N> node1[N_edges];
int<lower=1, upper=N> node2[N_edges];
real<lower=0.00001> scaling_factor;
vector [N] y;
matrix [N, K] X;
int <lower=1> S; // Spaces
int <lower=1> TJ; // Times
}
transformed data {
for (i in 1:N_edges) {
if (node1[i] >= node2[i])
reject("edge pair i=", i, ", node1[i] must be less than node2[i]");
}
}
parameters {
//distribution
real<lower=0> inverse_phi; //the variance parameter
//BYM2 model
vector[N] theta;
vector[N] sai;
real<lower=0> sigma;
real<lower=0, upper=1> rho;
//Reg model
vector[K] beta;
// Random walk
real<lower=0> sigmaRw;
real<lower=0> sigmaR;
}
transformed parameters {
vector[N] convolved_re;
vector<lower=0>[N] mu;
vector[N] gama; //rate parameter for the gamma distribution
vector [TJ] Rw1;
vector [TJ] R2;
convolved_re = sqrt(1 - rho) * theta + sqrt(rho / scaling_factor) * sai;
for (i in 1:S) {
for (j in 1:TJ){
mu[i*TJ-TJ+j]= exp(X[i*TJ-TJ+j, ] * beta + convolved_re[i] * sigma + Rw1[j] + R2[j] );
}
}
gama = rep_vector(inverse_phi, N) ./ mu;
}
model {
inverse_phi ~ exponential(1); // prior on inverse dispersion parameter
beta[1] ~ cauchy(0,10);
beta[2:] ~ cauchy(0,2.5);
theta ~ normal(0, 1);
sigma ~ normal(0, 1);
rho ~ beta(3, 3);
sigmaRw ~ gamma (1, 0.001);
Rw1[1] ~ uniform(0,1) ;
for (j in 2:TJ) {
Rw1[j] ~ normal(Rw1[j-1], 1/sigmaRw);
}
sigmaR ~ gamma (1, 0.001);
R2 ~ normal(0, 1/sigmaR);
target += -0.5 * dot_self(sai[node1] - sai[node2]);
sum(sai) ~ normal(0, 0.001 * N);
y ~ gamma(inverse_phi,gama);
}
the error is at line 49 which is a constraint on mu:
vector<lower=0>[N] mu;
this is the error (as the same):
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: mu[1] is nan, but must be greater than or equal to 0 (in 'model2d587ac1394c_spatialtemporal_gamma' at line 49)
Chain 1:
Chain 1: Initialization between (-2, 2) failed after 100 attempts.
Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."
[1] "error occurred during calling the sampler; sampling not done"
I see that mu has length N but the loops assign to the first S*TJ places. I don’t know if those are equal but the error is what I’d expect if S*TJ=155 and N>155.
Well, 155=31\times 5 so first NaN at index 156 looks a bit like j=6 is failing–but the loops are the other way around so 156 actually corresponds to i=20,j=4 and all i<20 must work for any j.
Great, many thanks @nhuurre. It really helped me find the problem. I had a mistake in specifying data in R before running the code, and set TJ=5! I modified the code to TJ=8, and also I changed constraint on lambda from 0.1 to:
vector<lower=0.000001>[N] lambda;
the error changed to this:
...
target = -inf
target = -inf
target = -inf
target = -inf
target = -inf
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: gama[193] is 6.93183e-007, but must be greater than or equal to 1e-006 (in 'model2d585aa36987_BESAG_temporal_RW' at line 75)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: validate transformed params: lambda[40] is 1.97203e-007, but must be greater than or equal to 1e-006 (in 'model2d585aa36987_BESAG_temporal_RW' at line 73)
Chain 1: target = -inf
target = -inf
target = -inf
target = -inf
target = -inf
target = -inf
target = -inf
...
target = -inf
target = -inf
target = -inf
target = -inf
target = -inf
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1:
Chain 1: Initialization between (-2, 2) failed after 100 attempts.
Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."
[1] "error occurred during calling the sampler; sampling not done"
Lambda, alpha and gama are parameters of the distribution. The problem seems to be the constraint on lambda and gama. The constraint for these parameters are: