Dear all,
I am working on estimating the probability of a binary endpoint that we measure at different timepoints. Once it is observed at one timepoint, it is not longer followed and measured.
I build a stan that I fit with rstan, but it is very very long to run. I tried to optimize some parts, but I need to do better because I need to call this estimation a lot.
I would like to have your opinion and your help on how I could gain time when fitting it ?
Thank you in advance
data {
int nb_seq;
int nb_adm;
int n_pat;
int nb_adm_pat[n_pat];
matrix[nb_seq,nb_adm] sequences;
matrix[n_pat, nb_adm] seq_donnees;
int y[n_pat,nb_adm];
real doseref;
vector[9] prior;
}
transformed data {
matrix[nb_seq,nb_adm] norm_cum_sequences;
matrix[nb_seq,nb_adm] norm_sequences;
for(k in 1:nb_seq){
norm_cum_sequences[k,] = cumulative_sum(sequences[k,]/doseref);
norm_sequences[k,] = sequences[k,]/doseref;
}
}
parameters {
vector[3] param;
}
transformed parameters {
vector[3] exp_param;
vector[nb_seq] ptox_seq;
exp_param = exp(param);
for(k in 1:nb_seq){
vector[nb_adm] ptox_j;
vector[nb_adm] ptox_cum_j;
vector[nb_adm] temp_cumprod;
real exp_calcul;
exp_calcul = exp(param[1] - exp_param[2] + exp_param[3]*log(norm_sequences[k,1]));
ptox_j[1] = exp_calcul / (1 + exp_calcul);
ptox_cum_j[1] = ptox_j[1];
temp_cumprod[1] = 1-ptox_j[1];
for(j in 2:nb_adm){
exp_calcul = exp(param[1] - exp_param[2]*log(norm_cum_sequences[k,j-1]) + exp_param[3]*log(norm_sequences[k,j]));
ptox_j[j] = exp_calcul / (1 + exp_calcul);
temp_cumprod[j] = temp_cumprod[j-1] * (1-ptox_j[j]);
ptox_cum_j[j] = ptox_j[1] + sum(ptox_j[2:j] .* temp_cumprod[1:(j-1)]);
}
ptox_seq[k] = ptox_cum_j[nb_adm];
}
}
model {
vector[3] mu;
matrix[3,3] Sigma;
mu[1] = prior[1];
mu[2] = prior[2];
mu[3] = prior[3];
Sigma[1,1] = prior[4]*prior[4];
Sigma[2,2] = prior[5]*prior[5];
Sigma[3,3] = prior[6]*prior[6];
Sigma[1,2] = prior[7]*prior[4]*prior[5];
Sigma[2,1] = Sigma[1,2];
Sigma[1,3] = prior[8]*prior[4]*prior[6];
Sigma[3,1] = Sigma[1,3];
Sigma[2,3] = prior[9]*prior[5]*prior[6];
Sigma[3,2] = Sigma[2,3];
param ~ multi_normal(mu, Sigma);
for(i in 1:n_pat){
for(j in 1:nb_adm_pat[i]){
real ptox;
real exp_calcul;
if(j==1){
exp_calcul = exp(param[1] - exp_param[2] + exp_param[3]*log(seq_donnees[i,1]/doseref));
}
else if(j>1){
exp_calcul = exp(param[1] - exp_param[2]*log(sum(seq_donnees[i,1:(j-1)])/doseref) + exp_param[3]*log(seq_donnees[i,j]/doseref));
}
ptox = exp_calcul / (1 + exp_calcul);
y[i,j] ~ bernoulli(ptox);
}
}
}