stan-dev / example-models

Example models for Stan
http://mc-stan.org/
772 stars 479 forks source link

Add case study on transmission models. #173

Closed charlesm93 closed 4 years ago

charlesm93 commented 4 years ago

Add source files for case study on disease transmission models. Complements this PR on Stan's website.

ericnovik commented 4 years ago

Cool! The ODE code would be a little easier to read this way:

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

      real S = y[1];
      real I = y[2];
      real R = y[3];

      real beta = theta[1];
      real gamma = theta[2];

      real dS_dt = -beta * I * S;
      real dI_dt =  beta * I * S - gamma * I;
      real dR_dt =  gamma * I;

      return {dS_dt, dI_dt, dR_dt};
  }
}
charlesm93 commented 4 years ago

thanks, @ericnovik. Duely noted.

bob-carpenter commented 4 years ago

I'd suggest not repeating components. This also makes it clear that one of the three components is redundant given that they sum to a constant. And it's a bit more efficient by eliminating a redundant compountation. I really don't know why the SIR models aren't typically written like this:

real[] sir(real t, real[] y, real[] theta, real[] x_r, int[] x_i) {
      real S = y[1];
      real I = y[2];
      real R = y[3];
      real beta = theta[1];
      real gamma = theta[2];
      real dS_dt = -beta * I * S;
      real dR_dt =  gamma * I;
      real dI_dt =  -(dS_dt + dR_dt);
      return {dS_dt, dI_dt, dR_dt};
  }
elizavetasemenova commented 4 years ago

The advantage of keeping the compartment explicit real dR_dt = gamma * I; is that it allows adding vital dynamics, i.e. adding births and deaths, in a more straightforward way. These modifications would break the total volume of the flow and compartment sizes wouldn't sum up to a constant anymore. Since this is a case study and not a model aimed for production, probably the didactic angle is preferable.

bob-carpenter commented 4 years ago

@elizavetasemenova: Thanks for the feedback! That line real dR_dt = gamma * I is in all the suggested models. Were you referring to not making dI_dt just the negation of the sum of the other compartments or something else? That would make sense with birth and deaths making the total number of subjects vary over time.

I was just suggesting a way to make it easier to understand the basic SIR model by reusing terms, which I now understand might be counterproductive in a larger context. I'd still suggest reusing terms even if the I component isn't so simple.