Introduction à la modélisation en épidémiologie - l'exemple du Covid-19

Support du TD d'écologie des pathogènes du 8 juin 2020.

Frédéric Hamelin, le 10 mai 2020.

Introduction à la modélisation en épidémiologie - l'exemple du Covid-19Avertissements Modèle épidémiologiqueTéléchargement et traitement des donnéesSimulations et ajustement du modèle aux donnéesRésolution numérique du modèle épidémiologiqueModèle d'observationMaximisation de la vraisemblanceRésultats"Prédictions" et discussion des limites du modèles RéférencesRemerciements

Avertissements

Mauvais R

L'auteur de ce TD est débutant en R. Les lignes de codes sont proposées à titre indicatif et correspondent à une utilisation maladroite de R. L'un des objectifs de cette version bêta du TD est de recueillir par mail vos propositions d'améliorations, tant sur le code que sur les figures, ainsi que sur le contenu du TD bien sûr.

Contre-modèle

Nous allons ajuster un modèle simple à des données réelles concernant l'épidémie de Covid-19 en France. Il s'agit d'un exercice à visée strictement pédagogique. La critique du modèle fait partie de l'exercice. Des études de modélisation pour la recherche appliquée au Covid-19 sont listées en fin de document.

Modèle épidémiologique

Dans l'épidémie de Covid-19, les individus infectés peuvent être symptomatiques ou asymptomatiques. Les données sur le nombre d'individus infectés dans la population sont peu fiables car elles dépendent des tests éventuellement réalisés. Les individus symptomatiques peuvent faire une forme sévère de la maladie et être hospitalisés. Les données concernant les individus hospitalisés sont fiables en principe.

Nous ferons un certain nombre d'hypothèses volontairement grossières pour simplifier l'étude en première approximation. Ces hypothèses incluent :

la population française est spatialement homogène au sein de l'hexagone,
les hommes et les femmes sont égaux vis-à-vis du virus,
la structure en âge de la population peut être ignorée pour modéliser l'épidémie,
les individus hospitalisés ne transmettent plus l'infection,
la guérison confère une immunité permanente.

Ces hypothèses sont plus que grossières mais permettent néanmoins de commencer à travailler. Techniquement, cela nous permettra de limiter le nombre de paramètres à estimer. Autrement, le modèle risquerait d'être sur-paramétré au regard des données et des connaissances dont on dispose sur le virus.

Nous définissons les variables

$S(t)$ $t$ ,
$I(t)$ $t$ ,
$H(t)$ $t$ ,
$R(t)$ $t$ ,

et les paramètres

$\beta$ $\times$ probabilité d'infection),
$\rho$ : taux de "guérison" par unité de temps (inverse du temps moyen avant guérison ou décès),
$\alpha$ : taux d'hospitalisation par unité de temps (inverse du temps moyen avant hospitalisation),
$\gamma$ : taux de guérison ou décès à l'hôpital (inverse du temps moyen avant sortie de l'hôpital).

Nous considérons le modèle compartimental suivant :

\begin{eqnarray*} \text{variation du nombre d'individus sensibles} &=& - \text{ nouvelles infections}\\ \text{variation du nombre d'individus infectés} &=& + \text{ nouvelles infections } - \text{guérisons} - \text{hospitalisations}\\ \text{variation du nombre d'individus hospitalisés} &=& + \text{ hospitalisations } - \text{sorties suite à guérisons ou décès}\\ \text{variation du nombre d'individus retirés} &=& + \text{ guérisons } + \text{ sorties suite à guérisons ou décès} \end{eqnarray*}

ce qui se traduit mathématiquement par le système d'équations différentielles suivant :

\begin{eqnarray*} \frac{\mathrm{d} S}{\mathrm{d} t} &=& -\beta \frac{S}{N} I\,,\\ \frac{\mathrm{d} I}{\mathrm{d} t} &=& +\beta \frac{S}{N} I - (\alpha+\rho) I\,,\\ \frac{\mathrm{d} H}{\mathrm{d} t} &=& \alpha I - \gamma H\,,\\ \frac{\mathrm{d} R}{\mathrm{d} t} &=& \rho I+ \gamma H\,. \end{eqnarray*}

La taille de la population est définie comme

N= S+I+H+R=\mbox{constante}\,.

$S/N$ $\beta SI/N$ .

Le nombre de reproduction de base du pathogène est défini comme

\mathcal R_0 = \frac{\beta}{\alpha+\rho}\frac{S(0)}{N}\,.

C'est le nombre d'infections secondaires générées par un individu infecté dans la population initiale :

$\mathcal R_0>1$ une vague épidémique se forme.
$\mathcal R_0<1$ l'épidémie s'éteint progressivement.

La probabilité d'être hospitalisé suite à l'infection est

p = \frac{\alpha}{\alpha+\rho}\,.

$t$ $H(t)$ $t$ (par unité de temps) :

A(t)=\alpha I(t)\,.

Le modèle épidémiologique est à base d'équations différentielles ordinaires (ordinary differential equations). Les simulations se baseront sur la fonction ode de la librairie deSolve de R à importer en début de script.


library(deSolve) # import de la fonction ode

Dans un souci d'hygiène algorithmique, nous commencerons toujours par faire table rase des exécutions précédentes pour travailler sur une une base propre à chaque nouvelle exécution du programme.


xxxxxxxxxx
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

Téléchargement et traitement des données

Les données sont disponibles sur le site suivant :

https://www.data.gouv.fr/fr/datasets/donnees-hospitalieres-relatives-a-lepidemie-de-covid-19/

Nous baserons notre travail sur les données hospitalières suivantes (à la date du TD) :

donnees-hospitalieres-covid19-2020-05-10-19h00.csv
donnees-hospitalieres-nouveaux-covid19-2020-05-10-19h00.csv

Téléchargez les fichiers .csv et placez les dans le répertoire de travail donné par la commande getwd().

Commençons par importer ces données dans R :


xxxxxxxxxx
dataH = read.table("./donnees-hospitalieres-covid19-2020-05-10-19h00.csv", header=TRUE, sep=";");# Nombre d'individus hospitalisés (H) - à partir du 19 mars
dataA = read.table("./donnees-hospitalieres-nouveaux-covid19-2020-05-10-19h00.csv", header=TRUE, sep=";");# Nombre d'admissions (A) à l'hôpital par jour

La deuxième colonne distingue les hommes et les femmes par des 1 et des 2. Comme le modèle n'est pas structuré selon ces catégories, nous ne conserverons que les lignes 0 qui agrègent hommes et femmes :


xxxxxxxxxx
# Nous ne distinguerons pas les hommes et les femmes
Z=which(dataH[,2]==0);dataH=dataH[Z,];

Pour ce qui suit, nous ne conserverons que les informations relatives aux hospitalisations, décès et guérisons, par départements :


xxxxxxxxxx
dataHH=dataH[,c(1,4)];# Département et nombre de personnes hospitalisées
dataAA=dataA[,c(1,3)];# Département et nombre d'admissions par jour

Comme le modèle est spatialement homogène, nous agrégerons les données à l'échelle de la France :


xxxxxxxxxx
n=95;# Nombre de département en France métropolitaine 
v=c(1:19,21:95);# Vecteur des départements Corse exclue

Pour récupérer la longueur des séries temporelles dans ces fichiers, on peut faire la gymnastique suivante :


xxxxxxxxxx
ZH=which(as.integer(as.character(dataH[,1]))==1);
LH=length(dataH[ZH,2]);#longueur de la série H 
ZA=which(as.integer(as.character(dataA[,1]))==1);
LA=length(dataA[ZA,2]);#longueur de la série A

ce qui nous permet de définir les matrices suivantes, que nous remplissons de 0 dans un premier temps :


xxxxxxxxxx
HH=matrix(0,n,LH);# Matrice des hospitalisés
AA=matrix(0,n,LA);# Matrice des admissions

Nous pouvons ensuite construire ces matrices de la façon suivante :


xxxxxxxxxx
for (i in v){
    ZH=which(as.integer(as.character(dataHH[,1]))==i);
    HH[i,]=dataHH[ZH,2];
    ZA=which(as.integer(as.character(dataAA[,1]))==i);  
    AA[i,]=dataAA[ZA,2];
}

ce qui permet finalement de construire deux vecteurs qui contiennent uniquement les données qui nous intéressent :


xxxxxxxxxx
H=matrix(0,1,LH); 
for (i in 1:LH) H[i]=sum(HH[,i]);# Hospitalisés
A=matrix(0,1,LA); 
for (i in 1:LA) A[i]=sum(AA[,i]);# Admissions

Visualisation des données à partir du 19 mars :


xxxxxxxxxx
plot(1:LH,H,xlab="Temps écoulé depuis le 19 mars (en jours)",
     ylab="Nombre de personnes hospitalisées",col="blue")


xxxxxxxxxx
plot(1:LA,A,xlab="Temps écoulé depuis le 19 mars (en jours)",
     ylab="Nombre d'admissions à l'hôpital (par jour)",col="red")
# ajout d'une courbe de tendance
panel.smooth(1:LA,A,col="red",col.smooth = "black")

Comme le confinement a débuté le 17 mars, la tendance croissante des deux premières semaines suit probablement une dynamique de transition suite à la mise en place du confinement.

Comme le modèle épidémiologique suppose que les paramètres épidémiologiques sont constants, en particulier le taux de transmission du virus, nous commencerons les analyses 15 jours après le début du confinement, soit à partir du 1er avril :


xxxxxxxxxx
T0=13;# Décalage du point de départ 13 jours après le 19 mars (le 1er avril) 
H=H[T0:length(H)];A=A[T0:length(A)];# Troncations des données
LH=length(H);LA=length(A);# Mise à jour des longeurs des séries temporelles

Visualisation des données tronquées à partir du 1er avril :


xxxxxxxxxx
plot(1:LH,H,xlab="Temps écoulé depuis le 1er avril (en jours)",
     ylab="Nombre de personnes hospitalisées",col="blue")


xxxxxxxxxx
plot(1:LA,A,xlab="Temps écoulé depuis le 1er avril (en jours)",
     ylab="Nombre d'admissions à l'hôpital (par jour)",col="red")

Simulations et ajustement du modèle aux données

Résolution numérique du modèle épidémiologique

Commençons par définir les paramètres du modèle épidémiologiques. Les valeurs données aux paramètres à estimer sont de grossières estimations initiales (ordres de grandeurs approximatifs).

Pour que le problème d'estimation des paramètres soit bien posé (c'est-à-dire pour que le modèle ne soit pas sur-paramétré par rapport aux données utilisées), Sallet, Iggidr et Rapaport (2020) ont montré qu'il faut à minima fixer la valeur d'un paramètre. La période infectieuse a été estimée à 10.91 jours en moyenne (avec un écart-type de 3.95 jours) en Chine (You et al 2020). Nous retiendrons la valeur moyenne.

durant le confinement $R(t)-R(0)$ $R(0)$ $R(0)$ , comme vous pourrez le vérifier.


xxxxxxxxxx
N=64e6;                     # Population hexagonale approximative
rho=1/10.91;                # Temps moyen avant guérison ou décès : 10.91 jours
gamma=1/14;                 # Durée moyenne de l'hospitalisation : 14 jours
R_0=1;                      # Nombre de reproduction de base en confinement
beta=R_0*rho;               # Taux de transmission (calcul très approximatif)
p=0.1;                      # Probabilité d'être hospitalisé suite à l'infection
alpha=p*rho/(1-p);          # Taux d'hospitalisation
P0=c(beta,alpha,gamma);     # Vecteur des paramètres

Définissons ensuite les variables du modèle et les conditions initiales (au 1er avril) :


xxxxxxxxxx
H0=H[1];        # Le nombre de personnes hospitalisées au 1er avril
I0=A[1]/alpha;  # Les admissions correspondent à A(t) = alpha*I(t)
R0=N*0.01;      # C'est la grande inconnue : mettons 1% d'immunisés 
S0=N-I0-H0-R0;  # La taille de la population sensible au 1er avril
X0=c(S0,I0,H0); # Vecteur d'état. Pas besoin de simuler R=N-(S+I+H).

Nous aurons besoin d'un vecteur contenant les dates sur lesquelles comparer modèle et données :


xxxxxxxxxx
t=0:(LH-1);     #Vecteur temps

Définissons la fonction SIHR qui prend en arguments 3 vecteurs : temps, variables d'états, et paramètres :


x
SIHR = function(t, X, P){
  beta=P[1];            # Le taux de transmission
  alpha=P[2];           # Le taux d'hospitalisation
  gamma=P[3];           # Le taux de sortie d'hôpital
    
  S=X[1];I=X[2];H=X[3]; # Le vecteur d'état X contient: S, I, et H
  y = beta*S*I/N;       # Le nombre de nouvelles infections par jour
  dS = -y;              # On exprime dS/dt = - beta*S*I
  dI = y-(alpha+rho)*I; # On exprime dI/dt = beta*S*I - (alpha+rho)*I
  dH = alpha*I -gamma*H;# On exprime dH/dt = alpha*I - gamma*H
      
  dX=c(dS,dI,dH);       # Renvoie dX/dt tel que demandé par la fonction ode
  return(list(dX));
}

$\alpha,\beta,\gamma$ $I(0)$ dans la population) qui sont les plus vraisemblables d'après les données. La vraisemblance d'un jeu de paramètres et de conditions initiales à estimer est la probabilité d'observer les données sachant ces paramètres et conditions initiales. Pour pouvoir calculer cette probabilité, il faut modéliser le processus d'observation de façon stochastique (avec variables aléatoires).

Modèle d'observation

distribution de Poisson $t$ $A(t)$ $t$ $H(t)$ $\gamma H(t)$ qui devraient être modélisées de la sorte ; nous pourrons également tester cette variante ensuite).

$\theta$ le vecteur des paramètres et conditions initiales à estimer :

\theta=\{\alpha,\beta,\gamma,S(0),I(0),H(0)\}\,.

Sous l'hypothèse (discutable) que les observations sont indépendantes conditionnellement au modèle épidémiologique (ce qui veut dire "sachant la dynamique prédite par le modèle théorique"), la vraisemblance (likelihood) peut s'écrire :

\mathcal{L}\left(\theta\right)=\prod_{t=0}^T \text{Pr}(\text{observer A admissions le jour }t| A(t))\times\text{Pr}(\text{observer H hospitalisations le jour }t|H(t))\,.

$\theta$ $\mathcal L(\theta)$ $\theta$ log-likelihood $\log\mathcal{L}(\theta)$ qui est plus communément utilisée en pratique.

Maximisation de la vraisemblance

Construisons une fonction qui calcule la log-vraisemblance d'un jeu de paramètres et de conditions initiales. Cette fonction prend en argument le vecteur des paramètres et des conditions initiales à estimer :


xxxxxxxxxx
logLike=function(theta){
    P=theta[1:3];       # Les paramètres beta, alpha, et gamma
    X0=theta[4:6];      # Mise à jour des conditions initiales 
    X=ode(X0,t,SIHR,P); # Résolution du système d'EDO (modèle SIHR)
    h=X[,4];                # Hospitalisation théoriques : H(t)
    a=P[2]*tail(X[,3],-1);  # Admissions théoriques : alpha*I(t)
    
    LLH=dpois(H,h,log=T);# Probabilité d'observer H (loi de Poisson)
    LLA=dpois(A,a,log=T);# Probabilité d'observer A (Poisson)
    LL=sum(c(LLH,LLA));# Log transforme produit des probas en somme
    return(LL);        # Renvoie la log-vraisemblance (likelihood)
}

Avant d'appeler les fonctions ci-dessus, on définit le vecteurs des paramètres et conditions initiales :


xxxxxxxxxx
theta0=c(P0,X0);        # Concatène paramètres et conditions initiales

Utilisons la fonction optim pour trouver les paramètres et conditions initiales qui maximisent la vraisemblance du modèle.


xxxxxxxxxx
opt=optim(theta0,logLike,control=list(fnscale=-1));# Maximise logLike

Résultats

On récupère les paramètres et conditions initiales dont la vraisemblance est maximale :


xxxxxxxxxx
# Les paramètres optimaux
beta=opt$par[1];alpha=opt$par[2];gamma=opt$par[3];
# Les conditions initiales optimales
S0=opt$par[4];I0=opt$par[5];H0=opt$par[6];

On met à jour les vecteurs des paramètres et conditions initiales :


xxxxxxxxxx
X0=c(S0,I0,H0);                 # Vecteur des conditions initiales
P0=c(beta,alpha,gamma);         # Vecteur des paramètres mis à jour

On affiche trois quantités importantes :


xxxxxxxxxx
R_0=beta/(alpha+rho)*S0/N;      # Nombre de reproduction de base estimé
print(R_0);
p=alpha/(alpha+rho);print(p);   # Probabilité d'être hospitalisé 
print(1/gamma);                 # Temps moyen d'hospitalisation

On trouve :

$\mathcal R_0 \approx 0.62$
$p \approx 0.08$
$\frac{1}{\gamma} \approx 16.93$ jours

$\mathcal R_0$ $\mathcal R_0<1$ ).

On calcule la solution du modèle pour les paramètres et conditions initiales estimés :


xxxxxxxxxx
T=120;t=0:T;            # Mise à jour du vecteur temps
X=ode(X0,t,SIHR,P0);    # Calcul de la solution optimale

On compare visuellement la solution du modèle et les observations :


xxxxxxxxxx
# Affiche le nombre d'individus hospitalisés (données et modèle)
plot(0:(LH-1),H,xlab="Temps écoulé depuis le 1er avril (en jours)",
     ylab="Nombre de personnes hospitalisées",col="blue");      
lines(X[,1],X[,4]);


xxxxxxxxxx
plot(0:(LA-1),A,xlab="Temps écoulé depuis le 1er avril (en jours)",
     ylab="Nombre de nouvelles hospitalisations (par jour)",col="red");
lines(X[,1],alpha*X[,3]);

"Prédictions" et discussion des limites du modèles

"Prédire est un art difficile, surtout en ce qui concerne l'avenir." Néanmoins, juste pour prendre du recul sur la dynamique épidémique, nous pouvons extrapoler les dynamiques précédentes sur un horizon de 4 mois (ce qui correspondrait à la dynamique prédite si nous devions rester confinés jusque fin juillet) :

Enfin, on peut afficher la fraction de la population qui a été infectée par le virus durant le confinement (depuis le premier avril) et est donc immunisée ou en voie de l'être (sans décompter les décès) :

F(t)\approx\frac{N-S(t)-R(0)}{N}.


xxxxxxxxxx
F=(N-X[,2]-R0)/N;
plot(40,F[40],xlim=c(0,T),ylim=c(0,max(F)),
     xlab="Temps écoulé depuis le 1er avril (en jours)", 
     ylab="Fraction de la population infectée");
lines(X[,1],F,col="red")

durant le confinement $R(0)$ , qui est inaccessible par cette seule étude.

En comparaison, Roques et al (2020) estiment la proportion de la population infectée durant cette même période (du 1er avril au 10 mai) à 3.7% (95%-CI: 3.0-4.8%), sur la base de jeux de données différents (nombre de cas confirmés, tests, décès en EPHAD, etc.). Cependant, Salje et al (2020) estiment qu'environ 1% (entre 0,5% et 2%) de la population aurait été infectée sur la même période (Figure 3 ; panel E).

Références

Voici une liste de références choisies pour leurs qualités pédagogiques et scientifiques ou leur visibilité dans les cercles médiatiques et décisionnels en France :

Le site du "Groupe de modélisation de l’équipe ETE" (CNRS, IRD, Université de Montpellier).
Le site "COVID-19" de l'unité de Recherche BioSP (INRAE) et en particulier cette pré-publication.
La pré-publication d'une étude dirigée par l'Institut Pasteur avec Inserm, CNRS, Santé Publique France.
Le site "COVID-2019" de l'unité de Recherche EPIcx (Inserm, Sorbonne Université) et en particuliers les Rapport 9 et Rapport 10.
Le site du "Groupe Maladie Infectieuse" de l'unité de recherche REPERES (CHU Rennes et EHESP) et en particulier cette pré-publication.

Remerciements

Je remercie Gauthier Sallet, Alain Rapaport, Lionel Roques, Samuel Soubeyrand, Nik Cunniffe, Tewfik Sari, Ludovic Mailleret et Samuel Alizon pour divers retours qui m'ont aidé à construire ce TD.

Ce(tte) œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale 4.0 International.