<jupyter_start><jupyter_text>Sliding SIR Modelling for Different CountriesThis notebook contains the code to estimate effective reproduction index Rt for different countries, based on COVID-19 Data Repository by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University, which is available [here][CountryData].[CountryData]: https://github.com/CSSEGISandData/COVID-19  Initial setupLet's start by installing all required libraries...<jupyter_code>import sys
!{sys.executable} -m pip install --user --quiet pandas>=1.0.3 matplotlib scipy numpy

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import pandas as pd<jupyter_output><empty_output><jupyter_text>Loading data for all countriesData is loaded directly from [GitHub Repo][CountryData]. Since the repository is updated regularly with new actual data, we can obtain up-to-data estimations each time we run the code. For your convenience, the data files that we have used in this study are also stored in the `data` directory, so they can be loaded from there, should the current data become unavailable.<jupyter_code>infected_dataset_url = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"
recovered_dataset_url = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv"
deaths_dataset_url = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv"
countries_dataset_url = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/UID_ISO_FIPS_LookUp_Table.csv"
countries = pd.read_csv(countries_dataset_url)
infected_original = pd.read_csv(infected_dataset_url)
recovered_original = pd.read_csv(recovered_dataset_url)
deaths_original = pd.read_csv(deaths_dataset_url)
countries.head()<jupyter_output><empty_output><jupyter_text>Data on infected/recovered cases is located in three separate tables, and we also need the information about country population, which is located in a separate table. Since some countries are further broken down by territories (eg. US - by State/City), we also add those territories together (using `sum`) to get the overall numbers for the country:<jupyter_code>population = countries[countries['Province_State'].isnull()][['Country_Region','Population']].rename(columns={'Country_Region' : 'Country/Region'}).set_index('Country/Region')
infected = infected_original.groupby('Country/Region').sum().reset_index().set_index('Country/Region').join(population,on='Country/Region')
deaths = deaths_original.groupby('Country/Region').sum().reset_index().set_index('Country/Region').join(population,on='Country/Region')
recovered = recovered_original.groupby('Country/Region').sum().reset_index().set_index('Country/Region').join(population,on='Country/Region')
infected.head()<jupyter_output><empty_output><jupyter_text>Computing Parameters using Sliding SIR modelFirst, we define the function to solve SIR model with given starting values and parameters<jupyter_code># The SIR model differential equations.
def deriv(y, t, N, beta, gamma):
    S, I, R = y
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return dSdt, dIdt, dRdt

# Compute SIR model starting from given numbers of infected/removed ppl
def sir_model(infected,removed,N,beta,gamma,ndays):
    t = np.linspace(0,ndays,ndays)
    y0 = N-infected-removed,infected,removed
    ret = odeint(deriv, y0, t, args=(N, beta, gamma))
    return ret.T # S,I,R<jupyter_output><empty_output><jupyter_text>Function `model` computes the model given real values of newly infected people in `V` and removed cases in `R`. It computes SIR model for the given number of days, and then computes and returns loss function between real and computed data.`fit` then uses SciKit Learn optimization to find the best parameter $\beta$, given fixed $\gamma=$`the_gamma`<jupyter_code># Модель по заданному вектору заболевших V и удалённых R
def model(V,R,N,beta,gamma):
    S,I,R = sir_model(V[0],R[0],N,beta,gamma,len(V))
    dV = np.diff(V)
    dI = np.diff(I+R)
    return np.linalg.norm(dV-dI)

# Fix the gamma parameter
the_gamma = 1/30

# Подобрать параметры модели по векторам V и R
def fit(V,R,N):
    # res = minimize(lambda x:model(V,R,N,x[0],x[1]),x0=[0.5,1/20],method='powell')
    # return res.x[0],res.x[1]
    res = minimize(lambda x:model(V,R,N,x,the_gamma),x0=0.5,method='powell')
    return res.x,the_gamma<jupyter_output><empty_output><jupyter_text>Performing Analysis of a CountryHere we load the data for specified country into separate dataframe (so-called **country dataframe**). `make_frame` re-formats original country data into dataframe and performs running window averaging to smooth curves.<jupyter_code>def make_frame(country_name,smooth_window=3):
    f = pd.DataFrame([infected.loc[country_name],recovered.loc[country_name],deaths.loc[country_name]]).T
    population = f.iloc[-1,0]
    f = f.iloc[2:-1].reset_index()
    f.columns = ['Date','Infected','Recovered','Deaths']
    f['Removed'] = f['Recovered']+f['Deaths']
    f["Date"] = pd.to_datetime(f["Date"],format="%m/%d/%y")
    for x in ['Infected','Recovered','Deaths','Removed']:
        f[x+"_Av"] = f[x].rolling(window=smooth_window).mean()
    return population, f

def get_start_index(df):
    return df[df['Infected_Av']>1000].index[0]<jupyter_output><empty_output><jupyter_text>`compute_params` is the main method implementing Sliding SIR algorithm: it runs through all days from starting date and performs 7-day fitting of the SIR model using `fit` function above<jupyter_code>def compute_params(df,population, start_index, ndays=8):
    for i in range(start_index,len(df)-ndays):
        V = df['Infected_Av'][i:i+ndays].to_numpy()
        R = df['Removed_Av'][i:i+ndays].to_numpy()
        beta,gamma = fit(V,R,population)
        df.loc[i,'Beta'] = beta
        df.loc[i,'Gamma'] = gamma<jupyter_output><empty_output><jupyter_text>The method below puts everything together, it takes the country name, loads the data, and performs sliding SIR fitting and Rt calculations. It returns enriched country frame and population.<jupyter_code>def analyze(country_name,truncate_frame=True):
    population, df = make_frame(country_name)
    n = get_start_index(df)
    compute_params(df,population,n)
    df['Rt'] = df['Beta'] / df['Gamma']
    return population, df.iloc[n:] if truncate_frame else df

pop, df = analyze('Russia')

df.head()<jupyter_output><empty_output><jupyter_text>Below, we display the data for Russia. A couple of specific dates are highlighted as well.<jupyter_code>import matplotlib.dates as mdates
ax = df.plot('Date','Rt',linewidth=3,color='red',grid=True,ylim=(0,6),figsize=(12,8))
ax.axvline('2020-04-02', color='k', linestyle='--',label='Low quarantine measures')
ax.axvline('2020-04-15', color='k', linestyle='-.',label='Strict quarantine measures')
ax.xaxis.grid(True, which='minor')
ax.xaxis.set_major_locator(mdates.DayLocator(interval=7))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%d-%m'))
plt.gcf().autofmt_xdate()
plt.legend(fontsize=16)
def plot(population,df,ax=None):
    n = 0 # get_start_index(df)
    start_date = df.iloc[n,0]
    if ax is None:
        print("Population = {}, Infection start date (>100) = {}".format(population,start_date.date()))
    df['Days'] = df['Date'].apply(lambda x : (x-start_date).days)
    df['PI'] = df['Infected'] / population * 100
    df['Delta_Infected_Gr'] = (df['Infected'].diff()).clip(lower=0) / population * 1000000
    ax = ax or plt.gca()
    ax = df.plot('Days','Rt',ax=ax,legend=False,linewidth=3,color='red')
    ax.axhline(y = 1,linestyle='--', color='red')
    ax.set_ylim(0,8)
    ax.set_ylabel('Rt')
    df.plot(x='Days',y='Delta_Infected_Gr',secondary_y=True,kind='bar',ax=ax,legend=False)#,xticks=df['Days'][::10])
    ax.set_xlabel('Days',fontsize=14)
    ax.set_xticks(df['Days'][::10])
    ax.set_xticklabels(df['Days'][::10])
    ax.right_ax.set_ylabel('Infected (/day/1m)',fontsize=14)

plot(pop,df)<jupyter_output><empty_output><jupyter_text>Comparing Several CountriesWe define several contries as a grid, because it is easier to study them group by group. This grid is also used for plotting.We run computations for each country, and store its **country frame** into a dictionary.<jupyter_code>countries = [
    ['Spain','Italy', 'France', 'Germany'],
    ['Brazil', 'United Kingdom', 'US', 'Russia'],
    ['Sweden', 'Norway', 'Finland', 'Denmark'],        
    ['China', 'Japan', 'Korea, South', 'India']]

# Compute all countries
country_data = {}
for x in sum(countries,[]):
    print("Computing data for {}".format(x))
    country_data[x] = analyze(x)<jupyter_output>Computing data for Spain<jupyter_text>Producing a table of $R_0$ values for countries. It is as simple as taking the $R_t$ value at the index 0<jupyter_code>r0s = pd.DataFrame([ (c,k[1].iloc[0]['Rt']) for c,k in country_data.items() ],columns=['Country','R0'])
r0s<jupyter_output><empty_output><jupyter_text>Now let's produce multi-country plot.<jupyter_code>mx = max([len(x) for x in countries])
my = len(countries)
plt.rc('xtick',labelsize=14)
plt.rc('ytick',labelsize=14)
plt.rc('font',size=14)
fig,ax = plt.subplots(my,mx,figsize=(20,15),sharex=False)
for i in range(my):
    for j in range(mx):
        c = countries[i][j]
        pop,df = country_data[c]
        plot(pop,df,ax=ax[i,j])
        ax[i,j].set_xlim(0,90)
        ax[i,j].set_title(c,fontsize=20)
        ax[i,j].axes.xaxis.set_visible(True)
fig.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Comparing $R_t$ GraphsNow we want to plot the $R_t$ graphs on one plot. To do that, we convert the dictionary to one dataframe, with days as rows, and country names as columns:<jupyter_code>rtdf = [
    cd[1][['Days','Rt']].set_index('Days').rename(columns={'Rt':c }) for c,cd in country_data.items()]
from functools import reduce
rtdf = reduce(lambda df,f : df.join(f), rtdf)
rtdf
rtdf.plot(figsize=(13,9))<jupyter_output><empty_output><jupyter_text>It is impractival to see all graphs on one plot, so we plot only a few countries of interest.<jupyter_code>def plot_countries(l):
    rtdf = [cd[1][['Days','Rt']].set_index('Days').rename(columns={'Rt':c }) for c,cd in country_data.items() if c in l]
    rtdf = reduce(lambda df,f : df.join(f), rtdf)
    rtdf.plot(figsize=(10,7))

plot_countries(['Russia','US','Italy','China','Sweden'])<jupyter_output><empty_output><jupyter_text>Adding Apple Mobility IndexNow we want to see how people's mobility affects the epidemic spread. To do so, we will use Apple's mobility index. We use the data from Apple's site directly. For your convenience, the file that we used in our study is also provided in the `data` directory.<jupyter_code>apple_mobility_dataset_url = "https://covid19-static.cdn-apple.com/covid19-mobility-data/2010HotfixDev18/v3/en-us/applemobilitytrends-2020-06-14.csv"
apple_mobility = pd.read_csv(apple_mobility_dataset_url)
apple_mobility.head()<jupyter_output><empty_output><jupyter_text>In Apple mobility index, there are several different indices for each country: driving, walking, etc. There is also information about sub-regions. So in the following code we average different types of mobility, and group / average by regions, to get the average value of Apple Mobility index per country.<jupyter_code>apple_mobility['cou'] = np.where(np.logical_or(apple_mobility['geo_type']=='city',apple_mobility['geo_type']=='sub-region'),apple_mobility['country'],apple_mobility['region'])
apple = apple_mobility[apple_mobility['geo_type']!='county'].groupby('cou').mean().reset_index().set_index('cou')
apple.head()<jupyter_output><empty_output><jupyter_text>Because in Apple Mobility Index and in the previous dataset country names may differ, we introduce the mapping table. We only map those countries that we are going to analyze in the code below, if you find any other discrepancies - feel free to add them to the table.<jupyter_code>country_mapping = { 'US' : 'United States', 'China' : None, 'Korea, South' : 'Republic of Korea' }<jupyter_output><empty_output><jupyter_text>The following function adds Apple Mobility Index as another column into country dataframe. We also average out the Apple Mobility Index with 7-day sliding windows to get rid of weekly fluctuations.<jupyter_code>def add_apple_index(df,country_name):
    country_name = country_mapping.get(country_name,country_name)
    if country_name is None:
        df['AppleIndex'] = 0
        df['AppleIndexAv'] = 0
        return df
    start_date = df.iloc[0]['Date']
    tdf = apple.loc[country_name].T.reset_index().rename(columns={'index':'Date', country_name : 'AppleIndex'})
    tdf["Date"] = pd.to_datetime(tdf["Date"],format="%Y-%m-%d")
    tdf["AppleIndexAv"] = tdf["AppleIndex"].rolling(window=7).mean()
    tdf = tdf[tdf['Date']>=start_date].set_index('Date')
    return df.join(tdf,on='Date')<jupyter_output><empty_output><jupyter_text>For each country we have analyzed above, we now add Apple Mobility Index to the dataframe. The resulting dictionary of country names and enriched dataframes is called `country_datax`<jupyter_code>country_datax = { k : (d[0],add_apple_index(d[1],k)) for k,d in country_data.items() }
country_datax['Russia'][1].head()<jupyter_output><empty_output><jupyter_text>Now we can draw the graph of Rt and Mobility Index on one plot:<jupyter_code>mx = max([len(x) for x in countries])
my = len(countries)
plt.rc('xtick',labelsize=16)
plt.rc('ytick',labelsize=16)
plt.rc('font',size=16)
fig,ax = plt.subplots(my,mx,figsize=(20,15),sharex=False)
for i in range(my):
    for j in range(mx):
        c = countries[i][j]
        pop,df = country_datax[c]
        axe=df.plot('Days','Rt',ax=ax[i,j],linewidth=3,)
        axe.set_xlim(0, 90)
        axe.set_ylim(0, 10)
        axe.set_title(c,fontsize=24)
        ax3 = df.plot('Days','AppleIndexAv',secondary_y=True,ax=axe,legend=False,xticks=df['Days'][::10], linewidth=3,)
        ax[i,j].right_ax.set_ylim(0, 150)
fig.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>To see the relation more clearly, we compute the derivative $dR_t/dt$, and plot it together with Apple Mobility Index. Derivative is computed as simple first-order difference between daily values, and then averaged out using weekly average.<jupyter_code>mx = max([len(x) for x in countries])
my = len(countries)
plt.rc('xtick',labelsize=16)
plt.rc('ytick',labelsize=16)
plt.rc('font',size=16)
fig,ax = plt.subplots(my,mx,figsize=(20,15),sharex=False)
for i in range(my):
    for j in range(mx):
        c = countries[i][j]
        pop,df = country_datax[c]
        df["dRt/dt"] = df["Rt"].diff().rolling(window=7).mean()
        axe=df.plot('Days','dRt/dt',ax=ax[i,j],linewidth=3,)
        axe.set_xlim(0, 90)
        axe.set_ylim(-0.6, 0.1)
        axe.set_title(c,fontsize=24)
        ax3 = df.plot('Days','AppleIndexAv',secondary_y=True,ax=axe,legend=False,xticks=df['Days'][::10], linewidth=3,)
        ax[i,j].right_ax.set_ylim(0, 150)
        ax[i,j].set_title(c)
fig.tight_layout()
plt.show()
mx = max([len(x) for x in countries])
my = len(countries)
fig,ax = plt.subplots(my,mx,figsize=(20,15))
fig.tight_layout()
for i in range(my):
    for j in range(mx):
        c = countries[i][j]
        pop,df = country_datax[c]
        df["RtP"] = df["Rt"].diff()
        df["RtP_7"] = (df["RtP"]).rolling(window=7).mean()
        axe2=df.plot('Days','RtP_7',ax=ax[i,j])
        axe2.set_ylim(-0.6, 0.15)
        ax3 = df.plot('Days','AppleIndexAv',secondary_y=True,ax=axe2,legend=False,xticks=df['Days'][::10])
        ax[i,j].right_ax.set_ylim(0, 150)
        ax[i,j].set_title(c)

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14260651629072682, 'fraction_numerical': 0.016228070175438595, 'mean_word_length': 5.58457095709571, 'pattern_counts': {'":': 0, '<': 84, '<?xml version=': 0, '>': 88, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '515761', 'license': 'MIT', 'n_tokens_mistral': 5803, 'n_tokens_neox': 5514, 'n_words': 1268, 'path': 'notebooks/SlidingSIR.ipynb', 'repo_name': 'AzureAdvocateBit/SlidingSIR-1', 'size': '15960'}