COVID-19 Data Analysis

Stima infetti reali

Stima infetti reali con fast backcasting.

Max Pierini


ITALIA

Possiamo definire il tasso di rilevamento $\mathrm{DR}_t$ (Detection Rate) al tempo $t$ come la percentuale di eventi $n(\hat{\alpha}_t)$ stimati sugli eventi $n(\alpha_t)$ osservati in $t$

$$ \mathrm{DR}_t = \frac{ \sum_{i=0}^{t} n(\hat{\alpha_i}) }{ \sum_{i=0}^{t} n(\alpha_i) } $$

Un $\mathrm{DR}_t$ pari al 100% significa che tutti i casi vengono correttamente rilevati. Al contrario un $\mathrm{DR}_t$ minore del 100% significa che una parte dei casi, ovvero $1-\mathrm{DR}_t$, sfuggono al rilevamento.

L'attuale $\mathrm{DR}_t$ stimato è quindi pari a 36.8% (95% CI: 20.4%, 56.8%), ovvero si stima che il 63.2% dei casi sarebbe sfuggito al rilevamento (95% CI: 43.2%, 79.6%).

Dal grafico dei nuovi casi e del tasso di rilevamento, notiamo che la maggior parte dei mancati rilevamenti è avvenuta nella prima fase della curva epidemica, indicativamente tra Febbraio e Giugno. In seguito, verosimilmente grazie alle campagne di screening e contact-tracing, la differenza tra casi reali (stimati) e osservati risulta molto meno marcata.

METODO

Concetto di backcasting

Il metodo di fast backcasting di EpiData.it si basa sul backcasting per simulazione Monte Carlo di Phipps et al. 2020 ma implementato per esecuzione veloce senza uso di simulazioni.

Siano dati due eventi $\alpha$ e $\omega$.

Sia definita Alpha-Omega-Rate $\mathrm{A}\Omega\mathrm{R}$ la probabilità condizionata non nulla del verificarsi dell'evento $\omega$ in seguito all'osservazione dell'evento $\alpha$

$$ P(\omega | \alpha) = \mathrm{A}\Omega\mathrm{R} \;\,\; >0 $$

Sia posta pari a 1 la probabilità condizionata che si sia verificato $\alpha$ data l'osservazione di un evento $\omega$

$$ P(\alpha | \omega) = 1 $$

ovvero

  • l'evento $\omega$ ha probabilità $\mathrm{A}\Omega\mathrm{R}$ di accadere in seguito ad un evento $\alpha$

  • l'osservazione di un evento $\omega$ comporta sicuramente l'accadimento di un evento $\alpha$ nel passato

Pertanto, data l'osservazione di $n(\omega)$ eventi $\omega$, si può stimare che il numero $n(\alpha)$ di eventi $\alpha$ occorsi sia

$$ n(\alpha) = \frac{n(\omega)}{\mathrm{A}\Omega\mathrm{R}} $$

Sia data una distribuzione di probabilità $p(\omega_{t_0 + \tau})$ nel tempo $t$ del verificarsi dell'evento $\omega_{t_0+\tau}$ dopo $\tau$ unità di tempo dall'osservazione dell'evento $\alpha_{t_0}$ al tempo $t_0$

$$ p(\omega_{t_0 + \tau}) = \mathbf{PDF}[\omega]_{t_0 + \tau} $$

Noto il numero $n(\omega_t)$ di eventi $\omega$ nel tempo $t$, possiamo stimare che il numero $n(\alpha_{t_0})$ di eventi $\alpha$ occorsi in $t_0<t$ siano pari alla somma degli eventi $n(\omega_{t_0 + \tau})$ moltiplicati per la densità di probabilità $\mathbf{PDF}[\omega]_{t_0 + \tau}$ al tempo $t_0 + \tau$ divisa per la probabilità condizionata $p(\omega|\alpha)$

$$ n(\alpha_{t_0}) = \sum_{\tau=1} \frac{ n(\omega_{t_0 + \tau}) }{ \mathrm{A}\Omega\mathrm{R} } \mathbf{PDF}[\omega]_{t_0 + \tau} $$

Indicando con $t=T$ l'ultimo giorno di osservazione, possiamo affermare che soprattutto negli ultimi giorni della serie temporale, il numero di eventi $n(\alpha)$ sarà ovviamente sottostimato perché i futuri eventi $\omega$, di cui potrebbero essere causa, non sono ancora accaduti.

Supponiamo, ad esempio, di avere i seguenti dati:

Vogliamo stimare il numero di eventi $n(\alpha)$ con metodo di backcasting dagli eventi $\omega$ noti e sapendo che la probabilità di osservare $\omega$ in conseguenza di $\alpha$ nel tempo $t$ è una distribuzione Gamma con i seguenti parametri di shape e scale (scale è definita come 1/rate) e il seguente $\mathrm{A} \Omega \mathrm{R}$

AΩR and Gamma parameters
AΩR shape scale
value 6.00% 6.223439 3.462716

Notiamo la sottostima nell'ultima parte della serie temporale dovuta agli eventi $\omega$ non ancora occorsi e la "caduta a zero" dei nuovi eventi $\alpha$ stimati in $t=T$: per l'ultimo giorno della serie temporale non si ha alcun dato futuro da cui stimare i probabili eventi $\alpha$.

Correzione temporale

Nota la distribuzione di probabilità del verificarsi dell'evento $\omega$ in seguito ad $\alpha$ nel tempo $t$ e indicato con $t=T$ l'ultimo giorno di osservazione, la probabilità $p(\alpha|\omega)_{t_0}$ di aver osservato tutti gli eventi $\omega$ causati dagli eventi $\alpha$ nel tempo $[t_0...T]$ è esprimibile come la densità di probabilità cumulativa di $\omega$ al tempo $t=T-t_0$

$$ p(\alpha|\omega)_{t_0} = \mathbf{CDF}[\omega]_{T-t_0} $$

possiamo perciò correggere il numero di eventi $n(\alpha_{t_0})$ stimati dividendo per la probabilità $p(\alpha|\omega)_{t_0}$

$$ n(\alpha_{t_0}) = \sum_{\tau=1}^{T} \frac{ n(\omega_{t_0+\tau}) }{ \mathrm{A}\Omega\mathrm{R} } \frac{ \mathbf{PDF}[\omega]_{t_0+\tau} }{ \mathbf{CDF}[\omega]_{T-t_0} } $$

Notiamo che nonostante la correzione abbia evitato la "caduta a zero" dei nuovi eventi $\alpha$ e la sottostima, ora siamo di fronte ad una sovrastima: la correzione effettuata non è in grado di "prevedere" i futuri eventi $\omega$ ma solamente di stimare la probabilità di aver osservato in $t_0$ tutti gli eventi $\omega$ causati da $\alpha$.

Aggiustamento

Possiamo perciò ottenere un'aggiustamento $n(\alpha_{t_0})_{\mathrm{adj}}$ della sovrastima dovuta alla correzione, prendendo la media degli stimati grezzi $n(\alpha_{t_0})_{\mathrm{raw}}$ e corretti $n(\alpha_{t_0})_{\mathrm{cor}}$

$$ n(\alpha_{t_0})_{\mathrm{raw}} = \sum_{\tau=1}^{T} \frac{ n(\omega_{t_0+\tau}) }{ \mathrm{A}\Omega\mathrm{R} } \frac{ \mathbf{PDF}[\omega]_{t_0+\tau} }{ \mathbf{CDF}[\omega]_{T-t_0} } $$$$ n(\alpha_{t_0})_{\mathrm{cor}} = \sum_{\tau=1}^{T} \frac{ n(\omega_{t_0+\tau}) }{ \mathrm{A}\Omega\mathrm{R} } \mathbf{PDF}[\omega]_{t_0+\tau} $$$$ n(\alpha_{t_0})_{\mathrm{adj}} = \frac{ n(\alpha_{t_0})_{\mathrm{raw}} + n(\alpha_{t_0})_{\mathrm{cor}} }{2} $$

COVID-19 in Italia

Se definiamo

  • $\alpha$: infezione da SARS-CoV-2
  • $\omega$: decesso per COVID-19
  • $\mathrm{A}\Omega\mathrm{R}$: Infection Fatality Rate $\mathrm{IFR}$ per COVID-19
  • $\mathbf{PDF}[\omega]$: distribuzione di probabilità del tempo dall'infezione all'eventuale decesso per COVID-19

data l'osservazione dei decessi giornalieri per COVID-19 possiamo stimare le infezioni occorse precedentemente.

I dati dei decessi e dei casi osservati dal repository GitHub COVID-19 del Dipartimento di Protezione Civile.

L'IFR stimato è tratto da Imperial College Report 34 Brazeau, Verity et al. 2020 (tabella pag 11, High Income Country).

I parametri per la distribuzione dell'intervallo dall'infezione al decesso sono tratti da Phipps et al. 2020 e Flaxman et al. 2020.

Sappiamo dunque che:

IFR values (Imperial College)
IFR
2.5% 0.78%
mean 1.29%
97.5% 1.79%
infection-to-onset & onset-to-death values
infection-to-onset onset-to-death
2.5% 4.10 12.80
mean 5.55 16.00
97.5% 7.00 19.20

Da questi parametri e conoscendo il coefficiente di variazione cov delle distribuzioni Gamma degli intervalli del tempo infection-to-onset $io$ (periodo d'incubazione) e del tempo onset-to-death $od$ (dai sintomi al decesso), possiamo facilmente calcolare i parametri di shape e scale della distribuzione Gamma dell'intervallo infection-to-death $id$ (dall'infezione al decesso).

Sappiamo infatti che

$$ \mathbf{cov} = \frac{\sigma}{\mu} $$$$ shape = \frac{\mu^2}{\sigma^2} = \frac{1}{\mathbf{cov}^2} $$$$ scale = \mu \cdot \mathbf{cov}^2 $$

Da Phipps et al. 2020 e Flaxman et al. 2020 sappiamo che i coefficienti di variazione stimati per COVID-19 dell'intervallo infection-to-onset (periodo d'incubazione) $\mathbf{cov}_{io} = 0.86$ e dell'intervallo onset-to-death (dai sintomi al decesso) $\mathbf{cov}_{od} = 0.45$ da cui, nota la media $\mu$, possiamo ricavare la deviazione standard $\sigma$ dei due intervalli

$$ \sigma_{io} = \mu_{io} \cdot \mathbf{cov}_{io} $$$$ \sigma_{od} = \mu_{od} \cdot \mathbf{cov}_{od} $$

e le distribuzioni di probabilità dei due intervalli.

Da Stewart et al. 2006 sappiamo che media e deviazione standard della convoluzione di due distribuzioni Gamma possono essere stimate come

$$ \mu_{1,2} = \mu_1 + \mu_2 $$$$ \sigma_{1,2} = \sqrt{\sigma_1^2 + \sigma_2^2} $$

Pertanto la distribuzione $\mathbf{Gamma}_{id}$ dell'intervallo infection-to-death data dalla convoluzione delle distribuzioni $\mathbf{Gamma}_{io}$ infection-to-onset (periodo d'incubazione) e $\mathbf{Gamma}_{od}$ onset-to-death, sarà

$$ \mu_{id} = \mu_{io} + \mu_{od} $$$$ \sigma_{id} = \sqrt{\sigma_{io}^2 + \sigma_{od}^2} $$

da cui possiamo ricavare il coefficiente di variazione $\mathbf{cov}_{id}$ dell'intervallo infection-to-death

$$ \mathbf{cov}_{id} = \frac{\sigma_{id}}{\mu_{id}} $$

e infine i parametri shape e scale dell'intervallo infection-to-death

$$ shape_{id} = \frac{1}{\mathbf{cov}_{id}^2} $$$$ scale_{id} = \mu_{id} \cdot \mathbf{cov}_{id}^2 $$

Pertanto, la distribuzione $\mathbf{PDF}[\omega]$ per l'intervallo infection-to-death risulta

$$ \mathbf{Gamma} \left( \alpha = \frac{1}{\mathbf{cov}_{id}^2} \;,\; \beta = \mu_{id} \cdot \mathbf{cov}_{id}^2 \right) $$

Con i valori sopra definiti, otteniamo i seguenti parametri per la distribuzione Gamma dell'intervallo infection-to-death

infection-to-death
Gamma parameters
shape scale
2.5% 6.261966 2.698833
mean 6.223439 3.462716
97.5% 6.190279 4.232443

Calcolate le tre distribuzioni per l'intervallo infection-to-death

  • CI 2.5%
  • media
  • CI 97.5%

le combineremo con i tre valori di IFR

  • CI 2.5% di IFR
  • media di IFR
  • CI 97.5% di IFR

in modo da ottenere una matrice di nove possibili stime

Backcasting parameters matrix
id 2.5% id mean id 97.5%
IFR 2.5% 1 4 7
IFR mean 2 5 8
IFR 97.% 3 6 9

da cui stimeremo mediana e intervallo di confidenza al 95%.

Dal risultato del backcasting notiamo la differenza tra la sottostima dei valori grezzi, la sovrastima della correzione temporale e l'aggiustamento con metodo EpiData.it.

Confronto con Monte Carlo

Confrontiamo dunque i risultati dell'algoritmo di backcasting veloce implementato da EpiData.it con le simulazioni Monte Carlo dell'implementazione di Phipps et al. 2020.

IFR values (Phipps et al.)
IFR
2.5% 0.37%
mean 0.76%
97.5% 1.15%

La distribuzione Gamma dell'intervallo da infezione a decesso è il medesimo utilizzato sopra.

L'IFR generico usato nel lavoro di Phipps et al è invece molto minore rispetto a quello usato da noi e stimato dall'Imperial College per l'Italia: differenti nazioni hanno verosimilmente differenti IFR dovuti soprattutto alla distribuzione delle fascie d'età nella popolazione. L'Italia è nazione piuttosto "anziana" ed ha pertanto un IFR verosimilmente più elevato rispetto a nazioni più "giovani".

Per questo motivo, non confronteremo i risultati delle due stime con i casi osservati ma solamente tra loro per verificare che il metodo di stima di EpiData.it dia risultati sovrapponibili alle simulazioni Monte Carlo di Phipps et al. (casi dal 1 Gennaio 2020 al 14 Settembre 2020).

Si noti che l'implementazione di Phipps et al. con metodo Monte Carlo comporta 10'000 simulazioni che, in un normale personal computer, prendono qualche ora di tempo per una serie temporale come quella di COVID-19.

L'algoritmo implementato qui e visualizzabile su GitHub (come tutto il sito EpiData.it) all'indirizzo Real_infected è stato eseguito in:

Eseguito in 0:00:00.178361