… S区画の人は単位時間あたりβS(t)I(t)だけ減少する最終規模方程式：pi≒1-exp(-piR0) … 流行終息時におけるS区画の感染者の割合
：基本再生産数 … 1人の感染者が隔離・回復・死亡するまでの間に何人の人を感染させたかを表す値
SIRモデルにおけるS区画の微分方程式ではS区間は単位時間あたりにβS(t)I(t)だけ減少し、流行終息時は最終規模方程式によって規定される累積感染者の割合pi分が減少します。
… S区画の人は単位時間あたりβS(t)I(t)だけ減少する：基本再生産数 … 1人の感染者が隔離・回復・死亡するまでの間に何人の人を感染させたかを表す値
そこでモデルを単純にするためにβとγは流行期間中は一定と考え、最終的な感染者数(pi×S0)をあらためてS区画の初期値S0と考えます。 これは感染対策を何もせずに流行が進行し、最終的にS0全員が感染するモデルと解釈することができます。 そして累積感染者数をC(t)とし、これをC区画とすると、C区画の変化速度を表す微分方程式を次のように表すことができます。
S(t)=S0 - C(t)
… S区画の減少速度＝C区画の増加速度
この微分方程式を少し飛躍させ、C(t)の増加速度はI(t)の積分値であるC(t)自身に比例し、C(t)の最大値S0に対する増加可能非感染者の割合{S0-C(t)}/S0にも比例し、両者の比例定数は同じと考えると次のような微分方程式が成り立ちます。
：y0のロジット＝ロジット初期値
この微分方程式は、1883年にピエール＝フランソア・フェルフルスト(Pierre-François Verhulst)が人口予測のための関数として提唱した微分方程式と同じものであり、ロジスティック方程式と呼ばれています。 そしてこのロジスティック方程式の解がロジスティック関数です。 この関数の比例定数cはβと同じようなパラメータですが、βがI(t)とS(t)の比例定数であるのに対してcはC(t)と{S0-C(t)}/S0の比例定数であり、c/S0=c/Mがβに近いことになります。 (当館の「統計学入門・第9章第6節 ロジスティック曲線 (注1)」参照)
現実の感染で実際に観察できるのは、単位時間(通常は1日)あたりの感染者I(t)のうち発見されて隔離された感染者の数です。 それはγI(t)に近い数ですが、発見されずに回復したり、死亡したりする感染者もいるはずです。 そこで発見されて隔離される割合をαとすると、αI(t)≦γI(t)≦I(t)になります。 そこで上記のロジスティック方程式とロジスティック関数において、y=αC(t)、M=αS0、y0=αC(0)と置き直すと次のようになります。 これは実際の観察データに適用可能な簡便かつ実用的なモデルです。
：y0のロジット＝ロジット初期値このようにしてロジスティックモデルを導いた後で、インターネットを検索したら同じモデルを利用している論文を発見しました。 元々、感染症数理モデルは人口学的モデルを参考にしているので、以前からロジスティックモデルがよく使われるているようです。 (Christopher Y.Shen「Logistic growth modelling of COVID-19 proliferation in China and its international implications」)
ロジスティック関数のパラメータは関数に関して線形ではないので非線形最小2乗法を利用してパラメータを求めるのが普通であり、上記の論文でも非線形最小2乗法を用いています。 しかし非線形最小2乗法は適当な初期値を指定し、繰り返し計算によってパラメータを求めるややこしい手法です。 そのためデータの誤差が多かったり、データが最大値付近まで観測されておらず、途中までしかない時は正確な解を求めるのは困難です。 感染症対策を検討する時は流行の途中で将来の動向を予測したいことが多いので、これでは困ります。 (当館の「統計学入門・第14章第1節 コンパートメントモデル (注2)」参照)
そこで流行の途中までデータがあれば、まずまずの精度でパラメータを求められる微分回帰分析(Differential Regression Analysis)を利用してロジスティック関数のパラメータを推測することにしました。 微分回帰分析は国立長寿医療研究センターの岡田佑介先生が開発され、僕が数学的な理論付けをお手伝いした新しい統計手法です。 この手法は誤差が多かったり、データ数が少なかったりしても非線形関数の回帰分析を効率良く行うことができます。 微分回帰分析の詳しい説明は当館の「統計学入門・付録7 微分回帰分析」をご覧ください。
まず最初にクルーズ船「ダイヤモンド・プリンセス」のデータにロジスティックモデルを当てはめ、非線形最小2乗法と微分回帰分析の両方でパラメータを求めてみました。
左のグラフの青い曲線は非線形最小2乗法によって求めたロジスティック曲線であり、赤い曲線は微分回帰分析によって求めたロジスティック曲線です。 パラメータがよく似ているので、2本の曲線もよく似ていることが分かると思います。 右のグラフの赤い曲線は微分回帰分析によって求めた微分ロジスティック曲線であり、近似的な流行曲線になります。 この曲線は2月5日から10日後の2月15日にピークになり、その時のピーク値は64人です。
累積感染者C(t)の近似関数が求められれば、第2章と第3章で説明した計算式を用いて実効再生産数Rtや倍加時間Tdの近似値を理論的に求めることができます。
左のグラフの青い曲線はk値曲線であり、右のグラフの青い曲線は倍加時間曲線です。 k値は約0.34から始まり、微分ロジスティック曲線がピークになる2月15日に0になり、その後は負の値になっています。 (1+k)は1日あたりの再生産数を表すので、最初は1人の感染者が1日に1.34人を感染させていて、2、3日後から急に感染させる人数が減り、2月15日頃に1人だけになり、その後は1人より小さくなって流行が収束しました。
COVID-19の場合、平均感染期間は5日程度と言われています。 従って実効再生産数Rtで表すと、1月31日頃(2月5日の5日前)：Rt=1+0.34×5=2.7→2月10日頃(2月15日の5日前)：Rt=1→2月10日以後：Rt＜1と変化したことになります。 ただしこれは流行全体でβとγの値は変わらず、Rtの変化はS(t)の変化だけに依存していると仮定した時のものであり、実際のRtの変化を緩やかにしている可能性が高いと思います。
事実、1日あたりの感染者数を表す棒グラフを見ると、実際の流行曲線は左右非対称であり、流行が拡大するスピードと流行が収束するスピードは異なることがわかります。 そのため第5章で説明したように最初の感染者が発見される前はβが大きくてγが小さく、Rt=15くらいで、検疫・隔離処置が実施された後はβが小さくなってγが大きくなり、Rtは急激に小さくなり、流行終息時には流行全体の平均Rtが1.1程度になったと考えられています。
一方、倍加時間は約2日から始まり、微分ロジスティック曲線がピークになる少し前にようやく3日になり、2週間後の2月19日以後は急激に大きくなっています。 上の2つのグラフから、K値や倍加時間よりもk値やRtの方が反応が早く、流行の動向をより早く予測できることがわかります。 これは、k値の式とK値の式はほとんど同じものの、k値の方はyが2倍されているのでK値よりも速く変化することからもわかります。
補足3へ