Setting ODE integrator time step to be different than t

Modeling
1

Hi everyone,

I’m quite new to modeling and Stan, so this might be a stupid question, but I’m still on a learning curve!

I’m trying to model transmission of an infectious disease with yearly rates defining the flow between states (compartments). However, I’d like the ODE integrator (rk45) to solve the equations at a monthly time step to have smoother curves.

A particularity of my model is also that I’m fitting to 1 prevalence sample (1 data point).

I can’t find any examples of this, and when I try doing it on my own, I run into quite a few errors.

Here are some relevant parts of my Stan code (i’ve removed the paramters and equations etc to make it easier to read), with the integrator solving yearly for 40 yrs if I’m understanding and coding correctly:

functions {
  real[] dtb1(real t, real[] y, real[] theta, 
              real[] x_r, int[] x_i) {
              
  //define interger data//
    int n_cases = x_i[1];
    int n_samples = x_i[2];
    int N = x_i[3];
    
  //define real data//
    
  //define paramaters needed to solve model//
    
  //define compartments & init conditions //
    
  //define states & ODE//
    
    return {dS_dt, dE_dt, dL_dt, dI_dt, dR_dt};
  }
}

data {
  int<lower = 1> n_ts;        // number of times to solve for
  int<lower = 1> N;          // population size
  int<lower = 1> n_cases;     // observed cases
  int<lower = 1> n_samples;   // number of data points to fit to 
  int<lower = 1> n_comp;     // number of compartments

  real t0;                    //initial time point 0
  real ts[n_ts]; 

  int y[n_samples, 1];       // data, reported prev in ea sample


}

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

}

 parameters {
}

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

  real init[n_comp] = {1231, 524, 1057, 150, 2036};

  simplex [n_comp] 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 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 {
  // sampling distribution
for(i in 1:n_ts)  n_cases ~ binomial(N, output[i,4]);
}

generated quantities {
array[n_ts] real pred_cases;

for (i in 1:n_ts)
pred_cases[i] = binomial_rng(N, output[i,4]);

}


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

#data
n_cases<-c(47)
N<-1080
t<-1
ts<-seq(1, 40, by=1)
n_ts<-40
y <- array(c(I = n_cases), dim = c(1, 1))

stan_d = list(n_samples =1,
              N=N,
              n_comp = 5,
              n_ts= n_ts,
              n_cases=n_cases,
              ts=ts,
              y=y,
              f=f,
              r=r,
              t0 = 0)


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

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

Thank you in advance for any help!
Morgana

2

[edit: initially misunderstood question]

Just replace ts<-seq(1, 40, by=1) with ts<-seq(1, 40, by=1.0/12)

You’ll also need to code the sizes to be 480 rather than 40 (your n_ts I think).

3

yep, works perfectly.
thank you very much!

4

Just to be clear, your data is fit to yearly prevalence data, but you want to plot the model curves at smaller intervals, is that correct?

If that is the case, you need to use a time range compatible with the number of observations you have (n_ts as it seems in your code), but you will need another variable (say n_tp such that real tp[n_tp]) to integrate the ODE model for each accepted sample at a smaller interval.

You could just use single time range with a smaller interval and pick the time points where you have observations, but that would probably be more cumbersome.

5

Thank you for your response @caesoma,
This was something I’d thought of and initally tried, but I’m unsure about the implementation of what you suggest.
You would use your monthly time variable in your ode solver, and then model using a yearly time variable? or do everything post solving the ode with the smaller time step?

Thanks again.

6

I think the clearer way to write is to do it separately, in two separate blocks:

  1. in the model block you would solve the system for the current sample/proposal for the time range ts of type/size real ts[n_ts] (e.g. for n_ts = 40 years) with your y_hat = integrate_ode_rk45(dtb1, init, t0, ts, theta, x_r, x_i); (you could keep it in transformed parameters if you are want to save it for any post hoc additional calculation of the posterior/likelihood like WAIC, for instance).

  2. in the generated quantities block you do the same computation of the solution but with tp of type/size real tp[n_tp] where n_tp = 12*n_ts, for plotting purposes only.

The additional cost of a separate ODE solution is probably negligible, considering the number of times HMC needs to compute it and the gradient for each leapfrog step.

7

@caesoma Thank you, I was able to implement it and I think its doing exactly what I was looking for it to do!! the curves appear much smoother and the run time is pretty quick :)

1 Like