[PYTHON] Let's calculate the transition of the basic reproduction number of the new coronavirus by prefecture
Introduction
The center of infection for the new coronavirus infection (COVID-19) is already shifting from mainland China to Europe and the United States, and there are emergencies such as restricting free border movement within the EU and issuing a stay-at-home order. It's a situation. In Japan, the threat of the new coronavirus began to be recognized when the Diamond Princess arrived at Yokohama Port on February 3, 2020, the government announced its basic policy on February 25, and in Hokkaido on February 28. Although it was triggered by a state of emergency, concrete measures include a request to refrain from large-scale events from February 26, and temporary closure of elementary, junior high and high schools nationwide from March 2. About half a month has passed since then, and during that time, the Nikkei Stock Average has plummeted from the 22,000 yen level to the 17,000 yen level, and there is a dark cloud in the future. The effects of the new coronavirus infection have seriously affected not only health but also socio-economics, and when it will be resolved is of great concern. Therefore, in this article, I would like to present one method for calculating how the basic reproduction number, which is an index of the infectivity of the new coronavirus in Japan, has changed from the end of January to the beginning of March. I will. The calculation formula and code are a little long, so it may be better to look at the calculation result first.
What is the basic reproduction number?
The basic reproduction number seems to be very profound, [Research Institute for Mathematical Sciences, Kyoto University](https://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/140828/ If you look at 1 / 653-04.pdf) and Society for Mathematical Biology Letter, you can see the transition between generations. It doesn't seem simple, but in a nutshell,
- Expected number of secondary infections that a typical infected person in a population will reproduce during the entire period of infection
Seems to be defined as. Assuming that this number is R0, ** the condition for spreading infection is $ R_0> 1 $, and the condition for converging infection is $ R_0 <1 $ **.
British Prime Minister Boris Johnson explained at a meeting on March 16 that about 60% of the British people would control the peak until they gained herd immunity (Newsweek article. //www.newsweekjapan.jp/kimura/2020/03/post-74.php)), assuming that $ R_0 = 2.6 $, the herd immunity rate $ H $ in vaccines and autoimmunity is ,
(1-H)R_0 < 1
It seems that $ H> 1-1 / R_0 = 0.615 $ was assumed from the condition that satisfies. In other words, herd immunity relaxes the condition of $ R_0 $. However, the development of vaccines is expected to take more than a year, including clinical trials, and it is a bold strategy to acquire herd immunity by not controlling infection excessively.
Calculation model of basic reproduction number
Now, in order to calculate the basic reproduction number, I would like to consider a simple model based on the above definition and SEIR model. See the following figure.
The upper row shows the date, and the lower bar shows the infected population one day. The meaning of each symbol is
- E: Infectious disease is in the incubation period (Exposed)
- I: Infectious
- P: Positive test
- lp: Incubation period
- ip: Infection period
Represents. That is, the production of I to E on one day t is an infection from a population of I who was infected and developed before t, and the number of reproductions of t on that day is $ R_0 (t) / ip $. It is a model called. Expressing this as an expression,
\frac{1}{ip} R_0(t) \times \sum_{s=t+1}^{t+ip} P(s) = P(t+lp+ip), \\
\therefore R_0(t) = \frac{ ip \times P(t+lp+ip)} {\sum_{s=t+1}^{t+ip} P(s)}
It will be. Below, let's use this formula to calculate the time transition of the basic reproduction number $ R_0 (t) $ from the published data of test positives.
Original data
In this article, we use the csv data published in Map of the number of new coronavirus infections by prefecture (provided by Jag Japan Co., Ltd.). I was allowed to do it. The data of prefectures are collected in one, which makes it very easy to use.
Try to calculate with Python
Now, let's use Python to calculate the basic reproduction number R0 (t) by prefecture.
Prerequisites
As a prerequisite, use the following.
- lp: Incubation period = 5 [day]
- ip: Infection period = 8 [day]
- P in the above figure is the definitive diagnosis date.
There are various theories about lp and ip, but they are set as a guide. Also, due to the calculation formula, only $ R_0 (t) $ before lp + ip = 13 [day] can be calculated from the latest date of the data. Please note that on the graph, even if the value for the period from 13 days ago to the latest day is 0, it is meaningless.
code
First, import the library.
# coding: utf-8
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import datetime
import locale
Data by prefecture from csv provided in Map of the number of new coronavirus infections by prefecture (provided by Jag Japan Co., Ltd.) A function that reads and converts to a data frame.
def readCsvOfJapanPref(pref : None):
#Download from the URL below
# https://jag-japan.com/covid19map-readme/
fcsv = u'COVID-19.csv'
df = pd.read_csv(fcsv, header=0, encoding='utf8', parse_dates=[u'Confirmed date YYYYMMDD'])
#Fixed date,Extract only the prefectures where you receive a medical examination
df1 = df.loc[:,[u'Confirmed date YYYYMMDD',u'Consultation prefecture']]
df1.columns = ['date','pref']
#Extracted by consultation prefecture
if pref is not None: #If None, then the whole of Japan
df1 = df1[df1.pref == pref]
df1 = df1.loc[:,'date']
#Count on fixed date
df2 = pd.DataFrame( df1.value_counts() ).reset_index()
df2.columns = ['date','P']
df2 = df2.sort_values('date')
return df2
A function that defines a data frame for calculations. The column Ppre is for the denominator of the formula and Pat is for the numerator of the formula. Starting from January 24, 2020, we will create a calculation frame for days.
def makeCalcFrame(days):
t_1 = pd.Timestamp(2020,1,24) #Calculation start date
td = pd.Timedelta('1 days')
#
npd = [[t_1 + td * i, 0, 0, 0 ] for i in range(0,days)]
df1 = pd.DataFrame(npd)
df1.columns = ['date', 'Ppre','Pat', 'R0']
#
return df1
A function that combines prefectural data and calculation frames.
def mergeCalcFrame(df1, df2):
return pd.merge(df1, df2, on='date', how='left').fillna(0)
A function that calculates $ R_0 $ according to the formula. To simplify the calculation, I created a lambda expression that returns NaN when accessed outside the index range.
def calcR0(df, keys):
lp = keys['lp']
ip = keys['ip']
nrow = len(df)
getP = lambda s: df.loc[s, 'P'] if s < nrow else np.NaN
for t in range(nrow):
df.loc[t, 'Ppre'] = sum([ getP(s) for s in range(t+1, t + ip + 1)])
df.loc[t, 'Pat' ] = getP(t + lp + ip)
if df.loc[t, 'Ppre'] > 0:
df.loc[t, 'R0' ] = ip * df.loc[t, 'Pat'] / df.loc[t, 'Ppre']
else:
df.loc[t, 'R0' ] = np.NaN
return df
A function that graphs the results.
def showResult(df, title):
# R0=1 :Target for convergence
ptgt = pd.DataFrame([[df.iloc[0,0],1],[df.iloc[len(df)-1,0],1]])
ptgt.columns = ['date','target']
# show R0
plt.rcParams["font.size"] = 12
ax = df.plot(title=title,x='date',y='R0', figsize=(10,7))
ptgt.plot(x='date',y='target',style='r--',ax=ax)
ax.grid(True)
ax.set_ylim(0,)
plt.show()
This is the part that calculates for each prefecture.
def R0inJapanPref(pref, label):
keys = {'lp':5, 'ip':8 }
df1 = makeCalcFrame(60) # 60 days
df2 = readCsvOfJapanPref(pref)
df = mergeCalcFrame(df1, df2)
df = calcR0(df, keys)
showResult(df, 'COVID-19 R0 ({})'.format(label))
return df
preflist = [[None, 'Japan'], [u'Tokyo', 'Tokyo'],\
[u'Osaka', 'Osaka'], [u'Aichi prefecture', 'Aichi'],\
[u'Hokkaido', 'Hokkaido']]
dflist = [[R0inJapanPref(pref, label), label] for pref, label in preflist]
This is the part that displays the above calculation results together.
def showResult2(ax, df, label):
# show R0
plt.rcParams["font.size"] = 12
df1 = df.rename(columns={'R0':label})
df1.plot(x='date',y=label, ax=ax)
# R0=1
dfs = dflist[0][0]
ptgt = pd.DataFrame([[dfs.iloc[0,0],1],[dfs.iloc[len(dfs)-1,0],1]])
ptgt.columns = ['date','target']
ax = ptgt.plot(title='COVID-19 R0', x='date',y='target',style='r--', figsize=(10,8))
#
for df, label in dflist:
showResult2(ax, df, label)
#
ax.grid(True)
ax.set_ylim(0,10)
plt.show()
Calculation result
Now let's take a look at the calculation results.
Changes in the number of basic reproduction numbers in Japan as a whole
- There seems to be a large inflow of up to around 10 until February 10. This is probably due to the influence of Chinese tourists who came to Japan during the Chinese New Year (January 24th to January 30th).
- After that, although it once rises to around 3, it seems that it is staying around 1 (on average around 1.5).
Changes in the number of basic reproduction numbers in Tokyo
- It was a high value of up to 16 until around February 6, but it seems to be generally calm after February 7.
- However, there are occasional mountains of more than 3 from February 18th to February 23rd and March 1st, suggesting the possibility of clusters and inflows from overseas.
Changes in the number of basic reproduction numbers in Osaka Prefecture
- Up to 26 large peaks can be seen from February 17th to February 24th. This is probably a cluster outbreak.
- In fact, according to Takatsuki City Page, four live performances in Osaka City Participants in a live performance at the house from February 15th to February 24th have confirmed the outbreak of infected clusters.
Changes in the basic reproduction number in Aichi Prefecture
- Probably infected by overseas travelers by February 12th.
- After that, from February 18th to February 26th, you can see the mountains that change around 3. This seems to be a cluster that happened in gyms and welfare facilities, as reported.
- I am worried about cluster chaining because the peaks of the cluster are adjacent to each other.
Changes in the number of basic reproduction numbers in Hokkaido
- Up to 80 sharp peaks can be seen from February 8th to February 10th. This is probably a cluster that occurred during the Sapporo Snow Festival (February 4th to February 11th).
- Since February 16th, it has been around below 1, and it seems that it is heading toward convergence.
Compare the whole
- It is probable that there were infections mainly from overseas by February 18th, and sporadic outbreaks of clusters in Japan by March 1st.
- Since March 1st, it has generally settled down, and although it is still above the level of 1, it can be said that it is heading toward convergence.
Consideration
From the above, the following trends can be derived from the simulation regarding the estimation of the basic reproduction number R0 based on the infected person data of each prefecture in Japan.
- The inflow from overseas until mid-February has the highest R0, and then there is a tendency to sporadically jump due to the influence of what seems to be a domestic cluster.
- In other periods, especially in Hokkaido, we were able to keep the period below 1, which is generally in good condition.
- The cluster suspicion date of the calculation result and the reported cluster occurrence date (Osaka live house, Sapporo Snow Festival) match well, and there is a high possibility that the cluster occurrence can be guessed from the calculation model. (However, it seems that the findings of active epidemiological surveys are also effective.)
Furthermore ...
- From the simulation results, the basic reproduction number $ R_0 $ is not constant and fluctuates considerably.
- ** The main factors that increase the basic reproduction number $ R_0 $ are numerically supported by inflows from overseas and cluster generation. ** Measures to curb these seem to be restrictions on overseas travel and self-restraint of events and facilities that can cause clusters.
- The overall trend is that ** it has been converging since March 1st, and R0 is approaching 1 in Japan as a whole **, so it seems that the current measures are working.
- $ R_0 $ is the ratio of the number of positives to the number of positives, so even if the number of untested hidden infections is a certain number of positives, the effect will be minor.
Reference link
I referred to the following page. Map of the number of people infected with the new coronavirus by prefecture (provided by Jag Japan Co., Ltd.) Research Institute for Mathematical Sciences, Kyoto University Society for Mathematical Biology Letter Newsweek article SEIR model Takatsuki City Page
I haven't quoted it directly, but I'll add a very useful link. COVID-19 reports National Cluster Map