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
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.
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]$
- la media $\mu$ è pari al rapporto tra $\alpha$ e $\beta$ (ovvero tra $shape$ e $rate$)
- la deviazione standard $\sigma$ è pari al rapporto tra $\sqrt{\alpha}$ e $\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
$k_t = 300$ e $k_{t-1}=100$, ne ricaviamo $R[t] = 300 / 100 = 3$
$k_t = 100$ e $k_{t-1}=100$, ne ricaviamo $R[t] = 100 / 100 = 1$
$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
- $k_t = 1500$ e $k_{t-1}=1000$
- $k_t = 150$ e $k_{t-1}=100$
- $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).