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 = 2017 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.
OpenEdition suggests that you cite this post as follows:
Arthur Charpentier (April 22, 2020). Estimer la surmortalité. Freakonometrics. Retrieved September 14, 2024 from https://doi.org/10.58079/ovfr
Salut Arthur,
the first link (http://data.cquest.org/insee_deces) is dead, unfortunately… I think I found the same data here:
https://www.insee.fr/fr/statistiques/2117120
but the format is less friendly? Any idea where I could find the data set you used?
a+
salut Nicolas, yes, it is now on https://www.data.gouv.fr/en/datasets/fichier-des-personnes-decedees/
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
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))