直近一週間のデータで基本再生産数R₀を推定する方法

新型コロナウイルス（SARS-CoV-2）感染者がどこまで増加しうるのか、その大雑把な目安として使える数字に基本再生産数R₀と言うのがある。R₀<1であれば感染症はすぐに消えてしまうが、R₀>1であれば集団免疫を獲得するまで累積感染者数は増加していき、R₀に応じて累積感染者が定まる。

精緻にR₀を推定するのは専門的かつ労力が要る作業なのだが、教科書的なSIRモデルに準じて計算するのはそうでもない。推定モデルの導出までだったら、高校もしくは学部初年度の微分積分の範囲でいける。複雑なモデルは短期間のデータからは上手く推定できないし、今回は雑な方法で推定してみよう。

SIRモデルは、人口に対するt期の感受性保持者の比率（S_t）、感染者（I_t）の比率、免疫保持者（R_t）の比率を決定する微分方程式モデルで、以下の4式で表される。

(1)

(2)

(3)

(4)

βは（新規）感染率、γは回復率だ。

(2)式の両辺をI_tで割ると、もっとも基本的な微分方程式の問題になる。

短期での推定に用いるので、S_t=S₀と置いて、両辺を積分しよう。

(5)

C₀は積分定数。毎日の入院患者数などを使ってこの式を対数線形化して回帰分析をしてもよいが、退院者数のデータが見当たらない事も多い。また、初期のデータでは、最初に感染した人が早々に退院するなどして、感染者数ゼロが続く期間もあり得る。累積感染者数J_tから推定できるように整理しよう。

tで微分し、dS_t/dtに式(1)を代入したあとで、近似としてS_tをS₀に置き換える。

(6)

J_tは、0時点での累積数と、0時点からt時点までの累積数の増加に分けられることに注意すると、以下となる。

式(6)を代入し、右辺第2項の積分をする。

以下のように置いて整理しよう。

t期の累積感染者数J_tは以下のようになる。

J₀-C₁とexp((βS₀-γ)t)のtに対するオーダーから、J₀-C₁を無視して対数化して線形回帰をしてもよいのだが、一週間程度の短い期間で推定したいので、誤差項に正規分布を仮定して、最尤法を用いる^*1。累積感染者数J_tは都道府県などの報告、回復率γは平均9日間で半分が陰性になると言う中国の経験が使える。

平均して2週間ほどで回復することから以下を満たすγを推定してもよい。(1-γ)^tが0期に感染したがt期までに回復していない比率で、そのうちt期ではγだけ回復をするので、tをかけて積分すると回復日数の期待値が計算できる。

S₀は無症者・軽症者が多いことから、7割の感染者は露見していないとして1-累積感染者数/0.3/人口とした。もっとも、ほとんど1なので、現段階でこの暗数の仮定は推定結果に無影響。

こういうわけで、推定する未知のパラメーターはC₁とβと誤差項の標準偏差の3つだ。隔離措置など感染予防策の効果が入ってしまっているが、S₀を含めることで免疫保持者の影響をある程度はコントロールしているので、実効再生産数Rではないとする。J_tは、人口比に直していない累積感染者数をそのまま使っても、βの推定値に影響は無い。

SIRモデルを動かすにはγとβと初期値があれば十分だが、R₀が無いと話がしづらいので計算する。近似式のようなのがあるのだが、以下を用いる。なお、この式、積分の公式で代数的解-β/log(1 - γ)が出てくる。

これで潜伏期間を考慮し、先週、先々週の施策の結果を見ることができるようになった。どうも収束先が初期値に左右されるので、推定結果の当てあまりの確認が欠かせない分析方法になってしまったが。

ところで直近7日間の東京のデータからだと以下のように、近くの閑散とした神社に御参りに行きたくなるようなR₀が観察された。都民の9割ぐらいが感染しそうな勢いである。

*1小標本のためか推定結果が安定しないため、以下の式を対数化回帰して残差平方和が最小になるように、最適化アルゴリズムでRを定め、そのとき推定されたβを用いた。J₀-C₁をRと置いているわけだが、t=0を考えるとJ₀=R+C₁+εとなって、事後的に制約を満たすことになる。

なお、βS₀<γのときはC₁<0となるので、R-Jₜ=-C₁exp((βS₀-γ)t)を推定する。逓増から逓減に変わりそうなときは、βS₀とγの大小関係は事前には分からないので、C₁>0とC₁<0のそれぞれを仮定した推定を行い、残差平方和が小さい方を採用する。