Introduction
COVID-19のデータ（PCR陽性者数など）のデータを簡単にダウンロードして解析できるPythonパッケージ CovsirPhyを作成しています。
紹介記事：
今回はParameter estimation（SIR-F modelなどのパラメータ推定）のご紹介です。
英語版のドキュメントはCovsirPhy: COVID-19 analysis with phase-dependent SIRs, Kaggle: COVID-19 data with SIR modelをご参照ください。
1. 実行環境
CovsirPhyは下記方法でインストールできます！Python 3.7以上, もしくはGoogle Colaboratoryをご使用ください。
- 安定版：
pip install covsirphy --upgrade - 開発版：
pip install "git+https://github.com/lisphilar/covid19-sir.git#egg=covsirphy"
import covsirphy as cs
cs.__version__
# '2.8.2'
| 実行環境 | |
|---|---|
| OS | Windows Subsystem for Linux / Ubuntu |
| Python | version 3.8.5 |
2. 準備
本記事の表及びグラフは2020/9/12時点のデータを使用して作成しました。実データをCOVID-19 Data Hub1よりダウンロードするコード2はこちら：
data_loader = cs.DataLoader("input")
jhu_data = data_loader.jhu()
population_data = data_loader.population()
また、下記コードにより実データの確認とS-R trend analysis（感染拡大状況のトレンドの分析法）3を実行してください。
# 解析用クラスのインスタンス生成
snl = cs.Scenario(jhu_data, population_data, country="Japan")
# 実データの確認
snl.records(filename=None)
# S-R trend analysis
snl.trend(filename=None)
# Phase設定の確認
snl.summary()
※v2.13.0以降DataLoader.japan()の呼び出しが不要になりました。（2020/12/24追記）
実データのグラフ：
S-R trend analysis:
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. 実行例
S-R trend analysisにより、"Phase"（パラメータが一定となる期間）に分割することができました。本記事では、各"Phase"のデータ（例えば0th phaseであれば2020/2/6 - 2020/4/21のデータ）を用いてパラメータの値を推定します。
推定のしくみについては別記事を作成する予定です。optunaパッケージを使ってパラメータの値を提案、scipy.integrate.solve_ivpによって数値解を計算、実データとの誤差の少ないパラメータセットを選び出しています。
結果の見方は後述しますが、下記2行で実行・結果一覧の取得ができます。今回はSIR-F model4を使用しました。
# Parameter estimation with SIR-F model
snl.estimate(cs.SIRF)
# パラメータ一覧の取得
snl.summary()
# 標準出力の例（CPU数などによって処理時間は異なる）
# 詳細後述：最新のPhase = 6th phaseにてtauを含めた推定を行った後、tauを固定して0-5thのパラメータ推定
<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. パラメータ推定値
横に長いので順番に見ていきましょう。まずはパラメータの推定値です。
# cs.SIRF.PARAMETERS: SIR-F modelのパラメータ名リスト
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 model4:
\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$は「時間」を次元として有しています。この「時間」は連立常微分方程式$f'(T)=x(T)$を離散化した方程式$f(T+\Delta T) - f(T) = x(T) \Delta T$の$\Delta T$に依存しています。1日おきにデータを取得している以上$\Delta T$は1日（=1440 min）未満となりますが、具体的な値はわかりません。国や地域によって異なる可能性もあり、異なる国の有次元パラメータ$\alpha_2$, $\beta$, $\gamma$の比較を行うことは困難です。
そこで$\Delta T$に相当する変数$\tau$を設定し、有次元のパラメータを無次元化しました。
\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*}
このとき4、
\begin{align*}
& 0 \leq (x, y, z, w, \theta, \kappa, \rho, \sigma) \leq 1 \\
\end{align*}
同一国内では$\tau$の値が一定となるように、最新のPhaseデータを用いて$\tau$と無次元パラメータを推定した後、他のPhaseのパラメータ推定時には同じ値を$\tau$として使用するようにしています。今回は表の通り480 minとなりました。計算を単純化するため1日=1440 minの約数を使用するようにプログラミングしています。
では、各無次元パラメータの推移をグラフ化してみましょう。
\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*}
rhoの推移
Susceptible（感受性保持者）がInfected（感染者）と接触したとき、感染する確率 $\rho$の推移：
snl.history(target="rho", filename=None)
緊急事態宣言（区域限定2020/4/7, 全国2020/4/16, 全国で解除2020/5/25）の効果が4月後半に現れ、7月初旬まで維持されていた、と解釈しています。その後はね上がり、少しずつ低下してきています。詳細については議論が必要ですが、3密回避などの対策の効果が直接的に反映されるパラメータとなっています。
sigmaの推移
InfectedからRecovered（回復者）に移行する確率 $\sigma$の推移：
snl.history(target="sigma", filename=None)
4月後半に大きく上昇した後、上昇/下降を繰り返しながら上昇傾向にあります。医療提供体制、新薬の開発/供給状況が反映されるパラメータです。
kappaの推移
感染者の死亡率$\kappa$の推移：
snl.history(target="kappa", filename=None)
大きく変動しているように見えますが、値の絶対値が小さいため、ある程度一定に抑えられていると考えてよいのではないでしょうか（他国との比較による検証が必要）。医療体制を整え、新薬の十分な供給により$\kappa$を限りなく0に近づけていくことが必要です。
thetaの推移
確定診断を受けた感染者のうち、確定診断の時点で亡くなっていた感染者の割合$\theta$の推移：
snl.history(target="theta", filename=None)
感染初期においては医療提供体制自体が極めて逼迫していたため解釈は難しいですが、検査が遅れて適切な治療を行えなかった場合などに値が上昇すると考えられます。
5. 有次元パラメータの推移
有次元のパラメータについても表に含まれています。解釈しやすくするため、もともと無次元の$\alpha_1$を除き、逆数にして単位を[day]にしています。
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 |
有次元パラメータは無次元パラメータと同じ経過を示すため省略しますが、グラフは下記コードで取得できます。
# betaの場合 -> グラフ省略
snl.history(target="1/beta [day]", filename=None)
6. 実効再生産数の推移
SIR-F model4の基本/実効再生産数Rtは次の通り定義しています。
\begin{align*}
R_t = \rho (1 - \theta) (\sigma + \kappa)^{-1} = \beta (1 - \alpha_1) (\gamma + \alpha_2)^{-1}
\end{align*}
推移：
# 一覧
cols = ["Start", "End", "ODE", "tau", "Rt"]
snl.summary(columns=cols)
# グラフ
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 > 1$が感染拡大の1つの目安となるため、水平線$Rt=1$を表示させました。
7. パラメータの正確性
パラメータの正確性を測るRMSLE score, パラメータを推定するためにoptunaパッケージが提案したパラメータセットの数、実行時間についても表に含まれています。
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 |
RMSLE score (Root Mean Squared Log Error)5の定義式は次の通りです。0に近いほど実データをよく反映していると言えます。省略しますが、推定方法の検証自体は理論データ（SIR-F modelの式とパラメータセットの例から理論的に作成される患者数データ）を使って行っております。
\begin{align*}
& \sqrt{\cfrac{1}{n}\sum_{i=1}^{n}(log_{10}(A_{i} + 1) - log_{10}(P_{i} + 1))^2}
\end{align*}
$A$は実データ, $P$は予測値を示しています。$i=1, 2, 3, 4(=n)$のとき、$A_i$及び$P_i$はそれぞれ$S, I, R, F$の実データ/予測値です。
値だけではイメージがつきにくいためグラフ化してみました。まずはRMSLE値が最も大きい0th phaseについて。1番上のグラフは実データと予測値の差、2, 3, 4番目は変数ごとに実データと予測値の両方を表示しています。
snl.estimate_accuracy(phase="0th", filename=None)
誤差がある程度生じています。Scenario.separate()などを使って0th phaseを分割したほうが良いようです（別記事を作成予定）。
一方でRMSLE scoreが最も小さい6th phaseでは実データと予測値がよく重なっています。
snl.estimate_accuracy(phase="6th", filename=None)
8. あとがき
今回は各Phaseのパラメータを推定する方法についてご説明しました。流行開始から半年以上が経過し、パラメータの検証、今後のパラメータの予測、シナリオ分析が重要となる段階に来ていると思います。データサイエンスのテーマとしてCOVID-19の注目度が下がってきているようですが、データが蓄積されてきた分、ダッシュボードの作成に注力していた初期の頃より深い分析ができるようになってきました。
今回もお疲れさまでした！
-
Guidotti, E., Ardia, D., (2020), “COVID-19 Data Hub”, Journal of Open Source Software 5(51):2376, doi: 10.21105/joss.02376. ↩
-
Please refer to What’s the Difference Between RMSE and RMSLE? ↩