Error "Stan model 'anon_model' does not contain samples"

Modeling
, ,
1

Chains do not contain any samples

Hi everyone, I am new to mathematical modelling as well as Rstan and I’m trying to fit my first, pretty simplistic, compartmental model (for disease transmission using ODEs)
the model seems to compile just fine but when I try to fit the model to the data, I get an error message saying : “Stan model ‘anon_model’ does not contain samples.”

I’ve made a few adjustments to ts and the initial values, which seemed to help, but I continue to get the same error
when I look at the chains specifically, it says :
"Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: Exception: Exception: array[uni, …] index: accessing element out of range. index 2 out of range; expecting index to be between 1 and 1 (in ‘string’, line 7, column 4 to column 26) (in ‘string’, line 84, column 4 to column 68)
[1] “Error : Exception: Exception: array[uni, …] index: accessing element out of range. index 2 out of range; expecting index to be between 1 and 1 (in ‘string’, line 7, column 4 to column 26) (in ‘string’, line 84, column 4 to column 68)”

I’ve put my entire code below, if someone has any tips or notices where I might have gone wrong, I’d really appreciate the help!

#### MODEL FITTING WITH r STAN ####
cat(
  '
functions {
  real[] dtb1(real t, real[] y, real[] theta, 
              real[] x_r, int[] x_i) {
              
  //define interger data//
    int N = x_i[1];
    int n_latent = x_i[2];
    int n_infectious = x_i[3];
    int n_recov = x_i[4];
    int n_samples = x_i[5];
    int time_seed = x_i[6];
    
  //define real data//
    real f = x_r[1];         //proportion of fast progressing cases
    real r = x_r[2];         //proportion of successfully completed treatment 
    
  //define paramaters needed to solve model//
    real sigma = theta[1];     //transmission rate
    real mu = theta[2];       //all-cause mortality rate, with equal births
    real gamma = theta[3];    //disease progression rate latent fast
    real epsilon= theta[4];   //disease progression rate latent slow
    real delta = theta[5];    //case detection & treatment initiation rate
    real omega = theta[6];    //recovery rate (1/mean treatment duration time 2 yrs) 
    real chi = theta[7];      //TB specific mortality
    real pi = theta[8];       //spontaneous cure rate
    real rho = theta[9];      //relapse rate
    
  //define compartments & init conditions //
    real S = y[1];            //uninfected
    real E = y[2];           //latent fast progressing TB
    real L = y[3];           //latent slow progressing TB
    real I = y[4];            //active & infectious TB
    real R = y[5];            //recovered

  //define states & ODE//
    real dS_dt = -sigma*I*S/N - mu*S + mu*5000 + chi*I;
    real dE_dt = f*sigma*I*S/N - gamma*E - mu*E;
    real dL_dt = (1-f)*sigma*I*S/N - epsilon*L - mu*L;
    real dI_dt = gamma*E + epsilon*L - (mu+chi)*I - pi*I - r*delta*omega*I + rho*R;
    real dR_dt = r*delta*omega*I - rho*R  - mu*R;
    
    return {dS_dt, dE_dt, dL_dt, dI_dt, dR_dt};
  }
}

data {
  int<lower = 1> n_ts;        // number of days observed
  int<lower = 1> n_recov;     // recovered
  int<lower = 1> n_infectious;
  int<lower = 1> n_lf; 
  int<lower = 1> n_ls; 
  int<lower = 1> N;           // Total population size 
  int<lower = 1> n_samples;   // number of data points to fit to 
  int<lower = 1> n_comp;   // number of compartments
  int <lower = 0> time_seed;  // index when to fit the model to data
  
  real<lower=0> f;	      // prop fast progressing cases
  real<lower=0> r;	      // prop success complete Rx 
  real t0;                    //initial time point 0
  real ts[n_ts]; 
  int y[n_samples, n_comp];  // data, reported prev in ea sample


}

transformed data {
  real x_r[2] = {f, r};
  int x_i[1] = {N};

}

parameters {
  real<lower=0> sigma;
  real<lower=0> mu;
  real<lower=0> gamma;
  real<lower=0> epsilon;
  real<lower=0> delta;
  real<lower=0> omega;
  real<lower=0> chi;
  real<lower=0> pi;
  real<lower=0> rho;  
}

transformed parameters {
  real y_hat[n_ts, n_comp];         // Matrix to store the integrated ODE solutions
 
  real theta[9] = {sigma, mu, gamma, epsilon, delta, omega, chi, pi, rho};

  real init[n_comp] = {N - n_lf - n_ls - n_infectious - n_recov, n_lf, n_ls, n_infectious, n_recov};
  
  simplex [5] output[n_ts];


// Solve the ODE using integrate_ode_rk45
    y_hat = integrate_ode_rk45(dtb1, init, t0, ts, theta, x_r, x_i);
  

for(i in time_seed+1:n_ts) output[i] = ( to_vector(y_hat[i,1:5]) ) ./ sum( to_vector(y_hat[i,1:5])); // get the proportion in each category  
}

model {
  // Prior distribution
  sigma ~ normal(2, 1);
  mu ~ normal(0.0013, 0.05);
  gamma ~ normal(0.8, 0.1);
  epsilon ~ normal(0.05, 0.1);
  delta ~ normal(0.6, 0.1);
  chi ~ normal(0.16667, 0.05);
  rho ~ normal(0.05, 0.02);
  pi ~ normal(0.16667, 0.1);
  omega ~ normal(0.5, 0.2);
  
  
  // Likelihood
  for(i in time_seed+1:n_ts) y[i] ~ multinomial(output[i] );
}


', 
file = "dtb1_multinom.stan", sep="", fill=T)

#data
N <- 5000
n_lf<-round(5000*0.2606)
n_ls<-round(5000*0.0652)
n_infectious<-round(5000*0.0443)
n_recov<-round(5000*0.1337)
time_seed<-1000
n_ts<-2000
f=0.8
r=0.76   ##(6/41 discontinued rx=35/41*0.9Rx success=0.76)
y <- array(c(S = N - n_lf - n_ls - n_infectious - n_recov,
             E = n_lf,
             L = n_ls,
             I = n_infectious,
             R = n_recov), dim = c(1, 5))
ts <- ts$ts

stan_d = list(n_samples = 1,
              n_comp = 5,
              n_lf = n_lf,
              n_ls = n_ls,
              n_infectious = n_infectious,
              n_recov = n_recov,
              n_ts= n_ts,
              ts=ts,
              time_seed = time_seed,
              y=y,
              f=f,
              r=r,
              t0 = 0,
              N=N)


#compile stan model
model <- stan_model("dtb1_multinom.stan")

fit <- sampling(model,
                data = stan_d,
                iter = 2000,
                chains = 4, 
                seed = 0)

Thanks again in advance!

2

Sorry, @morg_ana — these involved modeling questions tend to scare people off.

The error messages are where to look—thanks for including. First, I would recommend putting your Stan code in its own file so that you can see the line numbers. It will also let you use strings (") and transposition (') in the same program.

Here’s the error you report:

The relevant code line is here:

    int n_latent = x_i[2];

What this is saying is that index 2 is out of range from x_i and saying the index must be between 1 and 1 which indicates x_i only has one element. This system function gets called with the integer data from the ODE solver:

y_hat = integrate_ode_rk45(dtb1, init, t0, ts, theta, x_r, x_i);

That’s the bug. You define x_i to only have one element here:

transformed data {
  real x_r[2] = {f, r};
  int x_i[1] = {N};
}
3

Thank you so much, Bob!
indeed that was where the issue was

and thank you for the suggestion of putting the stan code in a separate file, it really does help with deciphering error messages… might seem obvious but it’s extremely helpful to a beginner such as myself