2020/07/23

COVID-19 感染者数プロファイルの計算モデルと見方


COVID-19 感染者数プロファイルの計算モデルと見方

 
ロジスティック関数は,
   N = K/{1+exp[-r(t-T)]}                      1)
で表されます。

N はある時刻 t での個体数(感染者数の累計数)。
K は環境収容力とよばれ,ここでは最大となる感染個体数を表します。
r は内的自然増加率とよばれ,1個体がとりうる最大の増加率です。
T は変曲点(N K の半分となる,またはNの増加数が最大となる点)の時刻です。

ロジスティック関数は,関数の微分の形であるロジスティック方程式
   dN/d = r N (1-N/K)                        2)
の微分方程式の解です。ここでは,時刻を日の単位としているので,ロジスティック方程式の値は各時刻(日)における増加数を表しています。

ある集団について,第1波,第2波のように,あるいは,いくつかの小集団1,小集団2のように,複数の感染者数プロファイルから成っているとき,それらの合計として累計感染者数 NT を考えています。NT は,それぞれが式1)で表される感染者数プロファイル j の累計感染者数 Nj の和です。
   NT = Nj                            3)

解析では,報告されている日々の累計感染者数 No について, 累計感染者数 NT の計算値 Nc を最適化しています。最適化は,次式の Rw を最小にする非線形の最小二乗法によっています。
   Rw = w(No-Nc)2                       4)
ここでの w は個々(日々)の累計感染者数の確からしさを反映した重みです。この最小二乗法で, 各々の感染者数プロファイル j についての環境収容力 Kj,内的自然増加率 rj と変曲点の時刻 Tj の3つのパラメータが得られます。

ロジスティック関数モデルによる感染者数プロファイルを図に例として示します。計算は,東京都の第2波を想定し,環境収容力 K を 17,600,内的自然増加率 r を 0.08,変曲点の日付を8月1日としています。このブログでは確定日別の感染者数を用いており,実際に感染が起きたのは7日間ほど前と考えられ,確定日が8月1日ならば感染日は7月25日頃と思われます。

日本全体のモデルでは,内的自然増加率は東京都の場合とほぼ同じなので,プロファイルの概形もほぼ同じです。ただし,変曲点の日付は数日の遅れ,環境収容力は数倍となっています。


ロジスティック関数による計算モデルと見方 [図をクリックすると拡大

グラフの見方


感染確定日データの日別の感染者数の累計が,"累計obs"です。ただし,最新の値で割って,最大値が1となるようにした"累計obs'"をグラフにプロットしています。

累計obsに合致するようにロジスティック関数を最適化し,最適化した関数による計算値が"累計calc"です。この値を最新の累計obsで割った"累計obs'"と"累計calc'"をプロットしています。東京都の場合は第1波と第2波の合計で算出しており,最新の"累計obs'"は1です。

"日別obs"は日別の感染者数です。最適化した関数から計算される日別の感染者数が"日別calc"です。

最適化した関数から計算される内的自然増加率 r から
   RLe = τ r/{1+exp[r(t-T)]}                     5)
により計算される実効再生産数RLeが,"τ×増加率"です。ここでの τ (tau) は,感染者が感染させてしまう平均日数で,値は7を採用しています。初期の頃の"τ×増加率"(τr)に1を加えた数が基本再生産R0に対応すると考えられ,東京都の第1波では2,第2波では1.55程度です。よく使用される実効再生産数 R0e とここでの RLe との関係は R0e ≈ 2RLe-τr+1 です。

日別の感染者数から, τ ×(日増加数)÷(累計数)により見積もることができる"τ×平均"について,素のデータは曜日ごとのばらつきが大きいため,7日間の移動平均をとった値を"τ×平均"としています。第1波について"τ×平均1",第2波について"τ×平均2"としています。最新の3日間では7日間移動平均が適用できませんが,動向を把握するために,最新日は実際の値そのもの,前日では3日間の,前々日では5日間の移動平均を採用しています。そのため,最新日と前日の値の変動の幅は大きくなっています。

なお,東京都の確定日別データの最新日(発表日の前日)の数値は集計途中の値が発表され,集計が進んだときの値の約80%と,小さく見積もられています。最適化では,最新日の発表日ベースの値が後日の集計後の値に近いことから,この発表日ベースの値を確定日別の最新値の仮の値としています(前々日以後は確定日別の数値に更新します)。

"τ×平均"は関数モデルには依存しません。関数モデルが妥当ならば,"τ×平均""τ×増加率"に次第に合致するはずです。"τ×平均1"は第1波の"τ×増加率"によく沿っていて,"τ×平均2"は変化しながらも第2波の"τ×増加率"に追随しています。

"累計calc'""日別calc""τ×増加率"は日付を指定すれば計算できるので,数日後の値もプロットしています。

日別感染者数がピークに達するとき,"日別calc""τ×増加率"は変曲点に来ます。変曲点に来ると"τ×増加率"が初めの頃の値の1/2となります。"τ×増加率""τ×平均"が次第に小さくなって,半分となる時期が日別感染者数のピークです。このときの累計感染者数を2倍すると,最大値になります。

"日別calc"はピークを挟んでグラフでは左右対称となります(偶関数です)。ピークの前と後では日別感染者数,および,その累計値(こちらは奇関数)はほとんど同じ値になります。"τ×増加率"(奇関数)は累計値を日付を反転した形になっています。

0 件のコメント:

コメントを投稿