Model running extremely slow

Modeling
1

Hi all,

The following model takes 2.625 seconds, and it is taking really long to run.
I do not know if the problem belongs to the data.? (I have around 700 ids and 2085 obs), or the problem is the prior distributions? But I also do not think that the priors are very informative.

functions {
  real[] sldecay(real t, real[] y, real[] theta,  real[] x_r, int[] x_i) {
    real dydt[3];
    dydt[1] = -theta[1]*y[1];
    dydt[2] = -theta[2]*y[2];
    dydt[3] = theta[1]*y[1] + theta[2]*y[2] - theta[3]*y[3];
    
    return dydt;
  }
}
data {
  int<lower=1> N;
  int<lower=1> nsubjects;
  int<lower=1> subj_ids[N];
  real <lower=0> y[N];
  real ts[N];
  real t0;
}
transformed data {
    real x_r[0];
    int x_i[0]; 
}
parameters {
  real logps;
  real logpl;
  real logmuab;
  
  real logAB0;
  real logps0;
  real logpl0;
  real <lower= 0> sigma;
}
model {
  real new_mu[N,3];
  real mu[1,3];
  real eval_time[1];
  real theta[3];
  
  real inits[3];
  
  real ps = exp(logps);
  real pl = exp(logpl);
  real muab = exp(logmuab);
  real AB0 = exp(logAB0);
  real ps0 = exp(logps0);
  real pl0 = exp(logpl0);
  
  //prior distributions
  logps ~ normal(0, 1);
  logpl ~ normal(0, 1);
  logmuab ~ normal(0, 10);
  logAB0 ~ normal(5, 10);
  logps0 ~ normal(0, 1);
  logpl0 ~ normal(0, 1);
  sigma ~ normal(1, 10);
  
  for (obs_id in 1:N){
    inits[1] = AB0;
    inits[2] = inits[1]/ps;
    inits[3] = inits[1]/pl;
    
    theta[1] = ps;
    theta[2] = pl;
    theta[3] = muab;
     
    eval_time[1] = ts[obs_id];
    if (eval_time[1] == 0){
      new_mu[obs_id,1] = log(inits[1]) - pow(sigma, 2.0)/2;
      new_mu[obs_id,2] = inits[2];
      new_mu[obs_id,3] = inits[3];
      }
    else{
      mu = integrate_ode_rk45(sldecay, inits, t0, eval_time, theta, x_r, x_i);
      new_mu[obs_id,1] = log(mu[1,1]) - pow(sigma, 2.0)/2;
      new_mu[obs_id,2] = mu[1,2];
      new_mu[obs_id,3] = mu[1,3];
    }  
    //likelihood function for AB levels
    y[obs_id] ~ lognormal(new_mu[obs_id,1], sigma); 
  }
}
generated quantities {
  real z_pred[N];
  real sigma2 = sigma;
  for (t in 1:N){
   z_pred[t] = lognormal_rng(log(y[t] + 1e-5), sigma2); 
  }
}
2

you should integrate the entire time-series at once with the integrate call… and you should evaluate the likelihood in a vectorized manner rather than one-by-one. I hope this points you in the right direction.

3

Thanks! I do not really understand the call to integrate, cause that is already on the program in the function calling mu?

4

“integrate the entire time series at once” mean something like

model {
  real new_mu[N,3];
  real mu[N,3];
  real eval_time[1];
  real theta[3];
  
  real inits[3];
  
  real ps = exp(logps);
  real pl = exp(logpl);
  real muab = exp(logmuab);
  real AB0 = exp(logAB0);
  real ps0 = exp(logps0);
  real pl0 = exp(logpl0);
  
  //prior distributions
  logps ~ normal(0, 1);
  logpl ~ normal(0, 1);
  logmuab ~ normal(0, 10);
  logAB0 ~ normal(5, 10);
  logps0 ~ normal(0, 1);
  logpl0 ~ normal(0, 1);
  sigma ~ normal(1, 10);

  theta[1] = ps;
  theta[2] = pl;
  theta[3] = muab;
  inits[1] = AB0;
  inits[2] = inits[1]/ps;
  inits[3] = inits[1]/pl;

  mu[1] = inits;
  mu[2:] = integrate_ode_rk45(sldecay, inits, t0, ts[2:], theta, x_r, x_i);
  
  for (obs_id in 1:N){
    new_mu[obs_id,1] = log(mu[obs_id,1]) - pow(sigma, 2.0)/2;
    new_mu[obs_id,2] = mu[obs_id,2];
    new_mu[obs_id,3] = mu[obs_id,3];
    //likelihood function for AB levels
    y[obs_id] ~ lognormal(new_mu[obs_id,1], sigma); 
  }
}

You may have to sort the data ahead of time for that to work.

The model doesn’t use nsubjects or subj_ids. Are there parts missing?

5

Hi,

Thanks for the message. The subj_ids or subjects were used in the loop before, and will be used when I add random effects to the model for each individual.

6

Hi,

My data is already sorted (time is from 0 to 96), but when I run the model I get an initialization chain failed. What exactly should the data look like?

7

what does the error message say?

8

Initialization failed chain (-2, 2) error in ts 0 expecting to be higher than 2.

9

Not sure what that means. Maybe the first ts is 0 and the initial time t0 is 2?