玄関雑学の部屋雑学コーナー感染症数理モデル

補足3.対数ロジスティックモデル

実際の流行曲線が左右非対称であるのに合わせて、ロジスティックモデルを左右非対称にする方法を考えてみました。 ロジスティックモデルを左右非対称にする方法は色々ありますが、最も単純な方法は対数ロジスティックモデルにすることです。

ロジスティック曲線は正規分布の近似関数として利用されます。 そして正規分布において、説明変数xを対数変換してln(x)にすると左右非対称の対数正規分布になります。 それと同様にロジスティック関数の説明変数tを対数変換してln(t)にすると左右非対称な対数ロジスティック関数になり、対数正規分布の近似関数として利用できます。 (当館の「統計学入門・付録1 各種の確率分布」と「統計学入門・第10章第2節 各種のシグモイド曲線」参照)

対数ロジスティック関数: (t>0)
微分対数ロジスティック関数:
:y0のロジット=ロジット初期値
y:時点tまでの累積感染者数  M:累積感染者数最大値
c:比例定数  y0=C(0):ln(t)=0→t=1における初期感染者数

このモデルでは微分対数ロジスティック関数が微分方程式に相当します。 これをロジスティックモデルと同様に解釈すると、「累積感染者数yの増加速度vはy自身と、yの最大値Mに対する増加可能非感染者の割合{M-y}/Mに比例し、さらに経過時間tに反比例し、それらの比例定数cは同じ」という微分方程式になります。

vがtに反比例するということは増加速度が経過時間とともに遅くなる、つまり時間が経つほど感染速度が遅くなるということを表します。 でも普通はこんな現象は有り得ないので、時間が経つほど感染率が小さくなると解釈した方が合理的です。 本来のSIRモデルは感染時点からの経過時間である感染齢を考慮していて、感染率βと隔離率γが感染齢に依存して変化するというモデルです。 この対数ロジスティックモデルは、それをかなりラフに近似したモデルと考えられなくはありません。

前章と同様にクルーズ船「ダイヤモンド・プリンセス」のデータに対数ロジスティックモデルを当てはめ、微分回帰分析によってパラメータを求めると次のようになります。

DPのロジスティック曲線 DPの微分ロジスティック曲線
図 補足3.1 DPの対数ロジスティック曲線と微分対数ロジスティック曲線
ロジスティック関数:M=709.7  c=0.361018  y0=16.5233  重寄与率R2=0.998568
対数ロジスティック関数:M=728.024  c=3.51883  y0=0.210684  重寄与率R2=0.99639

上記のように対数ロジスティックモデルにするとMが少し大きくなり、流行のピークが10日後から8日後になってピークの値が64から69になります。 そして重寄与率を見ると、このデータの場合はロジスティック関数の方がわずかに当てはまりが良いようです。

対数ロジスティックモデルではk値、実効再生産数Rt、K値、倍加時間Tdは次のようになります。


2階微分対数ロジスティック関数:
⊿Rt=1 + k:1日あたりの再生産数
Rt=1 + Dk  D:平均感染期間

DPのk値曲線 DPの倍加時間曲線
図 補足3.2 DPのk値曲線と倍加時間曲線

上記の式を用いて対数ロジスティックモデルのk値曲線と倍加時間曲線を描いてみると、倍加時間曲線はロジスティックモデルとあまり変わりませんが、k値曲線はかなり異なります。 k値曲線を見ると、流行の最初は1日当たりの再生産数(1+k)=2.8、実効再生産数Rt=1+5×1.8=10から始まり、すぐに急激に小さくなって8日後の2月13日頃に(1+k)=Rt=1になり、その後は1より小さくなっています。 そしてロジスティックモデルよりもk値が大きいので、流行はゆっくりと収束しています。 実際のRtの変化はロジスティックモデルよりも対数ロジスティックモデルに近いので、実効再生産数に関しては対数ロジスティックモデルの方が近似が良いと考えられます。

「ダイヤモンド・プリンセス」の場合は流行が速やかに収束したので、ロジスティックモデルと対数ロジスティックモデルの違いは大きくありません。 しかし流行がゆっくりと収束しつつあるイタリアのデータに両方のモデルを当てはめると、両者の違いがもっと大きくなります。

ITの累積感染者数 ITの感染者数
図 補足3.3 イタリアの1万人当たりの累積感染者数と感染者数
ロジスティック関数:M=38.6046(実数=229427)  c=0.0885754   y0=0.127725(実数=759.064)   重寄与率R2=0.9983
対数ロジスティック関数:M=40.7785(実数=242346)  c=5.15631   y0=1.8×10-8(実数=0.0001)   重寄与率R2=0.999663

上の2つのグラフを見るとイタリアの場合は対数ロジスティックモデルの方が当てはまりが良く、流行が収束するにはもう少し時間がかかると予想されます。

イタリアのk値曲線 イタリアの倍加時間曲線
図 補足3.2 イタリアのk値曲線と倍加時間曲線

一方、k値曲線を見ると、ロジスティックモデルでは約0.09(1日あたりの再生産数=1.09、実効再生産数Rt=1.45)から始まり、ゆっくりと小さくなって、流行ピーク時(約64日後の4月4日)に0になり、その後もゆっくりと小さくなっています。 それに対して対数ロジスティックモデルでは約0.24(1日あたりの再生産数=1.24、実効再生産数Rt=2.2)から始まり、急に小さくなって流行ピーク時(約61日後の4月1日)に0になり、その後はゆっくりと小さくなり、約75日(4月15日)以後は約-0.05〜-0.04くらいの値を持続しています。

また倍加時間曲線はロジスティックモデルでは約8日から始まり、ゆっくりと大きくなっていますが、対数ロジスティックモデルでは最初の約1ヶ月間は3日以下になっています。 累積感染者数を見ると、最初の1ヶ月間で累積感染者が約0.3人/1万人(実数にして約2000人)になり、その後1ヶ月後には約17人/1万人(実数にして約10万人)に急増しています。 そのため最初の2ヶ月間は医療崩壊を起こすほどの爆発的患者急増状態だったと思われます。 このことを考えると、イタリアの場合はロジスティックモデルよりも対数ロジスティックモデルの方が実態に合っていると思われます。