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]$ è 1.431 mentre il minimo dell'intervallo di credibilità al 95% è 1.384 e il valore massimo 1.479.

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 1.953307 0.461315
1 Basilicata 1.043011 0.560379
2 Calabria 1.427702 0.246669
3 Campania 1.467443 0.147484
4 Emilia-Romagna 1.454793 0.166408
5 Friuli Venezia Giulia 1.673953 0.292189
6 Lazio 1.576819 0.154276
7 Liguria 1.129330 0.270237
8 Lombardia 1.235738 0.127724
9 Marche 1.655860 0.340031
10 Molise 2.600000 1.907878
11 P.A. Bolzano 1.547112 0.362862
12 P.A. Trento 1.549133 0.500726
13 Piemonte 1.264659 0.167508
14 Puglia 1.867692 0.283645
15 Sardegna 0.861635 0.389531
16 Sicilia 1.490428 0.153290
17 Toscana 1.235669 0.156480
18 Umbria 2.382075 0.560904
19 Valle d'Aosta 0.642857 0.654654
20 Veneto 1.333333 0.132578

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
2021-10-06 0.854197 0.033908 0.820288 0.888105
2021-10-07 0.877759 0.035111 0.842647 0.912870
2021-10-08 0.880025 0.035498 0.844526 0.915523
2021-10-09 0.918536 0.036801 0.881735 0.955337
2021-10-10 0.961858 0.038363 0.923495 1.000221
2021-10-11 0.971380 0.038654 0.932726 1.010034
2021-10-12 0.981022 0.038815 0.942207 1.019838
2021-10-13 1.058923 0.040849 1.018074 1.099772
2021-10-14 1.139608 0.042702 1.096906 1.182310
2021-10-15 1.225709 0.044659 1.181050 1.270368
2021-10-16 1.262226 0.045013 1.217213 1.307238
2021-10-17 1.323997 0.045892 1.278105 1.369890
2021-10-18 1.370978 0.046593 1.324385 1.417570
2021-10-19 1.431287 0.047335 1.383952 1.478622