Survival analysis: initial values rejected

Modeling
21

For what it’s worth, it seems that, after the two changes mentioned above (swapping the lines assigning to mu[j] and dL0_star[j] and adding lower bound to vector<lower=0>[T] dL0;), the problem with the initialization is with the gradient: Gradient evaluated at the initial value is not finite. I think this is caused by some of the Idt values being zeroes because some of the values in Y are zero, while the parameter of the Poisson distribution should be positive. Without actually thinking about what the model is doing, I don’t know what would be the correct solution to this problem (or if this even is the real issue).

(edit: ok, one final guess without understanding details of the model: it seems that Poisson with parameter -> 0 tends to point mass at 0, so based on this one could skip all likelihood terms where Y[i, j] and dN[i, j] are both zero.)

1 Like
22

I’ve found the paper that this model is most likely based on:

They provide the jags code at the end. It seems like a “litteral translation”

23

@danielw2904 Exactly! I also saw this paper, but did not know they provided the code. Thank you for posting it here! I tried to run this model in jags without any problems. So I was wondering if it is how Stan works that gives this type of error.

24

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;
real<lower=0> dL0[T];
}

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 
    dL0_star[j] <- r * (t[j + 1] - t[j]);
    mu[j] <- dL0_star[j] * c;   
    
    
}

}

model {

for(j in 1:1) {
    for(i in 1:Nsubj) {
        # Likelihood
        if (Y[i,j]==1) { //Exclude Idt[i,j]==0, since it gives no information about the parameter?
        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(’~/learning_R/stan_workshop_2016-master/code/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=T)

print(fit)