[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:

• SIR model
• SIR-F model
• Data loading
• S-R trend analysis

** 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'

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:

S-R trend analysis:

Réglage de phase:

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

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)

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)

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)

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)

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)

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())

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 ")

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)

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)

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)

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!