{ "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 }