Calcolo del numero di riproduzione di base Rt in Italia
G. Amati
18 September, 2020
EpiModel
Modelli di diffusione e calcolo di \(R_t\)
“Il numero di riproduzione netto \(R_t\)indica il numero medio di infezioni secondarie generate da una persona infetta ad una certa data ed è una grandezza fondamentale per capire l’andamento dell’epidemia, Se \(R_t\)ha un valore inferiore alla soglia critica di 1 il numero di nuove infezioni tenderà a decrescere tanto più velocemente quanto più è lontano dall’unità, Per contro, quanto più \(R_t\)supera 1 tanto più rapidamente aumenterà il numero dei contagi. Pertanto, un valore di Rt sopra la soglia, anche in presenza di un numero basso di casi, è un campanello di allarme sull’andamento epidemico. Il numero di riproduzione netto in un certo contesto geografico può essere stimato a partire dalla conoscenza della data di inizio sintomi dei casi, una volta nota la distribuzione dell’intervallo seriale (ovvero la distanza temporale fra la comparsa dei sintomi in una persona infettata e nei casi da essa generati) [1]”
Per calcolare \({R_t}\), ovvero il numero medio di persone infettate da parte di un singolo infetto al tempo t, occorre prima definire i modelli di diffusione della malattia. Esistono due possibili modelli per la simulazione della diffusione del virus. I modelli a compartimenti secondo i quali la popolazione viene distribuita in compartimenti separati, Deterministic compartmental models (DCM), e quelli con approccio micro, Stochastic individual contact models (ICM). I primi sono estremamente facili da comprendere ma soprattutto non richiedono risorse computazionali notevoli. I secondi invece sono basati su simulazione Monte Carlo e diventano proibitivi se il numero di nodi diventa dell’ordine del milione. I modelli con compartimenti descrive tre tipi di processi: il processo dei contatti, il processo relativo all diffusione dell’infezione (cioè relativo alla persona e tra le persone), e il processo relativo all’evoluzione demografica (partenze e arrivi, morti inclusi). I gruppi in cui viene suddivisa la popolazione sono in prima analisi quelli relativi all’evoluzione della malattia. Esistono i suscettibili \(S\), cioè le persone sane alle quali è possibile trasmettere la malattia, gli infetti \(I\)cioè quelli che hanno già contratto la malattia. Gli infetti si possono ulteriormente dividere in altri gruppi, ad esempio gli asintomatici \(A\)e quelli con sintomi \(Y\), poi i guariti \(R\)o i morti \(D\). Si possono estendere questi modelli anche a ulteriori compartimenti come i maschi $M $e femmine \(F\),gli anziani e i giovani, le città o le provincie. La relazione tra questi gruppi sono espresse da equazioni differenziali. Se ad esempio consideriamo il modello più semplice secondo il quale abbiamo solo i due gruppi dei suscettibili e degli infetti, allora possiamo considerare il modello \(S \rightarrow I\)seguente:
\[\displaystyle \frac{\partial I}{\partial t} = - \lambda\cdot S \]
\[\displaystyle \frac{\partial S}{\partial t} = - \lambda\cdot I \]
\[\displaystyle\frac{\partial S}{\partial t} + \frac{\partial I}{\partial t} = 0 \]
\[\lambda = \beta\cdot c \cdot\frac{I}{N} \]
\(\beta\)la probabilità di trasmissione per contatto
\(c\) il tasso di contatti per persona nell’unità di tempo
\(N\)la popolazione al tempo \(t\)
Il parametro \(c\)è anche chiamato act rate e denotato (act.rate nel modello EpiModel [2]).
All’inizio dell’infezione l’andamento dell’infezione è esponenziale (la popolazione è in pratica l’insieme dei suscettibili, cioè \(N\sim S\)). Infatti abbiamo:
\[\displaystyle \frac {\partial I} {\partial t}=-\frac {\partial S} {\partial t}\]
e quindi
\[\displaystyle\frac{\partial I}{\partial t}=\lambda\cdot I\]
cioè:
\[\displaystyle I\sim I_0\cdot e^{\lambda t} \]
dove
\[\displaystyle\lambda = \frac{R_0-1} {T_g}\]
\(T_g\)è anche chiamato il serial interval cioè il tempo nel quale il virus prima di morire in un ospite fa in tempo a infettare altre \(R_0-1\)persone quindi \(R_0-1\)è il numero dei nuovi infettati che si debbono aggiungere alla rete degli infettati [3]. Più semplicemente il serial interval è il tempo che intercorre tra il sorgere dei sintomi in un infetto (primario) e il sorgere dei sintomi in uno che lui ha infettato (secondario).
Il serial interval mediano è stato stimato intono ai 4.6 giorni (95% CrI: 3.5, 5.9) [4] e a 3.96 giorni [5]. Sono stime molto minori rispetto alle prime stime fatte che collocavano il serial interval intorno ai 7 giorni [6]. In effetti nel momento iniziale di sviluppo dell’epidemia, anche in Lombardia è stato calcolato un serial interval simile, cioè di 6.6 giorni.
Otteniamo dunque una relazione interessante che dal fitting della curva esponenziale di crescita ci fornisce il nostro \(R_0\)correlato al tempo \(DT\)di duplicazione della popolazione infetta:
\[\displaystyle I\sim I_0\cdot 2^{\frac{t} {DT}} \]
Questo è il modello logistico che è molto semplice e descrive bene l’evoluzione dell’epidemia al suo inizio. Ovviamente va sempre modificato in quanto è sempre \(R_0(t+\Delta t) < R(t)\)e dunque il fitting va fatto sui dati degli ultimi giorni.
\[\displaystyle DT\sim \frac{\ln(2)\cdot T_g}{R_0-1} \]
Il tempo di duplicazione dunque si può ottenere come \(DT=\frac{1}{a}\)dal coefficiente di regressione \(a\):
\[\log_2(I) \sim a\cdot t + b \]
dal quale otteniamo:
\[R_0 \sim 1+ \frac{ a\cdot\ln(2)}{T_g} \]
Nel caso che si prenda il modello SIR con popolazione variabile allora l’insieme degli infetti è \(I= N-(R+D)\)invece della cumulativa degli infettati \(I\)come nel modello più semplice SI, allora in fase di regressione del contagio \(a\)sarà negativo e quindi \(R_0<1\)(mentre nel modello SI è sempre >=1).
Il calcolo di \(R_t\)
Utilizzando il coefficiente di regressione \(a\)per l’equazione
\[\log_2(I) \sim a\cdot t + b \]
calcolato su un periodo di \(n\)giorni, che precedono il nostro istante \(t\)otteniamo una funzione \(a(t)\), il cui valore inverso, come abbiamo visto, fornisce il tempo di duplicazione \(DT\).
Il valore \(R_t\)è dunque dato da:
\[\displaystyle R_t \sim 1+ \frac{ a(t)\cdot\ln(2)}{T_g} \]
Per il Lazio, ad esempio, otteniamo il seguente grafico:
In blu i valori di \(R_t\)del Lazio. La curva in rosso è ottenuta mediante un filtro di Kalman sui valori del Lazio. L’ultimo valore è relativo al 3 settembre.
Possiamo mettere in unico grafico tutte le regioni (attenzione alla diversa percezione delle curve dovute a una variazione delle proporzioni di scala tra gli assi!): 
L’impressione è che le curve tra le varie regioni si somiglino. Cioè, anche se con maggiori o minori intensità, gli andamenti siano temporalmente simili. In effetti possiamo mettere in grafico anche l’\(R_t\)dell’Italia e notare un andamento simile:
La curva è ottenuta mediante un filtro di Kalman sui valori aggregati italiani. L’ultimo valore è relativo al 3 settembre 2020.
Il metodo bayesiano per il calcolo di Rt
Qui presentiamo due metodi di tipo bayesiano per calcolare i valori di \(R_t\)[7] e quello usato in [8] basato sul software EpiModel in R [9] per il calcolo di \(R_t\)utilizzato in Italia.
Il primo metodo è spiegato molto bene in [10], ed ha il codice disponibile in Python.
Dato il valore \(\lambda\), che ci siamo già calcolati precedentemente \(\lambda = \frac{R_0-1} {T_g}\), possiamo calcolarci la probabilità di osservare \(k_t= I_t-I_{t-1}\)nuovi casi mediante la distribuzione di Poisson (likelihood)
\[\displaystyle P(k_t|\lambda)= \frac{\lambda^{k_t} e^{ \lambda}} {I_t!}=P(k_t|R_t)\]
come prior prendiamo \(P(R_{t-1}|k_{t-1})\)e quindi possiamo calcolare il valore della distribuzione a posteriori. Applicando iterativamente Bayes su \(P(R_{t-1}|k_{t-1})\)otteniamo
\[\displaystyle P(R_{t}|k_{t})\propto P(R_{t_0})\cdot \prod_{x=t_0}^{x=t}\frac{\lambda^{k_x} e^{ \lambda}} {k_x!}\]
Ad esempio nel Lazio abbiamo: 
Mettendo tutte le regioni:
I valori in grigio sono delle altre regioni. La curva è ottenuta mediante un filtro di Kalman sui valori aggregati del Lazio (in blu). L’ultimo valore è relativo al 3 settembre 2020.
L’andamento in Italia è:
La curva Rt italiana è ottenuta con il metodo bayesiano e mediante un filtro di Kalman sui valori aggregati da tutte le regioni. L’ultimo valore è relativo al 3 settembre 2020.
Il secondo metodo Bayesiano utilizza la distribuzione Gamma per definire il serial interval (con valore medio del serial interval = 6.6 e deviazione standard =4.9) ottenuti dai dati (non disponibili in formato aperto!) della regione Lombardia. Come prior si utilizza ancora una funzione Gamma con due nuovi parametri. La likelihood è come nell’altro metodo ancora una Poisson.
Metodo Bayesiano di Cori et al. Valori calcolati con il software EpiModel in R (valori in grigio). Smoothing dei valori mediante un filtro di Kalman (linea rossa)
Conclusioni
Tutti e tre i metodi evidenziano un aumento del valore di \(R_t\)negli ultimi tre mesi, con \(R_0 >1\)(crescita esponenziale) nell’ultimo mese e mezzo.Il metodo Bayesiano è più sensibile alle oscillazioni, ma questo dipende anche dal periodo di osservazione per effettuare la stima.
![In nero il metodo con regressione. Blu il metodo Bayesiano [9], in rosso il metodo Bayesiano 10
Il 4 settembre \(R_t\)secondo il metodo Bayesiano è stato 1.06 e secondo il metodo di regressione è 1.27. La stima ufficiale è stata 1.18. (https://www.repubblica.it/cronaca/2020/09/04/news/coronavirus_indice_rt_a_1_18_aumentano_i_sintomatici-266271299/)
Oggi
$R_t=$1.11 (Bayesiano)
$R_t=$1.1881238 (regressione)
Bibliografia
[1] Istituto Superiore Sanità,https://www.epicentro.iss.it/coronavirus/bollettino/Bollettino-sorveglianza-integrata-COVID-19_11-agosto-2020.pdf
[2]Jenness S, Goodreau SM, Morris M. EpiModel: Mathematical Modeling of Infectious Disease. R package version 1.6.5. 2018
[3] Anderson R, May R. Infectious Diseases of Humans: Dynamics and Control. Oxford University Press; 1992.
[4] Serial interval of novel coronavirus (COVID-19) infections. Nishiura H, Linton NM, Akhmetzhanov AR. Int J Infect Dis. 2020 Mar 4;93:284-286.
[5] Serial Interval of COVID-19, Stephen G. Baum, MD reviewing Du Z et al. Emerg Infect Dis 19-Mar-2020
[6] Li Q, Guan X, Wu P, Wang X, Zhou L, Tong Y, et al. Early transmission dynamics in Wuhan, China, of novel coronavirus-infected pneumonia. N Engl J Med. 2020 Jan 29
[7] Bettencourt LMA, Ribeiro RM (2008) Real Time Bayesian Estimation of the Epidemic Potential of Emerging Infectious Diseases. PLoS ONE 3(5)
[8] COVID-19 working group,Epidemiological characteristics of COVID-19 cases in Italy and estimates of the reproductive numbers one month into the epidemic, doi:https://doi.org/10.1101/2020.04.08.20056861
[9] A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics, Anne Cori, Neil M. Ferguson, Christophe Fraser, Simon Cauchemez,American Journal of Epidemiology, Volume 178, Issue 9, 1 November 2013, Pages 1505–1512