Combien de personnes sont décèdées des suites du coronavirus depuis le début de l'épidémie ? Au 13 avril, le chiffre officiel est de 9 588 décès (en milieu hospitalier et en EHPAD), mais chacun se doute que le "vrai" chiffre est largement supérieur. Mais de combien ? Nous allons voir qu'il y a eu environ 15 400 décès additionnels causés par l'épidémie jusqu'au 13 avril, soit environ 1.6 fois le chiffre officiel.
Pour télécharger ce notebook : https://pierremarion23.github.io/notebooks/analyse_surmortalite_1.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 suivant : https://pierremarion23.github.io/notebooks/analyse_surmortalite_2.ipynb
Source de données : INSEE (voir fichier de préparation des données pour les détails)
Afin d'éclairer la mortalité due à l'épidémie du coronavirus, l'INSEE a récemment publié des données de mortalité des premiers mois des années 2018 à 2020, agrégées selon différentes dimensions (âge, sexe, département, jour de décès) [1]. Ces données ont été reprises dans les médias, voir par exemple un article du Monde, afin d'estimer la surmortalité en 2020 causée par l'épidémie. Néanmoins, une comparaison directe du nombre de décès en mars 2020 à celui en 2018 et 2019 comprend des biais importants, et ce pour deux raisons principales :
Nous allons essayons de corriger ces deux biais.
Attention, changer cette valeur avant d'exécuter le notebook (voir nb de préparation des données)
folder_treated_data = 'treated_data'
import os
import zipfile
import datetime
import sklearn
import bqplot
import pickle
import copy
import requests as rq
import numpy as np
import pandas as pd
import ipywidgets as ipw
from bqplot import pyplot as plt
from bqplot import DateScale, LinearScale, Scatter, Lines, Axis, Figure
from scipy import stats as scs
from numpy import random as npr
from sklearn.linear_model import LinearRegression
from bqplot.traits import convert_to_date
Dates de début et de fin de la période d'étude
date_debut = datetime.datetime(2020, 2, 15)
date_fin = datetime.datetime(2020, 4, 13)
period = [date_debut + datetime.timedelta(days=k) for k in range((date_fin-date_debut).days+1)]
Notre source de donnée principale est une 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, département de décès).
path = os.path.join(folder_treated_data, 'deces9820.pkl')
deces9820 = pd.read_pickle(path)
deces9820.tail()
total2020 = ((deces9820['deces_date'] >= date_debut) & (deces9820['deces_date'] <= date_fin)).sum()
total2019 = ((deces9820['deces_date'] >= date_debut-datetime.timedelta(days=366)) & (deces9820['deces_date'] <= date_fin-datetime.timedelta(weeks=52))).sum()
total2018 = ((deces9820['deces_date'] >= date_debut-datetime.timedelta(days=365+366)) & (deces9820['deces_date'] <= date_fin-datetime.timedelta(weeks=2*52))).sum()
print("Deces du {} jusqu'au {} : {}".format(date_debut, date_fin, total2020))
print("Deces du {} jusqu'au {} : {}".format(date_debut-datetime.timedelta(days=366), date_fin-datetime.timedelta(days=366), total2019))
print("Deces du {} jusqu'au {} : {}".format(date_debut-datetime.timedelta(days=365+366), date_fin-datetime.timedelta(days=366+365), total2018))
print(total2020 - 0.5*(total2019 + total2018))
En analysant directement les données de l'INSEE, on conclut donc à une augmentation des décès de 5 884 sur la période par rapport aux deux années précédentes. Or on sait qu'au jour du 13 avril, 9 588 personnes étaient décèdées du coronavirus en milieu hospitalier et EHPAD (chiffres officiels de Santé Publique France), donc on a largement sous-estimé l'impact de l'épidémie par cette méthode naïve.
Pour tenir compte de l'augmentation de la population, on va diviser le nombre de décès quotidiens par son espérance, donnée par des tables de mortalité et de population mises à jour chaque année. Commençons par préparer ces tables grâce aux données démographiques de l'INSEE. On agrège les décès par âge, sexe, et année de décès.
deces_sexe_annee_age = pd.pivot_table(deces9820,values='depdeces', index=['sexe', 'annee'], columns=['age'], aggfunc='count', fill_value=0)
deces_sexe_annee_age[100] = deces_sexe_annee_age.loc[:,[x for x in deces_sexe_annee_age.columns if x >= 100]].sum(axis=1)
deces_sexe_annee_age.drop(columns=[x for x in deces_sexe_annee_age.columns if x > 100], inplace=True)
deces_sexe_annee_age.head()
On introduit une deuxième source de données : la popluation par âge et par sexe de 1991 à 2020.
path = os.path.join(folder_treated_data, 'population_age_sexe_9820.pkl')
population_sexe_annee_age = pd.read_pickle(path)
population_sexe_annee_age.head()
plt.figure()
plt.plot(np.arange(20, 101), deces_sexe_annee_age.loc[(1, 2014)], labels=['Décès'])
plt.plot(np.arange(20, 101), population_sexe_annee_age.loc[(1, 2014)], colors=['green'], labels=['Population'])
plt.title("Population et nombre de décès dans l'année par âge au 1er janvier chez les hommes en 2014")
plt.legend()
plt.show()
En divisant le nombre de morts par la population, on obtient des probabilités de décès par année, sexe, et âge.
pourc_mortalite = deces_sexe_annee_age.loc[(slice(None), slice(1998, 2019)), :] / population_sexe_annee_age.loc[(slice(1, 2), slice(1998, 2019)), :]
pourc_mortalite.head()
On fait une moyenne glissante sur 5 ans, pour lisser les fluctuations aléatoires, ce qui nous donne des tables de décès par année, sexe, et âge, de 2000 à 2017.
pourc_mortalite_np = pourc_mortalite.to_numpy().reshape(2, 22, 81)
pourc_mortalite_np = np.hstack((np.zeros((2, 1, 81)), pourc_mortalite_np))
table_mortalite_av = np.cumsum(pourc_mortalite_np, axis=1)
# moyenne glissante entre 2000 et 2017
table_mortalite_av = ((table_mortalite_av[:,5:,:]-table_mortalite_av[:,:-5,:]) / 5)
Pour obtenir les tables de décès en 2018 et 2020, on projette de manière un peu brutale : on calcule la variation moyenne de probabilité de décès pour chaque âge et chaque sexe, entre 2014 et 2017, puis on suppose que cette variation continue de la même manière entre 2017 et 2020.
# Gain annuel moyen entre 2014 et 2017
gain_annuel = np.mean(np.diff(table_mortalite_av[:, -4:, :], axis=1), axis=1)
# projection sur 2018-2020
n_annee_proj = 3
table_mortalite_proj = np.zeros((2, n_annee_proj, 81))
for k in range(n_annee_proj):
table_mortalite_proj[:, k, :] = table_mortalite_av[:, -1, :] + (k+1) * gain_annuel
table_mortalite_np = np.hstack((table_mortalite_av, table_mortalite_proj))
table_mortalite = pd.DataFrame(table_mortalite_np.reshape(42,81))
table_mortalite.columns = pourc_mortalite.columns
table_mortalite.index = pd.MultiIndex.from_arrays([[1]*21 + [2]*21, np.hstack((np.arange(2000, 2021), np.arange(2000, 2021)))])
On obtient finalement une table de mortalité jusqu'en 2020.
table_mortalite.tail()
On peut visualiser ces tables de mortalité en traçant la fonction de survie (c'est-à -dire la probabilité d'atteindre un âge donné), en fonction du sexe et de l'année. On visualise le gain d'espérance de vie au cours des années, ainsi que la différence entre hommes et femmes.
fct_survie_H = np.cumprod(1-table_mortalite_np[0], axis=1)
fct_survie_F = np.cumprod(1-table_mortalite_np[1], axis=1)
plt.figure()
for k in range(0, 24, 4):
plt.plot(np.arange(20, 101), fct_survie_H[k], colors=[bqplot.CATEGORY10[k//4]], labels=['20{:02}'.format(k)], display_legend=True)
plt.plot(np.arange(20, 101), fct_survie_F[k], colors=[bqplot.CATEGORY10[k//4]])
plt.xlabel('Âge')
plt.ylabel('Valeur de la fonction de survie')
plt.title('Evolution de la fonction de survie des adultes sur 20 ans')
plt.show()
A partir des tables de mortalité et des tables de population, on calcule facilement l'espérance du nombre de décès annuel.
# espérance du nombre de morts annuels entre 2000 et 2020
morts_par_an_esp = np.einsum('ijk,ijk->j', table_mortalite_np, population_sexe_annee_age.loc[(slice(1, 2), slice(2000, 2020)), :].to_numpy().reshape(2, 21, 81))
plt.figure()
plt.plot(np.arange(2000, 2021), morts_par_an_esp)
plt.title('Espérance du nombre de décès annuels dans la population adulte')
plt.ylim(0, 640000)
plt.show(display_toolbar=False)
On constate que l'espérance du nombre de décès annuels, qui stagne autour de 550 000 entre 2000 et 2012 augmente rapidement dans la décennie 2010 pour atteindre 650 000 en 2020.
Traçons maintenant le nombre de décès journaliers en mars, pour les années entre 2000 et 2020, divisés par 1/365ème de l'espérance des décès l'année correspondante, afin de tenir compte de l'augmentation de la population. Grâce à cette normalisation, on s'attend donc à ce que ces courbent présentent des fluctuations aléatoires autour de 1, sauf la courbe de 2020 qui devrait exploser dans la deuxième moitié de mars.
decalage = (date_debut - datetime.datetime(2020, 1, 1)).days
duree_periode = (date_fin-date_debut).days + 1
n_training = 20
deces9820['delta_debut'] = deces9820['jour_annee'] - decalage
deces_par_jour = pd.pivot_table(deces9820[(deces9820['delta_debut'] >= 0) & (deces9820['delta_debut'] <= duree_periode-1) & (deces9820['annee'] >= 2000)], values='depdeces', index=['annee'], columns=['delta_debut'], aggfunc='count', fill_value=0)
deces_par_jour.tail()
deces_par_jour_np = deces_par_jour.to_numpy()
# Entre 2000 et 2019
y_0019 = 365 * deces_par_jour_np[:-1] / np.repeat(morts_par_an_esp[:-1], repeats=duree_periode).reshape(n_training, duree_periode)
y_20 = 365 * deces_par_jour_np[-1] / np.repeat(morts_par_an_esp[-1], repeats=duree_periode)
plt.figure()
x = pd.to_datetime(period)
for k in range(n_training):
plt.plot(x, y_0019[k], colors=[bqplot.CATEGORY10[k%10]])
plt.plot(x, y_20, colors=['red'], labels=['2020'], display_legend=True, marker='circle')
plt.xlabel("Période d'étude")
plt.ylabel('Ratio du nombre de décès')
plt.title("Nombre de décès journaliers rapporté à 1/365e de l'espérance des décès annuels, entre 2000 et 2020")
plt.ylim(min=0.7, max=1.5)
plt.show(display_toolbar=False)
On constate que les courbes sont tendanciellement au dessus de 1, en particulier pendant la première quinzaine de mars. Ceci est dû à la répartition inégale des décès dans l'année : il y a plus de décès journaliers en mars qu'en moyenne sur l'année, en particulier à cause de l'épidémie de grippe. Les quelques courbes présentant des ratios supérieurs à 1.2 correspondent à des épidémies de grippe très virulentes en mars, sauf 2020, où l'explosion en deuxième moitié de mois est due au coronavirus.
On remarque aussi que mars 2020 avait commencé avec une mortalité très faible (même historiquement faible autour du 7-8 mars). Cela est sans doute dû à une combinaison de facteurs : grippe clémente en 2020, diminution des morts accidentelles grâce au confinement, problèmes de remontée de données les 7 et 8 (qui sont un week-end).
Malgré cela, on fait l'hypothèse en première approximation que le nombre de décès pendant la première quinzaine de mars n'a quasiment pas été influencé par l'épidémie. Pour connaître le nombre de morts dus au coronavirus, on souhaiterait donc "prolonger" la courbe de début mars 2020 pour obtenir une distribution de courbes de mortalité dans la seconde moitié du mois de mars, en l'absence de coronavirus. Ensuite il suffira de calculer la différence entre la courbe réalisée et la moyenne de cette distribution, pour estimer le nombre de morts excédentaires.
Pour cela, on propose un modèle de régression linéaire bayésienne. On suppose donc que les décès au $i$-ème jour de mars l'année $t$ s'écrivent :
$$y_{t,i} = a_t i + b_t \exp(-(i-15)^2/128) + c_t + \epsilon_{t,i}$$où les $(\epsilon_{t,i})_i$ sont i.i.d. de loi normale de variance $\sigma_t^2$, les $(a_t, b_t, c_t)$ sont i.i.d. et suivent une loi normale (à 3 dimensions) également de variance $\sigma_t^2$, et enfin les $\sigma_t^2$ sont i.i.d. de loi inverse-gamma. Le rôle du coefficient $b_t$ est de modéliser une grippe tardive en mars (courbe en cloche centrée autour du 1er mars).
Pour justifier cette hypothèse de normalité sur $(a_t, b_t)$, on peut estimer les coefficients $a_t$ et $b_t$ pour $t$ entre 2000 et 2019. On commence par tracer chaque régression linéaire.
x1 = np.arange(duree_periode)
x2 = scs.norm.pdf(np.arange(duree_periode), loc=20, scale=8)
X = np.array([x1, x2])
x_sc = DateScale()
y_sc = LinearScale()
x = pd.to_datetime(period)
ax_x = Axis(label='X', scale=x_sc, grid_lines='solid', tick_values = pd.to_datetime([date_debut + datetime.timedelta(days=k) for k in range(0, (date_fin-date_debut).days+1, 15)]))
ax_y = Axis(label='Y', scale=y_sc, orientation='vertical', grid_lines='solid', max=1.4)
scatters = []
lines = []
coefs = []
intercepts = []
for k in range(n_training):
# Instantiate a Linear regression model
reg = LinearRegression()
# Fit to data
reg.fit(X.T, y_0019[k])
y_pred = reg.predict(X.T)
scatters.append(Scatter(x=x, y=y_0019[k], scales={'x': x_sc, 'y': y_sc}, default_size=32))
lines.append(Lines(x=x, y=y_pred, scales={'x': x_sc, 'y': y_sc}))
coefs.append(reg.coef_)
intercepts.append(reg.intercept_)
coefs = np.array(coefs)
figy=[]
for i in range(5):
figx=[]
for j in range(4):
figx.append(Figure(axes=[ax_x, ax_y], marks=[scatters[i*4+j], lines[i*4+j]], title=str(2000 + i*4+j), fig_margin = dict(top=30, bottom=20, left=20, right=20), layout={'width':"240px", 'height':'200px'}))
figy.append(ipw.HBox(figx))
display(ipw.VBox(figy))
Ensuite on trace le scatter plot des coefficients $a_t$, $b_t$ et $c_t$ de la régression linéaire, afin de vérifier que leur distribution est approximativement gaussienne.
fig = []
scatters = [Scatter(x=coefs[:,0], y=coefs[:,1], scales={'x': LinearScale(), 'y': LinearScale()}),
Scatter(x=coefs[:,1], y=intercepts, scales={'x': LinearScale(), 'y': LinearScale()}),
Scatter(x=coefs[:,0], y=intercepts, scales={'x': LinearScale(), 'y': LinearScale()})]
for scatter, title in zip(scatters, ['a versus b', 'b versus c', 'a versus c']):
ax_x = Axis(scale=x_sc, grid_lines='solid')
ax_y = Axis(scale=y_sc, orientation='vertical', grid_lines='solid')
fig.append(Figure(axes=[ax_x, ax_y], marks=[scatter], title=title, fig_margin = dict(top=30, bottom=10, left=20, right=20), layout={'width':"400px", 'height':'400px'}))
display(ipw.HBox(fig))
Pour les paramètres a priori, on utilise des estimateurs basés sur la moyenne des années 2000 à 2019, en enlevant les années de grippe virulente (ratio de décès le 1er mars supérieur à 1.15). On met à jour les paramètres a posteriori dans la fonction compute_post_param en utilisant les n premiers jours de mars 2020. Les formules sont basées sur celles données dans l'article Wikipedia sur la régression linéaire bayésienne.
n = -1
x0 = np.ones((1, duree_periode))
f = 30
y = np.mean(y_0019, axis=0)
T = y_0019.shape[0]
X = np.vstack((x0, x1, x2)).T
lambd_0 = np.matmul(X.T, X) / f
mu_0 = np.matmul(np.matmul(np.linalg.inv(lambd_0), X.T), y) / f
a_0 = duree_periode/(2*f)
b_0 = 0.5 * (np.sum(y*y) - f * np.sum(mu_0 * (lambd_0 @ mu_0)))
lambd_post = {}
mu_post = {}
a_post = {}
b_post = {}
def compute_post_param(n_start):
X_2020 = np.vstack((x0, x1, x2))[:,:n_start].T
lambd_post[n_start] = lambd_0 + np.matmul(X_2020.T, X_2020)
mu_post[n_start] = np.matmul(np.linalg.inv(lambd_post[n_start]), lambd_0 @ mu_0 + X_2020.T @ y_20[:n_start])
a_post[n_start] = a_0 +n_start/2
b_post[n_start] = b_0 + 0.5 * (np.sum(y_20[:n_start]*y_20[:n_start]) + np.sum(mu_0 * (lambd_0 @ mu_0)) - np.sum(mu_post[n_start] * (lambd_post[n_start] @ mu_post[n_start])))
def sample_post(N, n_start):
y_sims = np.zeros((N, duree_periode))
for k in range(N):
var = scs.invgamma.rvs(a_post[n_start], size=1) * b_post[n_start]
beta = npr.multivariate_normal(mu_post[n_start], var * np.linalg.inv(lambd_post[n_start]), size=1)
y_sims[k] = npr.multivariate_normal(np.dot(X, beta.T).flatten(), var * np.identity(duree_periode))
y_sims[k,:n_start] = y_20[:n_start]
return y_sims
n_starts = [14, 20, 29] # on fait la simulation pour différents jours de départ, ici le 14ème, 20ème, 28ème
N = 10**3
y_sims = {}
moys = {}
env_infs = {}
env_sups = {}
for n_start in n_starts:
compute_post_param(n_start)
y_sims[n_start] = sample_post(N, n_start)
moys[n_start] = np.mean(y_sims[n_start], axis=0)
env_infs[n_start] = np.sort(y_sims[n_start], axis=0)[N//10]
env_sups[n_start] = np.sort(y_sims[n_start], axis=0)[N-N//10]
On peut maintenant tracer l'intervalle de confiance (à 90%) de l'évolution à partir du 15 mars du ratio de décès journalier. On constate que 2018 est entièrement en dehors de l'IC (à cause de la grippe), tandis que 2020 sort de l'IC à partir du 16 mars.
plt.figure(legend_location='bottom-left')
x = pd.to_datetime(period)
for k in range(len(y_0019[:-2])):
plt.plot(x, y_0019[k], colors=['gray'], opacities=[0.2])
plt.plot(x, y_0019[-2], colors=['blue'], labels=['2018'], display_legend=True)
plt.plot(x, y_0019[-1], colors=['violet'], labels=['2019'], display_legend=True)
plt.plot(x, y_20, colors=['red'], labels=['2020'], display_legend=True)
#plt.plot(x, pred, colors=['green'], labels=['2020 non bay'], display_legend=True)
for n_start in [29]:
plt.plot(x[n_start-1:], env_infs[n_start][n_start-1:], colors=['black'], labels=['2020 sans covid (IC)'], line_style='dashed', display_legend=True)
plt.plot(x[n_start-1:], env_sups[n_start][n_start-1:], colors=['black'], line_style='dashed')
plt.plot(x[n_start-1:], moys[n_start][n_start-1:], colors=['black'], line_style='dashed')
plt.xlabel("Période d'étude")
plt.ylabel('Ratio de décès journalier')
plt.title("Evolution du ration de décès : 2020 sans covid")
plt.ylim(min=0.7, max=1.5)
plt.show(display_toolbar=False)
On finit en calculant le nombre de décès excédentaires : on constate qu'en espérance, avec notre modèle, il y a eu 15 400 décès supplémentaires pendant la période, soit 1.6 fois le nombre de décès constatés en milieu hospitalier et Ehpad au 13 avril (9 558). En mars, il y a eu 6 800 décès supplémentaires pendant la période, soit 1.8 fois le nombre de décès constatés en milieu hospitalier au 31 mars (3 523).
n_start = 29
moy_post = np.cumsum(moys[n_start])
sup_post = np.sort(np.cumsum(y_sims[n_start], axis=1), axis=0)[N-N//20]
inf_post = np.sort(np.cumsum(y_sims[n_start], axis=1), axis=0)[N//20]
deces_realises = np.cumsum(y_20)
deces_excedentaires = (pd.DataFrame([(deces_realises - moy_post), (deces_realises - sup_post), (deces_realises - inf_post)]) * morts_par_an_esp[-1] / 365).astype('int64').T
deces_excedentaires.index = period
deces_excedentaires = deces_excedentaires.loc[deces_excedentaires.index >= date_debut + datetime.timedelta(days=n_start)]
path = os.path.join(folder_treated_data, 'deces_hop_ehpad.pkl')
deces_rapportes = pd.read_pickle(path)
deces_rapportes = deces_rapportes.loc[deces_excedentaires.index]
plt.figure(legend_location='top-left')
plt.plot(deces_excedentaires.index, deces_rapportes['dc off'], colors=['green'], labels=['Valeurs rapportées'], display_legend=True)
plt.plot(deces_excedentaires.index, deces_excedentaires[0], colors=['black'], labels=['Valeurs estimées (moy)'], display_legend=True)
plt.plot(deces_excedentaires.index, deces_excedentaires[1], colors=['black'], line_style='dashed', labels=['Valeurs estimées (IC)'], display_legend=True)
plt.plot(deces_excedentaires.index, deces_excedentaires[2], colors=['black'], line_style='dashed')
plt.title("Nombre de décès excédentaires dus à l'épidémie")
plt.show()
Remarque 1 : la taille de cet intervalle de confiance semble importante. On peut la comparer à la variabilité usuelle du nombre de décès. Pour cela, on trace les ratios de décès historiques cumulés sur la période considérée pour 2000-2019, auquel on superpose notre intervalle de confiance pour 2020. On voit que sa taille est cohérente par rapport à la variabilité historique.
historic = np.sum(y_0019, axis=1)
plt.figure(legend_location='bottom-left')
plt.scatter(np.arange(2000, 2020), historic, labels=['historique'])
plt.hline(moy_post[-1], stroke_width=2, colors=['orangered'], labels=['moyenne 2020'], display_legend=True)
plt.hline(sup_post[-1], stroke_width=2, colors=['orangered'], line_style='dashed', labels=['IC 2020'], display_legend=True)
plt.hline(inf_post[-1], stroke_width=2, colors=['orangered'], line_style='dashed')
plt.ylim(min=49, max=69)
plt.title("Ratio de décès cumulés sur la période d'étude")
plt.show()
Remarque 2 : le jour du début de la divergence entre "2020 avec covid" et "2020 sans covid" est choisie de manière arbitraire, car elle semble correspondre au début de l'augmentation significative de la mortalité. Si on choisit des dates antérieures, on obtient des résultats similaires (voir plot ci-dessous), avec des intervalles de confiance plus larges.
plt.figure()
x = pd.to_datetime(period)
plt.plot(x, y_20, colors=['red'], labels=['Courbe réelle'])
for n_start in n_starts:
plt.plot(x[n_start-1:], env_infs[n_start][n_start-1:], colors=['gray'], labels=['Enveloppe inférieure'], line_style='dashed')
plt.plot(x[n_start-1:], moys[n_start][n_start-1:], colors=['black'])
plt.plot(x[n_start-1:], env_sups[n_start][n_start-1:], colors=['gray'], labels=['Enveloppe supérieure'], line_style='dashed')
plt.xlabel('Jour depuis le début de la période')
plt.ylabel('Ratio de décès')
plt.title('Intervalles de confiance partant de trois dates pendant la période')
plt.ylim(min=0.7, max=1.5)
plt.show(display_toolbar=False)
Pour l'illustration, on fit des processus gaussiens sur les données de 2000-2019.
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
from bqplot import LinearScale, Scatter, Lines, Axis, Figure
x_sc = DateScale()
y_sc = LinearScale()
ax_x = Axis(label='X', scale=x_sc, grid_lines='solid', tick_values = pd.to_datetime([date_debut + datetime.timedelta(days=k) for k in range(0, (date_fin-date_debut).days+1, 15)]))
ax_y = Axis(label='Y', scale=y_sc, orientation='vertical', grid_lines='solid')
scatters = []
lines = []
for k in range(20):
# Instantiate a Gaussian Process model
kernel = RBF(80, (4e1, 1e2))
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9, alpha=8*10**-4)
# Fit to data using Maximum Likelihood Estimation of the parameters
x = np.array(range(duree_periode)).reshape(-1, 1)
gp.fit(x, y_0019[k])
y_pred = gp.predict(x)
scatters.append(Scatter(x=pd.to_datetime(period), y=y_0019[k], scales={'x': x_sc, 'y': y_sc}, default_size=32))
lines.append(Lines(x=pd.to_datetime(period), y=y_pred, scales={'x': x_sc, 'y': y_sc}))
figy=[]
for i in range(5):
figx=[]
for j in range(4):
figx.append(Figure(axes=[ax_x, ax_y], marks=[scatters[i*4+j], lines[i*4+j]], title=str(2000 + i*4+j), fig_margin = dict(top=30, bottom=20, left=20, right=20), layout={'width':"240px",'height':'200px'}))
figy.append(ipw.HBox(figx))
display(ipw.VBox(figy, align_content = 'stretch'))