Dear all, I am trying to fit a spatial model where parameters exhibit regime change behavior. The model itself is rather long and the data is several MB, so I will post the model component.
- The Spatial component model is based on the paper by Morris et al 2019 (ICAR)
- The regime switching model is based on the explanation of Jim Savage Regime-switching models in Stan
The idea is to investigate whether through bi-weekly time components we can observe changes in the covariate coefficients while controlling for spatial effects.
When we run temporal-spatial model only, the model runs.
When we run the regime change model only, the model runs.
The problem is more than likely due to f[t] which evaluate to 0 at time period 1 due since it represents the likelihood and that of course goes to 0 which means that the probabilities of the transitional matrix are not being calculated. At run time I end up getting the message, that the likelihood can not be negative infinity which is not surprising due to what happens to f[t].
I have tried the log_sum_exp but did not help, however feel free to suggest it as perhaps I am using it wrong.
The changepoint model (the double looped version for K=3) is running
However I would like to have it work for the regime switching models as well. Any help is much appreciated.
functions {
real icar_normal_lpdf(row_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);
}
}
data {
int<lower=1> T;//number of time periods
int<lower=1> K;//number of regimes for Dirichlet
vector [K] iota;
// real<lower = 0> beta[K]; //priors for Dirichlet
int<lower=0> N;//number of spatial areas
int<lower=0> N_edges;
int<lower=1, upper=N> node1[N_edges]; // node1[i], node2[i] neighbors
int<lower=1, upper=N> node2[N_edges]; // node1[i] < node2[i]
int<lower=0> y[T,N]; // count outcomes
real<lower=0> x[3,N]; // covariates
// int <lower=0>
matrix [T,N] log_E;
real<lower=0> scaling_factor; // scales the variance of the spatial effects
}
transformed data {
real log_unif;
log_unif = -log(T);
}
parameters {
simplex[K] xiinit;
real logitrho;
vector<lower = 0, upper = 1>[4] p;
real<lower = -3, upper=3> rho1;
real<lower = -3, upper=3> rho2;
real<lower = -3, upper=3> rho3;
real<lower = -3, upper=3> rho4;
real<lower = -3, upper=3> rho5;
real<lower = -3, upper=3> rho6;
vector<lower = 0>[3] sigma;
row_vector[N] phi;
row_vector[N] theta;
real beta00;
real beta01;
real beta02;
real beta03;
}
transformed parameters {
matrix[3,N] alpha [T] ;
matrix[3,N] eta[T];
matrix[T, 3] xi;
real f[T];
real<lower=0,upper=1> rho = inv_logit(logitrho);
row_vector[N] convolved_re = sqrt(rho / scaling_factor) * phi + sqrt(1 - rho) * theta;
for(t in 12:T) {
if(t==12){
for(n in 1:N){
alpha[t][1,n]=exp(log_E[t,n]+beta00+convolved_re[n]*sigma[1]);//
alpha[t][2,n]=exp(log_E[t,n]+beta00+convolved_re[n]*sigma[1]);//
alpha[t][3,n]=exp(log_E[t,n]+beta00+convolved_re[n]*sigma[1]);//
eta[t][1,n] =(poisson_lpmf(y[t,n]| alpha[t][1,n]));
eta[t][2,n] =(poisson_lpmf(y[t,n]| alpha[t][2,n] ));
eta[t][3,n] =(poisson_lpmf(y[t,n]| alpha[t][3,n] ));
}
f[t] = p[1]*xiinit[1]*exp(sum(eta[t][1,1:N])) + // stay in state 1
(1 - p[1])*xiinit[2]*exp(sum(eta[t][2,1:N])) + // transition from 1 to 2
p[2]*(xiinit[2])*exp(sum(eta[t][2,1:N])) + // stay in state 2
p[3]*(xiinit[2])*exp(sum(eta[t][1,1:N])) + // transition from 2 to 1
(1-p[2]-p[3])*exp(sum(eta[t][3,1:N])) + //transition from 2 to 3
p[4]*(1-xiinit[1] -xiinit[2] )*exp(sum(eta[t][2,1:N])) + //transition from 3 to 2
(1-p[4])*(1-xiinit[1] -xiinit[2] )*exp(sum(eta[t][3,1:N])) ; //transition from 3 to 3
xi[t,1] = (p[1]*xiinit[1]*exp(sum(eta[t][1,1:N])) +p[3]*(xiinit[2])*exp(sum(eta[t][1,1:N])))/f[t];
xi[t,2] = ((1 - p[1])*xiinit[2]*exp(sum(eta[t][2,1:N]))+ p[2]*(xiinit[2])*exp(sum(eta[t][2,1:N])) + p[4]*(1-xiinit[1] -xiinit[2] )*exp(sum(eta[t][2,1:N]))) /f[t] ;
xi[t,3] = 1-xi[t,1]-xi[t,2];
}
else {
for(n in 1:N){
alpha[t][1,n] =exp(log_E[t,n]+beta01+convolved_re[n]*sigma[1]+rho1*(y[t-1,n])+rho4*(x[1,n]));//
alpha[t][2,n] =exp(log_E[t,n]+beta01+convolved_re[n]*sigma[1]+rho2*(y[t-1,n])+rho5*(x[1,n]));//
alpha[t][3,n] =exp(log_E[t,n]+beta01+convolved_re[n]*sigma[1]+rho3*(y[t-1,n])+rho6*(x[1,n]));//
eta[t][1,n] =(poisson_lpmf(y[t,n]| alpha[t][1,n]));
eta[t][2,n] =(poisson_lpmf(y[t,n]| alpha[t][2,n]));
eta[t][3,n] =(poisson_lpmf(y[t,n]| alpha[t][3,n]));
print(eta);
}
f[t] = p[1]*xi[t-1,1]*exp(sum(eta[t][1,1:N])) + // stay in state 1
(1 - p[1])*xi[t-1,1]*exp(sum(eta[t][2,1:N])) + // transition from 1 to 2
p[2]*xi[t-1,2]*exp(sum(eta[t][2,1:N])) + // stay in state 2
p[3]*(xi[t-1,2])*exp(sum(eta[t][1,1:N])) + // transition from 2 to 1
(1-p[2]-p[3])*xi[t-1,2]*exp(sum(eta[t][3,1:N])) + //transition from 2 to 3
p[4]*(1-xi[t-1,1] -xi[t-1,2] )*exp(sum(eta[t][2,1:N])) + //transition from 3 to 2
(1-p[4])*(1-xi[t-1,1] -xi[t-1,2])*exp(sum(eta[t][3,1:N])) ; //transition from 3 to 3
xi[t,1] = (p[1]*xi[t-1,1]*exp(sum(eta[t][1,1:N])) +p[3]*(xi[t-1,2])*exp(sum(eta[t][1,1:N])) )/f[t];
xi[t,2] = ((1 - p[1])*xi[t-1,2]*exp(sum(eta[t][2,1:N])) + p[2]*(xi[t-1,2])*exp(sum(eta[t][2,1:N])) + p[4]*(1-xi[t-1,1]-xi[t-1,2] )*exp(sum(eta[t][2,1:N])) ) /f[t] ;
xi[t,3] = 1-xi[t,1]-xi[t,2];
}
}
}
model {
// priors
p ~ beta(10, 2);
rho ~ normal(0, 1);
rho1 ~ normal(0, 1);
rho2 ~ normal(0, 1);
rho3 ~ normal(0, 1);
rho4 ~ normal(0, 1);
rho5 ~ normal(0, 1);
rho6 ~ normal(0, 1);
sigma ~ cauchy(0, 1);
xiinit~dirichlet(iota);
target += sum(log(f));
}