Hello all!

I’m trying to fit recurrent events with a Bayesian survival model with frailty term which specifies the baseline hazard as Weibull distribution. But currently I meet a problem as the title. My STAN code is shown as following.

data {

int<lower=0> N; // number of subject

int<lower=0> N_row; // number of rows in the dataset

int<lower=0> K; // number of observed events

int<lower=0> Cn; // number of censored observations

int ncov; // number of covariates

matrix[K,ncov] Xobs; // covariate matrix for observed event

matrix[Cn,ncov] Xcen; // covariate matrix for censored event

vector[K] tobs; // time for observed event

vector[Cn] tcen; // time for censored event

int<lower=1,upper=N> IDobs[K];

int<lower=1,upper=N> IDcen[Cn];

}

parameters {

real<lower = 0> nu; // nu = shape

real<lower = 0> theta; // lambda=scale

real<lower = 0> sigmaw; // variance parameter for frailty term

vector[N] w; // frailty term

real beta1;

real beta2;

}

model {

real lambda_Tmp1;

real lambda_obs;

real lambda_Tmp2;

real lambda_cen;

// prior

nu ~ gamma(0.01,0.01);

theta ~ normal(0,sqrt(1000));

beta1 ~ normal(0,sqrt(1000));

beta2 ~ normal(0,sqrt(1000));

w ~ gamma(1/sigmaw,1/sigmaw); // my guess: the part caused the error message

sigmaw ~ inv_gamma(0.01, 0.01);

// likelihood - observed times:

for (i in 1:K){

lambda_Tmp1 = w[IDobs[K]] + beta1*Xobs[i,1] + beta2*Xobs[i,2] + log(theta);

lambda_obs = exp(-lambda_Tmp1/nu);

target += weibull_lpdf(tobs[i] | nu, lambda_obs);

}

// likelihood - censored times:

for (i in 1:Cn){

lambda_Tmp2 = w[IDcen[Cn]] + beta1*Xcen[i,1] + beta2*Xcen[i,2] + log(theta);

lambda_cen = exp(-lambda_Tmp2/nu);

target += weibull_lccdf(tcen[i] | nu, lambda_cen);

}

}

I don’t understand why it gives error message as below.

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”

Even though I tried giving different initial values, it still doesn’t work. I changed the distribution of theta, it doesn’t work as well.(since I think the log(0) may be caused by theta.) However, interestingly, if I replaced the w ~ gamma(1/sigmaw,1/sigmaw) by w ~ normal(0,sqrt(sigmaw)). The whole code could work, but it will take around 8 to 12 hours to finish running (and a series of w values do not converge).

Could anyone help me with this question? Why does the gamma distribution of frailty term w cause the problem?

Thank you very much in advanced!