Estimer la surmortalité causée par l'épidémie de Covid-19 en France (2)¶
Combien de personnes sont décèdées des suites du coronavirus depuis le début de l'épidémie ? Au-delà du chiffre aggrégé sur la France entière, on peut aussi souhaiter connaître le nombre de décès par région, par âge, par sexe, etc. Pour cela, il nous faut stratifier les données selon ces différentes dimensions.
Téléchargements¶
Pour télécharger ce notebook : https://pierremarion23.github.io/notebooks/analyse_surmortalite_2.ipynb
Avant de l'exécuter, il faut télécharger et exécuter le noteboook de préparation des données : https://pierremarion23.github.io/notebooks/prep_donnees_covid.ipynb
Notebook précédent : https://pierremarion23.github.io/notebooks/analyse_surmortalite_1.ipynb
Source de données : INSEE (voir fichier de préparation des données pour les détails)
import os
import zipfile
import datetime
import sklearn
import bqplot
import pickle
import copy
import json
import requests as rq
import numpy as np
import pandas as pd
import ipywidgets as ipw
from ipywidgets import Layout
from bqplot import (Figure, LinearScale,ColorScale,
Color, Axis,GridHeatMap, HeatMap, ColorAxis, Label, OrdinalScale, Lines, Map, Mercator, PanZoom)
from IPython.display import display, Image
from sklearn.preprocessing import OneHotEncoder
from sklearn.linear_model import LinearRegression
from bqplot import pyplot as plt
from scipy import stats as scs
from numpy import random as npr
folder_treated_data = 'treated_data'
Notre source de donnée principale est à nouveau la table recensant tous les décès depuis 1998 jusqu'au 6 avril 2020, avec des informations sur chaque décès (sexe de la personne, âge, date du décès).
path = os.path.join(folder_treated_data, 'deces9820.pkl')
deces9820 = pd.read_pickle(path)
deces9820.tail()
Stratification des décès par classe d'âge et par région¶
Pour stratifier les données, on choisit le niveau de granularité suivant :
- en zone de décès : la région administrative (la Corse est fusionnée avec PACA car les effectifs sont trop petits autrement)
- en âge : 4 classes d'âge (20-39 ans, 40-59 ans, 60-74 ans, et 75 ans et plus)
- en date de décès : le mois
La granularité choisie répond à un équilibre entre le fait d'avoir des données suffisamment agrégées pour limiter les fluctuations statistiques, et suffisamment détaillées pour obtenir une information pertinente.
regions = {'Auvergne-Rhône-Alpes': ['01', '03', '07', '15', '26', '38', '42', '43', '63', '69', '73', '74'],
'Bourgogne-Franche-Comté': ['21', '25', '39', '58', '70', '71', '89', '90'],
'Bretagne': ['22', '29', '35', '56'],
'Centre-Val-de-Loire': ['18', '28', '36', '37', '41', '45'],
'Grand Est': ['08', '10', '51', '52', '54', '55', '57', '67', '68', '88'],
'Hauts-de-France': ['02', '59', '60', '62', '80'],
'Ile-de-France': ['75', '92', '93', '94', '95', '78', '91', '77'],
'Normandie': ['14', '27', '50', '61', '76'],
'Nouvelle-Aquitaine': ['16', '17', '19', '23', '24', '33', '40', '47', '64', '79', '86', '87'],
'Occitanie': ['09', '11', '12', '30', '31', '32', '34', '46', '48', '65', '66', '81', '82'],
'Pays de la Loire': ['44', '49', '53', '72', '85'],
"Provence-Alpes-Côte d'Azur-Corse": ['04', '05', '13', '83', '84', '2A', '2B']
# on enlève provisoirement le 06 à cause de pbs de remontées de donnée, cf https://www.insee.fr/fr/information/4470857
}
classes = {
'20-39': ['20 à 24 ans', '25 à 29 ans', '30 à 34 ans', '35 à 39 ans'],
'40-59': ['40 à 44 ans', '45 à 49 ans', '50 à 54 ans', '55 à 59 ans'],
'60-74': ['60 à 64 ans', '65 à 69 ans', '70 à 74 ans'],
'75+': ['75 à 79 ans', '80 à 84 ans', '85 à 89 ans', '90 à 94 ans', '95 ans et plus']
}
dep2region = {}
for region in regions:
for dep in regions[region]:
dep2region[dep] = region
age2class = {}
for c in classes:
for age in classes[c]:
age2class[age] = c
ages = list(age2class.keys())
def batch_age(a):
if pd.isna(a):
return a
return age2class[ages[min(int((a-20)//5), 15)]]
def batch_dep(d):
try:
return dep2region[d]
except:
return 'Autre'
deces9820.loc[:,'age'] = deces9820['age'].transform(batch_age)
deces9820.loc[:,'region'] = deces9820['depdeces'].transform(batch_dep)
#deces9820.drop(columns=['deces_annee'], inplace=True)
Un extrait de notre base de données de décès, avec les agrégations par classe d'âge et par région.
deces9820.tail()
On compte le nombre de décès par mois, classe d'âge, région et sexe.
deces_pivot = pd.pivot_table(deces9820[(deces9820["annee"] >= 1999) &(deces9820["region"] != "Autre")], values=["depdeces"], index=['age', 'sexe', 'region'], columns=['annee', 'month'], aggfunc='count', fill_value=0)
deces_pivot
deces_pivot_np = deces_pivot.to_numpy()
deces_pivot_np = deces_pivot_np[:,4:] # enlever les 4 premiers mois de 99, pour avoir un nombre de mois multiple de 12
Stratification de la population par classe d'âge et par région¶
On fait les mêmes opérations sur la table de population.
path = os.path.join(folder_treated_data, 'population_dep_classeage_sexe_9920.pkl')
pop9920 = pd.read_pickle(path)
def agg_region(index):
dep, sexe, annee = index
if dep == '06':
region = 'Autre'
else:
region = dep2region[dep]
return region, sexe, annee
pop9920_regions = pop9920.groupby(by=agg_region).sum()
pop9920_regions.index = pd.MultiIndex.from_tuples(pop9920_regions.index, names=('Région', 'Sexe', 'Année'))
pop9920_regions_classes = pop9920_regions.groupby(by=lambda a: age2class[a], axis=1).sum()
pop9920_regions_classes = pop9920_regions_classes.unstack().stack(-2)
pop9920_regions_classes = pop9920_regions_classes.iloc[8:].swaplevel(-1, -3).sort_index() # enlever région "Autre"
pop9920_regions_classes
pop9920_regions_classes_np = pop9920_regions_classes.to_numpy()
On calcule les taux de décès en divisant le nombre de morts par la population (comme dans le notebook précédent). On visualise ensuite l'évolution de ces taux au cours des années. Ils ont tendance à diminuer grâce à l'amélioration de l'espérance de vie. La composante saisonnière est beaucoup plus forte pour les personnes âgées que pour les personnes jeunes.
taux0020 = deces_pivot_np / np.repeat(pop9920_regions_classes_np[:,1:], 12).reshape(96, 21* 12)
i = 0
x = np.arange(1999.42, 2020.33, 1/12) # todo: better formula
for age in classes.keys():
for sexe in [1, 2]:
plt.figure(legend_location="bottom-left")
for j, reg in zip(range(12), list(regions.keys())):
plt.plot(x, taux0020[i*12+j, :-1], colors=[bqplot.CATEGORY10[j%10]], labels=[reg])
xax, yax = plt.axes()['x'], plt.axes()['y']
xax.tick_values = np.arange(2000, 2021, 1)
plt.ylim(min=0, max=taux0020[i*12:(i+1)*12, :-1].max())
plt.title('Taux de décès par région, sexe = {}, age = {}'.format(sexe, age))
plt.show()
i += 1
On peut enfin passer à la partie régression. Notre modèle est le suivant : le taux de mortalité au mois $m$, l'année $t$, pour la région $r$, la classe d'âge $a$ et le sexe $s$ s'écrit :
$$T_{m, t, a, s} = \alpha_{a} \rho_{r} \sigma_{s} \prod_{i=1}^{5}T_{m, t-i, a, s}^{ \beta_{i}}$$autrement dit comme une moyenne géométrique pondérée des taux de mortalité des cinq mois précédents, corrigé par des facteurs multiplicatifs pour la classe d'âge, la région et le sexe.
Concrètement, on passe au log, et on code chacune des caractéristiques (classe d'âge, région, sexe) en One Hot Encoding.
X0019 = taux0020.reshape(96, 21, 12)[:,:-1,-7:-2]
y0019 = taux0020.reshape(96, 21, 12)[:,:-1,-2]
X20 = np.log(taux0020.reshape(96, 21, 12)[:,-1,-7:-2])
X = np.log(np.rollaxis(X0019, 2).reshape(5, 96*20))
y = np.log(y0019.reshape(-1))
features = list(pop9920_regions_classes.index)
enc = OneHotEncoder(drop=['75+', 1, 'Auvergne-Rhône-Alpes'])
enc.fit(features)
one_hot_features = enc.transform(features)
one_hot = np.repeat(one_hot_features.toarray(), 20).reshape((96, 15, 20))
X = np.vstack([X, np.swapaxes(one_hot, 0, 1).reshape(15, 96*20)])
X20 = np.hstack([X20, one_hot_features.toarray()])
reg = LinearRegression()
reg.fit(X.T, y)
y_pred_hist = reg.predict(X.T).reshape(96, 20)
std_hist = np.std(y.reshape(96,20)-y_pred_hist, axis=1)
y_pred = reg.predict(X20)
ratio_type = (np.log(taux0020.reshape(96, 21, 12)[:,-1,-2]) - y_pred)/std_hist
deces_mars20 = np.exp(y_pred)*pop9920_regions_classes_np[:,-1]
On peut calculer le nombre de décès additionnels en mars. On trouve une valeur proche de celle trouvée dans le notebook précédent (c'est heureux !).
print('Décès additionnels en mars: ' + str(deces_pivot_np[:,-2].sum()-deces_mars20.sum()))
Mais l'intérêt de notre démarche était justement d'avoir des informations détaillées au niveau de la région, de l'âge, du sexe. Pour les visualiser, plutôt que de représenter directement la mortalité excédentaire, on visualise la déviation par rapport à notre prédiction, normalisée par l'écart-type empirique pour la strate correspondante. Par exemple, le premier graphe ci-dessous indique que la déviation de mortalité chez les hommes de 20-39 ans en Grand Est en mars 2020 atteint environ 3 fois la déviation standard dans cette population à cette période de l'année. En dessous de $\pm$1,5 déviation standard, on considère que la déviation est normale, et on l'écrête à zéro (c'est du "soft thresholding").
On constate des déviations positives importantes dans les régions touchées par l'épidémie :Ile-de-France, Grand Est et dans une moindre mesure Bourgogne-Franche-Comté, Hauts-de-France et Auvergne-Rhône-Alpes. Notons que ces déviations, si elles touchent en premier lieu la classe d'âge des 75 ans et plus, affectent aussi les autres classes d'âge, en particulier en Ile-de-France.
A l'inverse, des déviations négatives sont observées dans certaines régions peu touchées par l'épidémie, chez les personnes de moins de 74 ans. Ces déviations peuvent peut-être s'expliquer par le ralentissement de l'activité économique.
path = os.path.join(folder_treated_data, 'france.json')
with open(path) as data_file:
data = json.load(data_file)
regions2num = {'Auvergne-Rhône-Alpes': 12,
'Bourgogne-Franche-Comté': 2,
'Bretagne': 3,
'Centre-Val-de-Loire': 4,
'Grand Est': 0,
'Hauts-de-France': 7,
'Île-de-France': 11,
'Normandie': 8,
'Nouvelle-Aquitaine': 1,
'Occitanie': 6,
'Pays de la Loire': 9,
"Provence-Alpes-Côte d'Azur-Corse": 10
}
def soft_threshold(x, thres):
if x > - thres and x < thres:
return 0
return x
sc_geo = Mercator(scale_factor = 2700, center=(3.5, 48.8))
sc_c1 = ColorScale(scheme='Spectral', min=-7, max=7, reverse=True)
axis = ColorAxis(scale=sc_c1, side='right', orientation='vertical')
panzoom = PanZoom(allow_zoom=False)
i=0
figy=[]
for age in classes.keys():
figx=[]
for sexe in [1, 2]:
dic = {}
j = 0
for reg, num in regions2num.items():
dic[num] = soft_threshold(ratio_type[i*12+j], 1.5)
j += 1
map_styles = {'color': dic, 'scales': {'projection': sc_geo, 'color': sc_c1}, 'colors': {'default_color': 'Grey'}}
chloro_map = Map(map_data=data, **map_styles)
if sexe == 1:
fig = Figure(marks=[chloro_map], axes=[axis], title='Sexe={}, age={}'.format(sexe, age), layout=Layout(width='600px', height='500px'), interaction=panzoom)
else:
fig = Figure(marks=[chloro_map], title='Sexe={}, age={}'.format(sexe, age), layout=Layout(width='600px', height='500px'), interaction=panzoom)
figx.append(fig) #fig_margin = dict(top=30, bottom=10, left=20, right=20)
#if sexe == 1:
# figx.append(ipw.Image(value=open("ax_c2.png", "rb").read(), format='png', layout=Layout(max_height='300px')))
i += 1
figy.append(ipw.HBox(figx, layout=Layout(align_items='center')))
display(ipw.VBox(figy, align_content = 'stretch'))
On peut aussi représenter directement le nombre "brut" de décès excédentaire.
sc_c2 = ColorScale(scheme='BuPu', min=0, max=1000)
axis2 = ColorAxis(scale=sc_c2, side='right', orientation='vertical')
i=0
figy=[]
for age in classes.keys():
figx=[]
for sexe in [1, 2]:
dic = {}
j = 0
for reg, num in regions2num.items():
dic[num] = (deces_pivot_np[:,-2]-deces_mars20)[i*12+j]
j += 1
map_styles = {'color': dic, 'scales': {'projection': sc_geo, 'color': sc_c2}, 'colors': {'default_color': 'Grey'}}
chloro_map = Map(map_data=data, **map_styles)
if sexe == 1:
fig = Figure(marks=[chloro_map], axes=[axis2], title='Sexe={}, age={}'.format(sexe, age), layout=Layout(width='600px', height='500px'), interaction=panzoom,
fig_margin = dict(top=60, bottom=0, left=20, right=100))
else:
fig = Figure(marks=[chloro_map], title='Sexe={}, age={}'.format(sexe, age), layout=Layout(width='600px', height='500px'), interaction=panzoom)
figx.append(fig)
#if sexe == 1:
# figx.append(ipw.Image(value=open("ax_c.png", "rb").read(), format='png', layout=Layout(max_height='300px')))
i += 1
figy.append(ipw.HBox(figx, layout=Layout(align_items='center')))
display(ipw.VBox(figy, align_content = 'stretch'))
Transferts entre régions¶
L'analyse précédente est basée sur la région de décès des personnes. Or il est évident qu'on ne décède pas nécessairement dans la région où on est domicilié. On a donc potentiellement introduit un biais, en particulier à la lumière des transferts de patients ayant eu lieu pendant l'épidémie : il se peut qu'une partie de la mortalité apparaissant en Bretagne par exemple soit la conséquence du transfert de patient depuis l'Ile-de-France. Pour quantifier ce phénomène, on utilise les données exceptionnelles pour 2018-2020 mises à disposition par l'INSEE, avec en particulier la région de domiciliation des personnes décèdées. On va donc pouvoir comparer le nombre de personnes qui décèdent en dehors de leur région de domiciliation en temps normal et depuis le début de l'épidémie.
def to_region(dep):
try:
return dep2region[dep]
except:
return 'Outre-Mer'
path = open(os.path.join(folder_treated_data, "detail_deces1820.pkl"), "rb")
region_deces = pickle.load(path)
region_deces['regdeces'] = region_deces['depdeces'].transform(to_region)
region_deces['regdom'] = region_deces['depdom'].transform(to_region)
region_deces.tail()
def compute_transfert(date_debut, date_fin):
regdec = region_deces['regdeces'][(region_deces['deces_date'] >= date_debut) & (region_deces['deces_date'] <= date_fin)]
regdom = region_deces['regdom'][(region_deces['deces_date'] >= date_debut) & (region_deces['deces_date'] <= date_fin)]
transfert = pd.concat([regdec, regdom], axis=1)
transfert_pd = pd.pivot_table(transfert, index=[regdec], columns=[regdom], aggfunc='count', fill_value=0)
transfert_matrix = np.array(transfert_pd)[:,:13] / (date_fin - date_debut).days
for i in range(13):
transfert_matrix[i,i]=0
return transfert_matrix
def plot_heatmap(transfert_matrix, date_begin, date_end, show_c=True):
x = ['ARA', 'BFC', 'Bre.', 'CVL', 'GE', 'HdF', 'IdF', 'Nor.', 'NA', 'Occ.', 'OM', 'PdL', 'PACA']
x_sc, y_sc, col_sc = OrdinalScale(domain=x), OrdinalScale(domain=x), ColorScale(scheme='BuPu', min=0, max=3)
heat = GridHeatMap(column=x, row=x, color=transfert_matrix,
scales={'row': y_sc, 'column': x_sc, 'color': col_sc}, stroke='white')
# heat.display_format = '.2f'
heat.font_style={'font-size': '12px', 'fill':'black', 'font-weight': 'bold'}
ax_x = Axis(scale=x_sc, label='Lieu de résidence')
ax_y = Axis(scale=y_sc, orientation='vertical', grid_lines='none', label='Lieu de décès', label_offset='45px')
ax_c = ColorAxis(scale=col_sc, side='right', orientation='vertical')
if show_c:
axes = [ax_x, ax_y, ax_c]
else:
axes = [ax_x, ax_y]
fig = Figure(marks=[heat], axes=axes,
title='Transferts entre {} et {}'.format(date_begin.strftime("%d/%m/%Y"), date_end.strftime("%d/%m/%Y")),
layout=Layout(width='700px', height='600px'),
min_aspect_ratio=1, max_aspect_ratio=1, padding_y=0, background_style={'fill':'white'})
return fig
Le graphe représente le nombre moyen journalier de personnes qui décèdent dans une autre région que leur région de domicile. On observe que les transferts concernent de manière importante la région Ile-de-France, ce qui peut s'expliquer par le dynamisme économique de la région (beaucoup de voyageurs de passage) et par sa concentration importante en hopitaux. Ces transferts représentent au total environ 50 personnes par jour, ce qui est faible devant la mortalité totale (environ 1500 personnes par jour en France). Le biais introduit est donc raisonnable.
On trace le même graphe pendant l'épidémie. On n'observe pas de différence significative, ce qui nous amène à penser que les transferts de patients n'ont pas eu d'effet visible sur la mortalité trans-régionale.
date_begin = datetime.datetime(2018,1,1)
date_end = datetime.datetime(2020,2,29)
tm1 = compute_transfert(date_begin, date_end)
fig1 = plot_heatmap(tm1, date_begin, date_end)
date_begin = datetime.datetime(2020,3,1)
date_end = datetime.datetime(2020,4,13)
tm2 = compute_transfert(date_begin, date_end)
fig2 = plot_heatmap(tm2, date_begin, date_end, show_c=False)
display(ipw.HBox([fig1, fig2]))
#display(Image("heatmaps.png"))
