Introduction à la modélisation en épidémiologie : COVID-19 - 1ère flambée - croissance exponentielle des décès

Frédéric Hamelin, le 22 Mars 2021.

---

Dans ce premier TD, nous allons modéliser la croissance initiale de l'épidémie de COVID19 en France sur la base des données de décès uniquement.

Modèle épidémiologique

Les variables du modèle sont :

• $S(t)$ : le nombre de personnes sensibles au temps $t$,
• $I(t)$ : le nombre de personnes infectées et infectieuses,
• $R(t)$ : le nombre de personnes rétablies ou guéries,
• $D(t)$ : le nombre de personnes décédées.

La taille de la population est $N=S+I+R$.

Les paramètres du modèles sont

• $\beta$ : le taux de transmission de la maladie,
• $\gamma$ : le taux de guérison,
• $\alpha$ : le taux de mortalité due à la maladie.

Ce sont des "taux" par unité de temps.

Le modèle s'écrit : $$\begin{eqnarray} \frac{\mathrm{d}S}{\mathrm{d}t}&=&-\beta \frac{I}{N}S\,,\\ \frac{\mathrm{d}I}{\mathrm{d}t}&=&\beta \frac{I}{N}S - (\alpha+\gamma)I\,,\\ \frac{\mathrm{d}R}{\mathrm{d}t}&=&\gamma I\,,\\ \frac{\mathrm{d}D}{\mathrm{d}t}&=&\alpha I\,. \end{eqnarray}$$ Comme nous nous intéressons à la phase initiale de l'épidémie, nous faisons l'hypothèse que pour tout $t$ dans cette phase initiale, $I,R,D \ll S$, de sorte que $S\approx N$ (la taille de la population est approximativement constante durant cette période). Le modèle simplifié s'écrit : $$\begin{eqnarray} \frac{\mathrm{d}I}{\mathrm{d}t}&=&(\beta - \alpha-\gamma)I\,,\\ \frac{\mathrm{d}R}{\mathrm{d}t}&=&\gamma I\,,\\ \frac{\mathrm{d}D}{\mathrm{d}t}&=&\alpha I\,. \end{eqnarray}$$ Ainsi, $$I(t) = I(0)\exp\left(\left(\beta-\alpha-\gamma\right)t\right),$$ et $$\begin{eqnarray} D(t)&=&D(0)+\int_0^t\alpha I(\tau) \mathrm{d}\tau\,,\\ &=& D(0)+\alpha\int_0^t I(0)\exp\left(\left(\beta-\alpha-\gamma\right)\tau\right) \mathrm{d}\tau\,,\\ &=& D(0) + \frac{\alpha I(0)}{\beta-\alpha-\gamma}\left[\exp\left(\left(\beta-\alpha-\gamma\right)t\right)-1\right]\,. \end{eqnarray}$$ Soient $$r=\beta-\alpha-\gamma\,,\quad\mbox{et}\quad q=\alpha I(0)\,.$$ Le modèle s'écrit enfin $$D(t)=D(0)+\frac{q}{r}\left(\exp(rt)-1\right)\,.$$

Selon le modèle, le nombre de décès cumulés croît de façon exponentielle en fonction du temps.

Traitement des données

Commençons par nettoyer le plan de travail :

rm(list=ls())  # Efface les variables créées lors des exécutions précédentes
graphics.off() # Ferme les fenêtres ouvertes lors des exécutions précédentes

Nous aurons besoin de la librairie readr :

library(readr) # Nécessaire pour utiliser la fonction read_csv ci-dessous

Importons les Données relatives à l’épidémie de COVID-19 en France : vue d’ensemble :

data <- read_csv(url("https://www.data.gouv.fr/fr/datasets/r/d3a98a30-893f-47f7-96c5-2f4bcaaa0d71"))

Remplaçons les données manquantes par des zéros pour faciliter le traitement des données ensuite :

data[is.na(data)]=0

Nous sommes intéressés par les décès totaux (hôpital + EHPAD) :

Dobs=data$total_deces_hopital+data$total_deces_ehpad # décès cumulés

Affichons les décès en fonction du temps :

L=length(data$date) # longueur de la série de données
Tobs=seq(0,L-1,by=1) # Vecteur temps
plot(Tobs,Dobs,xlab="Temps écoulé depuis le 2 mars 2020 (en jours)",ylab="Nombre de personnes décédées",col="blue")

Restreignons la période étudiée au mois de mars 2020 (1ère flambée) :

Tmax=30 # Fin de la fenêtre temporelle (le 31 mars 2020)
w=seq(0,Tmax,by=1) # Fenêtre temporelle considérée

Tobs=Tobs[w]  # Troncation des données
Dobs=Dobs[w]

Affichons les décès cumulés sur cette période :

plot(Tobs,Dobs,xlab="Temps écoulé depuis le 2 mars 2020 (en jours)",ylab="Nombre de personnes décédées",col="blue")

La croissance est-elle exponentielle ? Si oui, le logarithme du nombre de décès devrait croître de façon linéaire en fonction du temps :

plot(Tobs,log(Dobs),xlab="Temps écoulé depuis le 2 mars 2020 (en jours)",ylab="Log(Nombre de décès cumulés)",col="blue")
reg=lm(log(Dobs)~Tobs) # Régression linéaire
intercept=reg$coefficients[[1]]
slope=reg$coefficients[[2]]
abline(reg)

Nous voyons que la croissance du logarithme est relativement linéaire sur la période considérée. L'hypothèse d'une croissance exponentielle du nombre de décès cumulés paraît raisonnable sur la période considérée (mars 2020).

Ajustement du modèle aux données

Au 2 mars ($t=0$), le nombre de décès cumulés est $D(0)=3$:

D0=Dobs[1] # Donnée initiale

Comme $D(0)$ est petit par rapport au nombre de décès cumulés (ce nombre dépasse 3000 à la fin du mois de mars), on peut considérer la relation suivante en première et grossière approximation : $$D(t)=D(0)+\frac{q}{r}\left(\exp(rt)-1\right))\approx \frac{q}{r}\exp(rt)\,.$$ Cela donne $$\log(D(t))\approx\log\left(\frac{q}{r}\right)+rt\,,$$ ce qui permet d'obtenir des estimations initiales des valeurs des paramètres $r$ et $q$ d'après la régression linéaire réalisée précédemment :

#estimations grossières initiales
r=slope             # r = beta - alpha - gamma
q=r*exp(intercept)  # q = alpha*I_0

Calculons la courbe issue du modèle pour ces valeurs de paramètres :

Tmodel=Tobs # on calcule la solution du modèle pour les dates d'observation
Dmodel=D0 + (q/r)*(exp(r*Tmodel)-1)

Comparons le modèle aux données :

plot(Tobs,Dobs,xlab="Temps écoulé depuis le 2 mars 2020 (en jours)",ylab="Nombre de décès cumulés",col="blue")
lines(Tmodel,Dmodel,col="red")

L'ajustement pourrait être meilleur. Appelons $\theta=\{r,q\}$ le vecteur des paramètres à optimiser :

theta=c(r,q) # Vecteur des paramètres à optimiser

Méthode des moindres carrés

Recherchons les valeurs des paramètres $r$ et $q$ qui minimisent la somme des carrés des écarts entre modèle et données (SCE) : $$\mbox{SCE} = \sum_{t=0}^T \left(D_{\small\mbox{model}}(t)-D_{\small\mbox{data}}(t)\right)^2\,.$$ Créons une fonction qui prend en entrée les paramètres à optimiser et renvoie en sortie la quantité à minimiser (SCE) :

SCE_Covid=function(theta){
  r=theta[1]    
  q=theta[2]

  Dmodel = D0 + (q/r)*(exp(r*Tmodel)-1) 

  ecarts = Dobs - Dmodel
  SCE = sum(ecarts^2)

  return(SCE)
}

Cherchons les valeurs des paramètres qui minimisent SCE :

opt=optim(theta,SCE_Covid) # La fonction optim minimise par défaut
r=opt$par[1]
q=opt$par[2]

Affichons la solution optimale au sens des moindres carrés :

Dmodel = D0 + (q/r)*(exp(r*Tmodel)-1) 
plot(Tobs,Dobs,xlab="Temps écoulé depuis le 2 mars 2020 (en jours)",ylab="Nombre de décès cumulés",col="blue")
lines(Tmodel,Dmodel,col="red")

L'ajustement est meilleur que celui obtenu en première et grossière approximation. Affichons les valeurs des paramètres obtenues :

print(r)
print(q)

On trouve $r \approx 0.17$ et $q \approx 4.11$.

Méthode du maximum de vraisemblance

Une méthode alternative consiste à maximiser la probabilité des observations sachant les paramètres. C'est ce qu'on appelle la vraisemblance (likelihood en anglais).

Nous faisons l'hypothèse que le nombre de décès cumulés observé au jour $t$ est tiré dans une loi de Poisson de moyenne $D(t)$ telle que donné par le modèle : $$D_{\small\mbox{data}}(t)\sim \mbox{Poisson}\left(D_{\small\mbox{model}}(t)\right)$$ En admettant que les observations sont indépendantes conditionnellement au modèle, la vraisemblance s'écrit : $$\mbox{Like}(\theta)=\prod_{t=0}^{T} \mbox{Prob}(\mbox{observe }D_{\small\mbox{data}}(t)|D_{\small\mbox{model}}(t))\,.$$ Rechercher les valeurs de $\theta$ qui maximisent la vraisemblance est équivalent à rechercher les valeurs de $\theta$ qui maximisent la log-vraisemblance, qui transforme le produit en somme : $$\mbox{logLike}(\theta)=\log(\mbox{Like}(\theta))=\sum_{t=0}^T \log\left (\mbox{Prob}(\mbox{observe }D_{\small\mbox{data}}(t)|D_{\small\mbox{model}}(t))\right )\,.$$ Créons une fonction qui prend en entrée les paramètres à optimiser et renvoie en sortie la quantité à minimiser (-logLike) :

LL_Covid=function(theta){ #Log-likelihood

  r=theta[1]
  q=theta[2]

  Dmodel = D0 + (q/r)*(exp(r*Tmodel)-1)

  probas = dpois(Dobs,Dmodel)
  LL=sum(log(probas))

  return(-LL)  #on renvoie l'opposé pour minimiser
}

Cherchons les valeurs des paramètres qui maximisent logLike :

opt=optim(theta,LL_Covid)
r=opt$par[1]
q=opt$par[2]

Affichons la solution optimale au sens du maximum de vraisemblance :

Dmodel = D0 + (q/r)*(exp(r*Tmodel)-1) 
plot(Tobs,Dobs,xlab="Temps écoulé depuis le 2 mars 2020 (en jours)",ylab="Nombre de décès cumulés",col="blue")
lines(Tmodel,Dmodel,col="red")

L'ajustement est comparable à celui obtenu via les moindres carrés. Affichons les valeurs des paramètres obtenues :

print(r)
print(q)

On trouve $r \approx 0.19$ et $q \approx 2.47$, ce qui diffère un peu des valeurs obtenues via les moindres carrés.

Interprétation des résultats

L'ajustement du modèle aux données a permis

• de confirmer le caractère exponentiel de la dynamique épidémique durant la première flambée, telle que mesurée par les décès cumulés en mars 2020,
• d'estimer les paramètres $r=\beta-\alpha-\gamma$ et $q=\alpha I(0)$.

En admettant que le taux de mortalité due à la maladie $(\alpha)$ est très petit devant le taux de guérison $(\gamma)$, nous avons l'approximation $r\approx \beta - \gamma$. Cela permet d'estimer la reproductivité du virus ($\mathcal{R}_0$), définie comme le nombre d'infections secondaires générées par un individu infecté dans une population initialement naïve : $$\mathcal{R}_0 = \frac{\beta}{\alpha+\gamma}\approx \frac{\beta}{\gamma} \approx 1 + \frac{r}{\gamma}\,.$$ Admettons que le temps de guérison est de 10 jours en moyenne, soit $\gamma=0.1$ par jour. On obtient :

• $\mathcal{R}_0\approx 2.74$ via la méthode des moindres carrés,
• $\mathcal{R}_0\approx 2.98$ via la méthode du maximum de vraisemblance.

Ces estimations sont cohérentes avec celles obtenues par différentes études (ex. Roques et al 2020).

Le taux de létalité du virus (la probabilité moyenne de mourir lorsqu'on est infecté) est $$p=\frac{\alpha}{\alpha+\gamma}\,.$$ Supposons $p=0.8\%$ (Roques et al 2020), ce qui donne $\alpha=p\gamma/(1-p)\approx p\gamma = .0008$. Nous obtenons

• $I(0)=q/\alpha\approx 5145$ via les moindres carrés,
• $I(0)=q/\alpha\approx 3092$ via le maximum de vraisemblance.

Le modèle permet d'estimer un ordre de grandeur du nombre total de personnes infectées et infectieuses au 3 mars, alors que le nombre total de cas confirmés à cette date s'élevait à 191 d'après les Données relatives à l’épidémie de COVID-19 en France : vue d’ensemble. D'après le modèle, on aurait approximativement $4000/200=20$ personnes infectées pour 1 cas confirmé au 2 mars 2020.