COVID-19. Modelli matematici

(Ultimo aggiornamento: 21 marzo 2020)

 

Non sono un epidemiologo, sono solo uno statistico piccolo piccolo con un minimo di abitudine a dati e modelli. I modelli vengono normalmente usati per formalizzare l'andamento passato di un fenomeno, al fine di poterne prevedere l'andamento futuro. Non sono un epidemiologo, ma ascolto gli epidemiologi quando dicono che in epidemiologia i modelli matematici non servono per fare previsioni, ma per valutare se le misure adottate per rallentare l'epidemia stanno funzionando (intervista di Per Luigi Lopalco a Wired). Vediamo quindi i modelli matematici più usati per fenomeni che crescono nel tempo, come possono essere concretamente elaborati, cosa ci possono dire.

NB: ne parlo perché possono dare una prima rozza idea di quello che succede, e anche perché ne stanno girando veramente tanti sui social e può forse essere utile spiegare in modo semplice come funzionano. Preferisco approcci diversi, a cui vorrei però arrivare con gradualità.

Per le elaborazioni userò R, uno strumento molto diffuso tra gli statistici ma scaricabile liberamente da chiunque (è disponibile un'ampia documentazione, anche in italiano). Trovate nel file modellimat.R il codice utilizzato per questa pagina.

Il modello esponenziale

È uno dei modelli più semplici. In un modello esponenziale un numero \(N\) cresce sempre più rapidamente perché la velocità della crescita è proporzionale a \(N\): \[\dfrac{dN}{dt}=rN\] dove \(r\) è un tasso di crescita. La soluzione dell'equazione è: \[N_t=N_0e^{rt}\] dove \(N_0\) è il valore di \(N\) quando la crescita inizia e \(N_t\) è il suo valore al tempo \(t\). Si tratta di trovare i parametri \(N_0\) e \(r\).

È relativamente semplice lavorare con modelli lineari, e si intende lineari nei parametri; ad esempio, il modello \(y=\beta_0+\beta_1x_1+\beta_2x^2\) non è lineare in sé, perché interviene il quadrato di \(x\), ma i parametri \(\beta_i\) non sono altro che fattori dei diversi termini. Nel modello esponenziale, invece, il parametro \(r\) è un esponente e in tali modelli, detti non lineari, è necessario usare metodi iterativi: si scelgono valori iniziali dei parametri e poi un opportuno algoritmo li cambia fino a trovare (se ci riesce) i valori che si adattano meglio ai dati.

La scelta dei valori iniziali è critica (partendo da valori inadeguati l'algoritmo "non converge", non riesce a trovare una soluzione) e tutt'altro che banale; in questo caso possiamo presumere che \(N_0\) sia noto, ma indovinare \(r\) non è così facile. Una delle tecniche più semplici si basa sul confronto visivo tra i dati e i valori che si otterrebbero sulla base di alcuni possibili parametri (il blog di Andrea Onofri propone un altro utile approccio). Proviamo a generare dati simulati coerenti con un modello esponenziale, usando \(N_0=1\):

> N0 <- 1
> r <- 0.7
> tt <- seq(0, 10, length = 50) 
> N <- N0 * exp(r * tt) + rnorm(50, 0, 50)
> df <- data.frame(tt, N)
      
Definiamo ora la formula del modello esponenziale e usiamo la funzione preview() del package nlstools per compiere un primo tentativo con r=0.5:

> library(nlstools)
> formulaExp <- as.formula(N ~ N0 * exp(r * tt))
> preview(formulaExp, data=df, start=list(N0=1, r=0.5))
      
Preview Preview Preview

 

Il risultato (grafico a sinistra) non è esaltante: i simboli rossi, determinati dai valori iniziali scelti, mostrano un andamento molto più blando rispetto ai simboli neri, che rappresentano i dati "veri". Proviamo con r=1, ma otteniamo l'effetto contrario (grafico al centro). Va meglio con una via di mezzo, r=0.75 (grafico a destra) e decidiamo quindi che i nostri valori iniziali sono \(N_0=1,r=0.75\).

Procediamo ora con la stima del modello, usando la funzione nls() di R:

> mod.exp <- nls(formulaExp, data=df, start=list(N0=1, r=0.75))
> mod.exp
[...]
    N0      r 
0.7161 0.7369 
[...]
      
Come si vede, abbiamo ottenuto valori dei parametri ragionevolmente vicini a quelli "veri", cioè a quelli usati per generare i dati. Si può apprezzare il risultato confrontando i dati (i punti) con un modello i cui parametri siano quelli trovati (curva rossa):

Ottenere valori esatti sarebbe impossibile, perché i dati non seguono pedissequamente il modello matematico e presentano comunque un andamento in parte erratico. Obiettivo del modello è spiegare quanto più possibile la variabilità dei dati.

Procedendo in modo analogo, si ottengono \(\hat{N}_0=16.45\) e \(\hat{r}=0.19\) per le vittime COVID-19 in Italia dal 21 febbraio al 19 marzo (\(N_0\) e \(r\) sono i parametri del modello, \(\hat{N}_0\) e \(\hat{r}\) indicano le loro stime):

Anche qui l'adattamente non è perfetto, ma finché i dati non cambiano tanto da rendere evidente l'opportunità di scegliere un modello diverso, si deve concludere che le misure di rallentamento non hanno ancora avuto un effetto appezzabile. Incombe così una caratteristica tipica dei modelli esponenziali: crescono sempre più rapidamente all'infinito.

Previsioni? Meglio cercare di fare qualcosa perché il modello esponenziale non funzioni più così bene...

Il modello logistico

Molti fenomeni crescono fino a un massimo, che viene detto capacità portante e può essere indicato con \(K\). Ad esempio, la statura di un individuo cresce da quando è neonato in poi, ma alla fine si attesta a quella che sarà la statura dell'individuo adulto. Per tenere conto di tali fenomeni, si può aggiungere un fattore \(1-N/K\) al modello esponenziale, un fattore che è prossimo a \(1\) quando \(N\) è \(0\) o comunque nettamente più piccolo di \(K\) (all'inizio si osserva una crescita esponenziale), ma poi diventa sempre più piccolo, fino ad annullarsi quando \(N=K\): \[\dfrac{dN}{dt}=rN\left(1-\dfrac{N}{K}\right)\] La soluzione dell'equazione è: \[N_t=\frac{KN_0}{(K-N_0)e^{-rt}+N_0}\]

Possiamo provare come prima a generare dati simulati e a stimare il modello (salto la fase di ricerca dei valori iniziali con preview, che non presenta problemi):

> K <- 1000
> N0 <- 1
> r <- 0.7
> tt <- seq(0, 30, length=50)
> N <- K * N0 / ((K - N0)*exp(-r * tt) + N0) + rnorm(50, 0, 50)
> df <- data.frame(tt, N)
> formulaLog <- as.formula(N ~ K * N0 / ((K - N0) * exp(-r * tt) + N0))
> mod.log <- nls(formulaLog, data=df, start=list(N0=1, K=1000, r=0.75))
> mod.log
[...]
       N0         K         r 
   1.5756 1013.1380    0.6576 
[...]
      
Il risultato è sicuramente soddisfacente:

Possiamo notare le caratteristiche del modello logistico osservando la curva rossa nel grafico. La curva cresce prima a un ritmo sempre più veloce, come un'esponenziale, ma poi, superato un punto di flesso che si trova più o meno dove \(t=10\), comincia a crescere sempre più lentamente, fino a raggiungere un asintoto orizzontale, un livello massimo che non può superare (il parametro \(K\) del modello, in questo caso pari a circa \(1013\)).

La curva è simmetrica rispetto al punto di flesso: tra \(t=10\) e \(t=20\) cresce sempre più lentamente tanto quanto cresce sempre più rapidamente tra \(t=0\) e \(t=10\).

Si tratta di una simmetria che a volte si osserva nel fenomeno oggetto di studio, ma a volte no. Sappiamo che nella provincia cinese di Hubei la COVID-19 sembra ormai sotto controllo e appare ragionevole applicare un modello logistico. Si ottengono \(\hat{N}_0=57.88,\hat{K}=3090.72\) e \(\hat{r}=0.16\):

Sembra accettabile, ma si notano due stranezze: l'accelerazione del ritmo di crescita dei dati reali sembra maggiore all'inizio (parte più in basso, il vero \(N_0\) è \(17\), poi raggiunge e supera la curva logistica), mentre la decelerazione sembra minore alla fine, tanto che i dati reali sembrano tendere a un asintoto più alto. Si può fare di meglio?

Il modello di Gompertz

In realtà vi sono più modelli di Gompertz, che differiscono sia per il numero che per il significato dei parametri. Vi sono anche altri modelli che condividono la caratteristica essenziale di un modello Gompertz: il ritmo della prima fase, fino al punto di flesso, è diverso dal ritmo della seconda fase, dopo il punto di flesso. Ad esempio, si può avere prima una rapida accelerazione, poi una decelerazione più lenta.

Un modello Gompertz "classico" è: \[N_t=Ke^{-e^{-r(t-T_f)}}\equiv K\exp(-\exp(-r(t-T_f)))\] dove \(K\) è ancora il valore massimo che \(N\) può assumere, \(r\) è un tasso di crescita e \(T_f\) è il tempo in cui si verifica il punto di flesso.

Si può procedere come già visto sia per la generazione di dati simulati (un utile esercizio), sia per la stima del modello su dati reali, usando:

> formulaGom <- as.formula(N ~ K * exp(-exp(-r * (tt - Tf))))
        

Si ottengono \(\hat{K}=3246, \hat{T}_f=21.43\) (tra l'11 e il 12 febbraio) e \(\hat{r}=0.097\):

Va meglio rispetto al modello logistico. La curva parte più vicina ai dati reali (si può calcolare \(\hat{N}_0=3246\exp(-\exp(-0.097(-21.43)=1\)) e, soprattutto, l'ultimo suo tratto li segue molto bene.

Confronto Italia-Cina

La vista di quello che è successo nella provincia di Hubei è incoraggiante: prima una crescita esponenziale, poi una lenta decelerazione e una curva che si avvicina sempre più a un massimo. Sta succedendo lo stesso anche da noi?

Il grafico proposto sopra, che mostra un chiaro andamento esponenziale per i deceduti in Italia, non sembra consentire una risposta affermativa. Possiamo forse sperare di non essere comunque troppo lontani dall'andamento di Hubei?

I dati relativi alla provincia di Hubei (59 milioni di abitanti) iniziano dal 22 gennaio, giorno in cui si contavano 17 decessi. In Italia (60 milioni di abitanti) si sono contati 17 decessi il 27 febbraio. Ho sovrapposto i dati di Hubei e italiani, facendo partire i primi dal 22 gennaio, i secondi dal 27 febbraio. Il risultato:

La partenza è simile, ma 10 giorni dopo il 27 febbraio la situazione in Italia è peggiorata notevolmente, tanto da costringere ad abbandonare qualsiasi ipotesi ottimistica, anche moderatamente ottimistica.

Siamo veramente messi male, e sarà opportuno chiedersi perché.

Letture utili

Florent Baty, Christian Ritz, Sandrine Charles, Martin Brutsche, Jean-Pierre Flandrois e Marie-Laure Delignette-Muller (2015), "A Toolbox for Nonlinear Regression in R: The Package nlstools", The Journal of Statistical Software, 66(5), http://dx.doi.org/10.18637/jss.v066.i05

Andrea Onofri (2019), "Some useful equations for nonlinear regression in R", https://www.statforbiology.com/articles/usefulequations/

Andrea Onofri (2020), "A collection of self-starters for nonlinear regression in R", https://www.statforbiology.com/2020/stat_nls_usefulfunctions/

Kathleen M. C. Tjørve e Even Tjørve (2017), "The use of Gompertz models in growth analyses, and new Gompertz-model approach: An addition to the Unified-Richards family", PLoS ONE, 12(6), https://doi.org/10.1371/journal.pone.0178691

A. Tsoularis e J. Wallace (2001), "Analysis of logistic growth models", Mathematical Biosciences, 179(1), 21-55, https://doi.org/10.1016/S0025-5564(02)00096-2