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.827 mentre il minimo dell'intervallo di credibilità al 95% è 0.811 e il valore massimo 0.843.

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.000000 0.165359
1 Basilicata 0.856716 0.154495
2 Calabria 0.682074 0.077230
3 Campania 0.828440 0.042903
4 Emilia-Romagna 0.880986 0.062234
5 Friuli Venezia Giulia 0.706144 0.132689
6 Lazio 0.837727 0.054919
7 Liguria 0.746237 0.105989
8 Lombardia 0.835409 0.041514
9 Marche 0.826995 0.111367
10 Molise 0.922222 0.378757
11 P.A. Bolzano 1.347826 0.310676
12 P.A. Trento 0.847170 0.211556
13 Piemonte 0.851031 0.060999
14 Puglia 0.790810 0.051142
15 Sardegna 0.627475 0.104269
16 Sicilia 0.876261 0.058637
17 Toscana 0.861879 0.063409
18 Umbria 0.857523 0.178806
19 Valle d'Aosta 0.706949 0.244545
20 Veneto 0.774836 0.059735

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, possiamo stimare $R[t]$ dell'intera serie temporale (nella prossima tabella e nel grafico, è stimato utilizzando le medie di settimane complete, da lunedì a domenica). Escludiamo la prima settimana perché non abbiamo una settimana precedente da confrontare e l'ultima settimana se incompleta.

alpha beta Rt mean Rt SD Rt lo Rt hi
data
2020-02-23 240.857143 NaN NaN NaN NaN NaN
2020-03-01 811.571429 240.857143 3.369514 0.118278 3.132958 3.606070
2020-03-08 2481.714286 811.571429 3.057912 0.061383 2.935146 3.180679
2020-03-15 4913.000000 2481.714286 1.979680 0.028244 1.923193 2.036167
2020-03-22 5507.285714 4913.000000 1.120962 0.015105 1.090752 1.151172
... ... ... ... ... ... ...
2021-04-04 14518.142857 19473.142857 0.745547 0.006188 0.733172 0.757922
2021-04-11 14340.714286 14518.142857 0.987779 0.008248 0.971282 1.004276
2021-04-18 13224.000000 14340.714286 0.922130 0.008019 0.906092 0.938167
2021-04-25 11730.285714 13224.000000 0.887045 0.008190 0.870665 0.903425
2021-05-02 5350.857143 11730.285714 0.456157 0.006236 0.443686 0.468629

63 rows × 6 columns