Hi, I am trying to do survival analysis, and saw this rstan code posted on Github, so I copied it into R, and used the provided data, but got the initial values rejection error. I have been working on it for 2 weeks, did a lot of googling in my efforts, tried to revise code many times, but still no luck.

Any suggestions are appreciated, thanks!

Here is the data https://github.com/fonnesbeck/stan_workshop_2016/blob/master/code/survival.data.r

the code: https://github.com/fonnesbeck/stan_workshop_2016/blob/master/code/survival.r

Initialization between (-2, 2) failed after 100 attempts.

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â€ť

```
library(ggplot2)
library(plyr)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
cox_model <-
"data {
int<lower=0> Nsubj;
int<lower=1> T;
vector[Nsubj] pscenter;
vector[Nsubj] hhcenter;
int<lower=0,upper=1> ncomact[Nsubj];
int<lower=0,upper=1> rleader[Nsubj];
int<lower=0,upper=1> dleader[Nsubj];
vector[Nsubj] inter1;
vector[Nsubj] inter2;
int<lower=0,upper=1> FAIL[Nsubj];
int<lower=0> obs_t[Nsubj];
int<lower=0> t[T+1];
}
transformed data {
int<lower=0> Y[Nsubj, T];
int<lower=0> dN[Nsubj, T];
# Set up data
for(i in 1:Nsubj) {
for(j in 1:T) {
# risk set = 1 if obs_t >= t
Y[i,j] <- int_step(obs_t[i] - t[j]);
# counting process
dN[i,j] <- Y[i, j] * int_step(t[j + 1] - obs_t[i]) * FAIL[i];
}
}
}
parameters {
vector[7] beta;
real<lower=0> c;
real<lower=0> r;
vector[T] dL0;
}
transformed parameters {
vector[T] mu;
matrix[Nsubj, T] Idt;
vector[T] dL0_star;
for(j in 1:T) {
for(i in 1:Nsubj) {
# Intensity
Idt[i, j] <- Y[i, j] * exp(beta[1]*pscenter[i] + beta[2]*hhcenter[i] + beta[3]*ncomact[i]
+ beta[4]*rleader[i] + beta[5]*dleader[i] + beta[6]*inter1[i] + beta[7]*inter2[i]) * dL0[j];
}
# prior mean hazard
mu[j] <- dL0_star[j] * c;
dL0_star[j] <- r * (t[j + 1] - t[j]);
}
}
model {
for(j in 1:T) {
for(i in 1:Nsubj) {
# Likelihood
dN[i, j] ~ poisson(Idt[i, j]);
}
dL0[j] ~ gamma(mu[j], c);
}
c ~ gamma(0.0001, 0.00001);
r ~ gamma(0.001, 0.0001);
beta ~ normal(0.0, 100000);
}"
data <- read_rdump('survival.data.r')
inits <- list(c1=list(beta=c(-.36,-.26,-.29,-.22,-.61,-9.73,-.23), c=0.01, r=0.01,
dL0=c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)))
fit <- stan(model_code=cox_model, data=data, init=inits, iter=2000, chains=1, verbose=F)
print(fit)
```

â€¦

Edit: Dr. Chris Foonesbeck replied me in the email, and he said he ran the code without any problem when he wrote it, so he suggested me to use an older version of Rstan. and thatâ€™s what I am doing now, I will report the result. I am curious why older version of Rstan works instead of the newest one.