Data-Science-For-Beginners/2-Working-With-Data/07-python/notebook-covidspread.ipynb

Estimation of COVID-19 Pandemic

Loading Data

We will use data on COVID-19 infected individuals, provided by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University. Dataset is available in this GitHub Repository.

In [124]:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (10,3) # make figures larger

We can load the most recent data directly from GitHub using pd.read_csv. If for some reason the data is not available, you can always use the copy available locally in the data folder - just uncomment the line below that defines base_url:

In [125]:

base_url = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/" # loading from Internet
# base_url = "../../data/COVID/" # loading from disk
infected_dataset_url = base_url + "time_series_covid19_confirmed_global.csv"
recovered_dataset_url = base_url + "time_series_covid19_recovered_global.csv"
deaths_dataset_url = base_url + "time_series_covid19_deaths_global.csv"
countries_dataset_url = base_url + "../UID_ISO_FIPS_LookUp_Table.csv"

Let's now load the data for infected individuals and see how the data looks like:

In [126]:

infected = pd.read_csv(infected_dataset_url)
infected.head()

Out [126]:

	Province/State	Country/Region	Lat	Long	1/22/20	1/23/20	1/24/20	1/25/20	1/26/20	1/27/20	...	8/20/21	8/21/21	8/22/21	8/23/21	8/24/21	8/25/21	8/26/21	8/27/21	8/28/21	8/29/21
0	NaN	Afghanistan	33.93911	67.709953	0	0	0	0	0	0	...	152448	152448	152448	152583	152660	152722	152822	152960	152960	152960
1	NaN	Albania	41.15330	20.168300	0	0	0	0	0	0	...	138132	138790	139324	139721	140521	141365	142253	143174	144079	144847
2	NaN	Algeria	28.03390	1.659600	0	0	0	0	0	0	...	190656	191171	191583	192089	192626	193171	193674	194186	194671	195162
3	NaN	Andorra	42.50630	1.521800	0	0	0	0	0	0	...	14988	14988	14988	15002	15003	15014	15016	15025	15025	15025
4	NaN	Angola	-11.20270	17.873900	0	0	0	0	0	0	...	45583	45817	45945	46076	46340	46539	46726	46929	47079	47168

5 rows × 590 columns

We can see that each row of the table defines the number of infected individuals for each country and/or province, and columns correspond to dates. Similar tables can be loaded for other data, such as number of recovered and number of deaths.

In [127]:

recovered = pd.read_csv(recovered_dataset_url)
deaths = pd.read_csv(deaths_dataset_url)

Making Sense of the Data

From the table above the role of province column is not clear. Let's see the different values that are present in Province/State column:

In [128]:

infected['Province/State'].value_counts()

Out [128]:

Australian Capital Territory    1
Xinjiang                        1
Martinique                      1
Guadeloupe                      1
French Polynesia                1
                               ..
Fujian                          1
Chongqing                       1
Beijing                         1
Anhui                           1
Turks and Caicos Islands        1
Name: Province/State, Length: 87, dtype: int64

From the names we can deduce that countries like Australia and China have more detailed breakdown by provinces. Let's look for information on China to see the example:

In [129]:

infected[infected['Country/Region']=='China']

Out [129]:

	Province/State	Country/Region	Lat	Long	1/22/20	1/23/20	1/24/20	1/25/20	1/26/20	1/27/20	...	8/20/21	8/21/21	8/22/21	8/23/21	8/24/21	8/25/21	8/26/21	8/27/21	8/28/21	8/29/21
58	Anhui	China	31.8257	117.2264	1	9	15	39	60	70	...	1008	1008	1008	1008	1008	1008	1008	1008	1008	1008
59	Beijing	China	40.1824	116.4142	14	22	36	41	68	80	...	1112	1113	1115	1115	1115	1115	1115	1115	1115	1115
60	Chongqing	China	30.0572	107.8740	6	9	27	57	75	110	...	603	603	603	603	603	603	603	603	603	603
61	Fujian	China	26.0789	117.9874	1	5	10	18	35	59	...	780	780	780	782	782	783	783	784	785	786
62	Gansu	China	35.7518	104.2861	0	2	2	4	7	14	...	199	199	199	199	199	199	199	199	199	199
63	Guangdong	China	23.3417	113.4244	26	32	53	78	111	151	...	3001	3007	3012	3020	3023	3032	3040	3043	3046	3055
64	Guangxi	China	23.8298	108.7881	2	5	23	23	36	46	...	289	289	289	289	289	289	289	289	289	289
65	Guizhou	China	26.8154	106.8748	1	3	3	4	5	7	...	147	147	147	147	147	147	147	147	147	147
66	Hainan	China	19.1959	109.7453	4	5	8	19	22	33	...	190	190	190	190	190	190	190	190	190	190
67	Hebei	China	39.5490	116.1306	1	1	2	8	13	18	...	1317	1317	1317	1317	1317	1317	1317	1317	1317	1317
68	Heilongjiang	China	47.8620	127.7615	0	2	4	9	15	21	...	1614	1614	1614	1614	1614	1614	1615	1615	1615	1615
69	Henan	China	37.8957	114.9042	5	5	9	32	83	128	...	1521	1522	1523	1524	1525	1525	1527	1528	1528	1528
70	Hong Kong	China	22.3000	114.2000	0	2	2	5	8	8	...	12049	12052	12057	12062	12069	12074	12077	12094	12100	12107
71	Hubei	China	30.9756	112.2707	444	444	549	761	1058	1423	...	68287	68289	68289	68289	68289	68289	68290	68290	68290	68290
72	Hunan	China	27.6104	111.7088	4	9	24	43	69	100	...	1181	1181	1181	1181	1181	1181	1181	1181	1181	1181
73	Inner Mongolia	China	44.0935	113.9448	0	0	1	7	7	11	...	412	412	412	412	412	412	412	412	412	412
74	Jiangsu	China	32.9711	119.4550	1	5	9	18	33	47	...	1583	1584	1584	1584	1586	1586	1587	1587	1589	1589
75	Jiangxi	China	27.6140	115.7221	2	7	18	18	36	72	...	937	937	937	937	937	937	937	937	937	937
76	Jilin	China	43.6661	126.1923	0	1	3	4	4	6	...	574	574	574	574	574	574	574	574	574	574
77	Liaoning	China	41.2956	122.6085	2	3	4	17	21	27	...	443	443	443	443	443	444	445	446	446	446
78	Macau	China	22.1667	113.5500	1	2	2	2	5	6	...	63	63	63	63	63	63	63	63	63	63
79	Ningxia	China	37.2692	106.1655	1	1	2	3	4	7	...	77	77	77	77	77	77	77	77	77	77
80	Qinghai	China	35.7452	95.9956	0	0	0	1	1	6	...	18	18	18	18	18	18	18	18	18	18
81	Shaanxi	China	35.1917	108.8701	0	3	5	15	22	35	...	668	668	668	669	669	669	669	669	669	669
82	Shandong	China	36.3427	118.1498	2	6	15	27	46	75	...	923	923	923	923	923	923	923	923	923	924
83	Shanghai	China	31.2020	121.4491	9	16	20	33	40	53	...	2420	2432	2436	2445	2451	2454	2462	2466	2471	2476
84	Shanxi	China	37.5777	112.2922	1	1	1	6	9	13	...	255	255	256	256	256	256	256	258	258	259
85	Sichuan	China	30.6171	102.7103	5	8	15	28	44	69	...	1181	1182	1183	1185	1185	1185	1186	1187	1188	1188
86	Tianjin	China	39.3054	117.3230	4	4	8	10	14	23	...	458	459	462	463	464	465	466	466	470	472
87	Tibet	China	31.6927	88.0924	0	0	0	0	0	0	...	1	1	1	1	1	1	1	1	1	1
88	Unknown	China	NaN	NaN	0	0	0	0	0	0	...	0	0	0	0	0	0	0	0	0	0
89	Xinjiang	China	41.1129	85.2401	0	2	2	3	4	5	...	980	980	980	980	980	980	980	980	980	980
90	Yunnan	China	24.9740	101.4870	1	2	5	11	16	26	...	1000	1007	1010	1014	1021	1031	1039	1047	1064	1067
91	Zhejiang	China	29.1832	120.0934	10	27	43	62	104	128	...	1420	1420	1421	1428	1428	1429	1429	1429	1429	1430

34 rows × 590 columns

Pre-processing the Data

We are not interested in breaking countries down to further territories, thus we would first get rid of this breakdown and add information on all territories together, to get info for the whole country. This can be done using groupby:

In [130]:

infected = infected.groupby('Country/Region').sum()
recovered = recovered.groupby('Country/Region').sum()
deaths = deaths.groupby('Country/Region').sum()

infected.head()

Out [130]:

	Lat	Long	1/22/20	1/23/20	1/24/20	1/25/20	1/26/20	1/27/20	1/28/20	1/29/20	...	8/20/21	8/21/21	8/22/21	8/23/21	8/24/21	8/25/21	8/26/21	8/27/21	8/28/21	8/29/21
Country/Region
Afghanistan	33.93911	67.709953	0	0	0	0	0	0	0	0	...	152448	152448	152448	152583	152660	152722	152822	152960	152960	152960
Albania	41.15330	20.168300	0	0	0	0	0	0	0	0	...	138132	138790	139324	139721	140521	141365	142253	143174	144079	144847
Algeria	28.03390	1.659600	0	0	0	0	0	0	0	0	...	190656	191171	191583	192089	192626	193171	193674	194186	194671	195162
Andorra	42.50630	1.521800	0	0	0	0	0	0	0	0	...	14988	14988	14988	15002	15003	15014	15016	15025	15025	15025
Angola	-11.20270	17.873900	0	0	0	0	0	0	0	0	...	45583	45817	45945	46076	46340	46539	46726	46929	47079	47168

5 rows × 588 columns

You can see that due to using groupby all DataFrames are now indexed by Country/Region. We can thus access the data for a specific country by using .loc:|

In [131]:

infected.loc['US'][3:].plot()
recovered.loc['US'][3:].plot()
plt.show()

Note

how we use [3:] to remove the first three elements of a sequence that contain non-date metadata (the Province/State and geolocation columns). We can also drop those three columns altogether:

In [132]:

infected.drop(columns=['Lat','Long','Province/State'],inplace=True)
recovered.drop(columns=['Lat','Long','Province/State'],inplace=True)
deaths.drop(columns=['Lat','Long','Province/State'],inplace=True)

Investigating the Data

Let's now switch to investigating a specific country. Let's create a frame that contains the data on infections indexed by date:

In [133]:

def mkframe(country):
    df = pd.DataFrame({ 'infected' : infected.loc[country],
                        'recovered' : recovered.loc[country],
                        'deaths' : deaths.loc[country]})
    df.index = pd.to_datetime(df.index)
    return df

df = mkframe('US')
df

Out [133]:

	infected	recovered	deaths
2020-01-22	1	0	0
2020-01-23	1	0	0
2020-01-24	2	0	0
2020-01-25	2	0	0
2020-01-26	5	0	0
...	...	...	...
2021-08-25	38223029	0	632272
2021-08-26	38384360	0	633564
2021-08-27	38707294	0	636720
2021-08-28	38760363	0	637254
2021-08-29	38796746	0	637531

586 rows × 3 columns

In [134]:

df.plot()
plt.show()

Now let's compute the number of new infected people each day. This will allow us to see the speed at which pandemic progresses. The easiest way to do it is to use diff:

In [135]:

df['ninfected'] = df['infected'].diff()
df['ninfected'].plot()
plt.show()

We can see high fluctuations in data. Let's look closer at one of the months:

In [136]:

df[(df.index.year==2020) & (df.index.month==7)]['ninfected'].plot()
plt.show()

It clearly looks like there are weekly fluctuations in data. Because we want to be able to see the trends, it makes sense to smooth out the curve by computing running average (i.e. for each day we will compute the average value of the previous several days):

In [137]:

df['ninfav'] = df['ninfected'].rolling(window=7).mean()
df['ninfav'].plot()
plt.show()

In order to be able to compare several countries, we might want to take the country's population into account, and compare the percentage of infected individuals with respect to country's population. In order to get country's population, let's load the dataset of countries:

In [138]:

countries = pd.read_csv(countries_dataset_url)
countries

Out [138]:

	UID	iso2	iso3	code3	FIPS	Admin2	Province_State	Country_Region	Lat	Long_	Combined_Key	Population
0	4	AF	AFG	4.0	NaN	NaN	NaN	Afghanistan	33.939110	67.709953	Afghanistan	38928341.0
1	8	AL	ALB	8.0	NaN	NaN	NaN	Albania	41.153300	20.168300	Albania	2877800.0
2	12	DZ	DZA	12.0	NaN	NaN	NaN	Algeria	28.033900	1.659600	Algeria	43851043.0
3	20	AD	AND	20.0	NaN	NaN	NaN	Andorra	42.506300	1.521800	Andorra	77265.0
4	24	AO	AGO	24.0	NaN	NaN	NaN	Angola	-11.202700	17.873900	Angola	32866268.0
...	...	...	...	...	...	...	...	...	...	...	...	...
4191	84056037	US	USA	840.0	56037.0	Sweetwater	Wyoming	US	41.659439	-108.882788	Sweetwater, Wyoming, US	42343.0
4192	84056039	US	USA	840.0	56039.0	Teton	Wyoming	US	43.935225	-110.589080	Teton, Wyoming, US	23464.0
4193	84056041	US	USA	840.0	56041.0	Uinta	Wyoming	US	41.287818	-110.547578	Uinta, Wyoming, US	20226.0
4194	84056043	US	USA	840.0	56043.0	Washakie	Wyoming	US	43.904516	-107.680187	Washakie, Wyoming, US	7805.0
4195	84056045	US	USA	840.0	56045.0	Weston	Wyoming	US	43.839612	-104.567488	Weston, Wyoming, US	6927.0

4196 rows × 12 columns

Because this dataset contains information on both countries and provinces, to get the population of the whole country we need to be a little bit clever:

In [139]:

countries[(countries['Country_Region']=='US') & countries['Province_State'].isna()]

Out [139]:

	UID	iso2	iso3	code3	FIPS	Admin2	Province_State	Country_Region	Lat	Long_	Combined_Key	Population
790	840	US	USA	840.0	NaN	NaN	NaN	US	40.0	-100.0	US	329466283.0

In [140]:

pop = countries[(countries['Country_Region']=='US') & countries['Province_State'].isna()]['Population'].iloc[0]
df['pinfected'] = df['infected']*100 / pop
df['pinfected'].plot(figsize=(10,3))
plt.show()

Computing R_t

To see how infectious is the disease, we look at the basic reproduction number R_0, which indicated the number of people that an infected person would further infect. When R_0 is more than 1, the epidemic is likely to spread.

R_0 is a property of the disease itself, and does not take into account some protective measures that people may take to slow down the pandemic. During the pandemic progression, we can estimate the reproduction number R_t at any given time t. It has been shown that this number can be roughly estimated by taking a window of 8 days, and computing R_t=\frac{I_{t-7}+I_{t-6}+I_{t-5}+I_{t-4}}{I_{t-3}+I_{t-2}+I_{t-1}+I_t} where I_t is the number of newly infected individuals on day t.

Let's compute R_t for our pandemic data. To do this, we will take a rolling window of 8 ninfected values, and apply the function to compute the ratio above:

In [141]:

df['Rt'] = df['ninfected'].rolling(8).apply(lambda x: x[4:].sum()/x[:4].sum())
df['Rt'].plot()
plt.show()

You can see that there are some gaps in the graph. Those can be caused by either NaN, if inf values being present in the dataset. inf may be caused by division by 0, and NaN can indicate missing data, or no data available to compute the result (like in the very beginning of our frame, where rolling window of width 8 is not yet available). To make the graph nicer, we need to fill those values using replace and fillna function.

Let's further look at the beginning of the pandemic. We will also limit the y-axis values to show only values below 6, in order to see better, and draw horizontal line at 1.

In [142]:

ax = df[df.index<"2020-05-01"]['Rt'].replace(np.inf,np.nan).fillna(method='pad').plot(figsize=(10,3))
ax.set_ylim([0,6])
ax.axhline(1,linestyle='--',color='red')
plt.show()

Another interesting indicator of the pandemic is the derivative, or daily difference in new cases. It allows us to see clearly when pandemic is increasing or declining.

In [143]:

df['ninfected'].diff().plot()
plt.show()

Given the fact that there are a lot of fluctuations in data caused by reporting, it makes sense to smooth the curve by running rolling average to get the overall picture. Let's again focus on the first months of the pandemic:

In [150]:

ax=df[df.index<"2020-06-01"]['ninfected'].diff().rolling(7).mean().plot()
ax.axhline(0,linestyle='-.',color='red')
plt.show()

Challenge

Now it is time for you to play more with the code and data! Here are a few suggestions you can experiment with:

• See the spread of the pandemic in different countries.
• Plot R_t graphs for several countries on one plot for comparison, or make several plots side-by-side
• See how the number of deaths and recoveries correlate with number of infected cases.
• Try to find out how long a typical disease lasts by visually correlating infection rate and deaths rate and looking for some anomalies. You may need to look at different countries to find that out.
• Calculate the fatality rate and how it changes over time. You may want to take into account the length of the disease in days to shift one time series before doing calculations

References

You may look at further studies of COVID epidemic spread in the following publications:

• Sliding SIR Model for Rt Estimation during COVID Pandemic, blog post by Dmitry Soshnikov
• T.Petrova, D.Soshnikov, A.Grunin. Estimation of Time-Dependent Reproduction Number for Global COVID-19 Outbreak. Preprints 2020, 2020060289 (doi: 10.20944/preprints202006.0289.v1)
• Code for the above paper on GitHub

In [ ]: