Priesemann-Group / covid19_inference_forecast

GNU General Public License v3.0
178 stars 53 forks source link

Lambda parameter in the SIR system #11

Closed truEngineer closed 4 years ago

truEngineer commented 4 years ago

Hello. Having studied your article and the corresponding code, I had questions about the λ parameter. The article says that λ is a spreading rate (infection rate, FIG.1 in article). Let's look at the code now

def SIR_model(λ, μ, S_begin, I_begin, N):
    new_I_0 = tt.zeros_like(I_begin)
    def next_day(λ, S_t, I_t, _):
        new_I_t = λ/N*I_t*S_t
        S_t = S_t - new_I_t
        I_t = I_t + new_I_t - μ * I_t
        return S_t, I_t, new_I_t
    outputs , _  = theano.scan(fn=next_day, sequences=[λ], 
                               outputs_info=[S_begin, I_begin, new_I_0])
    S_all, I_all, new_I_all = outputs
    return S_all, I_all, new_I_all
# fraction of people that are newly infected each day
λ = pm.Lognormal("λ", mu=np.log(0.4), sigma=0.5)

This code corresponds to a system of ODEs

The λ parameter here, according to the literature, is not infection rate (0<λ<1), but effective contact rate (λ>0)!

If λ is infection rate (Murray, p.320) or transmission rate (per capita), the system of equations should look like this:

Maybe I'm wrong, but in my opinion it greatly affects the interpretation of the model results.

This is also important when calculating R0. The basic reproduction number R0 is calculated as *infection_rateN/recovery_rate=effective_contact_rate/recovery_rate. In article R0 is calculated as infection_rate/recovery_rate** (0.41/0.12=3.4).

jdehning commented 4 years ago

Thanks for your interest, but I would say that it is more or less correct what we did, at least on how we calculate R0. Look at the Wikipedia article about the SIR model, that's basically the equations we implemented: https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model_without_vital_dynamics

I agree that our lambda is the effective contact rate of the first article you cited literature and not the infection rate of Murray. There seems to be different names for the parameter lambda, I found transmission rate/contact rate and infection rate.

truEngineer commented 4 years ago

Yes that's right. In Wikipedia (The SIR model without vital dynamics), R0 is calculated by the formula effective_contact_rate/recovery_rate=Tr/Tc, as I wrote earlier.

I see that you slightly corrected the article by deleting "infection rate" (FIG.1). Now everything is clear. λ parameter is effective contact rate in your case (or spreading rate in your article), R0 estimate is correct. Thanks for clarifying.

There is one more question. Informative prior for spreading rate λ is a log-normal distribution with mu=log(0.4) and sigma=0.5, informative prior for recovery rate μ is a log-normal distribution with mu=log(0.125) and sigma=0.2:

with pm.Model() as model:
    ...
    # fraction of people that are newly infected each day
    λ = pm.Lognormal("λ", mu=np.log(0.4), sigma=0.5)

    # fraction of people that recover each day, recovery rate mu
    μ = pm.Lognormal('μ', mu=np.log(1/8), sigma=0.2)
    ...

From your article "As median estimates, we obtain for the spreading rate λ = 0.41, µ = 0.12" (FIG. 1 C,E for the posterior distributions and the 95% credible intervals).

You write that "We chose informative priors based on available knowledge for λ, µ". What is this knowledge? If, for example, you need to customize your model for another country, then where to look for this information?

jdehning commented 4 years ago

Our interpretation, is that the μ shouldn't change from country to country, as it represents how long the infectious period is. The lambda should be chosen to approximately represent how strong the interventions/mitigation of the spread in a country is, but can the prior can also be chosen to be relatively uninformative. A median of 0.4 is about the unmitigated spread.

truEngineer commented 4 years ago

Thanks for clarifying