Estimer la surmortalité

Ce matin, Baptiste Coulmont m’envoyait un tweet avec un joli graphique, présentant le nombre quotidien de décès en France.

Comme les données sont en ligne, je me suis dit que je pourrais jouer un peu avec. Pour les plus curieux, on a la liste de tous les décès depuis… longtemps ! (on a plus de 25 millions de décès). La seule chose qui nous intéresse, c’est la date, alors on la récupère. Puis on compte, par jour, combien de décès on a eu,

D = read.csv("insee_deces.csv",header=TRUE)
Vecteur_Dates = as.character(D[,8])
TV = table(Vecteur_Dates)
B = data.frame(dateC=names(TV),nb=as.numeric(TV))
B$date=as.Date(B$dateC,"%Y-%m-%d")
B$year = format(B$date,"%Y")
B$month = format(B$date,"%m")
B$day = format(B$date,"%d")
B$deb=as.Date(paste(B$year,"01-01",sep="-"),"%Y-%m-%d")
B$fin=as.Date(paste(B$year,"12-31",sep="-"),"%Y-%m-%d")
B$dif=as.numeric(B$date-B$deb)/as.numeric(B$fin-B$deb)

Je crée ici une une variable qui me dit où je suis dans l’année en cours (0 pour le 1er janvier, et 1 pour pour 31 décembre). Ça me permet de contourner un petit problème technique : les années bissextiles. On peut ensuite représenter toutes les années,

plot(B2$dif[B2$year == "2020"],B2$NormCpte[B2$year == "2020"],type="l",xlim=0:1,ylim=c(1200,4000),col="white")
for(i in as.character(2000:2019)) lines(B2$dif[B2$year == i],B2$NormCpte[B2$year == i],col="light blue")

On veut ensuite avoir la tendance moyenne. Baptiste prenait la moyenne par jour (la courbe noire). Ici, on va lisser avec des splines

library(gam)
reg = gam(NormCpte~bs(dif,40),data=B2[B2$year != "2020",])
vx = seq(0,1,length=501)
vy = predict(reg,newdata=data.frame(dif=vx))
lines(vx,vy,col="blue",lwd=3)

ce qui donne la courbe bleue. C’est l’écart à cette courbe qui donne la surmortalité. Par exemple (prenons un exemple bien connu, et largement visible sur le graphique précédant) l’année 2003, et les quinze premiers jours d’août, correspondant (plus ou moins) à la période de canicule.

i = 2003
x = B2$dif[(B2$year == i)&((B2$month == "08"))&(B2$day %in% c(paste("0",1:9,sep=""),10:15))]
y = B2$NormCpte[(B2$year == i)&((B2$month == "08"))&(B2$day %in% c(paste("0",1:9,sep=""),10:15))]
yp = predict(reg,newdata=data.frame(dif=x))
e = y-yp
sum(e)
[1] 14294.39

Plusieurs chiffres sont mentionnés sur wikipédia, dont le rapport de l’INSERM qui annonçait 14,802 morts. On n’est pas trop loin… On peut visualiser cette surmortalité sur le graphique ci-dessous

plot(B2$dif[B2$year == "2020"],B2$NormCpte[B2$year == "2020"],type="l",col="white",xlim=0:1,ylim=c(1200,4000))
for(i in as.character(2000:2019)) lines(B2$dif[B2$year == i],B2$NormCpte[B2$year == i],col="light blue")
lines(vx,vy,col="blue")
lines(B2$dif[B2$year == 2003],B2$NormCpte[B2$year == 2003],col="red")
for(u in 1:length(x)) segments(x[u],y[u],x[u],yp[u],lwd=3,col="red")

Prenons un autre exemple qu’on voit bien sur le graphique de Baptiste, la grippe de 2017 (dite de 2016-2017). Si on prend les 45 premiers jours de 2017, on obtient les chiffres suivants

i = 2003
x=B2$dif[(B2$year == i)&((B2$month == "01")|
  ((B2$month == "02")&(B2$day %in% c(paste("0",1:9,sep=""),10:15))))]
y=B2$NormCpte[(B2$year == i)&((B2$month == "01")|
  ((B2$month == "02")&(B2$day %in% c(paste("0",1:9,sep=""),10:15))))]
yp=predict(reg,newdata=data.frame(dif=x))
yp = predict(reg,newdata=data.frame(dif=x))
e = y-yp
sum(e)
[1] 21177.33

Le rapport de Santé Publique France parle de 21,200 décès (pour la “surmortalité toute cause”) ce qui est là aussi comparable…

Maintenant, avant de conclure, et avant de me faire troller pendant quelques jours, je tiens à préciser que c’est juste un exercice de modélisation, et oui, j’ai exclu 2020. Je ne tiens pas à faire de comparaison avec la Covid-19, à qui 21,000 morts sont attribués, en France, ce soir. Si l’ordre de grandeur est le même, qu’on en me fasse pas dire qu’on est dans une situation comparable, et que la Covid-19 est finalement une simple grippe. Je ne me prononcerais pas sur ce point car je ne suis pas virologue. Mais surtout, on a 21,000 morts malgré des mesures incroyables de confinement. Je reviendrais probablement dans les jours qui viennent sur l’impact du contrôle par confinement, mais l’objet du court billet de ce soir était de voir comment on pouvait rapidement quantifier une surmortalité.

  • petite note technique : pour les 3 graphiques (mais pas pour les codes et les calculs) j’avais initialement normalisé les données. En effet, la population française a relativement pas mal changé en 20 ans, augmentant de 10%… Aussi, 1000 morts en 2003, ça n’est pas autant (relativement à la taille de la population) qu’en 2020.

14 thoughts on “Estimer la surmortalité”

  1. Bonsoir Monsieur Charpentier,

    Merci pour cette démonstration, c’est intéressant de voir comment utiliser R, avec des données officielles, pour avoir des beaux graphiques, ils ne sont pas finalement si ‘moches’ que ça, comparativement à ceux fournis par python.

    En revenant au choc de mortalité, est ce qu’on a atteint les 15% prévu dans Solvabilité2? et quel serait son impact sur les branches de l’assurance vie?

    Cordialement

  2. Et pour l’obtention du premier graphique, en Stata (sacrilège sur ce blog ? 😉 ), normé de la taille de la population (c’est dans une base pop qu’il est facile de faire), corrigé des années bissextiles et avec une simple moyenne hebdomadaire, je propose (on ne sait jamais !) :

    *** Création du fichier issu des dernières stats INSEE “special Covid”
    insheet using 2020-04-17_deces_quotidiens_departement.csv, clear delim(“;”)
    keep if zone == “France”
    keep date_evenement total_deces_2020
    replace total_deces_2020 = “” if total_deces_2020 == “NA”
    destring total_deces_2020 , replace
    keep if !missing(total_deces_2020 )
    generate nbdc = total_deces_2020 in 1
    replace nbdc = total_deces_2020- total_deces_2020[_n-1] if _n>1
    generate date_d = date(date_evenement , “DMY”)
    format date_d %td
    drop date_evenement total_deces_2020
    save stat_covid, replace

    * Importation et traitement de la base des décès INSEE
    insheet using insee_deces.csv, clear
    generate date_d = date(date_deces , “YMD”)
    format date_d %td
    drop if date_d= td(01mar2020)
    append using stat_covid

    *** identification des jours de l’année et de l’année
    gen jour_d=doy(date_d)
    gen annee_d=year(date_d)

    *traitement des années bisextiles
    *les 1er mars et suivants seront tous codés pareils
    foreach j of numlist 1972(4)2020 {
    if `j’== 1972 {
    local exclu 1972
    }
    else {
    local exclu `exclu’ ,`j’
    }
    }
    macro list _exclu

    replace jour_d = jour_d+1 if !inlist(annee_d, `exclu’) & jour_d>58

    *Normalisation par la taille annuelle de la population
    merge m:1 annee_d using pop
    replace nbdc = nbdc / poptot *1000000
    drop poptot _merge

    *Passage au mode large pour générer une variable par an
    drop date_d
    reshape wide nbdc , i(jour_d ) j(annee_d)

    *Création des moyennes, min et max
    tsset jour_d, daily
    foreach mth in mean min max {
    egen nbdc`mth’=row`mth'(nbdc2000-nbdc2019)
    }

    *Etiquette des variables
    label variable nbdcmean “Moyenne annuelle 2000-2019”
    forvalues j=2000/2020 {
    label variable nbdc`j’ “Décès quotidiens en `j'”
    }

    * Graphique
    twoway (rarea nbdcmin nbdcmax jour_d) (line nbdcmean jour_d) (line nbdc2020 jour_d), xtitle(“”) xlabel(,format(%2.0f)) legend(order(1 “Intervalle de fluctuation 2000-2019” 2 3) col(1))

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.