{ "cells": [ { "cell_type": "markdown", "source": [ "# Calcola l'indice di riproduzione Rt e gli intervalli di confidenza degli ultimi 30 giorni\r\n", "\r\n", "*Sulla base del tasso di crescita λ e su una stima del tempo di generazione g si calcolano l'indice Rt e le fasce di confidenza al 68% degli ultimi 30 giorni. Gli esiti vengono riportati a video nella forma* \r\n", "\r\n", " ------------------------------\r\n", " data Rt 68% [ min, max ]\r\n", " ------------------------------" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "Carica i moduli necessari" ], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [ "import json\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`)." ], "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(\" :> \")" ], "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", " 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": [ "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": [ "La variabile `date1` definisce la data cui associare i valori giornalieri di Rt. Qui, convenzionalmente, associamo il valore dell'indice alla data associata ai file iniziali in modo da reneder più immediato il confronto con i valori presenti nel foglio `nuoviPositivi-IT.xlsx`. Per un confronto con [CovidStat](https://covid19.infn.it/) è sufficiente sottrarre 4 giorni." ], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [ "date1 = [datetime.date(2020, 3, 8) + datetime.timedelta(days=int(intervalloGiorni/2)) + datetime.timedelta(days=i)\r\n", " for i in range(numGiorni-intervalloGiorni-6)]" ], "outputs": [], "metadata": {} }, { "cell_type": "markdown", "source": [ "Si riportano a video i risultati degli ultimi 30 giorni con un arrotondamento a due cifre decimali sia per Rt che per i corrispondenti estremi degli intervalli di confidenza al 68%. " ], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [ "print()\r\n", "print(\" data Rt int.confidenza\")\r\n", "print(\"-----------------------------------------\")\r\n", "for i in range(30):\r\n", " if media14Giorni[-30:][i] >= 1:\r\n", " print(date1[-30:][i].strftime(\"%d/%m/%y\"), ' ', '{:.2f}'.format(media14Giorni[-30:][i]), \" 68% [\",\r\n", " '{:.2f}'.format(media14GiorniInf[-30:][i]), \",\", '{:.2f}'.format(media14GiorniSup[-30:][i]), \"]\")\r\n", " else:\r\n", " print(date1[-30:][i].strftime(\"%d/%m/%y\"), ' ', '{:.2f}'.format(media14Giorni[-30:][i]), \" 68% [\",\r\n", " '{:.2f}'.format(media14GiorniSup[-30:][i]), \",\", '{:.2f}'.format(media14GiorniInf[-30:][i]), \"]\")\r\n", "print()" ], "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" }, "kernelspec": { "name": "python3", "display_name": "Python 3.7.4 64-bit" }, "interpreter": { "hash": "0600588c3b5f4418cbe7b5ebc6825b479f3bc010269d8b60d75058cdd010adfe" } }, "nbformat": 4, "nbformat_minor": 2 }