indiceRiproduzioneRt.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "source": [
                "# Calcola l'indice di riproduzione <em>R<sub>t</sub></em> e ne rappresenta l'andamento\r\n",
                "\r\n",
                "*Sulla base del tasso di crescita &lambda; (notebook `tassoCrescitaLambda.ipynb`) e su una stima del tempo di generazione si calcola l'indice <em>R</em><sub>t</sub>. 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 <em>R<sub>t</sub></em>."
            ],
            "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`).<br>\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 &lambda; (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 <a href=\"#ref01\">[1]</a> e <a href=\"#ref02\">[2]</a> ($g=5.8$ giorni), si calcola il valore giornaliero di <em>R<sub>t</sub></em> 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 ( <a href=\"#ref01\">[1]</a> e <a href=\"#ref02\">[2]</a>), si ricalcolano i valori di <em>R<sub>t</sub></em> 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 <em>R<sub>t</sub></em>. 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",
                "<p id=\"ref01\">\r\n",
                "      [1] G. Bonifazi, L. Lista et al., <a href=\"https://doi.org/10.1140/epjp/s13360-021-01339-6\">A simplified estimate of the Effective Reproduction <em>R</em><sub>t</sub> Number using its relation with the doubling time and application to Italian COVID-19 data</a>, The European Physical Journal Plus (2021)\r\n",
                "    </p>\r\n",
                "\r\n",
                "<p id=\"ref02\">\r\n",
                "      [2] D. Cereda et al., <a href=\"https://arxiv.org/abs/2003.09320\">The early phase of the COVID-19 outbreak in Lombardy, Italy</a>, arXiv:2003.09320 (2020)\r\n",
                "</p>"
            ],
            "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
}