fbpx Un nuovo metodo per la stima del numero di riproduzione R(t) | Scienza in rete
Covid-19/

Un nuovo metodo per la stima del numero di riproduzione R(t)

Tempo di lettura: 15 mins

La stima del numero di riproduzione netto R(t)1 di un'infezione, che rappresenta il numero di infetti generati in media da un singolo caso durante il proprio periodo infettivo, può essere fatta per mezzo di diverse assunzioni, ad esempio riguardo il tempo che intercorre tra contagio e diagnosi, e adottando diversi metodi. In particolare, la Fondazione Bruno Kessler, per conto dell’Istituto Superiore di Sanità, utilizza il metodo proposto dall’Imperial College2. In un articolo pubblicato su Scienza in rete Roberto Battiston ha mostrato come sia possibile usare un metodo più intuitivo e semplice per la stima di R(t), basato sul modello epidemiologico SIR3.

I metodi statistici come quello dell'Imperial College utilizzano come dato base il numero dei nuovi casi giornalieri registrati CR(t), e il tempo tra contagio e comparsa dei sintomi, mentre il metodo basato sul modello SIR fa uso degli infetti registrati IR(t). Quindi i problemi che riguardano le curve del numero dei nuovi casi registrati e del numero di persone infette si riflettono sulla stima del numero di riproduzione R(t) ottenuta con entrambi questi metodi.

In particolare, il numero di nuovi casi giornalieri registrati CR(t) ha un andamento che dipende dal numero dei tamponi fatti, quindi risulta in generale sottostimato rispetto al numero reale di nuovi casi, che indichiamo con C(t). Tuttavia, per la stima di R(t) ci interessa l'andamento di CR(t), cioè il modo in cui aumenta e diminuisce, piuttosto che il suo valore giorno dopo giorno, come spiegheremo più avanti. Se per esempio in un certo periodo di tempo il numero di tamponi fatti al giorno tende ad aumentare, il valore di CR(t) tenderà a crescere più rapidamente rispetto alla crescita del numero reale di nuovi casi e questo porterà a una sovrastima del valore di R(t). Al contrario, se in un certo periodo di tempo il numero di tamponi tende a ridursi progressivamente, avremo una decrescita di CR(t) che non fotografa l’andamento reale della diffusione della malattia. Nella figura qui sotto vediamo il confronto tra la curva CR(t) e le curve del tasso di positività del numero di tamponi effettuati al giorno, che indichiamo con PT(t), e del numero di casi testati al giorno, che indichiamo con PC(t), nel caso dell’Italia tra il 1 agosto 2020 e il 3 gennaio 2021 (fonte dati: Regioni/PPAA-Ministero della Salute-ISS).

In particolare, notiamo dopo il picco pandemico una diminuzione più veloce dei nuovi casi registrati rispetto al tasso di positività dei tamponi e del numero dei casi testati al giorno, dovuta evidentemente alla diminuzione del numero di tamponi effettuati. Dopo il picco si nota anche una notevole differenza negli andamenti del tasso di positività rispetto ai tamponi PT(t) e rispetto ai casi testati PC(t), quest'ultima appare più livellata.

Questa differenza tra PT(t) e PC(t) si evidenzia nel grafico sottostante, dove è mostrato il rapporto tra il numero di tamponi giornalieri e il numero di casi testati al giorno. Si nota che dall’inizio di novembre questo rapporto è andato aumentando e ciò ha determinato una diminuzione più veloce di PT(t) rispetto a PC(t).

È quindi chiaro che ci dobbiamo aspettare delle differenze nelle stime di R(t) ottenute considerando queste tre diverse curve. Si può inoltre affermare che certamente il numero dei nuovi casi registrati CR(t) è la curva meno affidabile fra le tre, dal momento che esiste una evidente correlazione con i tamponi fatti. Inoltre, anche il tasso di positività dei tamponi giornalieri PT(t) presenta una notevole criticità poiché che il numero di tamponi effettuati per singolo caso testato evidentemente non è una costante. Infine, come ben spiegato in quest'articolo di Cesare Cislaghi, anche la curva del numero dei casi testati al giorno PC(t), pur apparendo come la più affidabile, non fotografa esattamente la diffusione della malattia, perché il suo valore dipende anche dalla strategia con cui si effettuano i test.

Analizziamo adesso le problematiche della curva relativa al numero degli attuali infetti registrati IR(t). Ogni giorno questo numero aumenta di una quantità pari al numero di nuovi casi registrati CR(t) e diminuisce di una quantità pari alla somma dei guariti e dei morti. Quindi nella parte in salita IR(t) risente dei problemi di CR(t), mentre in discesa risente di un’importante questione che adesso descriviamo. È un fatto noto che il sistema di tracciamento non riesce, in particolare nelle fasi di picco, a identificare tutti gli infetti. In quelle fasi tra i casi registrati ci sarà un eccesso di casi gravi, che necessitano di un maggiore tempo di ospedalizzazione o isolamento e il cui decorso è in media più lungo, sia che porti alla guarigione che alla morte. È ragionevole dunque aspettarsi che il tasso di decrescita di IR(t) sia in generale sovrastimato rispetto a quello degli infetti reali. La principale conseguenza di ciò è che la curva di IR(t) risulterà eccessivamente ritardata rispetto a CR(t), come vediamo bene nella figura qui sotto, dove notiamo una differenza tra 15 e 20 giorni tra i picchi di CR(t) e IR(t) mentre sulla base della stima del tempo medio dell’infezione4 ci aspetteremmo un ritardo non superiore a 10 giorni (questo punto sarà ripreso nel seguito). Inoltre, sempre per lo stesso motivo, notiamo una velocità di discesa molto bassa di IR(t): questo perché i malati gravi per COVID-19 necessitano di tempi lunghi sia per guarire che per morire.

È quindi ragionevole pensare che la stima di R(t) risenta degli stessi problemi se questa viene fatta sulla base delle curve CR(t) o IR(t).

La nostra proposta è quindi di stimare R(t) a partire dai tassi di positività dei tamponi PT(t) e dei casi testati PC(t) ritenuti più affidabili. Nel seguito descriviamo il nuovo metodo dal punto di vista matematico, introducendo prima brevemente il metodo di cui ha parlato Roberto Battiston. Nell'ultima sezione mostriamo i risultati del nuovo metodo.

Descrizione del nuovo metodo

Per prima cosa ricaviamo l’equazione per R(t) ottenuta con il metodo descritto da Battiston partendo dal modello SIR. Quello che ci serve è l’equazione differenziale per il numero degli infetti I(t). Questa equazione ha un termine di crescita che dipende dal numero di nuovi casi C(t) e un termine di perdita che dipende dal fatto che i malati guariscono o muoiono:

\begin{equation} \frac{dI(t)}{dt} = C(t) - \frac{1}{\tau}I(t) \end{equation}

Il numero di nuovi casi C(t) è proporzionale al numero di suscettibili S(t) tramite due fattori. Il primo è il numero di contatti effettivi che ciascun suscettibile ha nell’unità di tempo; chiamiamo questo fattore β e lo definiamo così β=R(t)/τ, dove R(t) è proprio il numero di riproduzione netto che vogliamo stimare e τ è il tempo medio dell’infezione. Il secondo fattore è la probabilità di aver incontrato un infetto, ovvero il rapporto tra il numero degli infetti I(t) e il numero totale di individui della popolazione N(t). L'equazione differenziale che governa l'andamento degli infetti è dunque:

\begin{equation} \frac{dI(t)}{dt} = \frac{R(t)}{\tau}\frac{I(t)}{N(t)} S(t) - \frac{1}{\tau} I(t) \end{equation}

Facendo l’approssimazione S(t)/N(t)$\approx$1 otteniamo la relazione seguente:

\begin{equation} \frac{dI(t)}{dt}\approx\frac{R(t)}{\tau}I(t) - \frac{1}{\tau}I(t)= \frac{R(t)-1}{\tau} I(t) \end{equation}

L'equazione che abbiamo appena scritto ci permette di trovare il valore di R(t) a partire dall'andamento del numero di infetti I(t):

\begin{equation} R(t) = 1+ \tau \frac{1}{I(t)}\frac{dI(t)}{dt} = 1 + \tau\frac{d}{dt}ln\left(I(t)\right) \end{equation}

Questa stima di R(t) fornisce risultati molto simili a quelli ottenuti con il metodo proposto da Imperial College2. Idealmente, nelle equazioni scritte C(t) e I(t) sarebbero rispettivamente i nuovi casi reali e gli infetti realmente circolanti. In pratica, però, nell’implementazione del metodo viene usato il valore degli infetti registrati IR(t) e quindi in generale si ottiene un valore di R(t) diverso da quello reale.

Veniamo adesso al nuovo metodo che proponiamo per la stima di R(t). Facciamo subito notare che siamo autorizzati a sostituire nell’equazione per R(t) la funzione I(t) con la funzione k I(t), dove k è una costante qualsiasi. Questo è possibile perché nell'equazione per R(t) compare solo la derivata del logaritmo di I(t) e questa è uguale alla derivata del logaritmo di k I(t). Questa identità è ciò che ci consente di usare il nostro metodo:

\begin{equation} R(t)=1 + \tau \frac{d}{dt}ln\left(I(t)\right) = 1 + \tau\frac{d}{dt}ln\left(k\:I(t)\right) \end{equation}

Sfruttare questa semplice proprietà matematica ci torna utile per il motivo seguente. È ragionevole pensare che il tasso di positività, sia esso quello relativo al numero di tamponi che quello relativo al numero di casti testati, sia direttamente proporzionale al numero di nuovi casi giornalieri tramite un fattore che non sappiamo stimare ma che possiamo assumere costante nel tempo. Questo si traduce nella relazione seguente:

\begin{equation} P(t) = k\,C(t) \end{equation}

dove, come abbiamo anticipato, possiamo pensare di sostituire al tasso di positività P(t) sia PT(t) che PC(t). Introduciamo quindi una nuova funzione X(t) proporzionale al numero di infetti I(t) tramite la costante k:

\begin{equation} X(t)=k\,I(t) \end{equation}

Possiamo facilmente trovare l'equazione che governa l'andamento di X(t) nel tempo partendo da quella per il numero di infetti I(t)

\begin{equation} \frac{dI(t)}{dt}=C(t)-\frac{1}{\tau}I(t)\rightarrow k\,\frac{dI(t)}{dt}=k\,C(t)-\frac{1}{\tau}k\,I(t) \rightarrow \frac{dX(t)}{dt}=P(t)-\frac{1}{\tau}X(t) \end{equation}

Il metodo che proponiamo si riduce all'impiego delle due equazioni seguenti:

\begin{equation} \begin{cases} \frac{dX(t)}{dt}&=P(t)-\frac{1}{\tau}X(t)\\ R(t)&= 1+ \tau\:\frac{d}{dt}\left(ln(X(t))\right) \end{cases} \end{equation}

la prima ci permette di calcolare, risolvendo un'equazione differenziale nel modo che descriveremo fra poco, la funzione ausiliaria X(t) a partire dalla curva relativa al tasso di positività dei tamponi o dei casi testati, a seconda che sostituiamo a P(t) la PT(t) oppure PC(t). Una volta trovata X(t) possiamo sostituirla nella seconda equazione e ricavare l'andamento nel tempo del numero di riproduzione netto dell'epidemia, cioè R(t). Nel seguito applichiamo questo metodo ai dati che descrivono l'andamento dell'epidemia in Italia dal 1 agosto 2020 al 3 gennaio 202.

Risultato

Riportiamo sia i risultati ottenuti per R(t) che per X(t), perché è interessante sia il confronto tra gli R(t) ottenuti con metodi diversi che il confronto tra gli infetti registrati IR(t) e la funzione X(t) che dovrebbe essere proporzionale agli infetti reali.

Per ottenere la funzione X(t) si risolve numericamente l’equazione differenziale che abbiamo riportato nella sezione precedente. Per fare questo è necessario imporre una condizione iniziale non nota per il valore X(0). Data la forma dell’equazione, appare ragionevole imporre come condizione iniziale X(0)=τP(0), considerando che la soluzione tenderà a stabilizzarsi sul valore giusto con un tempo caratteristico τ. In appendice abbiamo riportato anche le soluzioni per valori diversi della condizione iniziale per mettere in evidenza questo tipo di comportamento. Si può affermare che non c’è una dipendenza critica dalla condizione iniziale vista la linearità della parte omogenea dell’equazione. Per ottenere la funzione X(t) abbiamo utilizzato un valore di τ=9 giorni4. Sempre in appendice viene analizzato come il risultato finale dipende dalla scelta di τ.

Nella figura qui sotto vediamo il confronto tra gli infetti registrati IR(t) e le funzioni XT(t) e XC(t) ottenute rispettivamente utilizzando le curve curve PT(t) e PC(t) riportate nella prima figura. Dal confronto emergono i difetti della curva IR(t) descritti nell’introduzione. Si nota un ritardo della IR(t) rispetto alle X(t) dovuto al fatto che evidentemente il tempo di guarigione o di morte degli infetti registrati è significativamente maggiore (circa il doppio) del tempo medio dell’infezione τ. Si nota anche come il livellamento della curva PC(t) dopo il picco impedisca alla curva XC(t) di decadere come la XT(t).

Nella grafico sottostante vediamo invece il confronto tra RS(t), RT(t) e RC(t) che sono rispettivamente la curva di R(t) ottenuta con il metodo descritto da Roberto Battiston qui partendo dalla curva IR(t) e le curve di R(t) ottenute con il nuovo metodo utilizzando i tassi di positività dei tamponi PT(t) e dei casi testati PC(t).

Ovviamente il confronto tra RS(t), RT(t) e RC(t) è strettamente correlato a quello tra IR(t), XT(t) e XC(t). Notiamo in generale un ritardo di RS(t) rispetto a RT(t) e RC(t). Notiamo che RT(t) e RC(t) raggiungono un primo picco importante già nella prima parte di agosto, mentre sulla RS(t) il picco è ritardato. Si nota una discesa ritardata di RS(t) rispetto a RT(t) e RC(t) dopo il picco di fine ottobre e anche un valore di RC(t) più vicino a 1 nell’ultimo periodo dovuto al livellamento della curva XC(t) dopo il picco. Infine, notiamo che il fatto che sul picco maggiore RS(t) è maggiore di RT(t) e RC(t) e ciò è una conseguenza del fatto che abbiamo avuto un picco di casi registrati guidato dal forte aumento di tamponi.

Per ottenere curve di R(t) più lisce è stato fatto uno smoothing delle P(t) prima di trovare la soluzione X(t) e anche della funzione IR(t). Questo è importante per evitare di vedere le oscillazioni dovute al sistema di tracciamento dei casi che tipicamente mostra un minimo il lunedì e quindi è caratterizzato da una struttura con oscillazione periodica di una settimana. Lo smoothing consiste in una media mobile settimanale con pesi gaussiani.

Un altro confronto interessante è quello tra il metodo che qui proponiamo e la stima del valore di R(t) regionale che la Fondazione Bruno Kessler effettua per conto dell'Istituto Superiore di Sanità effettua settimanalmente secondo il metodo descritto qui (questa stima è estratta dai report settimanali del Monitoraggio Fase 2 pubblicati sul sito del Ministero della salute). Riportiamo qui sotto il grafico relativo al caso della Lombardia tra il 1 giugno e il 4 gennaio 2020.

[commento - credo che l'R(t) stimato dall'ISS sia una media dei valori giornalieri ottenuti tramite metodo imperial college sui 14 giorni precedenti]

Conclusioni

Il metodo proposto serve a stimare il valore di R(t) evitando i difetti della curve CR(t) e IR(t) che sono dovuti a due fattori:

  • la curva CR(t) risente del numero di tamponi fatti e quindi ha un andamento che non dipende esclusivamente dalle variazioni nella diffusione della malattia;
  • la curva IR(t) risente sia del contributo di CR(t) sia del fatto che il tempo medio di perdita dei guariti e dei morti è notevolmente più lungo del tempo medio dell’infezione a causa del fatto che soprattutto nei picchi pandemici abbiamo un sistema di tracciamento in crisi che rileva prevalentemente solo i casi gravi che evidentemente hanno un tempo di permanenza maggiore.

Dal confronto tra RS(t), RT(t) e RC(t) appare un generale ritardo della RS(t) rispetto alle altre due curve dovuto al lungo tempo di permanenza dei contagiati, inoltre, si nota una sovrastima di RS(t) sul picco dovuta evidentemente al forte aumento di tamponi tipico delle fasi di forte aumento dei casi precedenti al picco. Appare quindi ragionevole prendere in considerazione un metodo per la stima di RT(t) che non parta dalle curve CR(t) e IR(t), che sono inaffidabili, ma che viene implementato partendo dai tassi di positività PT(t) e PC(t) che sono certamente più affidabili

Ringraziamenti: ringrazio il Prof. Carlo La Vecchia per i suggerimenti e il supporto.

Appendice

In questa appendice affrontiamo tre argomenti che riguardano le proprietà matematiche del metodo. In particolare tratteremo la dipendenza dalla condizione iniziale X(0), la dipendenza di R(t) dal parametro τ e infine analizziamo l’effetto introdotto da τ per capire meglio il ritardo tra la curva degli infetti registrati IR(t) e la curva X(t) che abbiamo introdotto per stimare il numero di riproduzione. Per comodità riportiamo qui sotto le equazioni più importanti che abbiamo scritto nell'articolo.

\begin{align} \textstyle \frac{dI(t)}{dt} &= \textstyle C(t) - \frac{1}{\tau}I(t) \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(1)\\ \textstyle R(t) &=\textstyle 1 + \tau\frac{d}{dt}ln\left(I(t)\right)\;\;\;\;\;\;\;\;\;\;\;\;\;(2)\\ &\\ \textstyle P(t) &=\textstyle k C(t),\;\; \textstyle X(t) = kI(t)\;\;\;\;\;\;\;\;(3)\\ &\\ \textstyle\frac{dX(t)}{dt}&=\textstyle P(t)-\frac{1}{\tau}\:X(t)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\,(4)\\ \textstyle R(t)&= \textstyle 1+ \tau\:\frac{d}{dt}\left(ln(X(t))\right)\;\;\;\;\;\;\;\;\;\;(5)\\ \end{align}

Dipendenza dalla condizione iniziale: nella figura sottostante vediamo le soluzioni ottenute per R(t) per scelte diverse della condizione iniziale per la X(0). Come si vede grazie al termine di perdita dell’equazione, la memoria della condizione iniziale viene persa nel tempo (con tempo caratteristico τ) e tutte le soluzioni convergono. Quindi per serie temporali lunghe il problema della condizione iniziale non sussiste. Una scelta ragionevole è quella di porre come condizione iniziale X(0)=τP(0) all’interno di un intervallo temporale a bassa crescita/decrescita in modo che la derivata temporale di X(t) abbiamo un effetto trascurabile.

Dipendenza di R(t) dal parametro τ: R(t) ha sostanzialmente una dipendenza lineare da τ. Per convincercene possiamo risolvere l'equazione differenziale che governa l'andamento di X(t), l'equazione (4), e calcolare la derivata del suo logaritmo per diversi valori di τ. Il risultato, mostrato nel grafico sottostante nel caso in cui si calcoli X(t) a partire dal tasso di positività dei tamponi, ci dice che questa derivata è debolmente dipendente da τ quasi ovunque. Ora R(t) dipende τ attraverso il secondo termine dell'equazione (5), che è uguale al prodotto tra τ e la derivata del logaritmo di X(t). Visto che, come abbiamo appena dimostrato, questa derivata è indipendente da τ, possiamo concludere che R(t) è lineare in τ con buona approssimazione e quindi è facilmente prevedibile l’effetto dell’incertezza sulla stima di τ sulla stima finale di R(t).

Effetto di τ su X(t) (e quindi anche su IR(t)): se si esegue tra trasformata di Fourier sull’equazione per X(t) (equazione 4) si ottiene subito l’equazione seguente

\begin{equation} i\omega\tilde{X}(\omega)=\tilde{P}(\omega)-\frac{1}{\tau}\tilde{X}(\omega) \end{equation}

dove i è il numero immaginario di Eulero. La relazione precedente può essere riscritta come

\begin{equation} \tilde{X}(\omega)=\frac{\tau}{1+i\omega\tau}\tilde{P}(\omega) \end{equation}

che ci dice che la funzione di trasferimento è quella di un filtro passa basso del primo ordine con una frequenza di taglio f=1/2πτ. Questo filtro introduce uno smoothing dell’andamento a causa di questa frequenza di taglio con conseguente ritardo temporale. Nel caso della IR(t) questo effetto di taglio è più pesante perché il tempo τ effettivo risulta maggiore a causa del fatto che il sistema di tracciamento individua prevalentemente casi più gravi con tempi di permanenza più lunghi.

Note
1Nel simbolo "R(t)" viene reso esplicito il fatto che il numero di riproduzione R dell'epidemia cambia con il tempo (t è il simbolo che rappresenta la variabile tempo). In tutto il resto dell'articolo viene usata la stessa convenzione: quando una certa quantità cambia nel tempo il suo simbolo viene fatto seguire dall'espressione "(t)". Così il numero di nuovi casi giornalieri registrati sarà abbreviato con Cr(t) poiché ogni giorno questo numero cambia, e così via per tutte le altre quantità.
2 Cori A, Ferguson NM, Fraser C, Cauchemez S. A new framework and software to estimate time-varying reproduction numbers during epidemics. Am J Epidemiol. 2013;178(9):1505-1512. doi:10.1093/aje/kwt133
3 Kermack, W. O., McKendrick, A. G. (1927). "A Contribution to the Mathematical Theory of Epidemics". Proceedings of the Royal Society A. 115 (772): 700–721. doi:10.1098/rspa.1927.0118.
4 D Cereda, M Tirani, F Rovida, V Demicheli, M Ajelli, P Poletti, F Trentini, G Guzzetta, V Marziano, A Barone, M Magoni, S Deandrea, G Diurno, M Lombardo, M Faccini, A Pan, R Bruno, E Pariani, G Grasselli, A Piatti, M Gramegna, F Baldanti, A Melegaro, S Merler. The early phase of the COVID-19 outbreak in Lombardy, Italy. ArXiv:2003.09320v1, 2020.

Aiuta Scienza in Rete a crescere. Il lavoro della redazione, soprattutto in questi momenti di emergenza, è enorme. Attualmente il giornale è interamente sostenuto dall'Editore Zadig, che non ricava alcun utile da questa attività, se non il piacere di fare giornalismo scientifico rigoroso, tempestivo e indipendente. Con il tuo contributo possiamo garantire un futuro a Scienza in Rete.

E' possibile inviare i contributi attraverso Paypal cliccando sul pulsante qui sopra. Questa forma di pagamento è garantita da Paypal.

Oppure attraverso bonifico bancario (IBAN: IT78X0311101614000000002939 intestato a Zadig srl - UBI SCPA - Agenzia di Milano, Piazzale Susa 2)

altri articoli

Pollution and Covid. Two vague clues don't make an evidence

In these days, newspapers and television programs (and the web, of course) are giving space to a statement by the Italian Society of Environmental Medicine (SIMA) announcing important discoveries on the link between airborne particulate matter and Coronavirus, even describing them as important for the decisions to be taken in the coming weeks.