Introduction à la modélisation en épidémiologie : COVID-19 - 1ère flambée - estimation du taux de létalité

Marie Joigneau, Sabine Lobligeois, Louise Belamy, Maimouna Diarra (M1 Agro Rennes), le 14 mars 2021.

Edité par Frédéric Hamelin et Marine Dorand (M1 MODE), le 23 avril 2021.

On se base sur l’étude de l’article suivant : Roques et al. (2020) Using Early Data to Estimate the Actual Infection Fatality Ratio from COVID-19 in France. Le but de notre étude est de reproduire les résultats de l’article avec les données de morts actualisées et publiées par Santé Publique France.

Modèle épidémiologique

Dans ce TD, nous allons modéliser la croissance initiale de l'épidémie de COVID19 en France sur la base des données de tests dans un premier temps. Les données concernant les décès seront utilisées dans un second temps pour estimer le taux de létalité de la maladie au début de l'année 2020 en France.

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.

plot0

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}$ Pour simplifier l'étude, nous faisons l'hypothèse que le taux de mortalité due à la maladie est très faible et négligeable d'un point de vue épidémiologique devant le taux de guérison $(\alpha\ll \gamma)$ . Cela permet de décomposer le modèle en deux parties, la première étant indépendante des décès : $\begin{eqnarray} \frac{\mathrm{d}S}{\mathrm{d}t}&=&-\beta \frac{I}{N}S\,,\\ \frac{\mathrm{d}I}{\mathrm{d}t}&\approx&\beta \frac{I}{N}S -\gamma I\,,\\ \frac{\mathrm{d}R}{\mathrm{d}t}&=&\gamma I\,. \end{eqnarray}$ Elle sera ajustée aux données de tests. La seconde partie du modèle sera ajustée aux données de décès : $\begin{eqnarray} \frac{\mathrm{d}D}{\mathrm{d}t}&=&\alpha I\,. \end{eqnarray}$

Modèle d'observation

L'ajustement du modèle se fera par maximisation de la vraisemblance. Cette méthode consiste à maximiser la probabilité des observations sachant les paramètres. Le modèle d'observation définit la probabilité qu'un test soit positif au temps $t$ d'après le modèle épidémiologique : $p(t)=\frac{\sigma I(t)}{I(t)+\kappa S(t)}\,,$ où $\sigma $ est la sensibilité des tests PCR (supposée connue) et $\kappa$ est probabilité relative pour les individus sensibles de se faire tester comparativement aux infectés (à estimer). Aussi, le nombre de tests positifs le jour $t$ peut être modélisé comme un tirage dans une loi Binomiale de paramètres $p(t)$ et $n(t)$ (le nombre de tests réalisés le jour $t$ ) : le nombre de nouveaux cas observés le jour $t$ est $c(t)\sim \mbox{Binomial}(n(t),p(t))$ En admettant que les observations sont indépendantes conditionnellement au modèle, la vraisemblance (likelihood en anglais) s'écrit : $\mbox{Like}(\theta)=\prod_{t=0}^{T} \mbox{Prob}(\mbox{observe }c(t)|\mbox{ model predicts }p(t))\,.$ Soit $\theta = \{ \beta, \kappa \}$ le vecteur des paramètres à estimer. 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 }c(t)|\mbox{ model predicts } p(t) )\right)\,.$ Nous travaillerons avec la log-vraisemblance.

Traitement des données

Commençons par nettoyer l'environnement 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

On charge les librairies suivantes :

library(deSolve) # pour résoudre le système d'équations différentielles (modèle SIRD)
library(R.matlab) # pour lire les fichiers de données au format matlab
library(downloader) # pour faciliter le téléchargement des données
library(zoo) # pour la moyenne mobile
library(readr) # pour utiliser la fonction read_csv

On télécharge l'archive contenant les données utilisées par Roques et al 2020 :

download("https://www.mdpi.com/2079-7737/9/5/97/s1", dest="suppMat.zip", mode="wb")

On décompresse l'archive :

unzip("suppMat.zip", exdir = "suppMat")

On récupère les données relatives aux tests (du 1er janvier au 17 mars) :

dataT = readMat("./suppMat/data_smooth5.mat")

TP=as.integer(dataT$DATA[1,]) # nombre de tests positifs cumulés
TT=as.integer(dataT$DATA[2,]) # nombre total de tests cumulés

LT=length(TP) # Longueur de la série temporelle

On récupère les données relatives aux décès (à partir du 2 mars) :

dataD <- read_csv(url("https://www.data.gouv.fr/fr/datasets/r/d3a98a30-893f-47f7-96c5-2f4bcaaa0d71"))
dataD[is.na(dataD)]=0 # on remplace les données manquantes par des zéros avant de sommer
DC=dataD$total_deces_hopital+dataD$total_deces_ehpad # nombre de décès cumulés

On crée un vecteur des décès cumulés du 1er janvier au 17 mars :

# Les données de décès sont disponible à partir du 2 mars
n=31+29+1 # nombre de jours en janvier + février 2020 + 1 mars
zeros=rep(0,n) # nombre de zéros à ajouter pour janvier/février
m=17-1 # la fenêtre de l'étude s'arrête au 17 mars
DC = c(zeros,DC[1:m]) # les décès cumulés du 1er janvier au 17 mars

On crée un vecteur pour le taux de positivité cumulé :

PP=TP/TT # taux de positivité cumulé

On crée un vecteur temps :

temps=seq.Date(from=as.Date("2020-01-01"), to=as.Date("2020-03-17"), by="days")

On affiche les tests cumulés en fonction du temps :

plot(temps,TT,xlab="Date (Annee 2020)",ylab="Nombre de tests",col="black")
points(temps,TP,col="red")
legend("topleft",c("Nombre de tests positifs","Nombre de tests realises"), fill=c("red","black"))
title("Resultats cumules des tests realises de covid 19")

On affiche la proportion de tests positifs en fonction du temps :

plot(temps,PP,xlab="Date (Annee 2020)",ylab="Proportion de tests positifs",col="purple")
title("Proportion de tests positifs en fonction du temps")

On affiche les décès cumulés en fonction du temps :

plot(temps,DC,xlab="Date (Annee 2020)",ylab="Nombre de morts",col="orange")
title("Morts cumules du covid 19")

On retrouve les données non-cumulées en utilisant la fonction diff :

tt=diff(TT) # nombre quotidien de tests réalisés
tp=diff(TP) # nombre quotidien de tests positifs
dc=diff(DC) # nombre quotidien de décès
lt=length(tt) # longueur de la série temporelle (l=L-1)

On utilise donc un intervalle de temps qui commence le 2 janvier 2020 :

temps2=seq.Date(from=as.Date("2020-01-02"), to=as.Date("2020-03-17"), by="days")

On affiche les tests quotidiens en fonction du temps :

plot(temps2,tt,xlab="Date (Annee 2020)",ylab="Nombre de tests quotidiens",col="black")
points(temps2,tp,col="red")
legend("topleft",c("Nombre de tests positifs","Nombre de tests realises"), fill=c("red",'black'))
title("Resultats des tests quotidiens realises de covid 19")

On affiche le taux de positivité quotidien :

pp=tp/tt # Proportion quotidienne de tests positifs
plot(temps2,pp,xlab="Date (Annee 2020)",ylab="Proportion de tests quotidiens positifs",col="purple")
title("Proportion de tests positifs en fonction du temps")

On affiche les décès quotidiens en fonction du temps :

plot(temps2,dc,xlab="Date (Annee 2020)",ylab="Nombre de morts",col="orange")
title("Morts quotidiens du covid 19")

Ajustement du modèle aux données

Les paramètres du modèle sont :

N=67e6; # Taille totale de la pop française
gamma=1/10; # "taux" de guerison par jour (1/gamma est la periode infectieuse)
R_0=3; # Reproductivite de la maladie
beta=gamma*R_0; # R0=beta/alpha (estimation initiale à optimiser)

sigma=0.7;# Sensibilite des tests PCR
# Probabilité relative pour les sensibles de se faire tester comparativement aux infectés
kappa=5e-3# (estimation initiale à optimiser)

Les conditions initiales du modèle sont :

I0=1; # On part d'un individu infecté. Ceci est fixé pour la premiere flambée.
R0=0; # On suppose que personne n'est immunisé dans la population initialement.
S0=N-I0-R0; # Nombre initial de sensibles (tout le monde sauf le premier infecté)
X0=c(S0,I0); # Vecteur d'etat

Créons une fonction qui permet de simuler le modèle SIR :

SIR = function(t, X, P){

  beta=P[1] # beta est le premier et unique élément du vecteur des paramètres P
            # gamma est un paramètre dont la valeur est fixée de façon globale 
  S=X[1]    # S est le premier élément du vecteur d'état X
  I=X[2]    # I est le second élément du vecteur d'état X

  y=beta*S*I/N # y est le nombre de nouvelles infections par unité de temps
  dS = -y # dS/dt : variation du nombre d'individus sains
  dI = +y - gamma*I # dI/dt : variation du nombre d'individus infectés

  dX=c(dS,dI)
  return(list(dX)) # la liste dX contient dS et dI les derivées de S et I
}

Créons une fonction qui permet de calculer la vraisemblance associée à un vecteur de paramètres à optimiser $\theta = \{ \beta, \kappa \}$ .

logLike = function(theta){
    # theta est le vecteur des paramètres
    beta=theta[1]   # beta est le premier élément du vecteur theta
    kappa=theta[2]  # kappa est le second élément du vecteur theta

    X=ode(X0,t,SIR,beta)
    # la fonction ode renvoie les solutions du modele SIR avec beta pour paramètre
    # X0 est le vecteur des conditions initiales
    # t est le vecteur des dates d'observations
    # La fonction ode renvoie un matrice X de la forme :
    # 1ère colonne : t, 2ème colonne : S(t), 3ème colonne : I(t)

    # Proportion théorique de tests quotidiens positifs (information issue du modele)
    p=sigma*X[,3]/(X[,3]+kappa*X[,2])

    # On met des zeros avant t0 (en pratique 1e-12 pour eviter log(0) plus bas)
    p=c(rep(1e-12,t0-1),p) 

    # Probabilite d observer "tp" positifs sachant le nombre de tests "tt" 
    # et "p(t)" telle que donnée par le modèle (loi binomiale de parametre "pt")
    L=dbinom(tp, tt, p, log=TRUE) # Vecteur des vraisemblances (likelihood) 
    #Log = TRUE car on travaille sur la log-vraisemblance

    LL=sum(L) # Le log transforme le produit des probabilités en somme

    return(LL) # renvoie la Log-vraisemblance associée au jeu de paramètres theta
}

En plus des paramètres $\beta$ et $\kappa$ , nous souhaitons estimer la date initiale de l'épidémie, $t_0$ , telle que $I(t_0)=I_0=1$ . Comme $t_0$ est un entier, nous le traitons séparément via une boucle d'optimisation qui parcours tous les $t_0$ possibles du 2 janvier au 1er février.

V=-1e6 # Valeur initiale de la vraisemblance (à maximiser)

for (t0 in 1:31){ # Boucle d'optimisation sur les t0

  t=t0:lt # Vecteur temps : dates sur lesquelles on simule le modèle
  theta0=c(beta,kappa) # Vecteur des parametres a estimer (estimations initiales)

  opt=optim(theta0,logLike,control=list(fnscale=-1));
  # Part d'un point de depart theta0 et suit la plus forte pente (le gradient de la log-vraisemblance en fonction de theta) jusqu'à maximiser la log-vraisemblance (localement) 
  # Selon le point theta0 de départ que l'on donne le maximum local renvoyé peut varier

  if (opt$value>V){ #Si le t0 considéré augmente la vraisemblance,
      V=opt$value   #on sauvegarde la vraisemblance ainsi que t0, beta et gamma.
      t0opt=t0
      beta=opt$par[1]
      kappa=opt$par[2]
  } 
}

On affiche les paramètres optimaux obtenus :

print(t0opt)
print(kappa)
print(beta)

On trouve

$t_0 = 22$ jours (soit le 23 janvier),
$\kappa \approx .0004$ ,
$\beta \approx 0.3$ , ce qui 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 3\,,$

puisque $\gamma=0.1$ par jour.

Ces estimations sont cohérentes avec celles obtenues par Roques et al 2020.

t=t0opt:lt # nouveau vecteur temps 
X=ode(X0,t,SIR,beta) # simulation du modèle SIR

La proportion théorique de tests quotidiens positifs est :

p=sigma*X[,3]/(X[,3]+kappa*X[,2]) # p(t) = sigma I(t)/(I(t) + kappa S(t))

Les tests positifs cumulés observés sont :

Sigma1=cumsum(tp[t]) # on reprend les notations de l'article Roques et al (2020)

Tests positifs cumulés théoriques sont :

Sigma2=cumsum(tt[t]*p) # on reprend les notations de l'article Roques et al (2020)

On définit la première flambée entre le 22 janvier et le 17 mars 2020 :

temps3=seq.Date(from=as.Date("2020-01-23"), to=as.Date("2020-03-17"), by="days")

On affiche le nombre de tests positifs cumulés en fonction du temps :

plot(temps3,Sigma1,xlab="Date (Annee 2020)",ylab="Nombre de tests positifs cumules",col="black")
lines(temps3,Sigma2, col='red')
legend("topleft",c("D'apres les donnees","D'apres le modele"),fill=c('black','red'))
title("Comparaison des tests positifs cumules des donnees et du modele")

On retrouve le graphe de l’article Roques et al 2020 (Fig. 1) donc l’article est bien reproductible.

On affiche le nombre de cas infectieux $I(t)$ et infectés $(N-S(t))$ en fonction du temps :

plot(temps3,N-X[,2],col="red",type="l",xlab="Date (Annee 2020)",ylab="Nombre de personnes")
lines(temps3,X[,3],col="purple")
lines(temps3,Sigma2)
points(temps3,Sigma1)
legend("topleft",c("Nombre de personnes infectees","Nombre de personnes infectieuses","Nombre de personnes testees positives d'apres les donnees (points)","d'apres le modele (courbe)"),fill=c('red','purple','black','black'),cex=0.8)
title("Evolution de la pandemie en France lors de la premiere flambee")

On retrouve approximativement la Figure 2 de l’article Roques et al 2020.

On affiche le taux de positivité quotidien :

plot(temps3,pp[t],xlab="Date (Année 2020)",ylab="Proportion de tests quotidiens positifs")
lines(temps3,p)
legend("topleft",c("D'après les données (points)","D'après le modèle (courbe)"),)
title("Evolution de la proportion de tests positifs dans le temps")

Estimation du taux de létalité

Le taux de létalité moyen (IFR pour infection fatality ratio) est la proportion de morts parmi les personnes infectées ou la probabilité de mourir des suites de l'infection : $\mbox{IFR}(t)=\frac{\alpha(t)}{\alpha(t)+\gamma}\,.$ Nous ne supposons pas qu'il est constant au cours de l'épidémie.

On estime $\alpha(t)$ (le taux de mortalité due à la maladie) via la relation suivante : $\alpha(t)=\frac{1}{I(t)}\frac{\mathrm{d} D(t)}{\mathrm{d} t}\,,$ où $\mathrm{d}D(t)/\mathrm{d}t$ est le nombre de décès quotidien.

dc=dc[-(1:(t0opt-1))] # toncation des données de décès à partir de t0
alpha=dc/X[,3] # estimation de alpha(t) = [dD(t)/dt]/I(t)

IFR=alpha/(gamma+alpha) # taux de létalité (infection fatality ratio)

On affiche le taux de létalité (IFR) brut :

plot(temps3, IFR, type ="l", xlab="Date (Annee 2020)", ylab="IFR")

On réalise une moyenne de l’IFR sur les valeurs non nulles :

w=which(IFR!=0)
print(mean(IFR[w]))

On obtient un IFR de 0.67%. Cette valeur est comparable à celles estimées par Roques et al 2020.

On lisse la sortie en faisant une moyenne mobile sur 5 jours :

IFRlisse <- rollmean(IFR, k=5, fill=NA) #fill=NA remplit les valeurs manquantes en raison de la fonction rollmean par des NA

plot(temps3, IFRlisse, type ="l", xlab="Date (Année 2020)", ylab="IFR",ylim=c(0,0.01))
points(temps3, IFR )
title("Evolution de l'IFR (ratio de létalité) lors de la première flambée")
legend("topleft",c("D'après l'IFR brut (points)","D'après IFR lissé (courbe)"))

Nous obtenons des valeurs d’IFR du même ordre de grandeur (0.8%) et de tendance comparable à l'article de Roques et al 2020.

Keys	Action
`?`	Open this help
`n`	Next page
`p`	Previous page
`s`	Search