感染症数理モデル-補足4

ロジスティックモデルと対数ロジスティックモデルは簡便なモデルなので、マルチタイプ型流行モデルを比較的簡単に構築できます。最も簡単なものは、次のように複数のロジスティック関数または対数ロジスティック関数を単純合計したモデルです。

ロジスティックモデル：

　(t_j=t-t_0j)
対数ロジスティックモデル：

　(t_j=t-t_0j＞0)
j番目のモデルにおいて　y_j：時点t_jまでの累積感染者数　　M_j：累積感染者数最大値
c_j：比例定数　　y₀：t_j=0(t=t_0j)における初期感染者数

このマルチタイプ型流行モデルを構築するには、例えば日本を北海道、本州、四国、九州、沖縄という5つの地域に分けて、各地域ごとのデータを用いて5つのモデルのパラメータを求め、最終的にそれを合計して日本全体のモデルにするという方法が考えられます。また同じ地域で流行を何度も繰り返した時には、薬物動態の区画モデルで用いられる皮むき法(peeling、stripping)とか残差法(method of residual)と呼ばれる方法で複数のモデルのパラメータを求め、それらを合計するという方法が考えられます。 (当館の「統計学入門・第14章第2節内服モデル」参照)

補足1で説明した日本のデータを用いて、皮むき法の手順を具体的に説明しましょう。説明を単純にするために、ロジスティック関数を用いたマルチタイプ型流行モデルについて説明します。

図補足4.1 ロジスティック関数を用いたマルチタイプ型流行モデル

最初に、実際のデータを用いて上図右側のような1日ごとの感染者数(この場合はPCR陽性者数/1万人)のグラフを描きます。そして1番目の流行と2番目の流行の分岐点、つまり2番目の流行が始まったと思われるポイントを探します。この場合は69日後(2020年3月24日)頃に、第2の流行が始まったと思われます。そこで0日(1月15日)〜68日(3月23日)までのデータに1番目のロジスティックモデルを当てはめ、微分回帰分析によってパラメータを求めます。その結果は次のようになりました。

1番目のロジスティックモデル：M₁=0.127026(実数1613.22) 　　c₁=0.112938
y₀₁=0.000143868(実数1.82712) 　　t₀₁=0　　重寄与率R²=0.999298

この結果から、1番目の流行は約60日後(3月15日)にピークになり、約120日後(5月14日)に累積PCR陽性者数約0.13人/1万人(実数にして約1,600人)に収束すると予想されます。その様子を描いたものが、上の2つのグラフで「1番目のロジスティックモデル」とラベルを付けた青い破線の曲線です。

次にこのロジスティック関数を利用して、69日目(3月24日)から最後のデータが観測された100日目(4月24日)までの累積PCR陽性者数予想値を求めます。そして69日目から100日目までの実際の累積PCR陽性者数から、その予想値を引きます。上の左側のグラフで言えば、69日以後の実際のデータを表す「×」の値から青い破線の値を引き、「×」を下にずらすわけです。これが皮をむく作業です。

69日以後のデータの皮をむいたら、そのデータに2番目のロジスティックモデルを当てはめ、微分回帰分析によってパラメータを求めます。その結果は次のようになりました。

2番目のロジスティックモデル：M₂=1.07655(実数13672.2) 　　c₂=0.167199
y₀₁=0.0307952(実数391.099) 　　t₀₂=69　　重寄与率R²=0.998611

この結果から、2番目の流行は約69+21=90日後(4月15日)にピークになり、約120+21=141日後(6月4日)に累積PCR陽性者数約1.1人/1万人(実数にして約14,000人)に収束すると予想されます。その様子を描いたものが、上の2つのグラフで「2番目のロジスティックモデル」とラベルを付けた青い破線の曲線です。

そして0日目から150日目について2つのロジスティック関数を合計すると、全体の流行は約90日後(4月15日)にピークになり、約141日後(6月4日)に累積PCR陽性者数約1.2人/1万人(実数にして約15,000人)に収束すると予想されます。その様子を描いたものが、上の2つのグラフで「マルチタイプ型流行モデル」とラベルを付けた赤い実線の曲線です。

ロジスティックモデルを用いたマルチタイプ型流行モデルのk値、実効再生産数R_t、K値、倍加時間T_dは次のようにして求められます。

⊿R_t=1 + k：1日あたりの再生産数
R_t=1 + Dk　　D：平均感染期間

上図左側のグラフのように、マルチタイプ型流行モデルのk値はピークが複数になります。そして右側のグラフのように、倍加時間はk値のピークから7日ほど遅れて極小値になります。

一方、対数ロジスティックモデルを用いると、マルチタイプ型流行モデルのk値、実効再生産数R_t、K値、倍加時間T_dは次のように少し複雑な式になります。補足1で説明したk値と倍加時間のグラフは、この式を用いて計算した結果をグラフにしたものです。

⊿R_t=1 + k：1日あたりの再生産数
R_t=1 + Dk　　D：平均感染期間

これらのグラフの解釈は補足1を見てください。そしてこれらのグラフのk値曲線と倍加時間曲線をロジスティックモデルを用いた時のものと比較すると、k値のピーク値と倍加時間の極小値がより極端になり、時期が少し前になっていることがわかります。これらのことから、やはり対数ロジスティックモデルの方が実態に合っていると思われます。

補足4.マルチタイプ型流行モデル