How can we compute the one dimensional integrals of survival function in RSTAN?

Dear all,
I’m new to Stan. The following paragraph is my first stan code for the hierarchical survival model. But I had trouble introducing the attached one-dimensional integral used to calculate the survival function part of likelihood.
Can anyone help me to compute the attached Integral in RSTAN?

Thank you in advance!

My stan code is as below (N=35000, n_grid=350, province=31, num=12).

fit ← stan(“final model.stan”,data = da,iter=iter, warmup=warmup, cores =4,chains = 1, thin = 10, init = “random”,control=list(max_treedepth=15))
Difinition of survival function.pdf (192.9 KB)

“”
data {
int <lower=0> province;
int <lower=0> n_grid;
int <lower=0> N;
int <lower=0> num;
int <lower=0> number;
int <lower =0 , upper = province> county[N];
vector <lower=0> [province] p;
row_vector[n_grid] kernel[province];
matrix [N,number] respon;
}
parameters{
vector [num] beta;
real <lower=0> sigma;
vector [n_grid] x;
}

transformed parameters{
vector [num] time_ratio;
vector [N] lambdaa;
real alpha = 1 / sigma;
vector [n_grid] exp_x;
vector[province] a; // a saved to output because it’s a toplevel transformed parameter
time_ratio = exp(beta);
exp_x = exp(x);
lambdaa = exp (((-1 * respon[,7:17] * beta) / sigma));
// }
{ // z not part of output
vector[province] landa;
for (k in 1 : province) {
landa[k] = kernel[k] * exp_x * p[k];
}
a = landa / sum(landa); // assign a
}
}
model{
vector [N] f;
vector [N] s;
real log_alpha;
log_alpha = log (alpha);
target += normal_lpdf(x|0,5);
target += normal_lpdf(beta | 0, 1000);
target += normal_lpdf(log_alpha | 0, 5);
target += -1 * log_alpha;
for (i in 1:N) {
f[i]= a[county[i]]((lambdaa [i]alpharespon[i,4]^(alpha-1))/((1+(lambdaa [i]respon[i,4]^alpha))^2));
// How can we introduce the integral related to the calculation of the mentioned survival function in the RSTAN?
s[i]= a[county[i]]
(1- ???);
target += respon[i,6] * (log(f[i])-log(s[i])) + log(s[i]);
}
}
generated quantities{
vector [N] log_lik;
vector [N] s;
vector [N] f;
real log_alpha;
log_alpha = log (alpha);
for (i in 1:N){
f[i]= a[county[i]]
((lambdaa [i]alpharespon[i,4]^(alpha-1))/((1+(lambdaa [i]respon[i,4]^alpha))^2));
// How can we introduce the integral related to the calculation of the mentioned survival function in the RSTAN?
s[i]= a[county[i]]
(1- ???);
log_lik [i]= respon[i,6] * (log(f[i])-log(s[i])) + log(s[i]) - log_alpha;
}
}
“”

Hi @Isa_Nazar,
I went ahead and re-formatted your model a bit to make it easier to read. You can put triple back-ticks around code like this to get syntax highlighting and to preserve indentation.

```stan
stan model goes here
.```

As posted, a couple bits of your model got italicized, but I made my best guess. To get to your question, Have you had a look at the description of integrate_1d in the Stan manual? I have just learned how to use it for something like a cumulative hazard function. The syntax looked a little intimidating, but the example of how to call the integrator is helpful.

1 Like

Thank you so much for the re-formatting of my code as well as your quick reply. Actually, I had looked at the Stan manual but I couldn’t use it for my survival function. Can you share with me your codes for cumulative hazard function or any other similar codes, please?

Thank you!

I’m happy to share the code I have so far. I’ve learned that when developing a model in stan, it’s important to start simple and build up to the more complex final model.

The idea here is to model a vector of data (rt) as having a distribution proportional to a normal cumulative density between two limits. To do that, the likelihood needs to be normalized by the definite integral of the cdf between the limits (so the likelihood integrates to one).

functions {
    real integrand(real x, real xc, real[] theta, real[] x_r, int[] x_i) {
        return normal_cdf(x, theta[1], theta[2]);
    }
}
data {
    int n_trial;
    real lb;
    real<lower=lb> ub;
    vector<lower=lb, upper=ub>[n_trial] rt; //relative to trial start
}
parameters {
    // parameters of the hazard function
    real mu_antic;
    real<lower=0> sig_antic;
}
transformed parameters {
    real log_integ = log(integrate_1d(integrand, lb, ub, {mu_antic, sig_antic}, {0.0}, {0}));
}
model {
    mu_antic ~ normal(1.35, 0.4);
    sig_antic ~ normal(0, 1);
    for (i in 1:n_trial) {
        target += normal_lcdf(rt[i] | mu_antic, sig_antic) - log_integ; // the issue is that this normal_cdf doesn't integrate to 1. We need a normalizing value.
    }
}
1 Like

Thanks for your help, I greatly appreciate it. I try to write my codes for computing one-dimensional integrals like your codes.

Kind regards,
Eisa

1 Like