CoDatMo Liverpool & UNINOVE models (slow ODE implementation and Trapezoidal Solver)

Modeling
,
45

Hm, okay everyone, forget every other model that has ever been mentioned in this thread. With the typical overreach I present to you the only SEIR model you’ll ever need. (I think) it allows any type of priors, regularization, ansatz or stochasticity for the infection rate, while converging perfectly and fast while probably being much more exact than any of the other models. Tagging the relevant people: @storopoli, @caesoma and @s.maskell and for good measure @bbbales2.

In less than two minutes we get > 4k N_eff for any of 485+ instead of 63 parameters.

So, while I don’t believe its applicable to many problems, here it fits perfectly. It’s essentially what @wds15 said:

The key idea is just that anything downstream of E1 is linear. Instead of parametrizing in terms of some inital condition and betas, we just parametrize in terms of new infections per day. Then, once per leapfrog iteration we can compute a fixed transition matrix using matrix_exp and just iteratively (per day) compute a single matrix vector product.

For sake of simplicity I don’t pool the daily deaths in this model, but this and anything else is perfectly possible and should not impact convergence at all ™.

Enjoy:

data {
  int no_days;
  int population;
  //observed new deaths per day
  int new_deaths[no_days];
  real beta_regularization;
  real likelihood;
}
parameters {
  /*
    The elements of a simplex sum to one.
    The first no_days entries represent daily new infections as a fraction
    of the total population while the last entry represents the proportion that
    remains susceptible at the end.
  */
  simplex[no_days+1] unit_dS;
  real<lower=0> dL;
  real<lower=0> dI;
  real<lower=0> dT;
  real<lower=0, upper=1> omega;
  real<lower=0> reciprocal_phi_deaths;
}
transformed parameters {
  vector[no_days] daily_infections = population * unit_dS[:no_days];
  vector[no_days] daily_deaths;
  vector[no_days] beta;
  if(likelihood){
    vector[7] state = [
        0, 0,
        0, 0,
        0, 0,
        0
    ]';
    matrix[7, 7] transition_matrix = matrix_exp([
    //[E1   ,E2   ,I1   ,I2         ,T1   ,T2   ,D]
      [-2/dL,0    ,0    ,0          ,0    ,0    ,0],//E1
      [+2/dL,-2/dL,0    ,0          ,0    ,0    ,0],//E2
      [0    ,+2/dL,-2/dI,0          ,0    ,0    ,0],//I1
      [0    ,0    ,+2/dI,-2/dI      ,0    ,0    ,0],//I2
      [0    ,0    ,0    ,+2/dI*omega,-2/dT,0    ,0],//T1
      [0    ,0    ,0    ,0          ,+2/dT,-2/dT,0],//T2
      [0    ,0    ,0    ,0          ,0    ,+2/dT,0]//D
    ]);
    real S = population;
    real last_D;
    for(i in 1:no_days){
      last_D = state[7];
      S -= daily_infections[i];
      state[1] += daily_infections[i];
      state = transition_matrix * state;
      daily_deaths[i] = state[7] - last_D;
      beta[i] = daily_infections[i] * population / (S*(state[3]+state[4]));
    }
  }
}
model {
  //One possible regularization
  if(beta_regularization){
    unit_dS[2:no_days] ~ lognormal(log(unit_dS[:no_days-1]), beta_regularization);
  }
  //This imposes a very wide prior on the proportion of still susceptible people!
  unit_dS[no_days+1] ~ uniform(0,1);
  dL ~ normal(4., .2);
  dI ~ normal(3.06, .21);
  dT ~ normal(16, .71);
  omega ~ beta(100, 9803);
  reciprocal_phi_deaths ~ exponential(5);
  if(likelihood){
    new_deaths ~ neg_binomial_2(daily_deaths, 1/reciprocal_phi_deaths);
  }
}

Small caveat: I don’t know how exact matrix_exp is, but probably quite exact.

Edit: also, unit_S is very inappropriately named, I did not know how simplex works.

Edit2: updated the model.

Edit3: fixed an off-by one error in beta[i] = daily_infections[i] * population / (S*(state[3]+state[4])); ← this is correct.

8 Likes