2020年5月3日日曜日

自分で計算できる実効再生産数

新型コロナウイルス感染症専門家会議が5月1日に出した提言・状況分析における実効再生産数のグラフについてあれこれいわれている。こういう場合は自分で確かめておよその感じをつかみたい。いやいやこんな緊急事態に素人が世情を惑わせるよけいな計算をしてはいけないという声も,かつてよけいなことをしまくったその分野の素人(別の分野の専門家)から聞こえてくるのである。ここでComputational Thinkingを標榜するのであれば,計算の自由(Freedom of Computation)を宣言しておきたい。もちろんそのためには最低限データとアルゴリズムを明らかにしておく必要がある(専門家会議はそうしていない)。

①使用するデータ:日本の新規感染数の日次統計(WHOに報告されたもの)WHOのSituation Reportsのデータ(確定感染者数 Confirmed)は次のようになった。ただし,2月22日から5月1日まで(WHOが公表している日付を基準としたもの)の70点のデータである。
{12,27,12,13,7,22,24,20,9,15,14,16,33,32,59,48,33,26,54,52,55,41,64,34,15,15,44,77,46,50,43,39,65,98,96,112,194,173,87,225,206,233,303,351,383,252,351,511,579,658,743,507,390,455,482,585,628,566,390,367,378,423,469,441,353,203,191,276,236,193}

②使用するモデル:単純なSIRモデルを仮定し,感染が人口の数%を越えて拡がっていないものとし,未感染者数(感受性保持者 Susceptible $S(t)$)が全人口($N$ 定数)とほぼ等しいとする。なお,感染者(Infected $I(t)$)と未感染者の接触による1人1日当りの感染率 $\beta$は対策効果を含めて時間の関数$\beta(t)$とした。さらに,感染期間(日)を$\alpha$として,実効再生産数を$R_t=\alpha \beta(t)$で定義する。このとき,感染者数$I(t)$は次の微分方程式を満足し,$R_t$は$I(t)$とその時間微分から求まる。
\begin{equation}
\begin{aligned}
\dfrac{d I(t) }{dt} = \beta(t) S(t) I(t) / N - I(t) / \alpha \approx \dfrac{ R_t - 1}{ \alpha} I(t) \\
\therefore R_t = 1 + \alpha  \dfrac{d I(t) }{dt} / I(t)
\end{aligned}
\end{equation}
③計算方法:誤差はあるけれど時間の単位を1日とする差分式に直して,エクセルで計算する。もとの新規感染数データのままではゆらぎが大きいので,5日移動平均を求めて$I(t)$とする。これから中心差分($I(t+1)-I(t-1))/2$)で1日当りの新規感染数の変化分を求める。さらに,このゆらぎを緩和するためにこの5日移動平均をもとめて$d I(t)$とする。これから②の式を用いて$R_t$を求めた。

④結果:図のとおりである。
図1 全国の新規感染数の推移(2/22-5/1)

図2 全国の実効再生産数の推移(2/24-4/29)

図2には社会的な事象を書き加えているが,WHOへの報告公表時点とは時間的なずれ(1週間〜10日)があることに注意する。→(注)12日ぐらいか(P. S. 参照)

⑤ 結論:1ヶ月前に1.5〜2.0近くまであった$R_t$が現在0.5〜1.0の範囲にまで落ちてきたことがみてとれ,これは専門家会議の結果とおおよそ近いものである。しかし,細かく見ればそのピークの位置や高さは異なっている。とくに,2月末から3月中旬までの結果は待ったく違う。そもそも出発点となる新規感染数のデータもWHOに報告されたものは,1週間の大きな周期構造を持っており,牧野さんのいうように(ちょっと意味が違うかもしれないけれど),専門家会議が確定感染者数(Confirmed)からどのように処理して新規感染数の推定値を導いているのかがはっきりしないのでもやもやが残る。

P. S. 高橋健太郎さん(@kentarotakahash)によれば,「潜伏期間+発症〜報告までの日数の平均を12日と考えてみたら、かなり納得できました。15日目あたりから一度、下降してRe=1を切りますが、2/24+15-12ですので、2月27日の休校要請の効果に見えます。その下降のピークは3月4日頃。一週間で緩んで、再び上昇が始まった。3月10日頃にRe=2を越えるピーク。」だそうだ。またさらに,「その後、Re=1.5以上の状態が続きますが、三連休の人手の反省があり、3月25日の都知事緊急記者会見を経て、3月27日頃から下降に転じる。3月31日にはRe=1を切り、専門家会議の4月1日には1を下回ったという分析と合致します。」
そして,「その後一度、上昇して、1を越えますが、4月5日頃から本格的な下降が始まる。これは緊急事態宣言が出るというムードの先取り。が、4月14日あたりから再び上昇の兆しというところでしょうか。」なるほどなのだった。


[2]新型コロナウイルス感染症対策専門家会議の見解等その2(牧野淳一郎,2020.05.12)
[3]感染症数理モデル(北海道大学医学統計学教室のSqquential SEIRモデル)
[4]Rt-COVID-19 Japan (都道府県別新型コロナウイルスの実効再生産数)
[5]山中伸弥による新型コロナウイルス情報発信
[6]A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics(Cori et al.)
[7]新型コロナ対策専門家会議が判断の拠り所にしている『実効再生産数・倍加時間』の算出方法に関する考察(@makirin1230  2020.05.06)

0 件のコメント: