Calcolo dell'indice di riproduzione \(R_t\) del Covid-19

Introduzione

La pandemia dovuta al virus SARS-CoV-2 e alle sue varianti continua a condizionare le nostre vite, sia nei contatti personali in ambienti ristretti e chiusi ma soprattutto là dove vi sia il rischio della presenza contemporanea di più persone o di affollamenti. La disponibilità di vaccini efficaci ha comunque fornito uno strumento fondamentale al contenimento della sua diffusione e a livello individuale con il completamento del ciclo vaccinale, una significativa e indubbia protezione nei riguardi dei sintomi e delle conseguenze di una seria patologia qual è quella di Covid-19. La ancora non completa copertura vaccinale favorisce comunque l'insorgere di mutazioni con maggiori capacità di replicazione e trasmissibilità e, non ultimo, potenzialmente in grado di ridurre l'efficacia degli attuali vaccini e portare ad una maggiore gravità della malattia. Per garantire a ciascun individuo, alla comunità sociale e alle sue istituzioni una significativa ripartenza è pertanto ancora necessario aderire a tutte quelle misure sanitarie e comportamentali che limitino la diffusione di nuove varianti.
Nonostante questi condizionamenti confidiamo che nella scuola, verso cui si indirizza questo lavoro, si possa affermare una stabile ripresa in presenza delle attività scolastiche e didattiche. Solo così e con il recupero della normalità dei rapporti interpersonali che sempre ha caratterizzato questa comunità si potranno garantire agli studenti significativi percorsi educativi e validi apprendimenti.

Tra i numerosi indicatori del rischio epidemico certamente uno dei più importanti è l'indice di riproduzione o indice di contagio o anche tasso di contagiosità \(R_t\). Nel nostro scritto intendiamo presentare, sulla base di un semplice modello matematico, una procedura di calcolo per tale indice in grado di fornire valori compatibili con quelli diffusi settimanalmente da siti istituzionali quali l'Istituto Superiore di Sanità e Epicentro.iss. In particolare si intende riprodurre quanto proposto giornalmente dal sito CovidStat nelle pagine riguardanti questo indicatore e, per quanto riguarda il metodo di calcolo, utilizzare la modalità più semplice tra quelle proposte dal medesimo gruppo e pubblicate nel giornale The European Physical Journal Plus (2021) [1]. Con ciò si vuole offrire a studenti delle classi terminali utili spunti multidisciplinari di approfondimento delle proprie conoscenze e competenze, occasioni peraltro già presenti nel blog orientato alla didattica nel sito CovidStat stesso.
I requisiti a tale scopo necessari sono abbastanza circoscritti in quanto includono solamente le nozioni collegate a

  • progressioni geometriche,
  • funzione esponenziale,
  • regressione lineare,
  • media aritmetica mobile.

Circa gli strumenti software che utilizziamo per l'elaborazione, questi si limitano ad un comune foglio di calcolo mentre, per sottolineare l'importanza in ambito scientifico di strumenti di elaborazione più avanzati, associamo a tale foglio di calcolo cinque notebook di Jupyter nei quali, in linguaggio Python, sviluppiamo, commentandolo, l'intero processo di calcolo. In particolare in questo lavoro intendiamo:

  • comprendere i significati degli indici di riproduzione \(R_0\) e \(R_t\) tramite una seppur breve indagine teorica affiancata a delle esemplificazioni.
  • Sviluppare una semplice procedura di calcolo che permetta di ridurre l'impatto delle variazioni settimanali così da
  • estrarre una rappresentazione attendibile dell'andamento effettivo della pandemia.

In base ad alcune ipotesi su una diffusione di tipo esponenziale nel breve periodo della pandemia, giungiamo a

  • definire un modello per l'indice di riproduzione per poi
  • dedurre una sua espressione teorica in funzione dei parametri a tal fine introdotti.
  • La stima di questi parametri permette infine di
  • calcolare l'indice \(R_t\) e
  • rappresentare il suo andamento temporale.

Indice

Due indici di riproduzione: \(R_0\) e \(R_t\)

L'indice \(R_0\)

In questa sezione presentiamo un modello esemplificativo di una possibile evoluzione temporale che, pur nella sua schematicità, ci permette di cogliere alcuni importanti aspetti delle grandezze coinvolte.
Nella figura 1 è rappresentata l'evoluzione temporale di una popolazione che a partire da un capostipite si sviluppa ad ogni generazione raddoppiando il proprio numero.

Evoluzione di una popolazione con R0=2
Fig. 1. Evoluzione di una popolazione con \(R_0=2\).

In altre parole, in un certo istante iniziale un soggetto infetto (il cosiddetto paziente zero) dà inizio alla diffusione di un virus e a seguito dei suoi contatti con altri soggetti, ne permette la propagazione. Questo virus è caratterizzato da una infettività tale che, nel nostro esempio, dà origine ad ogni generazione ad un raddoppio degli infetti. L'intervallo di tempo che separa le coppie formate dall'individuo che contagia con il contagiato, in sostanza l'intervallo temporale che distingue una generazione di infetti con la successiva, è evidentemente un parametro che può variare e che indichiamo con \(g\) (più avanti la definizione formale). Se indichiamo con \(I_t\) il numero degli infetti della quarta generazione e risaliamo indietro nel tempo di generazione in generazione possiamo esprimere il numero degli infetti come

\begin{equation} I_t=2 I_{t-g}, \qquad I_{t-g}=2 I_{t-2g}, \qquad I_{t-2g}=2 I_{t-3g}, \qquad I_{t-3g}=2 I_{t-4g},\qquad I_{t-4g}=I_0=1, \end{equation}

Dalla forma ricorsiva di tale successione possiamo risalire a quella esplicita collegando la generazione al tempo \(t\) con gli individui che l'hanno generata e sintetizzarla come

\begin{equation} \eqalign{I_t&=2 I_{t-g}\\[3pt] &=2 (2I_{t-2g})=2^2 I_{t-2g}\\[3pt] &=2^2(2I_{t-3g})=2^3 I_{t-3g}\\[3pt] &=2^3(2I_{t-4g})=2^4 I_0=2^4.\cr} \end{equation}

Questo risultato mostra come il meccanismo proposto comporti una crescita esponenziale con base pari a 2: potremmo quindi identificare il coefficiente di riproduzione \(R_0\) di tale virus con il valore \(R_0=2\) e riassumere il legame tra il valore iniziale e quello al tempo \(t\) come

\begin{equation} I_t=(R_0)^4I_{t-4g}.\tag{1}\label{eq01} \end{equation}

La successiva figura 2 rappresenta ancora l'evoluzione di quattro generazioni di una popolazione di infetti questa volta descritta dalla successione

Evoluzione di una popolazione con R0=1/2
Fig. 2. Evoluzione di una popolazione con \(R_0=1/2\).
\begin{equation} I_t=\frac{1}{2}I_{t-g}, \qquad I_{t-g}=\frac{1}{2} I_{t-2g}, \qquad I_{t-2g}=\frac{1}{2} I_{t-3g}, \qquad I_{t-3g}=\frac{1}{2}I_{t-4g},\qquad I_{t-4g}=I_0=16 \end{equation}

e, riportandoci alla forma analitica a partire dalla precedente ricorsiva, abbiamo

\begin{equation} \eqalign{I_t&=\frac{1}{2} I_{t-g}\cr &=\frac{1}{2} \biggl(\frac{1}{2} I_{t-2g}\biggr)=\frac{1}{2^2} I_{t-2g}\cr &=\frac{1}{2^2}\biggl(\frac{1}{2} I_{t-3g}\biggr)=\frac{1}{2^3} I_{t-3g}\cr &=\frac{1}{2^3}\biggl(\frac{1}{2} I_{t-4g}\biggr)=\frac{1}{2^4} I_{t-4g}\cr} \end{equation}

che sintetizziamo come

\begin{equation} I_t=(R_0)^4 I_{t-4g}, \qquad \hbox{con}\qquad R_0=\frac{1}{2}.\tag{2}\label{eq02} \end{equation}

Questi due esempi mostrano, come a partire da una relazione del tipo,

\begin{equation} \bbox[15px,#ffffcc, border:1px solid red]{ I_t=R_0 I_{t-g}\tag{3}\label{eq03} } \end{equation}

e cioè di una proporzionalità tra gli infetti di una generazione e quelli della successiva, si riscontri una crescita esponenziale se \(R_0>1\) mentre, al contrario, la popolazione si contrae esponenzialmente se \(R_0<1\). Seppure in un contesto esemplificativo ciò rappresenta il principale significato del parametro \(R_0\) detto anche numero di riproduzione di base: pertanto il

Numero di riproduzione di base \(R_0\)
rappresenta il numero di individui in una popolazione suscettibile che vengono contagiati da un singolo soggetto infetto nel periodo in cui questo può diffondere l'infezione.

Il suo valore maggiore o minore di uno informa sulla diffusione o sulla contrazione di una popolazione caratterizzata da una qualche proprietà. L'altro parametro coinvolto nei due esempi riguarda il parametro \(g\): a questo si dà il nome di tempo di generazione o generation time e a sua volta costituisce una importante caratteristica dei virus in quanto

Tempo di generazione \(g\)
rappresenta l'intervallo di tempo che intercorre tra una generazione di infetti e la successiva.

Da \(R_0\) all'indice \(R_t\)

Le relazioni \eqref{eq01} e \eqref{eq02} collegano il numero degli infetti \(I_t\) con gli infetti di quattro generazioni precedenti ma, partendo dalla relazione ricorsiva \eqref{eq03}, possiamo facilmente esprimere \(I_t\) come l'esito di una successione di \(N\) cicli di riproduzione. Difatti procedendo a ritroso da una generazione alla precedente discende formalmente

\begin{equation} \eqalign{I_t&=R_0 I_{t-g}\\[3pt] &=R_0 (R_0 I_{t-2g})=R_0^2 I_{t-2g}\\[3pt] &=R_0^2(R_0I_{t-3g})=R_0^3 I_{t-3g}\\[3pt] &\vdots\cr &=R_0^N I_{t-Ng}\cr} \end{equation}

per cui nell'arco temporale \(\Delta t=(t-Ng)-t=Ng\) ossia di \(N\) generazioni otteniamo ancora un legame di tipo esponenziale della forma

\begin{equation} \bbox[15px,#ffffcc, border:1px solid red]{ I_t=R_0^N I_{t-Ng}.\tag{4}\label{eq04} } \end{equation}

A questo punto va osservato come non sia possibile conoscere il valore istantaneo degli infetti \(I_t\) mentre possiamo solo stimare i suoi valori medi supposto che si siano eliminati o almeno ridotti tutti quei fattori e approssimazioni che sono indipendenti dalla diffusione del virus ma che possono comunque condizionare il numero rilevato di infetti. Tra questi fattori, come vedremo nella sezione successiva, appare evidente come la variazione del numero dei nuovi positivi sia modulata da un ciclo settimanale. Inoltre, pur supponendo che il valore di \(R_0\) dipenda almeno nel breve periodo dalle caratteristiche intrinseche del virus, il numero degli infetti delle generazioni successive certamente dipende pure

  • dal tempo di generazione \(g\) che, a prescindere dal nostro schema dove lo si è implicitamente ritenuto costante, questo non lo è necessariamente. Nel seguito potremo assegnare a tale grandezza solo una stima approssimativa. Inoltre, e ovviamente, \(I_t\) dipende
  • da tutte quelle misure di prevenzione che vengono attuate per limitarne l'aumento e quindi dall'impatto di tutti i conseguenti provvedimenti di contenimento presi a tal fine.

Se quindi il numero di riproduzione \(R_0\) si può intendere come una caratteristica a priori della diffusione di un virus quando questo non sia ostacolato da interventi mirati a ridurne la diffusione, con l'indice \(R_t\) si intende invece misurare il numero di riproduzione ma a posteriori considerando la situazione di fatto che si va delineando nell'evoluzione della pandemia stessa: è quindi una funzione che evolve del tempo e dipendente dalla serie temporale dei dati che si prendono in esame e dai metodi che si utilizzano per il suo calcolo. In altre parole, se sull'indice \(R_0\) non si può intervenire, possiamo invece agire sull'\(R_t\) per ridurlo a valori inferiori ad 1 ostacolando il contagio con tutti quegli accorgimenti come il distanziamento sociale, l'uso delle mascherine, il lavaggio delle mani che, oramai da tempo, ben conosciamo e, non ultimo, il sottoporsi a vaccinazione.
In base a tali considerazioni possiamo concludere che l'indice \(R_t\)

Numero di riproduzione effettivo \(R_t\) o tasso di contagio
rappresenta il numero medio di individui in una popolazione suscettibile che vengono contagiati da un singolo soggetto infetto nel periodo in cui questo può diffondere l'infezione e a seguito di tutte le misure attuate per rallentare il contagio.

Modifichiamo pertanto la relazione \eqref{eq04} nella più realistica

\begin{equation} \bbox[15px,#ffffcc, border:1px solid red]{ E(I_t)=R_t^N I_{t-Ng}.\tag{5}\label{eq05} } \end{equation}

dove l'indice \(R_0\) viene sostituito dall'indice di riproduzione effettivo \(R_t\) mentre il numero di infetti \(I_t\) viene sostituito da \(E(I_t)\) (expected value) e cioè dal numero medio di infetti.

Modulazione settimanale dei nuovi positivi

Il Dipartimento della Protezione Civile pubblica giornalmente alle ore 17 i dati nazionali e regionali provenienti da tutte le regioni allo scopo di fornire il quadro aggiornato dell'evoluzione della pandemia di COVID-19 [2]. Poiché in questo lavoro ci limitiamo allo studio dell'andamento nazionale e a quello delle singole regioni, i file di nostro interesse si riducono ai soli

Una volta prelevati e solo per evitare la sovrapposizione con file precedenti, rinominiamo rispettivamente questi come datiNazionaliYYYYMMDD.csv e datiRegionaliYYYYMMDD.json dove in luogo di YYYYMMDD inseriamo la data nel formato ISO del giorno a cui si riferiscono.

Entrambi i file sono leggibili con qualsiasi editor ma, mentre nel file con estensione CSV (Comma-Separated Values) i numerosi campi sono riportati in forma tabellare e sono facilmente interpretabili e importabili in un foglio di calcolo, per il formato JSON (JavaScript Object Notation) dei dati regionali, pur in forma strutturata, si mostra necessario un interprete (parser) adeguato che, per noi, sarà costituito da un modulo Python. Per i nostri scopi sarà necessario il solo campo relativo al numero di individui rilevati infetti nelle ultime 24 ore ossia il campo dei nuovi_positivi: di questo campo studieremo innanzitutto l'evoluzione temporale.

Limitandoci ai dati nazionali (l'elaborazione di quelli regionali verrà svolta nei notebook di Jupyter), importiamo in un foglio di calcolo i soli nuovi_positivi (per i particolari tecnici si rimanda all'Appendice) e diamo di essi una prima rappresentazione grafica estesa ad un intervallo temporale di 2/3 mesi. La figura 3 mostra, sovrapposta ad una evidente variabilità, pure una altrettanto chiara variazione periodica di periodo pari ad una settimana. Questa modulazione settimanale è manifestamente collegata al diverso numero di tamponi e quindi di nuovi positivi effettuati nei giorni centrali della settimana rispetto a quelli effettuati nei fine settimana. Ne segue che il numero dei nuovi_positivi risulta chiaramente condizionato da fattori che nulla hanno a che fare con la diffusione del contagio ma, al contrario, dall'organizzazione del lavoro.

Variazione settimanale dei nuovi positivi (Italia)
Fig. 3. Variazione settimanale dei nuovi positivi (Italia).

È pertanto necessario ridurre se non eliminare questo fattore di disturbo (con una terminologia audio chiamiamolo pure rumore) per poter cogliere con maggior chiarezza l'andamento pregresso della pandemia e individuarne eventualmente la tendenza. E poiché i dati da elaborare costituiscono sostanzialmente una serie temporale la funzione matematica adeguata a tale scopo (tecnicamente detto smoothing) consiste nell'eseguire una loro media mobile.
Rimandando per ora i particolari della implementazione ad una sezione successiva, forniamo qui un esempio schematico di questo semplice algoritmo: supponiamo quindi di disporre di un insieme \(A\) costituito da 7 valori numerici

\begin{equation} A=\{I_1, I_2, I_3, I_4, I_5, I_6, I_7\} \end{equation}

e decidiamo di eseguire una media mobile su 4 elementi. Ciò comporta la partizione dell'insieme \(A\) in sottoinsiemi di 4 elementi, ciascuno dei quali si ottiene dall'insieme di partenza con una traslazione (offset) di un solo elemento: nel nostro esempio otteniamo i 4 sottoinsiemi

\begin{equation} a_1=\{I_1, I_2, I_3, I_4\},\quad a_2=\{I_2,I_3,I_4,I_5\}, \quad a_3=\{I_3,I_4,I_5,I_6\},\quad a_4=\{I_4,I_5,I_6,I_7\} \end{equation}

e, in corrispondenza di ciascuno, calcoliamo la media aritmetica ottenendo i 4 valori \(\{b_1,b_2,b_3,b_4\}\)

\begin{equation} b_1=\frac{1}{4}(I_1+I_2+ I_3+ I_4),\quad b_2=\frac{1}{4}(I_2+I_3+I_4+I_5), \quad b_3=\frac{1}{4}(I_3+I_4+I_5+I_6),\quad b_4=\frac{1}{4}(I_4+I_5+I_6+I_7). \end{equation}

In generale, se disponiamo di \(n\) rilevazioni e decidiamo di calcolare medie mobili coinvolgenti \(m\) elementi con \(m< n\) allora disporremo di \(n-m+1\) medie mobili.
Nelle due figure che seguono (figg. 4 e 5) si rappresentano in due modi diversi i valori giornalieri dei nuovi_positivi così da mostrarne variabilità settimanale e, rispettivamente, dispersione. Da queste immagini appare evidente come l'applicazione di una media mobile estesa ad un intervallo di 14 giorni, se da un lato riduce la sensibilità ai cambiamenti che avvengono nell'arco di qualche giorno, dall'altro riesca a eliminare la periodicità settimanale assorbendone gli effetti sistematici e quindi fornendo una più chiara valutazione sull'andamento dei nuovi positivi: appaiono in tal modo evidenti le tre ondate di contagi comprese tra il marzo 2020 e il maggio 2021.

Andamento giornaliero dei nuovi positivi e media mobile (Italia)
Fig. 4. Andamento giornaliero dei nuovi positivi e media mobile (Italia).
Andamento giornaliero degli infetti, media mobile e dispersione (Italia)
Fig. 5. Andamento giornaliero degli infetti, media mobile e dispersione (Italia).

Il modello e i suoi parametri

Le osservazioni di carattere teorico sulla diffusione di una pandemia sono concordi nel riconoscere come, nelle sue fasi iniziali e in assenza di interventi limitatori, l'andamento temporale sia descritto da una curva di tipo esponenziale. Successivamente comunque, per il venir meno in una popolazione ovviamente finita di individui suscettibili, la medesima curva non può certo descrivere adeguatamente l'andamento temporale (si veda la pagina didattica relativa ai più semplici modelli epidemiologici [3]). Ciò nonostante facciamo qui l'ipotesi che

Ipotesi
l'andamento della diffusione (o della regressione) dell'epidemia sia, in media e in un arco temporale ristretto, di tipo esponenziale [1].

In base a tale ipotesi potremo quindi utilizzare la già discussa relazione \eqref{eq05}. Poiché questa mette in relazione la media della popolazione di infetti al tempo \(t\) con la popolazione di \(N\) generazioni precedenti mentre le rilevazioni di cui disponiamo sono giornaliere conviene rivedere tale relazione con una ridefinizione degli indici in modo da collegare due rilevazioni giornaliere separate da un intervallo temporale pari ad un multiplo \(N\) del tempo di generazione \(g\) ossia \(Ng\).

Decidiamo quindi di coinvolgere \(n\) rilevazioni giornaliere successive di nuovi infetti e poniamo in \eqref{eq05} \(h=t-Ng\) l'istante che dà origine alle \(N\) generazioni successive comprese nell'intervallo \(t-h=Ng\). In termini di giorni trascorsi, possiamo individuare il medesimo intervallo semplicemente come \((n-1)\) sottintendendo qui e nello schema di figura 6 (ma pure in quanto segue), l'unità \(giorni\).

Due modi per individuare un intervallo temporale
Fig. 6. Due modi per individuare un intervallo temporale.

Abbiamo pertanto le due espressioni equivalenti

\begin{equation} t-h=Ng\qquad \hbox{oppure}\qquad t-h=n-1\tag{6}\label{eq06} \end{equation}

dalle quali discende anche

\begin{equation} N=\frac{n-1}{g}.\tag{7}\label{eq07} \end{equation}

A seguito di questa ridefinizione di indici la relazione \eqref{eq05} assume la forma

\begin{equation} E(I_t)=R_t^N I_{t-Ng}\qquad\Longrightarrow \qquad E(I_t)=R_t^{(n-1)/g}\,I_h.\tag{8}\label{eq08} \end{equation}

Sfruttando l'ipotesi iniziale, supponiamo che nell'intervallo \(t-h\) la progressione (o regressione) giornaliera del numero medio dei nuovi positivi \(E(I_t)\) segua un andamento esponenziale per cui il nostro modello si esplicita nella funzione

\begin{equation} \bbox[15px,#ffffcc, border:1px solid red]{ I_t=A e^{\lambda t}}\tag{9}\label{eq09} \end{equation}

dove le grandezze \(A\) e \(\lambda\) costituiranno degli ulteriori parametri da determinare. In sostanza, l'ipotesi posta implica che scelto un adeguato intervallo temporale, la curva dei nuovi positivi si possa pensare come l'unione di singoli andamenti esponenziali.
Inserendo quest'ultima espressione in entrambi i membri della \eqref{eq08} (pure \(I_h\) rappresenta un valore medio), questa diviene

\begin{equation} E(I_t)=R_t^{(n-1)/g}\,I_h\qquad\Longrightarrow \qquad Ae^{\lambda t}=R_t^{(n-1)/g}\cdot A e^{\lambda h} \end{equation}

per cui, eliminato il fattore comune \(A\),

\begin{equation} e^{\lambda t}=R_t^{(n-1)/g}\cdot e^{\lambda h} \end{equation}

possiamo esplicitare a primo membro il termine dipendente dall'indice \(R_t\)

\begin{equation} R_t^{(n-1)/g} = e^{\lambda (t-h)}. \end{equation}

Considerando il logaritmo di entrambi i membri

\begin{equation} \frac{n-1}{g}\cdot \ln(R_t)= \lambda(t-h) \end{equation}

e quindi pure

\begin{equation} \ln(R_t)=\frac{\lambda g(t-h)}{n-1}, \end{equation}

in base alla \eqref{eq06} quest'ultima espressione si riduce alla semplice relazione

\begin{equation} \ln R_t= \lambda g. \end{equation}

Ritornando infine all'esponenziale il coefficiente di riproduzione effettivo \(R_t\) risulta descritto in termini di soli due parametri e dalla relazione esponenziale

\begin{equation} \bbox[15px,#ffffcc, border:1px solid red]{ R_t= e^{\lambda g}.\tag{10}\label{eq10} } \end{equation}

Questo importante risultato ci informa che per ottenere una stima dell'indice \(R_t\) dovremo innanzitutto conoscere

  • il tempo di generazione \(g\) e
  • il parametro \(\lambda\).

Per una stima del primo termine ci dovremo rifare ai valori reperibili in alcune pubblicazioni mentre potremo determinare il parametro \(\lambda\) eseguendo una regressione esponenziale basata sulle \(n\) rilevazioni giornaliere della funzione scelta come modello ossia l'espressione \eqref{eq09}.
Un importante significato di quest'ultimo parametro lo si coglie se, introdotto il parametro \(\tau\) espresso in giorni, riscriviamo la \eqref{eq09} nella base 2 come

\begin{equation} I_t=Ae^{\lambda t}=A 2^{t/\tau}.\tag{11}\label{eq11} \end{equation}

Poiché dopo un intervallo temporale pari a \(t=\tau\) il numero di infetti raddoppia

\begin{equation} I_\tau=A2^{\tau/\tau}=2A, \end{equation}

si giustifica per questo nuovo parametro il termine di tempo di raddoppio. D'altra parte è anche per \eqref{eq11}

\begin{equation} Ae^{\lambda \tau}=2A \end{equation}

dalla quale deduciamo che

\begin{equation} \lambda \tau=\ln 2\qquad\Longrightarrow\qquad \lambda=\frac{\ln 2}{\tau}.\tag{12}\label{eq12} \end{equation}

Quest'ultima relazione mostra come il parametro \(\lambda\) sia strettamente legato al tempo di raddoppio. Ne segue che una diminuzione del tempo di raddoppio \(\tau\) si traduce in un aumento inversamente proporzionale di \(\lambda\) indicando così una espansione nella diffusione del virus: viceversa si avrebbe invece una riduzione nella diffusione se \(\tau\) aumentasse. Sulla base di queste osservazioni associamo pertanto al parametro \(\lambda\) il termine di tasso di crescita del contagio.

Linearizzazione di una funzione esponenziale

La regressione lineare di un insieme di \(n\) coppie \((x_i,y_i)\) di dati sperimentali

\begin{equation} \{(x_i,y_i): i=1,2,\dots, n\} \end{equation}

è un procedimento numerico ben noto e viene applicato per individuare la retta \(y=a x+b\) che minimizza la grandezza \(s^2(a,b)\)

\begin{equation} s^2(a,b)=\sum_{i=1}^n d_i^2=\sum_i^n[y_i-(ax_i+b)]^2 \end{equation}

somma dei quadrati delle "distanze" \(d_i\) dei punti sperimentali dai punti della retta stessa con uguale ascissa (fig. 7). Questa grandezza dipende manifestamente dai parametri incogniti \(a\) e \(b\), ma imponendo su tale funzione la condizione di minimizzazione è possibile individuare la coppia di parametri in corrispondenza dei quali si ottiene il miglior adattamento della retta all'insieme dei dati rilevati.

Punti sperimentali e retta ottimale che li interpola
Fig. 7. Punti sperimentali e retta ottimale che li interpola.

La ricerca del minimo della funzione \(s^2\) si può svolgere sia in termini analitici che numerici: nel primo caso, conduce a espressioni finite per i parametri \(a\) e \(b\) caratterizzanti la retta ottimale (vedere, per es. [4]).
Nel nostro caso, disponendo comunque di una funzione esponenziale qual è la \eqref{eq09} sarebbe necessaria una regressione di tipo esponenziale ma, data la sua forma è possibile riportare quest'ultima nell'ambito di una relazione lineare se consideriamo il logaritmo di entrambi i suoi membri ossia

\begin{equation} \ln(I_t)=\ln(Ae^{\lambda t}) \qquad\hbox{che equivale alla}\qquad \ln(I_t)=\ln(A)+\lambda t.\tag{13}\label{eq13} \end{equation}

Posto quindi \(y_t=\ln(I_t)\), la relazione esponenziale è ricondotta ad una lineare tra la variabile indipendente \(t\) ma avente come variabile dipendente la \(y_t\).

\begin{equation} y_t=\ln(A)+\lambda t.\tag{14}\label{eq14} \end{equation}

Potremo quindi giungere ad una stima del parametro \(\lambda\) da inserire nella \eqref{eq10} determinando il solo coefficiente angolare di quest'ultima relazione lineare.

Elaborazione

Come accennato nell'introduzione, per il calcolo del fattore \(R_t\) utilizziamo due strumenti software ampiamente disponibili ossia un comune foglio di calcolo e alcuni notebook di Jupyter basati sul linguaggio Python. Pur con le ovvie diversità sintattiche realizziamo in entrambi lo stesso procedimento ma nel foglio di calcolo trattiamo solo i dati nazionali elaborandoli secondo il diagramma di flusso di figura 8.

Diagramma di flusso del foglio di calcolo
Fig. 8. Diagramma di flusso del foglio di calcolo.

Nei notebook invece, data la possibilità di integrare il codice Python con il testo esplicativo, descriviamo passo passo le singole istruzioni trattando sia i dati nazionali che quelli regionali.

Il foglio di calcolo e medie mobili

Dopo aver importato in Excel (oppure, in alternativa, LibreOffice Calc) il campo dei nuovi_positivi (Appendice) associamo a ciascun dato un indice numerico sequenziale che tenga conto dei giorni trascorsi dall'inizio delle rilevazioni (24 febbraio 2020) ed, eventualmente, per ridurre i valori numerici coinvolti nel calcolo dividiamo per 1000 il numero dei nuovi positivi. Nelle due colonne a fianco inseriamo la formula per il calcolo della media mobile in un arco temporale di 7 giorni e, rispettivamente, di 14 giorni. Come detto questa procedura filtra, riducendolo, l'effetto delle variazioni settimanali dovute ad un minor numero di tamponi nei fine settimana.
In corrispondenza delle due scelte su cui eseguire la media mobile otteniamo gli andamenti di figura 9 associati, come nelle figg. 4 e 5 ai valori rilevati giornalmente.

Andamenti settimanali e bisettimanali con media mobile
Fig. 9. Andamenti settimanali e bisettimanali con media mobile.

Innanzitutto notiamo come la media mobile su base settimanale sia più sensibile alle variazioni giornaliere dei dati rispetto alla media eseguita su un intervallo di 14 giorni. Inoltre, allo scopo di migliorare l'accordo visivo con i dati giornalieri applichiamo ad entrambe uno shift temporale associando il valore calcolato giornalmente alla data antecedente di 3 giorni per la media settimanale e di 7 per quella bisettimanale.
Per un confronto si veda la medesima rappresentazione grafica nella pagina web Italia - positivi giornalieri.

Calcolo dell'\(R_t\) nazionale

Per giungere alla stima del tasso di crescita \(\lambda\) presente nella \eqref{eq10} attraverso la regressione lineare rappresentata dalla \eqref{eq14}, va innanzitutto calcolato il valore giornaliero del logaritmo dei nuovi positivi ossia \(\ln(I_t)\) (colonna G del foglio di calcolo). Tramite la funzione predefinita pendenza(y_note;x_note) che fornisce il coefficiente angolare della retta di regressione dell'insieme di coppie (y_note;x_note), determiniamo il valore \(\lambda_1\) anch'esso giornaliero del parametro \(\lambda\) su un intervallo di 14 giorni. La figura 10 mostra come si sia scelto di riferire in modo assoluto l'intervallo temporale \([1,14]\) e in modo relativo i valori \(\ln(I_t)\).

Calcolo del parametro λ_1 con media mobile
Fig. 10. Calcolo del parametro \(\lambda_1\) con media mobile.

La copia all'intera colonna di \(\ln(I_t)\) della formula riportata nella figura 10 permette di ottenere i valori del tasso di crescita giornaliero \(\lambda_1\) che, convenzionalmente, associamo nel foglio al giorno finale dell'intervallo bisettimanale (ma potremmo associarlo ad uno qualsiasi dei giorni coinvolti). Il corrispondente grafico rappresentativo (fig. 11) mostra una evidente variabilità giornaliera di questa grandezza cosicché, con le stesse modalità di filtraggio applicate ai nuovi_positivi, procediamo al calcolo della sua media mobile estesa anche in questo caso ad un intervallo temporale di 14 giorni. In tal modo otteniamo nella medesima figura 11 la linea continua del tasso di crescita \(\lambda\) (vedere fig. 6, Italy - growth rate, in [1]).

Valori giornalieri del tasso di crescita e sua media mobile λ
Fig. 11. Valori giornalieri del tasso di crescita e sua media mobile λ.

In base alla relazione \eqref{eq10} il calcolo dell'indice di riproduzione \(R_t\) segue immediato non appena sia definito il valore del tempo di generazione \(g\). A tal fine e come anticipato, riprendiamo la stima di questo parametro come viene riportata nelle pubblicazioni [1] e [5]. Queste assegnano a \(g\) un valore medio pari a

\begin{equation} g=5.88\ giorni,\qquad\hbox{con una incertezza}\qquad \Delta g=1.88\ giorni.\tag{15}\label{eq15} \end{equation}

In corrispondenza di questo valore otteniamo l'andamento dell'indice \(R_t\) di figura (12) nella quale riportiamo pure i valori di \(R_t\) associati alle stime giornaliere di \(\lambda_1\) presenti in fig. 11. Come nel caso della media mobile dei nuovi positivi e allo scopo di ottenere un accordo visivo soddisfacente con i valori giornalieri, la curva rappresentativa dell'indice di riproduzione è stata anticipata di 7 giorni rispetto alla data di calcolo effettivo.

Andamento dell'indice di riproduzione Rt
Fig. 12. Andamento dell'indice di riproduzione \(R_t\).

Volendo considerare la non trascurabile incertezza \eqref{eq15} sul tempo di generazione e fissato l'intervallo di confidenza al 68%, eseguiamo il medesimo calcolo sostituendo a \(g\) i valori \(g-\Delta g=5.88-1.88\) e rispettivamente \(g+\Delta g=5.88+1.88\). Questo permette di individuare le fasce di confidenza la cui larghezza massima viene raggiunta in corrispondenza dei valori estremi di \(R_t\) mentre la loro larghezza tende ad annullarsi in intorni del valore unitario dell'indice (fig. 13).

Valori di Rt e fasce di attendibilità 68%
Fig. 13. Valori di \(R_t\) e fasce di attendibilità 68%.

Quale confronto conclusivo restringiamo la figura precedente alla seconda ondata pandemica considerando il periodo che va dal 19 settembre 2020 al 17 gennaio 2021. Riportiamo sia l'indice \(R_t\) che qui calcoliamo assieme ai valori corrispondenti e pubblicati nella pagina web Tasso di contagio \(R_t\) del sito CovidStat più volte citato. Seppure nel sito il calcolo dell'indice di contagio faccia uso di un approccio più generale la figura 14 mostra comunque un buon accordo essendo pari a 0,05 lo scarto massimo.

Confronto dell'indice Rt con i valori presenti in CovidStat
Fig. 14. Confronto dell'indice \(R_t\) con i valori presenti in CovidStat.

Ringraziamenti

Rivolgiamo un vivo ringraziamento al Gruppo di Lavoro CovidStat promosso in seno all'Istituto Nazionale di Fisica Nucleare per l'ottimo livello informativo raggiunto sulla diffusione dei dati della pandemia di Covid-19 in Italia. Il lavoro di tale gruppo è non solo utile alla comunità scientifica ma risulta di grande interesse anche ad un più ampio spettro di fruitori e, per le numerose occasioni didattiche ivi presenti, a tutti quegli studenti interessati ad approfondire in ambiti multidisciplinari le proprie conoscenze e competenze.

Appendice

Importare dati CSV in un foglio Excel

Avviare Excel, quindi selezionare Dati seguito da Carica dati esterni, da testo. Selezionata la virgola come separatore di campo, si ottengono tutti i dati disposti correttamente in forma tabellare. Selezionata la colonna dei nuovi_positivi con un copia e incolla la si riporta in un secondo foglio nel quale procedere all'elaborazione.

Importare dati CSV in un foglio LibreOffice Calc

Avviare Calc, quindi File/Apri e selezionato tra i tipi di file quello corrispondente a Fogli elettronici, nella finestra che si apre e tra le opzioni di separazione, selezionare la virgola controllando la corretta suddivisione in colonne nella parte inferiore della medesima finestra. Infine OK e selezionata la colonna dei nuovi_positivi con un copia e incolla la si riporta in un secondo foglio nel quale procedere all'elaborazione.

Le piccole frecce con sfondo azzurro accanto alla didascalia di alcune immagini aprono una pagina con la medesima figura ingrandita.

Materiali didattici

I supporti didattici collegati a questa pagina sono costituiti da

  • un foglio di calcolo Excel, nuoviPositivi-IT.xlsx dove seguendo il diagramma di flusso di figura 8, riportiamo formule, calcoli e rappresentazioni grafiche associati a tutti i parametri del modello.

Seguono quindi cinque notebook di Jupiter ciascuno dei quali ha per obiettivi,

  • prelevaRinomina.ipynb preleva dal repository del Dipartimento della Protezione Civile DPC Covid-19 i file di dati necessari ai successivi quattro notebook e li rinomina coerentemente alla data di prelievo;
  • medieMobiliNuoviPositivi.ipynb sulla base dei due file prelevati con il precedente notebook, rappresenta i dati giornalieri dei nuovi positivi ed esegue, di questi, due medie mobili di 7 e 14 giorni;
  • tassoCrescitaLambda.ipynb calcola il tasso di crescita λ tramite una regressione lineare, esegue una sua media mobile e ne rappresenta l'andamento;
  • indiceRiproduzioneRt.ipynb calcola l'indice di riproduzione sulla base di una stima del tempo di generazione e del parametro λ, esegue una sua media mobile e rappresenta i risultati assieme, eventualmente, agli intervalli di confidenza;
  • indiceRt-30giorni.ipynb calcola l'indice di riproduzione e gli intervalli di confidenza degli ultimi 30 giorni e li riporta in forma numerica a video.

Tutti questi materiali sono distribuiti nel file compresso

Preleva NBindiceRt.zip (305 kB)

Per una consultazione preventiva di questi notebook o per una sperimentazione online in forme interattive, consultare in GitHub la pagina

Referenze

[1] G. Bonifazi, L. Lista et al., A simplified estimate of the Effective Reproduction Rt Number using its relation with the doubling time and application to Italian COVID-19 data, The European Physical Journal Plus (2021)

[2] Repository Covid-19 del Dipartimento della Protezione Civile, DPC, Covid-19

[3] L. Roi, Semplici modelli matematici di diffusione dei virus, MATEpristem (2020)

[4] Regressione lineare, https://it.wikipedia.org/wiki/Regressione_lineare

[5] D. Cereda et al., The early phase of the COVID-19 outbreak in Lombardy, Italy, arXiv:2003.09320 (2020)

Siti istituzionali e non, consultati: