Initialization Error in Stochastic Epidemic SEIR Model

Hello,
My group members and I are trying to replicate a paper which uses MCMC for an assignment. The paper uses stochastic discrete‐time SEIR model to calculate the movement of population between compartment at each time point. Since this is our first time dealing with Stan or Bayesian modelling, we are having some difficulties in understanding where we are going wrong. The paper that we are trying to replicate for our assignment is “https://onlinelibrary.wiley.com/doi/full/10.1111/j.1541-0420.2006.00609.x”.
Our Code:

library(rstan)
data1 <- ebola_kikwit_1995
Cases_Sym_onset <- data1$onset
Cases_Removed <- data1$death
n_days <- length(Cases_Sym_onset)
ncases <- sum(Cases_Sym_onset)
tcontrol=130
N = 5364500
ST= N-ncases
m= ST-5364500
m #  final size of the epidemic

S <- vector(mode="numeric", length = 192)
S[1]=5364500 #  initial conditions for S
E <- vector(mode="numeric", length = 192)
E[1]= 1  #  initial conditions for E
I <- vector(mode="numeric", length = 192)
I[1]= 0 #  initial conditions for I
R <- vector(mode="numeric", length = 192)
R[1]= 0 #  initial conditions for R
Cases_Infected <- vector(mode="numeric", length = 192)
Cases_Infected[1] <- abs(m) #  initial conditions

DataSEIR <- list(n_days = n_days, S0 = S , E0 = E, I0 = I, R0 = R,
                 N = N, tcon = tcontrol, C = Cases_Sym_onset,
                 D = Cases_Removed, B0 = Cases_Infected)

data {
  int<lower=1> n_days;// total # of time point 
  int S0[n_days]; //of susceptible individuals at the start of the epidemic 
  int E0[n_days]; //of exposed individuals at the start of the epidemic
  int I0[n_days]; //of infectious individuals at the start of the epidemic
  int R0[n_days]; //of removed individuals at the start of the epidemic
  int N;          //population size
  int C[n_days];  // empirical data
  int D[n_days];  // empirical data
  int tcon;       // time point where control measures are introduced 
  int B0[n_days]; // initial conditions for unobserved entity. 
}
parameters {
  real<lower=0,upper=1> beta;
  real<lower=1/10.5,upper=1/3.5> gamma;
  real<lower=1/21.0,upper=1> rho;
  real<lower=0> q;
}
transformed parameters{
  real theta[4];
  theta[1]= beta;
  theta[2]= gamma;
  theta[3]= rho;
  theta[4]= q;
}
model{
  int S[n_days]; // # of susceptible individuals
  int E[n_days]; // # of exposed individuals
  int I[n_days]; // # of infectious individuals
  int R[n_days]; // # of removed individuals
  real P[n_days]; //  probability of leaving S compartment
  int B[n_days]; // unobserved entity. # of susceptitble people who become infected
  real Beta[n_days];
  real pc; //probability of leaving E compartment
  real pr; //probability of leaving I compartment
  
  S[1]= S0[1] - B0[1];
  E[1]= E0[1] + B0[1] - C[1];
  I[1]= I0[1] + C[1] - D[1];
  R[1]= N - S[1] - E[1] - I[1];
  
  beta ~ gamma(2,10); // prior
  gamma ~  gamma(2,14); // prior
  rho ~  gamma(2,10); // prior
  q ~ gamma(2,10); // prior
  
  for(t in 2:n_days){
      if(t >= tcon){
        Beta[t]= beta*exp(-q*(t-tcon));
      }
      else{
        Beta[t]= beta;
      }
      P[t]= 1 - exp(-(Beta[t]*I[t]/N));
      pc = 1 - exp(-rho);
      pr = 1- exp(-gamma);
      S[t]= S[t-1] - B[t-1];
      E[t]= E[t-1] + B[t-1] - C[t-1];
      I[t]= I[t-1] + C[t-1] - D[t-1];
      R[t] = N - S[t] - E[t] - I[t];
      
      B[t] ~ binomial(S[t],P[t]); //  
      C[t] ~ binomial(E[t],pc); // # of cases by date of symptom onset
      D[t] ~ binomial(I[t],pr); // # of cases removed from infectious class
    }
}",
"stan_model1.stan")
stan_model1 <- "stan_model1.stan"

parameter= c("beta","gamma","rho","q")
ini= function(){list(beta = 0.2, gamma = 0.143, rho=0.2, q = 0.2)}
fit_seir <- stan(file=stan_model1,
                 data = DataSEIR,
                 pars = parameter,
                 init = ini,
                 iter = 1000,
                 chains = 1)

The values used in the model are according to the paper.
The error that we get is

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: binomial_lpmf: Successes variable is -2147483648, but must be in the interval [0, -2142118856]  (in 'model4adc26af40d5_stan_model1' at line 63)

Chain 1: 
Chain 1: Initialization between (-2, 2) failed after 100 attempts. 
Chain 1:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."
[1] "error occurred during calling the sampler; sampling not done."

Going through posts with similar errors, we do understand that the prior is not initialized or drawn from default rather than the one specified but no matter what changes we make, we get the same error.
If anymore can point us where we are going wrong or, if we are entirely wrong in the modelling would be much appreciated.
Thank you for reading the post.
Best regards.

Update:
I added the initialisation of B in the model as suggested by few comments but I still get the same error.

write( "

data {
  int<lower=1> n_days;// total # of time point 
  int S0[n_days]; //of susceptible individuals at the start of the epidemic 
  int E0[n_days]; //of exposed individuals at the start of the epidemic
  int I0[n_days]; //of infectious individuals at the start of the epidemic
  int R0[n_days]; //of removed individuals at the start of the epidemic
  int N;          //population size
  int C[n_days];  // empirical data
  int D[n_days];  // empirical data
  int tcon;       // time point where control measures are introduced 
  int B0[n_days]; // initial conditions for unobserved entity. 
}
parameters {
  real<lower=0,upper=1> beta;
  real<lower=1/10.5,upper=1/3.5> gamma;
  real<lower=1/21.0,upper=1> rho;
  real<lower=0> q;
}
transformed parameters{
  real theta[4];
  theta[1]= beta;
  theta[2]= gamma;
  theta[3]= rho;
  theta[4]= q;
}
model{
  int S[n_days]; // # of susceptible individuals
  int E[n_days]; // # of exposed individuals
  int I[n_days]; // # of infectious individuals
  int R[n_days]; // # of removed individuals
  real P[n_days]; //  probability of leaving S compartment
  int B[n_days]; // unobserved entity. # of susceptitble people who become infected
  real Beta[n_days];
  real pc; //probability of leaving E compartment
  real pr; //probability of leaving I compartment
  
  B[1]= B0[1];
  S[1]= S0[1] - B0[1];
  E[1]= E0[1] + B0[1] - C[1];
  I[1]= I0[1] + C[1] - D[1];
  R[1]= N - S[1] - E[1] - I[1];
  
  beta ~ gamma(2,10); // prior
  gamma ~  gamma(2,14); // prior
  rho ~  gamma(2,10); // prior
  q ~ gamma(2,10); // prior
  
  for(t in 2:n_days){
    if(t >= tcon){
      Beta[t]= beta*exp(-q*(t-tcon));
    }
    else{
      Beta[t]= beta;
    }
    
    S[t]= S[t-1] - B[t-1];
    E[t]= E[t-1] + B[t-1] - C[t-1];
    I[t]= I[t-1] + C[t-1] - D[t-1];
    R[t] = N - S[t] - E[t] - I[t];
    
    P[t]= 1 - exp(-(Beta[t]*I[t]/N));
    pc = 1 - exp(-rho);
    pr = 1- exp(-gamma);
    
    B[t] ~ binomial(S[t],P[t]); //  
    C[t] ~ binomial(E[t],pc); // # of cases by date of symptom onset
    D[t] ~ binomial(I[t],pr); // # of cases removed from infectious class
  }
}", 
"stan_model1.stan")

Error

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: binomial_lpmf: Successes variable is -2147483648, but must be in the interval [0, 5363916]  (in 'model2d9c5c0523ff_stan_model1' at line 65)

Chain 1: 
Chain 1: Initialization between (-2, 2) failed after 100 attempts. 
Chain 1:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."
[1] "error occurred during calling the sampler; sampling not done"
2 Likes

I got the same error and also changed my initial values multiple times…

Just to clarify: the empirical data comes from the R library “outbreaks” and our main issue is regarding the sampling of the unobserved/latent variable B and it’s initialization.

This line seems to suggest that the data input is incorrect in some way. Probably too big an integer.

This value is sampled from B[t] ~ binomial(S[t],P[t]) and we cannot figure out why it is negative, while both S[t] and P[t] are positive

Can B or S ever hold a value larger than 2,147,483,647? The largest value an integer can hold in stan is 2,147,483,647, so adding one to that causes rollover to the largest negative value of -2,147,483,648.

2 Likes

You didn’t initialize your B[1], so S[2] is assigned with garbage.

3 Likes

S is initialised with the value 5364500. But even if we change it to smth like 500, we get the same mistake. B should be in the range of 0 to S as well.

We did initialize B in the DataSEIR, where B0=Cases_Infected and Cases_Infected[1]=m (m is the number specified by the authors of the paper). We also tried specifying B[1] explicitly in the model, got the same result

I think you only initialise B0, not B. You do need to specify B[1] as @yizhang said, but still, the rest of the array remains unassigned and thus equal to -2147483648 (default for int array). I think this explains the error you’re getting.

I also think there might be a confusion on the meaning of :

B[t] ~ binomial(S[t],P[t]);

This doesn’t sample values from the binomial to assign to B, which is why B remains unassigned. It only increments the log joint probability used by the Stan sampler (see here for more Info) by the log probability of having B[t] successes given S[t] and P[t]. B[2] for instance being -2147483648, it triggers the error you’re seeing.

1 Like

Thanks a lot for your comment @LeoGrin. We indeed thought that by

we are sampling values from the binomial distribution to assign to B. As the authors of the paper we are trying to replicate suggest, sampling should be done in 2 stages:

image

We initialised B (also added B[1] values as was suggested), and initialised our parameter vector. Now we are confused regarding how to update B from B | C, D, theta( our parameter vector). Maybe you could point us in the right direction of how to fix the code and actually assign values to B?

In a.stan file you define a model through its log posterior density. The lines with ~ are doing that. Then you can fit it with some of the MCMC algorithms available in stan. You should not try to implement an inference method or your own MCMC algorithm in a stan file. If you need random binomial numbers in a stan file, then you would use B[t] = binomial_rng(... but the rng functions can only be used in the transformed data and generated quantities blocks, because if you would use them anywhere else, you are most likely doing what you should not be doing.

So, in conclusion, I think you are trying to implement an MCMC method from that paper, and in Stan you cannot implement custom MCMC methods, and the .stan file is only supposed to define a model, not inference algorithm.

3 Likes

It’s perfectly ok to do what you’re doing, which is similar to manually discretizing a stochastic differential equation model, like in Stan example of volatility model. The problem is, as @LeoGrin pointed out, sampling is not equivalent to assignment, and your local variables in Model block remains undefined. What needs to be done is to move B/C/D to Parameters so their samples could be used in next time step.

1 Like

Inference for that model with Stan is probably possible, but if the intention was to replicate the original paper, as the original poster said, then that’s not possible because they seem to use different kind of MCMC than Stan.

I think @jtimonen is right you cannot do this problem in Stan. This is because B is a discrete parameter and Stan cannot sample discrete parameters. Also, according to the paper, you need a very specific proposal to sample B because it is subject to many constraints.

One thing you can do is assume a constant latent period to simplify the problem. The latent period for ebola is about 8-10 days. So, assume each individual shows symptoms 9 days after infection. Then you can construct B. Once B is known you can sample the rest of the parameters in Stan using the likelihood they have (8).

Only alternative I can see is to try to marginalize out B (I don’t know if that is possible) or use a different MCMC program that allows custom samplers (like Nimble).

1 Like

@jtimonen & @dirkdouwes is right. One will have to marginalize B out. I was mistaken not noticing it’s Binomial.

I was able to fit the model by replacing the binomial with a continuous approximation. I make no claims of correctness or appropriateness. Program is below.

functions{
  real cont_binomial_lpdf(real x, real n, real theta){
    real ans = lchoose(n, x)  + x*log(theta) + (n-x)*log1m(theta);
    if(x > n) return(negative_infinity());
    return(ans);
  }
}
data {
  int<lower=1> n_days;// total # of time point 
  int S0[n_days]; //of susceptible individuals at the start of the epidemic 
  int E0[n_days]; //of exposed individuals at the start of the epidemic
  int I0[n_days]; //of infectious individuals at the start of the epidemic
  int R0[n_days]; //of removed individuals at the start of the epidemic
  int N;          //population size
  int C[n_days];  // empirical data
  int D[n_days];  // empirical data
  int tcon;       // time point where control measures are introduced 
  int B0[n_days]; // initial conditions for unobserved entity. 
}
parameters {
  real<lower=0,upper=1> beta;
  real<lower=1/10.5,upper=1/3.5> gamma;
  real<lower=1/21.0,upper=1> rho;
  real<lower=0> q;
  real<lower=0> B[n_days]; // unobserved entity. # of susceptitble people who become infected
}
transformed parameters{
  real<lower=0> S[n_days]; // # of susceptible individuals
  real<lower=0> E[n_days]; // # of exposed individuals
  real<lower=0> I[n_days]; // # of infectious individuals
  real R[n_days]; // # of removed individuals
  real<lower=0> P[n_days]; //  probability of leaving S compartment
  real<lower=0> Beta[n_days];
  real<lower=0, upper=1> pc; //probability of leaving E compartment
  real<lower=0, upper=1> pr; //probability of leaving I compartment
  
  pc = 1 - exp(-rho);
  pr = 1- exp(-gamma);
  
  S[1] = S0[1] - B0[1];
  E[1] = E0[1] + B0[1] - C[1];
  I[1] = I0[1] + C[1] - D[1];
  R[1] = N - S[1] - E[1] - I[1];
  Beta[1] =  beta*exp(q*tcon); // dubious, please check
  P[1] = 1 - exp(-(Beta[1]*I[1]/N));
  
  for(t in 2:n_days){
    if(t >= tcon){
      Beta[t] = beta*exp(-q*(t-tcon));
    }
    else{
      Beta[t] = beta;
    }
    S[t]= S[t-1] - B[t-1];
    E[t]= E[t-1] + B[t-1] - C[t-1];
    I[t]= I[t-1] + C[t-1] - D[t-1];
    P[t] = 1 - exp(-(Beta[t]*I[t]/N));
    R[t] = N - S[t] - E[t] - I[t];
  }
  // print(R);
}
model{
  beta ~ gamma(2,10); // prior
  gamma ~  gamma(2,14); // prior
  rho ~  gamma(2,10); // prior
  q ~ gamma(2,10); // prior
  for(t in 2:n_days){
    B[t] ~ cont_binomial(S[t], P[t]); //  
    C[t] ~ cont_binomial(E[t], pc); // # of cases by date of symptom onset
    D[t] ~ cont_binomial(I[t], pr); // # of cases removed from infectious class
  }
}

3 Likes

This is what I was suggesting but as @jtimonen pointed out even if B can be sampled this is still not a replicate of the paper’s algorithm, neither is @Anastasiia_Malyk 's code example.

1 Like

Absolutely correct. Depending on exactly why the replication was being attempted though, might be useful. Hopefully it’ll be useful to anyone searching for “stochastic SEIR” on the forums as well.

1 Like