COVID-19 Data Analysis

Rt semplice

Stima di Rt con... carta e penna.

Max Pierini


Carta e penna

Abbiamo urgente bisogno di stimare velocemente il numero di riproduzione di COVID-19 in Italia senza uso di computer e complesse simulazioni?

Bene. Armiamoci di carta, penna e una calcolatrice.

Da una fonte accertata (ad esempio repository GitHub del Diparimento Protezione Civile) otteniamo i valori dei nuovi positivi giornalieri delle ultime due settimane.

Calcoliamo le medie di ultima e penultima settimana, che chiameremo rispettivamente $\alpha$ e $\beta$.

Stimiamo il valore medio del numero di riproduzione $R[t]_\mu$ come rapporto tra $\alpha$ e $\beta$

E stimiamo un'approssimazione dell'intervallo di credibilità al 95% come

Per cui il valore medio di $R[t]$ è 0.730 mentre il minimo dell'intervallo di credibilità al 95% è 0.721 e il valore massimo 0.739.

Possiamo disegnare sul foglio la media, l'intervallo di credibilità e il valore di riferimento 1

Regioni

Allo stesso modo possiamo stimare $R[t]$ e l'intervallo di credibilità al 95% delle singole regioni italiane.

regione Rt CI
0 Abruzzo 0.670266 0.045708
1 Basilicata 0.696474 0.084195
2 Calabria 0.755767 0.050822
3 Campania 0.717665 0.026439
4 Emilia-Romagna 0.733616 0.031301
5 Friuli Venezia Giulia 0.680556 0.064815
6 Lazio 0.777012 0.030057
7 Liguria 0.757370 0.059263
8 Lombardia 0.682642 0.022484
9 Marche 0.715809 0.050147
10 Molise 0.636935 0.100979
11 P.A. Bolzano 0.711751 0.097766
12 P.A. Trento 0.643087 0.098233
13 Piemonte 0.727291 0.036799
14 Puglia 0.714645 0.033820
15 Sardegna 0.862414 0.053961
16 Sicilia 0.911485 0.040387
17 Toscana 0.744358 0.037927
18 Umbria 0.785994 0.067328
19 Valle d'Aosta 0.731818 0.215801
20 Veneto 0.637056 0.027028

Metodo

Tra le molte funzioni di distribuzione di probabilità, la distribuzione Gamma (vedi WikiPedia) in particolare è definita come

$$ f(x; \alpha, \beta) = \frac{\beta^\alpha x^{\alpha-1} e^{-\beta x}}{\Gamma(\alpha)} $$

dove $\Gamma$ è la funzione Gamma

$$ \Gamma(\alpha) = \int_{0}^{\infty} t^{\alpha - 1} e^{-t} dt $$

e per numeri naturali è riducibile a

$$ \Gamma(\alpha) = (\alpha - 1)! $$

Complessità matematica a parte, che non è scopo di questa trattazione, la distribuzione di probabilità Gamma è parametrizzata sui parametri $\alpha$ e $\beta$, che sono chiamati anche $shape$ e $rate$

$$ f(x; \alpha, \beta) = \mathbf{Gamma}(\alpha, \beta) $$

ed ha due proprietà interessanti che useremo per la stima di $R[t]$

  1. la media $\mu$ è pari al rapporto tra $\alpha$ e $\beta$ (ovvero tra $shape$ e $rate$)
$$ \mu = \frac{\alpha}{\beta} $$
  1. la deviazione standard $\sigma$ è pari al rapporto tra $\sqrt{\alpha}$ e $\beta$
$$ \sigma = \frac{\sqrt{\alpha}}{\beta} $$

Il numero di riproduzione effettivo $R[t]$ può essere molto semplicemente definito come il rapporto tra i nuovi positivi $k_{t}$ osservati al tempo $t$ e i nuovi positivi $k_{t-1}$ osservati al tempo precedente $t-1$ (vedi Wallinga-Teunis 2004 e Wallinga-Lipsitch 2007)

$$ R[t] = \frac{k_{t}}{k_{t-1}} $$

Dunque

  • se $k_t = k_{t-1}$ ovvero se i positivi osservati in $t$ e in $t-1$ sono uguali, il loro rapporto è 1 ed $R[t]$ è 1: l'epidemia è in una fase incerta

  • se $k_t > k_{t-1}$ ovvero se i positivi osservati in $t$ sono maggiori che in $t-1$, il loro rapporto è maggiore di 1 ed $R[t]$ è maggiore di 1: l'epidemia è in fase di crescita esponenziale

  • se $k_t < k_{t-1}$ ovvero se i positivi osservati in $t$ sono minori che in $t-1$, il loro rapporto è minore di 1 ed $R[t]$ è minore di 1: l'epidemia è in fase di risoluzione

Vediamo degli esempi

  1. $k_t = 300$ e $k_{t-1}=100$, ne ricaviamo $R[t] = 300 / 100 = 3$

  2. $k_t = 100$ e $k_{t-1}=100$, ne ricaviamo $R[t] = 100 / 100 = 1$

  3. $k_t = 50$ e $k_{t-1}=100$, ne ricaviamo $R[t] = 50 / 100 = 0.5$

Chiediamoci però, con quale livello di credibilità possiamo affermare che $R[t]$ è 3 nel primo esempio, 1 nel secondo esempio e 0.5 nel terzo?

Ovvero, quanto siamo sicuri dei valori stimati?

Per questo, possiamo usare le proprietà della distribuzione Gamma (vedi Cori et al. 2013).

Osserviamo per esempio cosa accade nei seguenti tre esempi

  1. $k_t = 1500$ e $k_{t-1}=1000$
  2. $k_t = 150$ e $k_{t-1}=100$
  3. $k_t = 15$ e $k_{t-1}=10$

È evidente che in tutti i tre casi $R[t]=1.5$ ma, sapendo che nella funzione Gamma la media $\mu=\alpha/\beta$, possiamo utilizzarla ponendo $\alpha=k_t$ e $\beta=k_{t-1}$

$$ R_t = \mathbf{Gamma}(\alpha = k_t, \beta = k_{t-1}) $$

Come vediamo, la media $\mu$ corrisponde infatti a 1.5 in tutti i tre casi, ma molto diversa è la distribuzione di probabilità: minori sono $\alpha$ e $\beta$ maggiore è la deviazione standard.

Infatti sappiamo che per la distribuzione Gamma

$$ \sigma = \frac{\sqrt{\alpha}}{\beta} $$

In particolare notiamo che per $\alpha=15$ e $\beta=10$, la distribuzione di probabilità comprende 1 che è il valore di riferimento per determinare se l'epidemia sia in fase di crescita o di decadimento: in questo caso dunque, l'intervallo di credibilità non esclude 1 ovvero non siamo "sufficientemente certi" che il valore reale di $R[t]$ sia maggiore di 1.

Come vediamo, al diminuire di $\alpha$ e $\beta$ la deviazione standard $\sigma$ aumenta ovvero aumenta il "livello di incertezza" del risultato medio.

Per approssimare l'intervallo di credibilità* al 95%, complesso da calcolare a mano, possiamo definire come limiti

$$ R_{t,\mathrm{lo}} = \mu - 2 \sigma $$$$ R_{t,\mathrm{hi}} = \mu + 2 \sigma $$

*: si parla di intervallo di credibilità e non di intervallo di confidenza perché la stima di $R[t]$ si basa su metodo bayesiano e non frequentista; il significato dal punto di vista pratico è simile

Ma come definiamo i tempi $t$ e $t-1$?

L'intervallo tra un contagio e il contagio secondario da esso indotto si chiama intervallo seriale (o intervallo di generazione) ed è stimato con complessi metodi a partire dall'osservazione dei contagi e dei loro contatti.

Il problema è che l'intervallo seriale non è fisso ma è anch'esso una distribuzione di probabilità (solitamente Gamma o Weibull) con media stimata di circa 7 giorni per COVID-19 (vedi Cereda-Tirani 2020).

Per semplificare dunque la stima di $R[t]$ per COVID-19 con i dati reali italiani, possiamo confrontare i dati medi settimanali (7 giorni, pari circa alla media dell'intervallo seriale).

Pertanto, in riferimento al giorno $t$, la media dell'ultima settimana risulterà

$$ k_{t} = \frac{\sum_{i=0}^{6}k_{t-i}}{7} $$

e la media della penultima settimana

$$ k_{t-1} = \frac{\sum_{i=7}^{13}k_{t-i}}{7} $$

che corrispondono ai parametri $\alpha$ e $\beta$ della distribuzione Gamma

$$ R[t] = \mathbf{Gamma}\left( \alpha=\frac{\sum_{i=0}^{6}k_{t-i}}{7} \;,\; \beta=\frac{\sum_{i=7}^{13}k_{t-i}}{7} \right) $$

Procedendo così a ritroso nel tempo, giorno per giorno, possiamo stimare $R[t]$ dell'intera serie temporale. Assegneremo il valore di Rt al giorno centrale tra le due finestre di 7 giorni (l'ultima stima verrà pertanto assegnata a 7 giorni fa).

Rt mean Rt C.I. Rt lo Rt hi
data
2022-05-02 0.807594 0.007938 0.799656 0.815532
2022-05-03 0.725056 0.007203 0.717854 0.732259
2022-05-04 0.796013 0.007975 0.788038 0.803989
2022-05-05 0.819498 0.008345 0.811153 0.827843
2022-05-06 0.841058 0.008652 0.832405 0.849710
2022-05-07 0.862687 0.008951 0.853736 0.871638
2022-05-08 0.879647 0.009192 0.870455 0.888839
2022-05-09 0.872899 0.009184 0.863715 0.882083
2022-05-10 0.850908 0.009164 0.841744 0.860071
2022-05-11 0.823070 0.009089 0.813981 0.832160
2022-05-12 0.816958 0.009204 0.807754 0.826162
2022-05-13 0.788603 0.009136 0.779468 0.797739
2022-05-14 0.755795 0.009020 0.746775 0.764815
2022-05-15 0.729804 0.008927 0.720877 0.738730