Add death rate to the model

giulianobelinassi commented 4 years ago

This commit adds the death rate to the SIR model by decomposing

\gamma = a + b

then

ds/dr = - \betas(t)i(t) di/dr = \betas(t)i(t) - (a + b)*i(t) dr/dt = ai(t) dk/dt = bi(t)

where 'a' is the recovery rate and 'b' is the death rate.

Please see line 224. I don't know how to add l3 correctly into the equation correctly.

See #13

Signed-off-by: Giuliano Belinassi giuliano.belinassi@usp.br

giulianobelinassi commented 4 years ago

I just noticed that you used

return alpha * l1 + (1 - alpha) l2

To balance which value you want to give priority to.

To generalize for 3 variables, it is possible to use the convex closure. One possibility is

u = 0.1
v = 0.3
w = 1 - u - v
return u*l1 + v+l2 + w*l3

But I am not sure if this is the best value for this simulation.

bolinches commented 4 years ago

My personal view is I would rather not include death estimation. Besides the pollution of data on or if SIR is acceptable to this it can also send some wrong message. Just my 2 cents.

bolinches commented 4 years ago

you might find interesting https://github.com/ImperialCollegeLondon/covid19model they model deaths and how the different measures affect the count.

gasilva commented 4 years ago

My personal view is I would rather not include death estimation. Besides the pollution of data on or if SIR is acceptable to this it can also send some wrong message. Just my 2 cents.

But the deads with coronavirus is the most accurate data because they are tested after death. The number of cases or recovery are not at all. There is a lot of under notification.

I validated the model for China, Italy, UK and Canada. After I applied to Brazil. See attached.

In order to match China I had to make recovered initial condition negative. If the model would have a delay, it could be simulated without problems.

giulianobelinassi commented 4 years ago

Hi, Guilherme

That is very interesting!

Yesterday I implemented what I think is a better way to calculate the death rate. I start from the same equation about the death rate 'dk/dt', but then I use algebraic manipulation to find how to decompose it into 'cured' and 'death'.

It is implemented in this branch. I will provide a better documentation soon.

Thank you, Giuliano.

giulianobelinassi commented 4 years ago

Hi,

I have just written a document explaining the model. Check this branch.

Thank you, Giuliano.

gasilva commented 4 years ago

Great ...do you have plots for various countries? @giulianobelinassi

If you need to use the country with provinces, I implemented a function to sum province data under the same country. You can used it in place of remove_provinces.

see it here!

def sumCases_province(input_file, output_file):
    with open(input_file, "r") as read_obj, open(output_file,'w',newline='') as write_obj:
        csv_reader = reader(read_obj)
        csv_writer = writer(write_obj)

        lines=[]
        for line in csv_reader:
            lines.append(line)    

        i=0
        ix=0
        for i in range(0,len(lines[:])-1):
            if lines[i][1]==lines[i+1][1]:
                if ix==0:
                    ix=i
                lines[ix][4:] = np.asfarray(lines[ix][4:],float)+np.asfarray(lines[i+1][4:] ,float)
            else:
                if not ix==0:
                    lines[ix][0]=""
                    csv_writer.writerow(lines[ix])
                    ix=0
                else:
                    csv_writer.writerow(lines[i])
            i+=1

You will need to import some modules:

from csv import reader
from csv import writer

gasilva commented 4 years ago

@giulianobelinassi here are some improvements possible:

Use scipy.integrate.odeint to solve the ODE system. The results may differ from solve_ivp, which is common for nonlinear equations, but the code is very fast because is FORTRAN compiled. The solver is used by several people to solve Lotka Volterra equations (predator-prey type)

def lossOdeint(point, data, recovered, death, s_0, i_0, r_0, d_0):
    size = len(data)
    beta, a, b = point
    def SIR(y,t):
        return [-beta*y[0]*y[1], beta*y[0]*y[1]-(a+b)*y[1], a*y[1], b*y[1]]
    y0=[s_0,i_0,r_0,d_0]
    tspan=np.arange(0, size, 1)
    res=odeint(SIR,y0,tspan)
    l1 = np.sqrt(np.mean((res[:,1]- data)**2))
    l2 = np.sqrt(np.mean((res[:,2]- recovered)**2))
    l3 = np.sqrt(np.mean((res[:,3] - death)**2))
    #weight for cases
    u = 0.25
    #weight for recovered
    v = 0.02 ##Brazil France 0.04 US 0.02 China 0.01 (it has a lag in recoveries) Others 0.15
    #weight for deaths
    w = 1 - u - v - z
    return u*l1 + v*l2 + w*l3

Use different delays/lags in equations of recovered and deaths due to incubation, hospital time. China data shows a very clear delay in recovery for example. The optimizer must find the $\tau_{rec}$ and $\tau_{dead}$ in order to start the curves in the right time.
Implement an improvement in infected people because Covid19 has a lot of undetectable cases. See that at See paper here! A Time-dependent SIR model for COVID-19 with Undetectable Infected Persons from Cornell University.

giulianobelinassi commented 4 years ago

Hi, Guilherme.

Thank you for all this info. I will take a look closely

Could you provide me the parameters you used to generate the China graphic? (S_0, I_0, R_0, D_0). I see that your 'recovered' starts at a negative value. I am trying to reproduce it with my other implementation.

Thank you, Giuliano.

gasilva commented 4 years ago

country1="China"

    if country1=="Brazil":
        date="3/3/20"
        s0=25000
        i0=27
        r0=-35
        k0=-35

    if country1=="China":
        date="1/22/20"
        s0=170000
        i0=1200
        r0=-80000
        k0=200

    if country1=="Italy":
        date="1/31/20"
        s0=160000
        i0=23
        r0=15
        k0=100

    if country1=="France":
        date="2/25/20"
        s0=95e3
        i0=250
        r0=-75
        k0=105

    if country1=="United Kingdom":
        date="2/25/20"
        s0=80000
        i0=22
        r0=-5
        k0=-50

    if country1=="US":
        date="2/25/20"
        s0=470000
        i0=10
        r0=-50
        k0=50

gasilva commented 4 years ago

you might find interesting https://github.com/ImperialCollegeLondon/covid19model they model deaths and how the different measures affect the count.

I tried to run their model in R but I could not. A lot of errors.

giulianobelinassi commented 4 years ago

Hi, Guilherme

Please, checkout to branch death_rate_2

git checkout death_rate_2
git pull

Giuliano.

gasilva commented 4 years ago

Hi, Guilherme

Please, checkout to branch death_rate_2
git checkout death_rate_2
git pull
Giuliano.

Thanks

itsisthatis commented 4 years ago

@gasilva Can I get the whole code for generating these graphs ??

gasilva commented 4 years ago

@gasilva Can I get the whole code for generating these graphs ??

I will take a look and try to find it...a lot of codes here

Mamitiana2 commented 4 years ago

Hi every one, i'm a student, i want to fit beta an gamma in my model SIR to an data. Can any one help me?

Lewuathe / COVID19-SIR

Add death rate to the model #14