[CovsirPhy] Package Python COVID-19 pour l'analyse des données: estimation des paramètres

Introduction

Nous créons un package Python CovsirPhy qui vous permet de télécharger et d'analyser facilement les données COVID-19 (comme le nombre de PCR positifs).

Article d'introduction:

** Cette fois, nous aimerions introduire l'estimation des paramètres (estimation des paramètres comme le modèle SIR-F). ** **

La version anglaise du document est CovsirPhy: COVID-19 analysis with phase-depend SIRs, [Kaggle: COVID-19 data with SIR model]( Veuillez vous référer à https://www.kaggle.com/lisphilar/covid-19-data-with-sir-model).

1. Environnement d'exécution

CovsirPhy peut être installé par la méthode suivante! Veuillez utiliser Python 3.7 ou supérieur, ou Google Colaboratory.

--Version stable: pip install covsirphy --upgrade --Version de développement: pip install" git + https://github.com/lisphilar/covid19-sir.git#egg=covsirphy "

import covsirphy as cs
cs.__version__
# '2.8.2'
Environnement d'exécution
OS Windows Subsystem for Linux / Ubuntu
Python version 3.8.5

2. Préparation

Les tableaux et graphiques de cet article ont été créés à partir des données du 12/09/2020. Cliquez ici pour le code [^ 2] pour télécharger les données réelles depuis COVID-19 Data Hub [^ 1]:

data_loader = cs.DataLoader("input")
jhu_data = data_loader.jhu()
population_data = data_loader.population()

Utilisez également le code ci-dessous pour vérifier les données réelles et effectuer une analyse de tendance S-R [^ 3].

[^ 3]: [CovsirPhy] COVID-19 Python Package for Data Analysis: S-R trend analysis

# (Optional)Acquisition de données auprès du ministère de la Santé, du Travail et des Affaires sociales
japan_data = data_loader.japan()
jhu_data.replace(japan_data)
print(japan_data.citation)
#Génération d'instance de classe pour l'analyse
snl = cs.Scenario(jhu_data, population_data, country="Japan")
#Confirmation des données réelles
snl.records(filename=None)
# S-R trend analysis
snl.trend(filename=None)
#Vérifier les paramètres de phase
snl.summary()

Graphique de données réel: records.jpg

S-R trend analysis: trend.jpg

Réglage de phase:

Type Start End Population
0th Past 06Feb2020 21Apr2020 126529100
1st Past 22Apr2020 04Jul2020 126529100
2nd Past 05Jul2020 23Jul2020 126529100
3rd Past 24Jul2020 01Aug2020 126529100
4th Past 02Aug2020 14Aug2020 126529100
5th Past 15Aug2020 29Aug2020 126529100
6th Past 30Aug2020 12Sep2020 126529100

3. Exemple d'exécution

Par analyse de tendance S-R, il était possible de diviser en "Phase" (période où le paramètre est constant). Dans cet article, la valeur du paramètre est estimée à partir des données de chaque «Phase» (par exemple, les données du 2020/2/6 --2020 / 4/21 dans la 0ème phase).

J'écrirai un autre article sur le mécanisme d'estimation. Les valeurs des paramètres sont proposées en utilisant le package ʻoptuna, la solution numérique est calculée par scipy.integrate.solve_ivp`, et le jeu de paramètres avec peu d'erreur à partir des données réelles est sélectionné.

Comment lire le résultat sera décrit plus tard, mais vous pouvez exécuter et obtenir la liste des résultats dans les deux lignes suivantes. Cette fois, j'ai utilisé le modèle SIR-F [^ 4].

[^ 4]: [CovsirPhy] Package Python COVID-19 pour l'analyse des données: modèle SIR-F

# Parameter estimation with SIR-F model
snl.estimate(cs.SIRF)
#Obtenir la liste des paramètres
snl.summary()
#Exemple de sortie standard (le temps de traitement varie en fonction du nombre de processeurs, etc.)
#Détails ci-dessous: Dernière phase=Après avoir estimé incluant tau dans la 6ème phase, fixez tau et 0-Estimation du 5e paramètre
<SIR-F model: parameter estimation>
Running optimization with 8 CPUs...
        6th phase (30Aug2020 - 12Sep2020): finished  704 trials in 0 min 25 sec
        5th phase (15Aug2020 - 29Aug2020): finished  965 trials in 1 min  0 sec
        3rd phase (24Jul2020 - 01Aug2020): finished  965 trials in 1 min  0 sec
        1st phase (22Apr2020 - 04Jul2020): finished  913 trials in 1 min  0 sec
        4th phase (02Aug2020 - 14Aug2020): finished  969 trials in 1 min  0 sec
        0th phase (06Feb2020 - 21Apr2020): finished  853 trials in 1 min  0 sec
        2nd phase (05Jul2020 - 23Jul2020): finished  964 trials in 1 min  0 sec
Completed optimization. Total: 1 min 27 sec
Type Start End Population ODE Rt theta kappa rho sigma tau 1/alpha2 [day] 1/gamma [day] alpha1 [-] 1/beta [day] RMSLE Trials Runtime
0th Past 06Feb2020 21Apr2020 126529100 SIR-F 5.54 0.0258495 0.0002422 0.0322916 0.00543343 480 1376 61 0.026 10 1.17429 804 1 min 0 sec
1st Past 22Apr2020 04Jul2020 126529100 SIR-F 0.41 0.0730748 0.000267108 0.0118168 0.0264994 480 1247 12 0.073 28 1.11459 861 1 min 0 sec
2nd Past 05Jul2020 23Jul2020 126529100 SIR-F 2.01 0.000344333 7.92419e-05 0.0467789 0.023201 480 4206 14 0 7 0.0331522 910 1 min 0 sec
3rd Past 24Jul2020 01Aug2020 126529100 SIR-F 1.75 0.00169155 4.05087e-05 0.0459332 0.0260965 480 8228 12 0.002 7 0.0201773 923 1 min 0 sec
4th Past 02Aug2020 14Aug2020 126529100 SIR-F 1.46 0.000634554 0.000116581 0.0325815 0.0221345 480 2859 15 0.001 10 0.0751473 909 1 min 0 sec
5th Past 15Aug2020 29Aug2020 126529100 SIR-F 0.8 0.00244294 9.30884e-05 0.0272693 0.0337857 480 3580 9 0.002 12 0.0420563 907 1 min 0 sec
6th Past 30Aug2020 12Sep2020 126529100 SIR-F 0.69 5.36037e-05 0.000467824 0.0219379 0.0312686 480 712 10 0 15 0.0132161 724 0 min 25 sec

4. Estimations des paramètres

C'est long sur le côté, alors regardons-le dans l'ordre. La première est la valeur estimée du paramètre.

# cs.SIRF.PARAMETERS: SIR-Liste des noms de paramètres de modèle F
cols = ["Start", "End", "ODE", "tau", *cs.SIRF.PARAMETERS]
snl.summary(columns=cols)
Start End ODE tau theta kappa rho sigma
0th 06Feb2020 21Apr2020 SIR-F 480 0.0258495 0.0002422 0.0322916 0.00543343
1st 22Apr2020 04Jul2020 SIR-F 480 0.0730748 0.000267108 0.0118168 0.0264994
2nd 05Jul2020 23Jul2020 SIR-F 480 0.000344333 7.92419e-05 0.0467789 0.023201
3rd 24Jul2020 01Aug2020 SIR-F 480 0.00169155 4.05087e-05 0.0459332 0.0260965
4th 02Aug2020 14Aug2020 SIR-F 480 0.000634554 0.000116581 0.0325815 0.0221345
5th 15Aug2020 29Aug2020 SIR-F 480 0.000644851 0.000383424 0.0274104 0.0337156
6th 30Aug2020 12Sep2020 SIR-F 480 5.36037e-05 0.000467824 0.0219379 0.0312686

SIR-F model[^4]:

\begin{align*}
\mathrm{S} \overset{\beta I}{\longrightarrow} \mathrm{S}^\ast \overset{\alpha_1}{\longrightarrow}\ & \mathrm{F}    \\
\mathrm{S}^\ast \overset{1 - \alpha_1}{\longrightarrow}\ & \mathrm{I} \overset{\gamma}{\longrightarrow} \mathrm{R}    \\
& \mathrm{I} \overset{\alpha_2}{\longrightarrow} \mathrm{F}    \\
\end{align*}

$ \ Alpha_2 $, $ \ beta $, $ \ gamma $ a "temps" comme dimension. Ce "temps" est $ f (T + \ Delta T) --f (T) = x (T) \ Delta T $, qui est une équation discrète d'équations différentielles normales simultanées $ f '(T) = x (T) $. Cela dépend de \ Delta T $. Puisque $ \ Delta T $ est acquis tous les deux jours, c'est moins d'un jour (= 1440 min), mais je ne connais pas la valeur spécifique. Cela peut être différent selon le pays ou la région, et il est difficile de comparer les paramètres dimensionnels $ \ alpha_2 $, $ \ beta $, $ \ gamma $ de différents pays.

Par conséquent, nous définissons une variable $ \ tau $ correspondant à $ \ Delta T $ pour rendre les paramètres dimensionnels non dimensionnels.

\begin{align*}
(S, I, R, F) = & N \times (x, y, z, w)    \\
(T, \alpha_1, \alpha_2, \beta, \gamma) = & (\tau t, \theta, \tau^{-1}\kappa, \tau^{-1}\rho, \tau^{-1}\sigma)    \\
1 \leq \tau & \leq 1440    \\
\end{align*}

À ce moment [^ 4],

\begin{align*}
& 0 \leq (x, y, z, w, \theta, \kappa, \rho, \sigma) \leq 1  \\
\end{align*}

Après avoir estimé $ \ tau $ et les paramètres non dimensionnels en utilisant les dernières données de Phase afin que la valeur de $ \ tau $ soit constante dans le même pays, la même valeur est définie sur $ \ tau $ lors de l'estimation des paramètres des autres Phases. J'essaye de l'utiliser comme. Cette fois, c'était 480 minutes comme indiqué dans le tableau. Je programme pour utiliser une fraction de 1440 minutes par jour pour simplifier le calcul.

Maintenant, représentons la transition de chaque paramètre sans dimension.

\begin{align*}
& \frac{\mathrm{d}x}{\mathrm{d}t}= - \rho x y  \\
& \frac{\mathrm{d}y}{\mathrm{d}t}= \rho (1-\theta) x y - (\sigma + \kappa) y  \\
& \frac{\mathrm{d}z}{\mathrm{d}t}= \sigma y  \\
& \frac{\mathrm{d}w}{\mathrm{d}t}= \rho \theta x y + \kappa y  \\
\end{align*}

transition de rho

Lorsque Susceptible entre en contact avec Infected, la probabilité d'être infecté $ \ rho $:

snl.history(target="rho", filename=None)

rho.jpg

Il est interprété que l'effet de la déclaration d'urgence (zone limitée 2020/4/7, national 2020/4/16, annulation nationale 2020/5/25) est apparu dans la seconde moitié d'avril et a été maintenu jusqu'au début du mois de juillet. Je vais. Après cela, il a bondi et a progressivement diminué. Les détails doivent être discutés, mais c'est un paramètre qui reflète directement les effets de mesures telles que 3 évitement dense.

Transition de sigma

Transition de la probabilité $ \ sigma $ de transition d'Infecté à Récupéré:

snl.history(target="sigma", filename=None)

sigma.jpg

Après une forte hausse dans la seconde moitié du mois d'avril, il a suivi une tendance à la hausse tout en répétant des hausses / baisses. Ce paramètre reflète le système de fourniture de soins médicaux et l'état de développement / fourniture de nouveaux médicaments.

Transition de kappa

Changements dans le taux de mortalité des personnes infectées $ \ kappa $:

snl.history(target="kappa", filename=None)

kappa.jpg

Il semble qu'elle fluctue fortement, mais comme la valeur absolue de la valeur est faible, on peut considérer qu'elle est maintenue dans une certaine mesure constante (une vérification par comparaison avec d'autres pays est nécessaire). Il faut mettre en place un système médical et amener $ \ kappa $ le plus près possible de 0 avec un approvisionnement suffisant en nouveaux médicaments.

Transition de thêta

Pourcentage de personnes infectées ayant reçu un diagnostic définitif qui sont décédées au moment du diagnostic définitif $ \ theta $ Transition:

snl.history(target="theta", filename=None)

theta.jpg

Aux premiers stades de l'infection, le système de fourniture de soins médicaux lui-même était extrêmement serré, il est donc difficile à interpréter, mais on pense que la valeur augmentera si l'examen est retardé et qu'un traitement approprié ne peut être effectué.

5. Transition des paramètres dimensionnels

Les paramètres dimensionnels sont également inclus dans le tableau. Pour faciliter l'interprétation, sauf pour le $ \ alpha_1 $ initialement sans dimension, il est inversé et l'unité est [jour].

cols = ["Start", "End", "ODE", "tau", *cs.SIRF.DAY_PARAMETERS]
fh.write(snl.summary(columns=cols).to_markdown())
Start End ODE tau alpha1 [-] 1/alpha2 [day] 1/beta [day] 1/gamma [day]
0th 06Feb2020 21Apr2020 SIR-F 480 0.026 1376 10 61
1st 22Apr2020 04Jul2020 SIR-F 480 0.073 1247 28 12
2nd 05Jul2020 23Jul2020 SIR-F 480 0 4206 7 14
3rd 24Jul2020 01Aug2020 SIR-F 480 0.002 8228 7 12
4th 02Aug2020 14Aug2020 SIR-F 480 0.001 2859 10 15
5th 15Aug2020 29Aug2020 SIR-F 480 0.002 3580 12 9
6th 30Aug2020 12Sep2020 SIR-F 480 0 712 15 10

Le paramètre dimensionnel est omis car il montre la même progression que le paramètre non dimensionnel, mais le graphique peut être obtenu avec le code suivant.

#Pour la version bêta->Graphique omis
snl.history(target="1/beta [day]", filename=None)

6. Modifications du nombre de reproductions effectives

Le nombre de reproduction de base / effective Rt du modèle SIR-F [^ 4] est défini comme suit.

\begin{align*}
R_t = \rho (1 - \theta) (\sigma + \kappa)^{-1} = \beta (1 - \alpha_1) (\gamma + \alpha_2)^{-1}
\end{align*}

Transition:

#liste
cols = ["Start", "End", "ODE", "tau", "Rt"]
snl.summary(columns=cols)
#Graphique
snl.history(target="Rt", filename="rt.jpg ")
Start End ODE Rt
0th 06Feb2020 21Apr2020 SIR-F 5.54
1st 22Apr2020 04Jul2020 SIR-F 0.41
2nd 05Jul2020 23Jul2020 SIR-F 2.01
3rd 24Jul2020 01Aug2020 SIR-F 1.75
4th 02Aug2020 14Aug2020 SIR-F 1.46
5th 15Aug2020 29Aug2020 SIR-F 0.8
6th 30Aug2020 12Sep2020 SIR-F 0.69

rt.jpg

Puisque $ Rt> 1 $ est un indicateur de la propagation de l'infection, la ligne horizontale $ Rt = 1 $ est affichée.

7. Précision des paramètres

Le tableau comprend également le score RMSLE, qui mesure la précision des paramètres, le nombre de jeux de paramètres proposés par le package ʻoptuna` pour estimer les paramètres, et le temps d'exécution.

cols = ["Start", "End", "RMSLE", "Trials", "Runtime"]
snl.summary(columns=cols)
Start End RMSLE Trials Runtime
0th 06Feb2020 21Apr2020 1.17429 690 1 min 1 sec
1st 22Apr2020 04Jul2020 1.11459 764 1 min 0 sec
2nd 05Jul2020 23Jul2020 0.0331522 810 1 min 1 sec
3rd 24Jul2020 01Aug2020 0.0201773 816 1 min 1 sec
4th 02Aug2020 14Aug2020 0.0751473 808 1 min 0 sec
5th 15Aug2020 29Aug2020 0.0420563 804 1 min 0 sec
6th 30Aug2020 12Sep2020 0.0132161 658 0 min 25 sec

La formule de définition du score RMSLE (Root Mean Squared Log Error) [^ 5] est la suivante. On peut dire que plus il est proche de 0, meilleure est la réflexion des données réelles. Bien que omise, la vérification de la méthode d'estimation elle-même est effectuée à l'aide de données théoriques (données sur le nombre de patients théoriquement créées à partir de la formule du modèle SIR-F et de l'exemple de jeu de paramètres).

\begin{align*}
& \sqrt{\cfrac{1}{n}\sum_{i=1}^{n}(log_{10}(A_{i} + 1) - log_{10}(P_{i} + 1))^2}
\end{align*}

$ A $ sont les données réelles et $ P $ est la valeur prévue. Lorsque $ i = 1, 2, 3, 4 (= n) $, $ A_i $ et $ P_i $ sont les données réelles / les valeurs prédites de $ S, I, R, F $, respectivement.

Comme il est difficile d'obtenir une image avec juste la valeur, j'ai fait un graphique. Tout d'abord, à propos de la phase 0, qui a la plus grande valeur RMSLE. Le graphique du haut montre la différence entre les données réelles et la valeur prédite, et les deuxième, troisième et quatrième montrent à la fois les données réelles et la valeur prédite pour chaque variable.

snl.estimate_accuracy(phase="0th", filename=None)

accuracy_0th.jpg

Il y a une erreur. Il semble préférable de scinder la phase 0 en utilisant Scenario.separate () etc.

En revanche, dans la 6ème phase, qui a le score RMSLE le plus bas, les données réelles et la valeur prédite se chevauchent souvent.

snl.estimate_accuracy(phase="6th", filename=None)

accuracy_6th.jpg

8. Post-scriptum

Cette fois, j'ai expliqué comment estimer les paramètres de chaque phase. Plus de six mois se sont écoulés depuis le début de l'épidémie, et je pense que nous sommes à un stade où la vérification des paramètres, la prédiction des paramètres futurs et l'analyse de scénarios sont importants. Il semble que l'attention du COVID-19 en tant que thème de la science des données diminue, mais au fur et à mesure que les données s'accumulent, il est devenu possible d'effectuer une analyse plus approfondie que dans les premiers jours lorsque nous nous sommes concentrés sur la création de tableaux de bord. C'était.

Merci pour votre travail acharné cette fois aussi!

Recommended Posts

[CovsirPhy] Package Python COVID-19 pour l'analyse des données: estimation des paramètres
[CovsirPhy] Package Python COVID-19 pour l'analyse des données: chargement des données
[CovsirPhy] Package Python COVID-19 pour l'analyse de données: modèle SIR-F
[CovsirPhy] Package Python COVID-19 pour l'analyse des données: analyse des tendances S-R
[CovsirPhy] Package Python COVID-19 pour l'analyse des données: modèle SIR
Python pour l'analyse des données Chapitre 4
Python pour l'analyse des données Chapitre 2
Python pour l'analyse des données Chapitre 3
Modèle de prétraitement pour l'analyse des données (Python)
Outil de visualisation Python pour le travail d'analyse de données
Analyse de données python
Analyse de données avec python 2
Analyse de données avec Python
Analysons les données Covid-19 (Corona) en utilisant Python [Pour les débutants]
[Pour les débutants] Comment étudier le test d'analyse de données Python3
Mon conteneur d'analyse de données python
[Python] Notes sur l'analyse des données
Notes d'apprentissage sur l'analyse des données Python
Analyse de données à l'aide de pandas python
Conseils et précautions lors de l'analyse des données
Lequel dois-je étudier, R ou Python, pour l'analyse des données?
<Python> Construisez un serveur dédié pour l'analyse des données Jupyter Notebook
Implémentation de l'estimation des paramètres HMM en python
Cours Python pour la science des données - techniques utiles
Pratique de l'analyse de données par Python et pandas (Tokyo COVID-19 data edition)
Analyse de données pour améliorer POG 3 ~ Analyse de régression ~
Formatage des données pour les graphiques Python / couleur
Analyse de données à partir de python (visualisation de données 1)
Analyse de données à partir de python (visualisation de données 2)
Créer un Ubuntu de démarrage USB avec un environnement Python pour l'analyse des données
Résumé du livre électronique Python utile pour l'analyse de données gratuite
Techniques Python détaillées requises pour la mise en forme des données (1)
[Python] Première analyse de données / apprentissage automatique (Kaggle)
Analyse de données à partir de python (pré-traitement des données-apprentissage automatique)
Comment utiliser "deque" pour les données Python
Techniques Python détaillées requises pour la mise en forme des données (2)
J'ai suivi une formation à l'analyse de données Python à distance
Préparation à l'examen d'analyse de données certifié Python 3 Engineer
JupyterLab Basic Setting 2 pour l'analyse des données (pip)
Configuration de base de JupyterLab pour l'analyse des données (pip)
Analyse des données en Python Résumé des sources que les débutants devraient d'abord consulter
Analyse des données pour améliorer POG 2 ~ Analyse avec le notebook jupyter ~
Modèle Python qui effectue une analyse des journaux à une vitesse explosive
Préparer un environnement de langage de programmation pour l'analyse des données
Formation préalable à l'examen d'analyse des données de certification d'ingénieur Python 3
Introduction à la modélisation statistique pour l'analyse des données
[Python] Analyse de données, pratique du machine learning (Kaggle) -Prétraitement des données-
Comment utiliser les outils d'analyse de données pour les débutants
Afficher la bougie de données FX (forex) en Python
Analyse de données en Python: une note sur line_profiler
[Python] Flux du scraping Web à l'analyse des données
Environnement enregistré pour l'analyse des données avec Python
Astro: modules / fonctions python fréquemment utilisés pour l'analyse
[Mis à jour de temps en temps] Mémos Python souvent utilisés pour l'analyse des données [Division N, etc.]