How to include a changing coefficient in ODEs in rstan

Modeling
,
1

Hi all, I’m new to stan. I now have trouble in including a temperature-denpenent coefficient in stan. Below is my code. In this ODE, the coefficent “b” can be calculated based on a given temperature. However, I only know how to build a stan model by using a fixed value, but I have no idea on how to include a varying coefficent into stan.

functions {
  real[] seir(real t, real[] y, real[] theta,
             real[] x_r, int[] x_i) {

      real N = x_i[1];
      
      real beta = theta[1];
      real gamma = theta[2];
      real a = theta[3];
      real i0 = theta[4];
      real e0 = theta[5];
      real b;
      
      real init[4] = {N - i0 - e0, e0, i0, 0}; 
      real S = y[1] + init[1];
      real E = y[2] + init[2];
      real I = y[3] + init[3];
      real R = y[4] + init[4];
      
      real dS_dt = -b * beta * I * S / N;
      real dE_dt =  b * beta * I * S / N - a * E;
      real dI_dt = a * E - gamma * I;
      real dR_dt =  gamma * I;
      
      return {dS_dt, dE_dt, dI_dt, dR_dt};
  }
}
data {
  int<lower=1> n_days;
  real t0;
  real y0[4];
  real ts[n_days];
  int N;
  int cases[n_days];
  real b[n_days];
}
transformed data {
  real x_r[0];
  int x_i[1] = { N };
}
parameters {
  real<lower=0> gamma;
  real<lower=0> beta;
  real<lower=0> a;
  real<lower=0> phi_inv;
  real<lower=0> i0; 
  real<lower=0> e0;
}
transformed parameters{
  real y[n_days, 4];
  real phi = 1. / phi_inv;
  real theta[5] = {beta, gamma, a, i0, e0};
  y = integrate_ode_rk45(seir, rep_array(0.0, 4), t0, ts, theta, x_r, x_i);
}
model {
  beta ~ normal(0,1);
  gamma ~ normal(0,1);  
  a ~ normal(0,1);
  phi_inv ~ exponential(5);
  i0 ~ normal(0, 10);
  e0 ~ normal(0, 10);
  
  cases[1:(n_days)] ~ neg_binomial_2(y[,3], phi);
}
generated quantities {
  real M = beta / gamma;
}

My data looks like:
cases<-c(2,5,3,8,12,10,14, 20, 23…) (n=104)
b<-c(0.798, 0.752, 0.813, 0.806,0.774, 0.785,…) (n=104)

If I run the above code, I’ll get the error: “Exception: neg_binomial_2_lpmf: Location parameter[1] is nan, but must be > 0!” I know it must be something wrong with “b”, but I dont know how to correct it. Thank you so much if anyone can help me.

2

Hi Ruby,

Since b is not a parameter, you could pass b in x_r in the transformed data block. Then, inside your ode function, have a for statement which loops through time, combined with an if statement to assign to b the respective entry from x_r if the solver is currently at time point “t”.

Hope I got your question right.