Visualizes the data on a provincial level, focusing on the province of bolzano and models the infections of the next days for the province.
Disclaimer: Experts estimate that many of the infections remain unreported, and therefore the collected data does not necessarily give the full picture of the situation.
Covid-19 Outbreak Analysis Italian Provinces
import pandas as pd
from pandas.plotting import register_matplotlib_converters
import numpy as np
import seaborn as sns
sns.set(style="whitegrid")
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from datetime import datetime
%matplotlib inline
register_matplotlib_converters()
now = datetime.now()
print("Updated at (UTC Time):", now)
Updated at (UTC Time): 2020-03-19 19:31:06.408698
Data of italian provinces
Read data and preprocess
url = 'https://raw.githubusercontent.com/pcm-dpc/COVID-19/0c05e59ef9996edc20bd0d0b8ae65d4956dda383/dati-province/dpc-covid19-ita-province.csv'
df = pd.read_csv(url)
df = df[df['denominazione_provincia'] != 'In fase di definizione/aggiornamento']
print(df.size)
df['date'] = pd.to_datetime(df['data'])
26750
Data of Province Bolzano
df_bz = df[df["denominazione_provincia"] == 'Bolzano']
Plot the average cumulative cases per province in Italy
plt.figure(figsize=(10,5))
plot =sns.lineplot(df['date'], df['totale_casi'])
x_dates = df['date'].dt.strftime('%d-%m').sort_values().unique()
plot.set_xticklabels(x_dates, rotation=45)
plot.set_title("Average Cumulative Cases in Italy per Province")
plt.show()
Plot the total cumulative sum of cases in Italy
total_cases_sum = df.groupby('date')['totale_casi'].sum()
plt.figure(figsize=(10,5))
dates = np.array(total_cases_sum.index.to_pydatetime(), dtype=np.datetime64)
plot = sns.lineplot(dates, total_cases_sum)
x_dates = total_cases_sum.index.strftime('%d-%m')
plot.set_xticklabels(x_dates, rotation=45)
plot.set_title("Cumulative Cases in Italy")
plt.show()
Province Bolzano
plt.figure(figsize=(10,5))
plot = sns.lineplot(df_bz['date'], df_bz['totale_casi'])
plot =sns.lineplot(df['date'], df['totale_casi'])
x_dates = df['date'].dt.strftime('%d-%m').sort_values().unique()
plot.set_xticklabels(x_dates, rotation=45)
plot.set_title("Cumulative Cases in Province BZ")
plt.show()
New cases every day in Province Bolzano
cases_of_day = []
cases_before = 0
for index, row in df_bz.iterrows():
cases_of_day.append(row[9] - cases_before)
cases_before = row[9]
df_cases = df_bz.copy()
df_cases['cases_day'] = cases_of_day
fig, ax = plt.subplots(figsize=(10,6))
plot = sns.barplot(df_cases['date'], df_cases['cases_day'], palette="Blues")
x_dates = df_cases['date'].dt.strftime('%d-%m').unique()
plot.set_xticklabels(x_dates, rotation=45)
plot.set_title("New Cases Every Day in Province BZ")
Text(0.5, 1.0, 'New Cases Every Day in Province BZ')
Predicting cumulative cases in Province Bolzano
We’re going to use linear regression to predict the number of deaths for the next week. The growth is exponential, therefore we create a new column with the log of the infections, which we use to predict further cases.
import statsmodels.api as sm
# we will predict 2 days into the future
DAY_RANGE = range(1,df_bz.shape[0] + 2)
# create new DF with total cases and augment with log_infections and days since first case
df2 = pd.DataFrame()
days = [i for i in range(df_bz.shape[0])]
df2['date'] = df_bz['date']
df2['days'] = days
df2['infections'] = df_bz['totale_casi']
df2 = df2[df2.infections > 0]
df2['log_infections'] = np.log(df2['infections'])
df2.head()
| date | days | infections | log_infections | |
|---|---|---|---|---|
| 136 | 2020-02-25 18:00:00 | 1 | 1 | 0.0 |
| 264 | 2020-02-26 18:00:00 | 2 | 1 | 0.0 |
| 392 | 2020-02-27 18:00:00 | 3 | 1 | 0.0 |
| 520 | 2020-02-28 18:00:00 | 4 | 1 | 0.0 |
| 648 | 2020-02-29 17:00:00 | 5 | 1 | 0.0 |
Create a model using the days sice the first case and the log of the total infected count.
X = df2.days
X = sm.add_constant(X)
y = df2.log_infections
/home/kredde/anaconda3/lib/python3.7/site-packages/numpy/core/fromnumeric.py:2495: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.
return ptp(axis=axis, out=out, **kwargs)
model = sm.OLS(y, X)
res = model.fit()
We use the params of the fitted model to predict the growth:
a_exp = res.params.const
b_exp = res.params.days
y_pred = [np.exp(a_exp + b_exp * x) for x in DAY_RANGE]
Here we will do the same thing, but instead of fitting the model on the whole data, we will just use the past 4 days to get a recent trend of how the epidemic is evolving
df3 = df2[df2.days > df_bz.shape[0] - 5]
X = df3.days
X = sm.add_constant(X)
y = df3.log_infections
model = sm.OLS(y, X)
res = model.fit()
a = res.params.const
b = res.params.days
x_pred_last_four_days =[np.exp(a + b * x) for x in DAY_RANGE]
Same thing again, but now with just 2 days
df4 = df2[df2.days > df_bz.shape[0] - 5]
X = df4.days
X = sm.add_constant(X)
y = df4.log_infections
model = sm.OLS(y, X)
res = model.fit()
a = res.params.const
b = res.params.days
x_pred_last_two_days =[np.exp(a + b * x) for x in DAY_RANGE]
plt.figure(figsize=(10,5))
plot = sns.lineplot(DAY_RANGE, y_pred, label="Sim. Infections")
plot = sns.lineplot(DAY_RANGE, x_pred_last_four_days, label="Sim. Infections (Based on last 4 days)")
plot = sns.lineplot(DAY_RANGE, x_pred_last_two_days, label="Sim. Infections (Based on Last 2 days)")
plot = sns.lineplot(days, df_bz['totale_casi'], label="Confirmed Cases")
plt.ylabel('Cases', fontsize=14)
plt.xlabel('Days since first case', fontsize=14)
plot.set_title("Cumulative Cases in Province BZ")
Text(0.5, 1.0, 'Cumulative Cases in Province BZ')
Fitting a logistic model
from scipy.optimize import curve_fit, fsolve
import datetime as dt
def logistic_model(x,a,b,c):
return c/(1+np.exp(-(x-b)/a))
df_log = pd.read_csv(url, error_bad_lines=False)
df_log = df_log[df_log['denominazione_provincia'] == 'Bolzano']
df_log = df_log.loc[:,['data','totale_casi']]
FMT = '%Y-%m-%d %H:%M:%S'
date = df_log['data']
df_log['day_of_year'] = date.map(lambda x : (datetime.strptime(x, FMT) - datetime.strptime("2020-01-01 00:00:00", FMT)).days )
df_log.head()
| data | totale_casi | day_of_year | |
|---|---|---|---|
| 8 | 2020-02-24 18:00:00 | 0 | 54 |
| 136 | 2020-02-25 18:00:00 | 1 | 55 |
| 264 | 2020-02-26 18:00:00 | 1 | 56 |
| 392 | 2020-02-27 18:00:00 | 1 | 57 |
| 520 | 2020-02-28 18:00:00 | 1 | 58 |
x = list(df_log.iloc[:,2])
y = list(df_log.iloc[:,1])
fit = curve_fit(logistic_model,x,y,p0=[2,100,20000])
The result gives us three values:
a: the infection speed
b: day of year with maximum new infections
c: total number of recorded infected people at the end of the epidemic
a = fit[0][0]
b = fit[0][1]
c = fit[0][2]
print("a=", a, " b=", b, " c=", c)
a= 2.8623687755101614 b= 76.5093282568078 c= 683.5176105627821
B will be the inflection point, i.e. the point where new infections on a day will be less than on the previous day
date = dt.date(2020,1,1) + dt.timedelta(int(b))
print(date)
2020-03-17
We can compute the end of the epidemic according to the model by solving the equation of the model
sol = int(fsolve(lambda x : logistic_model(x,a,b,c) - int(c),b))
date = dt.date(2020,1,1) + dt.timedelta(sol)
print(date)
2020-04-07
Plotting the logistic curve
fig, ax = plt.subplots(figsize=(10,6))
# logistic model
x_pred_log = [logistic_model(x,a,b,c) for x in range(54,114)]
y_pred_log = [dt.date(2020,1,1) + dt.timedelta(i) for i in range(54,114)]
# exponential model
x_pred_infections = [np.exp(a_exp + b_exp * x) for x in range(30)]
y_pred_exp = [dt.date(2020,1,1) + dt.timedelta(i) for i in range(54,84)]
plot = sns.lineplot(df_bz['date'], df_bz['totale_casi'], label="Confirmed Infections")
plot = sns.lineplot(y_pred_log, x_pred_log, label="Logistic Model")
plot = sns.lineplot(y_pred_exp, x_pred_infections, label="Exponential Model")
plt.xticks(rotation=70)
ax.xaxis.set_major_locator(mdates.AutoDateLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter('%m.%d'))
plt.ylabel('Total Infections', fontsize=14)
plt.xlabel('Date', fontsize=14)
plot.set_title("Comparison of total infection models")
Text(0.5, 1.0, 'Comparison of total infection models')
