[CovsirPhy] COVID-19 Python-Paket für die Datenanalyse: Parameterschätzung

Introduction

Wir erstellen ein Python-Paket CovsirPhy, mit dem Sie COVID-19-Daten (z. B. die Anzahl der PCR-Positiven) einfach herunterladen und analysieren können.

Einführungsartikel:

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

** Dieses Mal möchten wir die Parameterschätzung einführen (Parameterschätzung wie das SIR-F-Modell). ** ** **

Die englische Version des Dokuments lautet CovsirPhy: COVID-19-Analyse mit phasenabhängigen SIRs, [Kaggle: COVID-19-Daten mit SIR-Modell]( Weitere Informationen finden Sie unter https://www.kaggle.com/lisphilar/covid-19-data-with-sir-model.

1. Ausführungsumgebung

CovsirPhy kann mit der folgenden Methode installiert werden! Bitte verwenden Sie Python 3.7 oder höher oder Google Colaboratory.

--Stabile Version: pip install covsirphy --upgrade

• Entwicklungsversion: pip install" git + https://github.com/lisphilar/covid19-sir.git#egg=covsirphy "

import covsirphy as cs
cs.__version__
# '2.8.2'

2. Vorbereitung

Die Tabellen und Grafiken in diesem Artikel wurden mit Daten vom 12.9.2020 erstellt. Klicken Sie hier für den Code [^ 2], um die tatsächlichen Daten von COVID-19 Data Hub [^ 1] herunterzuladen:

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

Verwenden Sie außerdem den folgenden Code, um die tatsächlichen Daten zu überprüfen und eine S-R-Trendanalyse durchzuführen [^ 3].

[^ 3]: [CovsirPhy] COVID-19 Python-Paket für die Datenanalyse: S-R-Trendanalyse

# (Optional)Datenerfassung beim Ministerium für Gesundheit, Arbeit und Soziales
japan_data = data_loader.japan()
jhu_data.replace(japan_data)
print(japan_data.citation)
#Instanzgenerierung einer Klasse zur Analyse
snl = cs.Scenario(jhu_data, population_data, country="Japan")
#Bestätigung der tatsächlichen Daten
snl.records(filename=None)
# S-R trend analysis
snl.trend(filename=None)
#Überprüfen Sie die Phaseneinstellungen
snl.summary()

Aktueller Datengraph:

S-R trend analysis:

Phaseneinstellung:

3. Ausführungsbeispiel

Durch S-R-Trendanalyse konnte in "Phase" (Zeitraum, in dem der Parameter konstant ist) unterteilt werden. In diesem Artikel wird der Wert des Parameters anhand der Daten jeder "Phase" geschätzt (z. B. der Daten von 2020/2/6 --2020 / 4/21 in der 0. Phase).

Ich werde einen weiteren Artikel über den Schätzmechanismus schreiben. Wir schlagen Parameterwerte mit dem Paket "optuna" vor, berechnen numerische Lösungen mit "scipy.integrate.solve_ivp" und wählen aus den tatsächlichen Daten einen Parametersatz mit geringem Fehler aus.

Das Lesen des Ergebnisses wird später beschrieben. Sie können die Ergebnisliste jedoch in den folgenden beiden Zeilen ausführen und abrufen. Diesmal habe ich das SIR-F-Modell [^ 4] verwendet.

[^ 4]: [CovsirPhy] COVID-19 Python-Paket für die Datenanalyse: SIR-F-Modell

# Parameter estimation with SIR-F model
snl.estimate(cs.SIRF)
#Parameterliste abrufen
snl.summary()

#Beispiel für eine Standardausgabe (die Verarbeitungszeit variiert je nach Anzahl der CPUs usw.)
#Details unten: Letzte Phase=Nach der Schätzung, einschließlich Tau in der 6. Phase, fixieren Sie Tau und 0-5. Parameterschätzung
<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. Parameterschätzungen

Es ist seitwärts lang, also schauen wir es uns der Reihe nach an. Erstens ist der geschätzte Wert des Parameters.

# cs.SIRF.PARAMETERS: SIR-F Liste der Modellparameternamen
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 $ hat "Zeit" als Dimension. Diese "Zeit" ist $ f (T + \ Delta T) - f (T) = x (T) \ Delta T $, was eine diskrete Gleichung simultaner normaler Differentialgleichungen $ f '(T) = x (T) $ ist. Es hängt von \ Delta T $ ab. Da $ \ Delta T $ jeden zweiten Tag erworben wird, ist es weniger als ein Tag (= 1440 min), aber ich kenne den spezifischen Wert nicht. Es kann je nach Land oder Region unterschiedlich sein, und es ist schwierig, die Dimensionsparameter $ \ alpha_2 $, $ \ beta $, $ \ gamma $ verschiedener Länder zu vergleichen.

Daher setzen wir eine Variable $ \ tau $, die $ \ Delta T $ entspricht, um die Dimensionsparameter nicht dimensional zu machen.

\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*}

Zu diesem Zeitpunkt [^ 4],

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

Nach dem Schätzen von $ \ tau $ und nicht-dimensionalen Parametern unter Verwendung der neuesten Phasendaten, sodass der Wert von $ \ tau $ im selben Land konstant ist, wird der gleiche Wert beim Schätzen der Parameter anderer Phasen auf $ \ tau $ gesetzt. Ich versuche es als zu verwenden. Diesmal waren es 480 Minuten, wie in der Tabelle gezeigt. Ich programmiere, um einen Bruchteil von 1440 Minuten pro Tag zu verwenden, um die Berechnung zu vereinfachen.

Lassen Sie uns nun den Übergang jedes dimensionslosen Parameters grafisch darstellen.

\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*}

Übergang von Rho

Wenn Susceptible mit Infected in Kontakt kommt, ist die Wahrscheinlichkeit einer Infektion $ \ rho $:

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

Es wird interpretiert, dass die Wirkung der Notfallerklärung (Flächenbegrenzung 2020/4/7, bundesweit 2020/4/16, landesweite Löschung 2020/5/25) in der zweiten Aprilhälfte auftrat und bis Anfang Juli beibehalten wurde. Ich werde. Danach sprang es auf und nahm allmählich ab. Details müssen diskutiert werden, aber es ist ein Parameter, der die Auswirkungen von Maßnahmen wie 3 dichte Vermeidung direkt widerspiegelt.

Übergang von Sigma

Übergang von $ \ sigma $ Wahrscheinlichkeit des Übergangs von infiziert zu wiederhergestellt:

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

Nach einem starken Anstieg in der zweiten Aprilhälfte war der Aufwärtstrend rückläufig, während sich das Steigen / Fallen wiederholte. Dieser Parameter spiegelt das medizinische Versorgungssystem und den Entwicklungs- / Lieferstatus neuer Medikamente wider.

Übergang von Kappa

Veränderungen der Sterblichkeitsrate infizierter Menschen $ \ kappa $:

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

Es scheint, dass es stark schwankt, aber da der absolute Wert des Wertes klein ist, kann davon ausgegangen werden, dass er bis zu einem gewissen Grad konstant gehalten wird (eine Überprüfung durch Vergleich mit anderen Ländern ist erforderlich). Es ist notwendig, ein medizinisches System einzurichten und $ \ kappa $ mit einer ausreichenden Menge neuer Medikamente so nahe wie möglich an 0 zu bringen.

Übergang von Theta

Prozentsatz der infizierten Personen, die eine endgültige Diagnose erhalten haben und zum Zeitpunkt der endgültigen Diagnose verstorben sind $ \ theta $ Übergang:

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

In den frühen Stadien der Infektion war das medizinische Versorgungssystem selbst extrem eng, so dass es schwierig zu interpretieren ist, aber es wird angenommen, dass sich der Wert erhöht, wenn die Untersuchung verzögert wird und keine angemessene Behandlung durchgeführt werden kann.

5. Übergang von Maßparametern

Die Maßparameter sind ebenfalls in der Tabelle enthalten. Zur Erleichterung der Interpretation wird es mit Ausnahme des ursprünglich dimensionslosen $ \ alpha_1 $ invertiert und die Einheit ist [Tag].

cols = ["Start", "End", "ODE", "tau", *cs.SIRF.DAY_PARAMETERS]
fh.write(snl.summary(columns=cols).to_markdown())

Der Dimensionsparameter wird weggelassen, da er den gleichen Fortschritt wie der nicht-dimensionale Parameter zeigt, der Graph jedoch mit dem folgenden Code erhalten werden kann.

#Für die Beta->Grafik weggelassen
snl.history(target="1/beta [day]", filename=None)

6. Änderungen in der Anzahl der effektiven Reproduktionen

Die grundlegende / effektive Reproduktionsnummer Rt des SIR-F-Modells [^ 4] ist wie folgt definiert.

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

Überleitung:

#Aufführen
cols = ["Start", "End", "ODE", "tau", "Rt"]
snl.summary(columns=cols)
#Graph
snl.history(target="Rt", filename="rt.jpg ")

Da $ Rt> 1 $ ein Indikator für die Ausbreitung der Infektion ist, wird die horizontale Linie $ Rt = 1 $ angezeigt.

7. Parametergenauigkeit

Die Tabelle enthält auch den RMSLE-Score, der die Genauigkeit der Parameter, die Anzahl der vom "optuna" -Paket zur Schätzung der Parameter vorgeschlagenen Parametersätze und die Ausführungszeit misst.

cols = ["Start", "End", "RMSLE", "Trials", "Runtime"]
snl.summary(columns=cols)

Die Definitionsformel des RMSLE-Scores (Root Mean Squared Log Error) [^ 5] lautet wie folgt. Man kann sagen, je näher es an 0 liegt, desto besser werden die tatsächlichen Daten reflektiert. Obwohl weggelassen, wird die Überprüfung der Schätzmethode selbst unter Verwendung theoretischer Daten durchgeführt (Daten zur Anzahl der Patienten, die theoretisch aus der SIR-F-Modellformel und dem Beispiel eines Parametersatzes erstellt wurden).

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

$ A $ sind die tatsächlichen Daten und $ P $ ist der vorhergesagte Wert. Wenn $ i = 1, 2, 3, 4 (= n) $, $ A_i $ und $ P_i $ die tatsächlichen Daten / vorhergesagten Werte von $ S, I, R, F $ sind.

Da es schwierig ist, ein Bild nur mit dem Wert zu erhalten, habe ich ein Diagramm erstellt. Zunächst etwa die 0. Phase, die den größten RMSLE-Wert aufweist. Das obere Diagramm zeigt die Differenz zwischen den tatsächlichen Daten und dem vorhergesagten Wert, und das zweite, dritte und vierte Diagramm zeigen sowohl die tatsächlichen Daten als auch den vorhergesagten Wert für jede Variable.

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

Es liegt ein Fehler vor. Es scheint besser, die 0. Phase mit "Scenario.separate ()" usw. aufzuteilen (ein separater Artikel wird erstellt).

Andererseits überlappen sich in der 6. Phase, die den niedrigsten RMSLE-Wert aufweist, die tatsächlichen Daten und der vorhergesagte Wert häufig.

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

8. Nachtrag

Dieses Mal erklärte ich, wie die Parameter jeder Phase geschätzt werden. Mehr als ein halbes Jahr ist seit Beginn des Ausbruchs vergangen, und ich denke, wir befinden uns in einem Stadium, in dem die Überprüfung der Parameter, die Vorhersage zukünftiger Parameter und die Szenarioanalyse wichtig sind. Es scheint, dass die Aufmerksamkeit von COVID-19 als Thema der Datenwissenschaft abnimmt, aber da die Daten akkumuliert werden, ist es möglich geworden, tiefere Analysen durchzuführen als in den frühen Tagen, als wir uns auf die Erstellung von Dashboards konzentrierten. Es war.

Vielen Dank auch für Ihre harte Arbeit!