Issues with OpenBUGS and switch bug's codes to Stan

Modeling
1

Hello! I am new to Stan and I am trying to switch the bug’s code of my spatial survival model to Stan. In OpenBUGS, my model is syntactically correct, data loaded, model compiled, initial values generated and model initialized. But, during the updating time of the model (model is updating), I encounter the following error in OpenBUGS:
“Sorry, something went wrong in procedure Stack. Value in module GraphStack”.
To tackle this problem, I am trying to switch the following bug’s code to Stan. Can anyone help me?

"# bug's code
model {
for (k in 1:q) {
for (j in 1:ng) {
landa[k,j] <- kernel [k,j] * exp(x[j])
}
z[k]<- sum(landa[k,])*p[k]
a[k]<- z[k]/sum(z[])
for (i in 1:N) {
#temp: (logt-beta0-beta*x-frailty) / sigma1
temp[k,i]<- (log(time[i]+0.01)-beta0-beta [1]*sex [i]-beta [2]*prison [i]-beta [3]*(age[i]-40)/20-beta [4]*smear.diag[i]-beta[5]*history[i]-W[county[i]]) / sigma
#survival and density function of extreme value distribution
s0[k,i]<-exp(-exp(temp[k,i]))
f0[k,i]<-s0[k,i]*exp(temp[k,i])
#loglikelihood function
L[k,i]<- censor.status[i]*log(f0[k,i])-censor.status[i]* log(sigma)-censor.status[i]*log(time[i]+0.01)+log(s0[k,i])+log(a[k])-censor.status[i]*log(s0[k,i])
zeros[k,i]<-0
zeros[k,i] ~dloglik(L[k,i])
}
}
#prior
inversesigma~ dgamma (1,2)
sigma<- 1/inversesigma
beta0~dnorm(0,0.001)
for(m in 1:5) {beta[m]~ dnorm(0,0.001)}
for(j in 1:ng) {x[j]~dnorm(0,0.04)}
for(i in 1:14) {W[i]~dnorm(0,0.001)}
}
#Initial Value
list(beta0=0,inversesigma=0.10,beta=c(0,0,0,0,0),x=c(0,0,0,0,0,0,0,0))
#Data
list(kernel =structure(.Data = c( 0.5158, 0.5853, 0.6622, 0.7440, 0.5257, 0.5994, 0.6833, 0.7787, 0.7102, 0.7641, 0.7793, 0.7452, 0.7686, 0.8550, 0.8872, 0.8216, 0.9454, 0.8349, 0.7330,0.6431, 0.8957, 0.8156, 0.7227, 0.6367, 0.8717, 0.8607, 0.7817, 0.6941, 0.9494, 0.9234, 0.8098, 0.7102,0.5822, 0.6598, 0.7436, 0.8250, 0.5936, 0.6769, 0.7717, 0.8798, 0.6648,0.7547, 0.8505, 0.9158, 0.6696, 0.7627, 0.8669, 0.9571, 0.9309, 0.9412, 0.8256, 0.7240, 0.8666, 0.8711, 0.7959, 0.7076,0.4144,0.4629,0.5131,0.5622,0.4393,0.4951,0.5553,0.6175,0.3315,0.3694,0.4091,0.4492,0.3537,0.3971,0.4437,0.4924,0.4767,0.5392,0.6073,0.6786,0.4914,0.5594,0.6362,0.7220,0.4178,0.4711,0.5288,0.5892,0.4350,0.4939,0.5596,0.6318,0.7955,0.9026,0.9452,0.8401,0.7850,0.8785,0.9074,0.8259,0.7840,0.7758,0.7242,0.6555,0.8896,0.8715,0.7849,0.6946,0.5022,0.5515,0.5957,0.6276,0.5424,0.6042,0.6636,0.7105), .Dim=c(14, 8)), q=14, ng=8,p=c(0.92,0.87,0.86,0.86,0.88,0.87,0.85,0.94,0.92,0.87,0.82,0.83,0.80,0.91))
list(N=2490)
county[] sex[]  prison[] age[] smear.diag[] history[]  censor.status[] time[]

END "

Best regards,
Eisa

@maxbiostat edited this post for syntax highlighting.

2

I would take a look at the stan user guide (in particular, the transitioning from BUGS section for some concrete differences in the languages and some of the example models to get a handle on the language). There are enough differences that it would probably be a good idea to try to build a few simple Stan models first before tackling something more complex.

One thing to note is that Stan is generally vectorized, so many of the large for loops found in BUGS models are usually unnecessary. Secondly, a lot of function names are different and some distributions are parameterized differently (e.g., the normal distribution is in terms of mean and standard deviation, not mean and precision). Finally, you’ll probably want to read up on recommended prior distributions, as they tend to differ from BUGS standard practices.

2 Likes
3

Thank you so much for your email. Can you guide me in defining the attached model in stan, please?
Actually, I have trouble coding and executing the model and I do not know exactly how to consider the censoring variable (in the codes, censor. status variable indicates right censorship) in the definition of likelihood and model.
Now, I am writing the codes of my spatial survival model in Stan. The likelihood function of my model is attached (where Xj is the white noise and Xj ∼ N (0, sigma) ).
I defined the likelihood of the attached model for exponential distribution as follows in the Stan software, but

"
functions {

real newweibul_log (vector time,vector censor.status,vector age, matrix kernel,vector x,vector p, int n_g,int province, real beta0,real beta, real sigma, vector lam) {
matrix [province,n_g] lambda;
vector [province] z;
vector [province] a;
matrix [province,num_elements(time)] L;
real L_i;

for (k in 1:province) {
  for (j in 1:n_g) {
      lambda [k,j] = kernel [k,j] * exp (x[j]);
        }
        z[k] = sum(lambda[k,]) * p[k];
        a[k] = z[k]/sum(z);
    }
    for (i in 1:num_elements(time)){
      log(lam[i]) = beta0 + beta1 * age[i];
    }
for (k in 1:province) {
  for (i in 1:num_elements(time)) {
    L [k,i] = ((a[k]* lam [i] * exp( - lam [i] * time[i])) ^ censor.status[i]) * ((a[k]*exp(-lam [i]*time[i]))^ (1-censor.status[i]));
  }               
}
L_i = sum(log(L));
return L_i;

}
}
data {
int <lower=0> province;
int <lower=0> n_grid;
int <lower=0> N;
vector <lower=0> [province] p;
vector [N] time;
vector [N] age;
vector [n_grid] x;
matrix [province,n_grid] kernel;
//vector<lower=0, upper=1> [N] censor.status;
}
parameters{
real <lower=0> sigma;
real beta0;
real beta;
// real mu;
real lambdaa;
real a;
}
transformed parameters{
real time_ratio;
time_ratio = exp(beta);
}
model{
for (j in 1:n_grid){
x ~ normal(0,5);
}
beta0 ~ normal( 0, 100);
beta ~ normal( 0, 100);
// real alpha;
//real bettta;
//alpha < - 1;
//bettta < - 1;
sigma ~ gamma(1, 1);
//mu ~ normal( 0, 100);
//real alpha0;
// real bettta0;
// alpha0 < - 1;
// bettta0 < - 1;
lambdaa ~ gamma(1, 1);
for (i in 1:N) {
time[i] ~ newweibul ( beta0, beta, lambdaa[i], a[k],) ;
}
}
"
Regards,
Eisa
Custom likelihood of spatial survival model.pdf (231.0 KB)