Fitting SEIR using a binomial distribution on the incidence data

Modeling
1

Hi,

This is my first try with Stan and R. I’m trying to calibrate the contact rate in an SEIR model from incidence data. I assume the incidence data is binomially distributed, where p indicates the probability of a case being reported.

My problem is that the incidence variable is real vector, and the binomial function requires an integer vector.

Is there a way to make this work or should I change my approach?

functions {
  real[] SEIR(real t,
            real[] y, // stocks
            real[] params,
            real[] x_r,
            int[] x_i) {
              
      real dydt[4];
      real forceOfInfection;
      real population;
      real prevaleceOfInfection;
      real newInfections;
      real completingLatentPeriod;
      real recovery;
      
      population = y[1] + y[2] + y[3] + y[4];
      prevaleceOfInfection = y[3] / population;
      forceOfInfection = params[1] * prevaleceOfInfection * 0.005; //beta
      newInfections = y[1] * forceOfInfection;
      completingLatentPeriod = y[2] / 0.42857; //Tao
      recovery = y[3] / 1; //Mu
      
      dydt[1] = - newInfections;
      dydt[2] = newInfections - completingLatentPeriod;
      dydt[3] = completingLatentPeriod - recovery;
      dydt[4] = recovery;
      
      return dydt;
    }
}

data {
  int<lower = 1> n_obs; // Number of days sampled
  int<lower = 1> n_params; // Number of model parameters
  int<lower = 1> n_difeq; // Number of differential equations in the system
  int y[n_obs]; // The binomially distributed data
  real t0; // Initial time point (zero)
  real ts[n_obs + 1]; // Time points that were sampled
}

transformed data {
  real x_r[0];
  int x_i[0];
}

parameters {
  real<lower = 0> params[n_params]; // Model parameters
  real<lower = 0, upper = 1> p; 
}

transformed parameters{
  real y_hat[n_obs + 1, n_difeq]; // Output from the ODE solver
  real y0[n_difeq]; // Initial conditions
  real incidence[n_obs];
    

  y0[1] = 999;
  y0[2] = 1;
  y0[3] = 0;
  y0[4] = 0;
  
  y_hat = integrate_ode_rk45(SEIR, y0, t0, ts, params, x_r, x_i);
  
  for (i in 0:n_obs) {
    incidence[i] = y_hat[i + 1, 3] - y_hat[i, 4] + y_hat[i + 1, 4]; 
  }
  
}

model {
  p ~ uniform(0, 1);
  params[1] ~ uniform(0, 500); 
  y ~ binomial(incidence, p); 
}

R code:

library(rstan)

stan_d = list(n_obs = 32,
              n_params = 1,
              n_difeq = 4,
              y = c(1,1,2,2,4,7,7,12,22,28,41,49,59,57,60,51,45,27,19,19,8,5,6,2,3,1,1,0,0,0,0,0),
              t0 = 0,
              ts = 0:31)


# Test / debug the model:
test = stan("SEIR.stan",
            data = stan_d,
            chains = 1, iter = 10)
1 Like
2

What did you intend when using binomial PMF where the size is a function of a solution to ODEs, rather than a fixed integer (array)? My guess is that you would have to use the normal approximation to a binomial.

3

I’m trying to reproduce an example in chapter 5 of this book. There, the authors assume the incidence is a binomial counting. They use the MCMCpack library. My goal is to get the same result using Hamiltonian MC.

There is no way to round up the incidence and convert it to an integer vector?

4

Can you write down the model so we can try and help you? There might be a few tricks we can pull, but I really need to see the maths.

5

That would lead to discontinuities in the posterior density and HMC would not work well.

1 Like
6

Well, it doesn’t exactly make sense to model a discrete parameter–the number of people infected–using differential equations.
If we assume that the exact (integer) number of infections is Poisson distributed with the average equal to the number predicted by the ODE then the problem simplifies a lot. The binomial becomes a Poisson distribution and you have

y ~ poisson(incidence*p);

This is also an appropriate approximation to the binomial if incidence is large and p is small.

1 Like
7

I agree with @nhuurre. What can be done is, if you measure n_t individuals at time t, then you can model how many would be infected using a \text{binomial}(n_t, I(t)), where I(t) is the scaled solution of a SIR (resp. SEIR) model. Something along the lines of this: https://jrmihalj.github.io/estimating-transmission-by-fitting-mechanistic-models-in-Stan/

2 Likes
8

Thanks for all the comments. Based on them, I’ll rethink my approach. I’m leaning towards considering the incidence data as following the Poisson distribution. I’ll give a try, and in case of doubt, will ask.

Thanks again.

1 Like
9

Hi, following your advice, I could fit an SEIR model. Based on that, I tried to fit a disaggregated SEIR. It has 16 parameters. I ran 3 chains with 1000 iterations, 500 warm-up.

Could you help to interpret this output?

Please let me know if it’s necessary to show the code. My question is aimed toward getting some intuition when I see that kind of results. Thanks.

Rplot
Rplot1024×686 97.4 KB

10

The intuition is that the geometry of your model is difficult to explore and this is most likely cause by an unidentifiable or nearly unidentifiable model. Can you share the code and same sample data?

11

This is the stan code:

    functions {
       real[] SEIR_extended(real t,
       real[] y,
       real[] params,
       real[] x_r,
       int[] x_i) {
          real dydt[16];
          real IR1;
          real InR1;
          real RR1;
          real IR2;
          real InR2;
          real RR2;
          real IR3;
          real InR3;
          real RR3;
          real IR4;
          real InR4;
          real RR4;
          IR1 = y[1]*(params[2]*y[3]+params[1]*y[7]+params[3]*y[11]+params[4]*y[15]);
          InR1 = y[2]/0.2142857142;
          RR1 = y[3]/0.2142857142;
          IR2 = y[5]*(params[5]*y[3]+params[6]*y[7]+params[7]*y[11]+params[8]*y[15]);
          InR2 = y[6]/0.2142857142;
          RR2 = y[7]/0.2142857142;
          IR3 = y[9]*(params[9]*y[3]+params[10]*y[7]+params[11]*y[11]+params[12]*y[15]);
          InR3 = y[10]/0.2142857142;
          RR3 = y[11]/0.2142857142;
          IR4 = y[13]*(params[13]*y[3]+params[14]*y[7]+params[15]*y[11]+params[16]*y[15]);
          InR4 = y[14]/0.2142857142;
          RR4 = y[15]/0.2142857142;
          dydt[1] = -IR1;
          dydt[2] = IR1-InR1;
          dydt[3] = InR1-RR1;
          dydt[4] = RR1;
          dydt[5] = -IR2;
          dydt[6] = IR2-InR2;
          dydt[7] = InR2-RR2;
          dydt[8] = RR2;
          dydt[9] = -IR3;
          dydt[10] = IR3-InR3;
          dydt[11] = InR3-RR3;
          dydt[12] = RR3;
          dydt[13] = -IR4;
          dydt[14] = IR4-InR4;
          dydt[15] = InR4-RR4;
          dydt[16] = RR4;
          return dydt;
      }
    }
    data {
        int<lower = 1> n_obs; // Number of weeks sampled
        int<lower = 1> n_params; // Number of model parameters
        int<lower = 1> n_difeq; // Number of differential equations in the system
        int y1[n_obs];
        int y2[n_obs];
        int y3[n_obs];
        int y4[n_obs];
        real t0; // Initial time point (zero)
        real ts[n_obs]; // Time points that were sampled
    }
    transformed data {
       real x_r[0];
       int x_i[0];
    }
    parameters {
        real<lower = 0> params[n_params]; // Model parameters
    }

    transformed parameters {
      real y_hat[n_obs, n_difeq]; // Output from the ODE solver
      real y0[n_difeq]; // Initial conditions
      real incidence1[n_obs];
      real incidence2[n_obs];
      real incidence3[n_obs];
      real incidence4[n_obs];
      y0[1] = 949;
      y0[2] = 0;
      y0[3] = 0;
      y0[4] = 0;
      y0[5] = 1690;
      y0[6] = 0;
      y0[7] = 0;
      y0[8] = 0;
      y0[9] = 3466;
      y0[10] = 0;
      y0[11] = 1;
      y0[12] = 0;
      y0[13] = 1893;
      y0[14] = 0;
      y0[15] = 1;
      y0[16] = 0;
      y_hat = integrate_ode_rk45(SEIR_extended, y0, t0, ts, params, x_r, x_i);
      incidence1[1] =  y_hat[1, 3] + y_hat[1, 4] - y0[3] - y0[4];
      incidence2[1] =  y_hat[1, 7] + y_hat[1, 8] - y0[7] - y0[8];
      incidence3[1] =  y_hat[1, 11] + y_hat[1, 12] - y0[11] - y0[12];
      incidence4[1] =  y_hat[1, 15] + y_hat[1, 16] - y0[15] - y0[16];
      for (i in 1:n_obs-1) {
        incidence1[i + 1] = y_hat[i + 1, 3] + y_hat[i + 1, 4] - y_hat[i, 3] - y_hat[i, 4] + 0.000001;
        incidence2[i + 1] = y_hat[i + 1, 7] + y_hat[i + 1, 8] - y_hat[i, 7] - y_hat[i, 8] + 0.000001;
        incidence3[i + 1] = y_hat[i + 1, 11] + y_hat[i + 1, 12] - y_hat[i, 11] - y_hat[i, 12] + 0.000001;
        incidence4[i + 1] = y_hat[i + 1, 15] + y_hat[i + 1, 16] - y_hat[i, 15] - y_hat[i, 16] + 0.000001;
      }
    }
    model {
       params ~ normal(0, 2);
       y1     ~ poisson(incidence1);
       y2     ~ poisson(incidence2);
       y3     ~ poisson(incidence3);
       y4     ~ poisson(incidence4);
    }

incidence
incidence475×536 11.9 KB

12

OK, there’s a lot to unpack here, so let’s take baby steps. It might very well be that a model with this many parameters is not identifiable. What is the smallest model you can fit reliably? Another question is whether the Poisson approximation works well for all of the data trajectories. As discussed by @nhuurre above, the Poisson trick is an approximation that deteriorates fast when data departs from the large N small p regime.

1 Like
13

I could fit this SIR model:

    functions {
      real[] SIR(real t,
                real[] y, // stocks
                real[] params,
                real[] x_r,
                int[] x_i) {
                  
          real dydt[3];
          real population;
          real probability;
          real contactsPerInfected;
          real IR; // Infection rate
          real RR; // Recovery rate
          
          population          = y[1] + y[2] + y[3];
          probability         = y[1] / population;
          contactsPerInfected = y[2] * params[1]; // effective contacts
          
          IR                  = contactsPerInfected * probability;
          RR                  = y[2] * params[2]; // recovery proportion
          
          dydt[1] = -IR;
          dydt[2] = IR - RR;
          dydt[3] = RR;
          
          return dydt;
        }
    }

    data {
      int<lower = 1> n_obs; // Number of weeks sampled
      int<lower = 1> n_params; // Number of model parameters
      int<lower = 1> n_difeq; // Number of differential equations in the system
      int y[n_obs];
      int y2[n_obs];
      real t0; // Initial time point (zero)
      real ts[n_obs]; // Time points that were sampled
    }
    transformed data {
      real x_r[0];
      int x_i[0];
    }
    parameters {
      real<lower = 0> params[n_params]; // Model parameters
      real<lower = 0> S0;
    }
    transformed parameters {
      real y_hat[n_obs, n_difeq]; // Output from the ODE solver
      real y0[n_difeq]; // Initial conditions
      real incidence[n_obs];
      real recovery[n_obs];
      y0[1] = S0;
      y0[2] = 10;
      y0[3] = 0;
      y_hat = integrate_ode_rk45(SIR, y0, t0, ts, params, x_r, x_i);
      incidence[1] =  y_hat[1, 2] + y_hat[1, 3] - y0[2] - y0[3];
      recovery[1] = y_hat[1, 3] - y0[3];
      for (i in 1:n_obs-1) {
        incidence[i + 1] = y_hat[i + 1, 2] + y_hat[i + 1, 3] - y_hat[i, 2] - y_hat[i, 3] + 0.000001; 
        recovery[i + 1] = y_hat[i + 1, 3] - y_hat[i, 3] + 0.000001;
      }
    }
    model {
      params[1] ~ uniform(0, 500);
      params[2] ~ uniform(0, 20);
      S0 ~ normal(791, 250);;
      y ~ poisson(incidence);
      y2 ~ poisson(recovery);
    }
14

Is it possible for you to show the resulting Stan run? Also, how do the posterior predictive trajectories look like?