Représentation du R0 par départements

Cet exemple reprend des données de tests par gouvernement pour calculer le coefficient R0 de l’épidémie. La méthode de Cori est utilisée avec la formule proposée dans covidtracker :

\[R = \frac{\sum_{i=t-6}^{t}C_i}{\sum_{i=t-13}^{t-7}C_i}\]

\(C_i\) est le nombre de cas positifs du jour i. Cette méthode est implémentée dans le package R EpiEstim ou epyestim pour le langage python. Dans cet exemple, il sera calculé directement à partir des données.

Sources de données:

Récupération des données

import warnings
from pandas import Timedelta, DataFrame
from geopandas import GeoDataFrame
from matplotlib.cbook.deprecation import MatplotlibDeprecationWarning
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from aftercovid.data import (
    data_covid_france_departments_tests,
    data_france_departments)

case = data_covid_france_departments_tests(metropole=True)
case = case[case.cl_age90 == 0]
case.tail()

Out:

somewhereaftercovid_39_std/aftercovid/aftercovid/data/data_insee.py:76: DtypeWarning: Columns (0) have mixed types.Specify dtype option on import or set low_memory=False.
  df = read_csv_cache(cache, url, sep=';')
dep jour P T cl_age90 pop
592371 95 2021-11-20 382 8722 0 1248354.0
592382 95 2021-11-21 67 1509 0 1248354.0
592393 95 2021-11-22 592 9656 0 1248354.0
592404 95 2021-11-23 559 9097 0 1248354.0
592415 95 2021-11-24 601 10335 0 1248354.0


Quelques aggrégations, par département.

deps = case.groupby(["dep", "jour"], as_index=False).sum()
deps.tail(n=10)
dep jour P T cl_age90 pop
53846 95 2021-11-15 503 7359 0 1248354.0
53847 95 2021-11-16 320 6610 0 1248354.0
53848 95 2021-11-17 386 7409 0 1248354.0
53849 95 2021-11-18 372 6975 0 1248354.0
53850 95 2021-11-19 420 9791 0 1248354.0
53851 95 2021-11-20 382 8722 0 1248354.0
53852 95 2021-11-21 67 1509 0 1248354.0
53853 95 2021-11-22 592 9656 0 1248354.0
53854 95 2021-11-23 559 9097 0 1248354.0
53855 95 2021-11-24 601 10335 0 1248354.0


Sur tout le territoire.

france = case.groupby(["jour"], as_index=False).sum()
france.tail(n=10)
jour P T cl_age90 pop
551 2021-11-15 22717 433115 0 64897954.0
552 2021-11-16 18547 411500 0 64897954.0
553 2021-11-17 19542 384603 0 64897954.0
554 2021-11-18 20789 422225 0 64897954.0
555 2021-11-19 22618 539677 0 64897954.0
556 2021-11-20 16814 396185 0 64897954.0
557 2021-11-21 4036 72917 0 64897954.0
558 2021-11-22 35701 584015 0 64897954.0
559 2021-11-23 30623 574411 0 64897954.0
560 2021-11-24 32307 562650 0 64897954.0


Calcul du R

def compute_r(df, col_day='jour', col='P', last_day=None):
    if last_day is None:
        last_day = df.jour.max()
    v1 = last_day - Timedelta(days=6)
    v2 = last_day
    p1 = last_day - Timedelta(days=13)
    p2 = last_day - Timedelta(days=7)
    w1 = df[(df[col_day] >= p1) & (df[col_day] <= p2)]
    w2 = df[(df[col_day] >= v1) & (df[col_day] <= v2)]
    return w2[col].sum() / w1[col].sum()


compute_r(france)

Out:

1.7542945148679066

On regarde quelle tête ça a sur six mois. Ce n’est pas le code le plus efficace mais c’est rapide à écrire. Dans le cas idéal, il faudra s’assurer que toutes les dates sont présents et les compléter le cas échéants puis calculer l’estimateur sur des fenêtres glissantes.

last_day = france.jour.max()
obs = []
for i in range(0, 180):
    ld = last_day - Timedelta(days=i)
    obs.append({'jour': ld, 'R': compute_r(france, last_day=ld)})

gr = DataFrame(obs).sort_values("jour")
gr.tail()
jour R
4 2021-11-20 1.843932
3 2021-11-21 1.837378
2 2021-11-22 1.793727
1 2021-11-23 1.768046
0 2021-11-24 1.754295


Et on dessine.

gr.set_index('jour').plot(title="Evolution de R en Métropole")
Evolution de R en Métropole

Carte du R par département

On calcule les R par départements.

obs = []
for d in set(deps.dep):
    r = compute_r(deps[deps.dep == d])
    obs.append({'dep': d, 'R': r})

depdf = DataFrame(obs).sort_values("dep")
depdf.tail()
dep R
4 91 1.685424
41 92 1.856895
17 93 1.799901
7 94 1.786769
14 95 1.670201


On récupère les contours des départements français.

loc = data_france_departments(metropole=True)
loc.tail()
ancienne_re code_ancien code_depart code_region departement region geometry
96 Limousin 74.0 87 75.0 Haute-Vienne Nouvelle-Aquitaine POLYGON ((0.92278 45.94610, 0.92856 45.94808, ...
97 Auvergne 83.0 63 84.0 Puy-de-Dôme Auvergne-Rhône-Alpes POLYGON ((2.45492 45.76131, 2.45169 45.76137, ...
98 Basse-Normandie 25.0 14 28.0 Calvados Normandie MULTIPOLYGON (((-0.42990 48.86405, -0.42881 48...
99 Rhône-Alpes 82.0 07 84.0 Ardèche Auvergne-Rhône-Alpes POLYGON ((4.30746 44.98596, 4.30412 44.98837, ...
100 Midi-Pyrénées 73.0 32 76.0 Gers Occitanie POLYGON ((-0.24307 43.66404, -0.24328 43.66476...


On fusionne avec les R.

locdep = GeoDataFrame(depdf.merge(loc, left_on='dep', right_on='code_depart'))
locdep.tail()
dep R ancienne_re code_ancien code_depart code_region departement region geometry
91 91 1.685424 Île-de-France 11.0 91 11.0 Essonne Île-de-France POLYGON ((2.04261 48.62638, 2.04045 48.62754, ...
92 92 1.856895 Île-de-France 11.0 92 11.0 Hauts-de-Seine Île-de-France POLYGON ((2.22040 48.92061, 2.23114 48.92773, ...
93 93 1.799901 Île-de-France 11.0 93 11.0 Seine-Saint-Denis Île-de-France POLYGON ((2.45949 48.95505, 2.46072 48.95840, ...
94 94 1.786769 Île-de-France 11.0 94 11.0 Val-de-Marne Île-de-France POLYGON ((2.32906 48.81378, 2.33190 48.81701, ...
95 95 1.670201 Île-de-France 11.0 95 11.0 Val-d'Oise Île-de-France POLYGON ((2.43548 49.13411, 2.44084 49.13419, ...


Et on dessine. Les départements en vert sont ceux pour lequel le R est > 1.

with warnings.catch_warnings():
    warnings.filterwarnings("ignore", category=UserWarning)
    warnings.filterwarnings("ignore", category=MatplotlibDeprecationWarning)
    fig, axs = plt.subplots(
        1, 2, figsize=(16, 10),
        gridspec_kw={'width_ratios': [2, 1]})

    # métropole
    ax = axs[0]
    cax = make_axes_locatable(ax).append_axes("right", size="5%", pad=0.1)
    locdep.plot(
        column="R", ax=ax, edgecolor='black',
        legend=True, cax=cax, cmap="OrRd")
    if (locdep.R < 1).sum() > 0:
        locdep[locdep.R < 1].geometry.boundary.plot(
            color=None, edgecolor='g', linewidth=2, ax=ax, label="R<1")
    ax.set_title("R par départments de la métropole\n%r" % deps.jour.max())

    for _, row in locdep.iterrows():
        p = row['geometry'].representative_point()
        ax.annotate("%1.1f" % row['R'], xy=(p.x, p.y),
                    horizontalalignment='center', color="black", fontsize=8)

    ax.legend()

    # Paris et sa région
    idf = set(['75', '77', '78', '91', '92', '93', '94', '95'])
    ax = axs[1]
    locdep2 = locdep[locdep.dep.isin(idf)]
    cax = make_axes_locatable(ax).append_axes("right", size="5%", pad=0.1)
    locdep2.plot(
        column="R", ax=ax, edgecolor='black',
        legend=True, cax=cax, cmap="OrRd")
    ax.set_title("R par départments de la métropole\n%r" % deps.jour.max())

    for _, row in locdep2.iterrows():
        p = row['geometry'].representative_point()
        ax.annotate("%1.1f" % row['R'], xy=(p.x, p.y),
                    horizontalalignment='center', color="black", fontsize=8)

plt.show()
R par départments de la métropole Timestamp('2021-11-24 00:00:00'), R par départments de la métropole Timestamp('2021-11-24 00:00:00')

Total running time of the script: ( 0 minutes 20.203 seconds)

Gallery generated by Sphinx-Gallery