[PYTHON] I tried to predict the number of people infected with coronavirus in Japan by the method of the latest paper in China

at first

I am not an infectious disease expert, so please read it with an understanding.

The new type of corona pneumonia (Covid-19) that occurred in Wuhan, Hubei Province, China from December 2019 has spread to Japan and the number of infected people is increasing. I am interested in how the number of infected people in Japan will increase in the future. Therefore, I searched for papers on infection prediction models. Infection models have already been announced from all over the world, but parameters such as infection rate and quarantine rate have a great influence on the accuracy when predicting with the infection model. If these parameters differ from the actual values, the prediction will be far from the actual situation.

The other day, a Chinese researcher announced a model for predicting the number of people infected with coronavirus. COVID-19 in Japan: What could happen in the future?

In this paper, we predict that the number of infected people can be predicted accurately by applying the prediction model to various parts of China (Wuhan, Beijing, Shanghai ...), and that the number of infected people will increase in the future by applying it to Japan. I am. Also, in this paper, the infection rate and quarantine rate, which are the parameters necessary for prediction, were announced (although it would be correct to say that the parameters were adjusted to match the actual number of infected people).

Predictive model

The SIR and SEIR models are used in the infectious disease model, but this paper uses a slightly different model. We call it the statistical time delay dynamic model.

[Quote COVID-19 in Japan: What could happen in the future?]

The calculation formula is as follows. t: day I (t): Cumulative number of infected people J (t): Cumulative number of infected people (confirmed onset at hospital) G (t): Infected person who occurred at that time (not the cumulative value, the onset is not confirmed) I0 (t): Number of potential infected persons (infected but not confirmed or quarantined) 　　　　　　I0(t) = I(t) − J(t) − G(t) [Quote COVID-19 in Japan: What could happen in the future?]

Parameters

There are important parameters for using this model: incidence and hospitalization. As for the parameters, a Chinese example was described in this paper.

r Growth rate

Growth rate is the rate at which one person infects another. This was set as follows based on the figures in China in the paper. 　　r = 0.3

l Isolation rate

The quarantine rate is the rate at which infected people are quarantined. This was set as follows based on the figures in China in the paper. l = 0.1 (until 2020/2/28) l = 0.5 (from February 29, 2020)

Incidence rate

f2 (t) Probability of transition from infection to onset This was not mentioned in the paper. I have heard the news that 70-80% of people do not develop the disease even if they are infected. In addition, since it takes 14 days at the longest to develop the disease, we set it as follows. 　f2(t) = 0.2/14 × t (t < 14) 　f2(t) = 0.2 (t >= 14)

Hospitalization rate

f4 (t) Probability of transition from infection to hospitalization The rate of infection and hospitalization is completely unknown. Therefore, assuming that 1/4 of the affected people will be hospitalized, the settings are as follows. 　f4(t) = 0.05/14 × t (t < 14) 　f4(t) = 0.05 (t >= 14)

program

predict.py

from __future__ import print_function
import numpy as np
import pandas as pd

default_output  = 'predict.csv'

class  Corona():
    def __init__(self, max_day):
        self.r  = 0.3 #growth rate
        self.tl = 15  #Quarantine start date
        self.l1 = 0.1 #Isolation rate before the quarantine start date
        self.l2 = 0.5 #Isolation rate after the quarantine start date
        self.i  = np.zeros(max_day + 1)
        self.i0 = np.zeros(max_day + 1)
        self.j  = np.zeros(max_day + 1)
        self.g  = np.zeros(max_day + 1)

    def set_start(self):
        # self.j[0]  = 11 # 2020/1/30
        self.j[0]  = 21 # 2020/2/14
        self.i[0]  = self.j[0] * 10
        self.g[0]  = 0
        self.i0[0] = self.func_i0(0)

    def f2(self, t):
        #Incidence rate
        if t < 14:
            a = 0.2/14.0
            b = 0.0
            y = a * t + b
        else:
            y = 0.2
        return y

    def f4(self, t):
        #Hospitalization rate
        if t < 14:
            a = 0.05/14.0
            b = 0.0
            y = a * t + b
        else:
            y = 0.05
        return y

    def func_l(self, t):
        #Isolation rate
        if t < self.tl:
            return self.l1
        else:
            return self.l2

    def func_i(self, t):
        #Cumulative number of infected people
        # I(t + 1) = I(t) + r I0(t),
        new_i = self.i[t] + self.r * self.i0[t]
        return new_i

    def func_j(self, t):
        #Cumulative number of infected people(Confirmed at the hospital)
        # J(t + 1) = J(t) + r Σs<t f4(t - s) I0(s)
        sum1 = 0
        for s in range(t):
            sum1 += self.f4(t - s) * self.i0[s]
        new_j = self.j[t] + self.r * sum1
        return new_j

    def func_g(self, t):
        #Infected person who occurred momentarily(Not confirmed to be infected at the hospital)
        # G(t + 1) = G(t) + f2(t) Σs<t f2(t - s) I0(s) - Σs<t f4(t - s) I0(s).
        sum1 = 0
        sum2 = 0
        for s in range(t):
            sum1 += self.f2(t - s) * self.i0[s]
        for s in range(t):
            sum2 += self.f4(t - s) * self.i0[s]
        new_g = self.g[t] + self.func_l(t) * sum1 - self.func_l(t) * sum2
        return new_g

    def func_i0(self, t):
        #Number of potential infections(Infected but not confirmed or quarantined)
        # I0(t) := I(t) - J(t) - G(t)
        new_i0 = self.i[t] - self.j[t] - self.g[t]
        if new_i0 < 0.0:
            new_i0 = 0.0
        return new_i0

    def predict(self, day):
        #Initialization
        period = day + 1
        predict_data = np.zeros([period, 5])
        df_predict = pd.DataFrame(predict_data, columns=['day', 'I', 'J', 'G', 'I0'])
        self.set_start()

        #Forecast
        for i in range(period - 1):
            self.i[i+1]  = self.func_i(i)
            self.j[i+1]  = self.func_j(i)
            self.g[i+1]  = self.func_g(i)
            self.i0[i+1] = self.func_i0(i)

            df_predict.loc[i, 'day'] = i+1
            df_predict.loc[i, 'I']   = self.i[i+1]
            df_predict.loc[i, 'J']   = self.j[i+1]
            df_predict.loc[i, 'G']   = self.g[i+1]
            df_predict.loc[i, 'I0']  = self.i0[i+1]

        return df_predict

def main():
    corona = Corona(25)
    predict = corona.predict(25)
    predict.to_csv(default_output, index=False)
 
if __name__ == "__main__":
    main()

result

Number of infected people (those who have been confirmed to be infected at the hospital)

We obtained the number of infected people at the start of the simulation from the materials released daily in Ministry of Health, Labor and Welfare Press Release. Predicted from 2/15 based on 21 infected people on 2/14/2020. The forecast results are as follows.

The error is large on the way, but on 2/28, the error was only 5 people. (Maybe it happens ...)

Number of potential infections (those who have not been confirmed to be infected at the hospital)

The above number of infected people is only those who have been confirmed to be infected at the hospital. The simulation also calculates the number of potential infections. The result is shown in the figure below. Orange: People confirmed to be infected at the hospital (cumulative) Navy blue: Hidden infected persons (cumulative) whose infection has not been confirmed at the hospital There are far more hidden infected people than those who have been confirmed to be infected at the hospital, about 10 times as many.

Impressions

Is the Ministry of Health, Labor and Welfare also conducting such a simulation?

References

https://www.medrxiv.org/content/10.1101/2020.02.21.20026070v2

2020/3/3 change

Program fix

Others have written the Hatena Blog using the program in this post. 　　https://kibashiri.hatenablog.com/entry/2020/03/02/171223 On his blog, he pointed out a mistake in the program. It was certainly a sign error in the calculation of func_g (). Fixed the program and results.

Prediction result

It is predicted that the number of infected people after 2/29 will increase by nearly 100 every day. The Ministry of Health, Labor and Welfare has announced a new number of infected people (domestic cases-excluding charter flight returnees), so I compared it with the forecast results.

The predicted value is significantly different from the actual number of infected people, and it seems that the predicted performance after 2/29 was not good. I tried to predict the number of infected people in Japan with the same parameters as in China, but I found that the situation is different between Japan and China, and there is a limit to the prediction with the same parameters.

According to the media, the PCR testing capacity of Japan was reported to be 3800 people / day, but the actual number of people tested per day is 130, 178, 96, and 71, which is unexpectedly small.