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).
“”
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.
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?
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.
}
}