{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Calcola l'indice di riproduzione Rt e ne rappresenta l'andamento\r\n",
"\r\n",
"*Sulla base del tasso di crescita λ (notebook `tassoCrescitaLambda.ipynb`) e su una stima del tempo di generazione si calcola l'indice Rt. Data l'incertezza del tempo di generazione, si rappresentano le fasce di confidenza al 68%.*"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Carica i moduli necessari"
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"import json\r\n",
"import matplotlib.pyplot as plt\r\n",
"import numpy as np\r\n",
"from scipy.stats import linregress\r\n",
"import datetime"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"La funzione `partiziona()` partiziona `xDati` in un array di array. Ciascun elemento ha lunghezza `intervalloGiorni` e il loro numero è pari a `n-intervalloGiorni+1`. Segue la funzione per il calcolo di Rt."
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"def partiziona(xdati, intervalloGiorni):\r\n",
" numElem = len(xdati)\r\n",
" elementi = []\r\n",
" for i in range(intervalloGiorni, numElem + 1):\r\n",
" elementi = np.append(elementi, xdati[i-intervalloGiorni:i])\r\n",
" return np.reshape(elementi, (numElem - intervalloGiorni + 1, intervalloGiorni))\r\n",
"\r\n",
"def rT(g, l):\r\n",
" return np.exp(g*l)\r\n",
"######### Predisposizioni ###############################################################################################"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Elementi da inserire in input per avviare l'elaborazione dei dati nazionali o regionali. Questi devono essere stati preventivamente memorizzati nella medesima cartella e devono corrispondere alla data immessa (avviare eventualmente il notebook `prelevaRinomina.ipynb`).
\r\n",
"Si richiedono inoltre l'unità temporale da porre in ascissa e l'eventuale visualizzazione delle fasce di confidenza."
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"print()\r\n",
"print(\"Si intendono elaborare i nuovi positivi nazionali (i) o quelli di una regione (r)? \")\r\n",
"scelta = input(\" :> \")\r\n",
"print(\"Inserire la data nel formato (YYYYMMDD) \")\r\n",
"dataISO = input(\" :> \")\r\n",
"print(\"Si vuole l'asse temporale espresso in mesi (m) o in giorni trascorsi (g) dall'inizio pandemia? \")\r\n",
"unitaAssex = input(\" :> \")\r\n",
"print(\"Si vogliono visualizzare le fasce di errore (s/n)? \")\r\n",
"fasce = input(\" :> \")"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Vengono composti i nomi dei file e, una volta letti i dati contenuti, questi vengono caricati nell'array `nuoviPositivi`."
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"if scelta == \"i\":\r\n",
" regioneScelta = \"ITALIA\"\r\n",
" nomeFile = 'datiNazionali' + dataISO\r\n",
" # lettura del file csv: la colonna 8 corrisponde al campo nuovi_positivi\r\n",
" nuoviPositivi = np.genfromtxt(nomeFile + '.csv', delimiter=',', skip_header=1, usecols=8, dtype=int)\r\n",
"else:\r\n",
" print(\"Inserire la regione \")\r\n",
" regioneScelta = input(\" :> \")\r\n",
" regioneScelta = regioneScelta.lower().capitalize()\r\n",
" nomeFile = 'datiRegionali' + dataISO\r\n",
" # lettura del file json e riportato l'array ad un array di Numpy\r\n",
" with open(nomeFile + '.json') as f:\r\n",
" datiGrezzi = json.load(f)\r\n",
" nuoviPositivi = []\r\n",
" for record in datiGrezzi:\r\n",
" if record['denominazione_regione'] == regioneScelta:\r\n",
" nuoviPositivi.append(record['nuovi_positivi'])\r\n",
" nuoviPositivi = np.array(nuoviPositivi)"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Dovendo calcolare il logaritmo dei `nuoviPositivi` nel caso vi sia un valore giornaliero nullo (per cui il suo logaritmo non avrebbe significato), lo si pone pari ad 1. L'array `giorniDelFit` rappresenta i giorni nei quali si eseguirà il fit lineare."
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"nuoviPositivi[nuoviPositivi == 0] = 1\r\n",
"logPositivi = np.log(nuoviPositivi)\r\n",
"numGiorni = len(nuoviPositivi)\r\n",
"intervalloGiorni = 14\r\n",
"# i giorni da 0 a 13 nel quale viene eseguito il fit\r\n",
"giorniDelFit = np.arange(0, intervalloGiorni, 1)"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Per ognuno degli intervalli in cui è stato partizionato l'insieme dei `logPositivi`, viene eseguita una regressione lineare in un intervallo temporale di 14 giorni. Della retta ottimale si determina la sola pendenza in quanto questo è il termine associato, nel modello teorico, al parametro λ (o *tasso di crescita*)."
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"valoriLambdaLin = []\r\n",
"for i in range(intervalloGiorni, numGiorni+1):\r\n",
" # fit lineare\r\n",
" esito = linregress(giorniDelFit, logPositivi[i-intervalloGiorni:i])\r\n",
" valoriLambdaLin = np.append(valoriLambdaLin, esito.slope)"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Assegnato al tempo di generazione il valore ripreso dalle pubblicazion1 [1] e [2] ($g=5.8$ giorni), si calcola il valore giornaliero di Rt e quindi allo scopo di ridurre i fattori di disturbo, viene eseguita una media mobile su 14 giorni. "
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"generationTime = 5.8\r\n",
"valoriMediDayRtLin = rT(generationTime, valoriLambdaLin)\r\n",
"media14Giorni = np.mean(partiziona(valoriMediDayRtLin, intervalloGiorni), axis=1)"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Assegnata l'incertezza $\\Delta g=1.88$ giorni al tempo di generazione ( [1] e [2]), si ricalcolano i valori di Rt a partire dalle nuove stime $g-\\Delta g=5.88-1.88$ e $g+\\Delta g=5.88+1.88$."
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"if fasce == 's':\r\n",
" generationTimeError = 1.88\r\n",
" valoriMediDayRtLinSup = rT(generationTime+generationTimeError, valoriLambdaLin)\r\n",
" valoriMediDayRtLinInf = rT(generationTime-generationTimeError, valoriLambdaLin)\r\n",
" media14GiorniSup = np.mean(partiziona(valoriMediDayRtLinSup, intervalloGiorni), axis=1)\r\n",
" media14GiorniInf = np.mean(partiziona(valoriMediDayRtLinInf, intervalloGiorni), axis=1)"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Predisposizioni per l'asse x associato ai giorni trascorsi dal 24/02/2020 e valore massimo per l'asse y associato all'indice Rt. Codifica italiana della data inserita inizialmente."
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"xValoriMediGiornalieri = np.arange(intervalloGiorni-1, numGiorni)\r\n",
"xMediaMobile = np.arange(intervalloGiorni+5, numGiorni-intervalloGiorni/2)\r\n",
"codificaItaliana = dataISO[6:] + '-' + dataISO[4:6] + '-' + dataISO[:4]\r\n",
"rtMassimo = 3"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Parte grafica comune e parte dipendente dalle scelte iniziali. Le variabili temporali `date1`, `date2` come le precedenti `xValoriMediGiornalieri` e `xMediaMobile`, tutte associate all'asse x, controllano la traslazione temporale dei valori in ordinata e sono scelte in modo da ottenere un accordo visivo soddisfacente con i dati giornalieri e in coerenza con il sito [CovidStat](https://covid19.infn.it/)."
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"# parte grafica comune\r\n",
"plt.rcParams['figure.figsize'] = [12, 6]\r\n",
"fig, ax = plt.subplots()\r\n",
"plt.ylim([0.4, rtMassimo])\r\n",
"y_major_ticks = np.arange(0.4, rtMassimo, 0.1)\r\n",
"ax.set_yticks(y_major_ticks)\r\n",
"ax.set_ylabel('indice di riproduzione Rt')\r\n",
"ax.grid(which='both', color='.85', linestyle='-', linewidth=1)\r\n",
"# parte grafica variabile\r\n",
"if fasce == \"s\": # assieme all'Rt si visualizzano le fasce di confidenza\r\n",
" ax.set_title(regioneScelta + \": Indice di riproduzione Rt e intervalli di confidenza 68%\")\r\n",
" if unitaAssex == \"m\":\r\n",
" ax.set_xlabel('date dal 24 febbraio 2020')\r\n",
" date2 = [datetime.date(2020, 3, 8) + datetime.timedelta(days=int(intervalloGiorni/2)) + datetime.timedelta(days=i)\r\n",
" for i in range(numGiorni-intervalloGiorni-12)]\r\n",
" ax.fill_between(date2, media14Giorni, media14GiorniSup, interpolate=True, facecolor='papayawhip')\r\n",
" ax.fill_between(date2, media14Giorni, media14GiorniInf, interpolate=True, facecolor='papayawhip')\r\n",
" ax.plot(date2, np.full(xMediaMobile.shape, 1), color='b')\r\n",
" ax.plot(date2, media14GiorniSup, color='lightskyblue', linewidth=2, label='media mobile Rt ($g+\\Delta g$)')\r\n",
" ax.plot(date2, media14GiorniInf, color='greenyellow', linewidth=2, label='media mobile Rt ($g-\\Delta g$)')\r\n",
" ax.plot(date2, media14Giorni, color='r', linewidth=2, label='media mobile Rt su 14 giorni', zorder=4)\r\n",
" date1 = [datetime.date(2020, 3, 1) + datetime.timedelta(days=int(intervalloGiorni/2)) + datetime.timedelta(days=i)\r\n",
" for i in range(numGiorni-intervalloGiorni+1)]\r\n",
" ax.scatter(date1, valoriMediDayRtLin, s=4, label='Rt giornaliero', zorder=3)\r\n",
" ax.text(datetime.date(2020, 11, 1), 0.8*rtMassimo, 'aggiornato il\\n' + codificaItaliana)\r\n",
" else:\r\n",
" ax.set_xlabel('giorni dal 24 febbraio 2020')\r\n",
" ax.fill_between(xMediaMobile, media14Giorni, media14GiorniSup, interpolate=True, facecolor='papayawhip')\r\n",
" ax.fill_between(xMediaMobile, media14Giorni, media14GiorniInf, interpolate=True, facecolor='papayawhip')\r\n",
" ax.plot(xMediaMobile, np.full(xMediaMobile.shape, 1), color='b')\r\n",
" ax.plot(xMediaMobile, media14GiorniSup, color='lightskyblue', linewidth=2, label='media mobile Rt ($g+\\Delta g$)')\r\n",
" ax.plot(xMediaMobile, media14GiorniInf, color='greenyellow', linewidth=2, label='media mobile Rt ($g-\\Delta g$)')\r\n",
" ax.plot(xMediaMobile, media14Giorni, color='r', linewidth=2, label='media mobile Rt su 14 giorni', zorder=4)\r\n",
" ax.scatter(xValoriMediGiornalieri, valoriMediDayRtLin, s=4, label='Rt giornaliero', zorder=3)\r\n",
" ax.text(140, 0.8*rtMassimo, 'aggiornato il\\n' + codificaItaliana)\r\n",
" x_major_ticks = np.arange(0, numGiorni+20, 20)\r\n",
" ax.set_xticks(x_major_ticks)\r\n",
"else: # visualizza solo l'andamento di Rt\r\n",
" ax.set_title(regioneScelta + \": Indice di riproduzione Rt\")\r\n",
" if unitaAssex == \"m\":\r\n",
" ax.set_xlabel('date dal 24 febbraio 2020')\r\n",
" date2 = [datetime.date(2020, 3, 8) + datetime.timedelta(days=int(intervalloGiorni/2)) + datetime.timedelta(days=i)\r\n",
" for i in range(numGiorni-intervalloGiorni-12)]\r\n",
" ax.plot(date2, np.full(xMediaMobile.shape, 1), color='b')\r\n",
" ax.plot(date2, media14Giorni, color='r', linewidth=2, label='media mobile Rt su 14 giorni', zorder=4)\r\n",
" date1 = [datetime.date(2020, 3, 1) + datetime.timedelta(days=int(intervalloGiorni/2)) + datetime.timedelta(days=i)\r\n",
" for i in range(numGiorni-intervalloGiorni+1)]\r\n",
" ax.scatter(date1, valoriMediDayRtLin, s=4, label='Rt giornaliero', zorder=3)\r\n",
" ax.text(datetime.date(2020, 11, 1), 0.8*rtMassimo, 'aggiornato il\\n' + codificaItaliana)\r\n",
" else:\r\n",
" ax.set_xlabel('giorni dal 24 febbraio 2020')\r\n",
" ax.plot(xMediaMobile, np.full(xMediaMobile.shape, 1), color='b')\r\n",
" ax.plot(xMediaMobile, media14Giorni, color='r', linewidth=2, label='media mobile Rt su 14 giorni', zorder=4)\r\n",
" ax.scatter(xValoriMediGiornalieri, valoriMediDayRtLin, s=4, label='Rt giornaliero', zorder=3)\r\n",
" ax.text(140, 0.8*rtMassimo, 'aggiornato il\\n' + codificaItaliana)\r\n",
" x_major_ticks = np.arange(0, numGiorni+20, 20)\r\n",
" ax.set_xticks(x_major_ticks)\r\n",
"\r\n",
"plt.legend()\r\n",
"plt.show()"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Referenze\r\n",
"\r\n",
"
\r\n", " [1] G. Bonifazi, L. Lista et al., A simplified estimate of the Effective Reproduction Rt Number using its relation with the doubling time and application to Italian COVID-19 data, The European Physical Journal Plus (2021)\r\n", "
\r\n", "\r\n", "\r\n", " [2] D. Cereda et al., The early phase of the COVID-19 outbreak in Lombardy, Italy, arXiv:2003.09320 (2020)\r\n", "
" ], "metadata": {} } ], "metadata": { "orig_nbformat": 4, "language_info": { "name": "python", "version": "3.7.4", "mimetype": "text/x-python", "codemirror_mode": { "name": "ipython", "version": 3 }, "pygments_lexer": "ipython3", "nbconvert_exporter": "python", "file_extension": ".py" }, "kernelspec": { "name": "python3", "display_name": "Python 3.7.4 64-bit" }, "interpreter": { "hash": "0600588c3b5f4418cbe7b5ebc6825b479f3bc010269d8b60d75058cdd010adfe" } }, "nbformat": 4, "nbformat_minor": 2 }