library(rstan)
model_code1="
functions{
//defined survival
vector log_s(vector t, real shape1, vector shape2){
vector[num_elements(t)] log_s;
for(i in 1:num_elements(t)){
log_s[i]=log(1-(1-(1+t[i])^(-shape2[i]) )^(shape1));
//log_s[i]=log(((1+t[i])^(scale[i]*shape) - ((1+t[i])^(scale[i]) -1)^shape )/((1 // +t[i])^(scale[i]*shape)));
}
return log_s;
}
//define log_ft
vector log_ft(vector t, real shape1, vector shape2){
vector[num_elements(t)] log_ft;
for(i in 1:num_elements(t)){
log_ft[i]=log((shape1 * shape2[i]) * (1-(1+t[i])^(shape2[i]))^(shape1-1) * (1+t[i])^-(shape2[i]+1));
}
return log_ft;
}
//define log hazard
vector log_h(vector t, real shape1, vector shape2){
vector[num_elements(t)] log_h;
vector[num_elements(t)] logft;
vector[num_elements(t)] logs;
logft=log_ft(t,shape1,shape2);
logs=log_s(t,shape1,shape2);
log_h=logft-logs;
return log_h;
}
//define the sampling distribution
real surv_expP_lpdf(vector t, vector d, real shape1 , vector shape2){
vector[num_elements(t)] log_lik;
real prob;
log_lik=d .* log_h(t,shape1,shape2)+log_s(t,shape1,shape2);
prob=sum(log_lik);
return prob;
}
}
//data block
data {
int N; // number of observations
vector<lower=0>[N] y; // observed times
vector<lower=0,upper=1>[N] censor;//censoring indicator (1=observed, 0=censored)
int M; // number of covariates
matrix[N, M] x; // matrix of covariates (with N rows and M columns)
}
parameters {
vector[M] beta; // Coefficients in the linear predictor (including intercept)
real<lower=0> shape1; // shape parameter
}
transformed parameters {
vector[N] linpred;
vector[N] shape2;
linpred = x * beta;
for (i in 1:N) {
shape2[i] = exp(linpred[i]);
}
}
model {
shape1 ~ uniform(0, 10);
beta ~ normal(0, 1);
y ~ surv_expP (censor, shape1, shape2);
}
generated quantities{
real dev;
dev=0;
dev=dev + (-2)*surv_expP_lpdf (y|censor,shape1,shape2);
}
"
data creation in R
y<-c(23,47,69,70,71,100,101,148,181,198,208,212,224,5,8,10,13,18,24,26,26,31,35,40,41,48,50,59,61,68,71,76,105,107,109,113,116,118,143,154,162,188,212,217,225)
x1<-c(rep(0,13), rep(1,32))
censor<-c(rep(1,3),rep(0,4),rep(1,2),rep(0,4),rep(1,18),rep(0,4),1,0,1,1,rep(0,6))
x <- cbind(1,x1)
N = nrow(x)
M = ncol(x)
event=censor
dat <- list( y=y, x=x, event=event, N=N, M=M)
M1<-stan(model_code=model_code1,data=dat,init=list(list(beta=c(-0.6,0.5))),iter=1000,chains=1)
I find this error!!!
Log probability evaluates to log(0), i.e. negative infinity.
Stan can’t start sampling from this initial value.
Rejecting initial value:
Log probability evaluates to log(0), i.e. negative infinity.
Stan can’t start sampling from this initial value.
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”