K値は実効再生産数の見積もりに役立ちます

K値は実効再生産数の見積もりに役立ちます

The K value yields reasonable estimates for effective reproduction numbers

(As of June 29, 2020)

先に，K値（末尾の引用を参照）は"日感染者数の7日間移動平均"×(7/累計感染者数），あるいは，"7日間の日感染者数の和"/(累計感染者数）に等しいことを記しました。

COVID-19 K値は日感染者数の7日間移動平均に対応する

https://ysatow.blogspot.com/2020/06/covid-19-k7.html

ここでは，私のブログで取り上げてきたロジステック関数との関連から，K値についてさらに考察してみます。考察をまとめると，以下のようになります。

1. 感染者数へのロジステック関数の最適化により，実効再生産数を見積もることができます。

2. K値は，感染のピークの前後のグラフでは，直線的になります。

3. K値は，感染の実効再生産数の見積もりに役立ちます。

4. K値は，最大となる感染個体数に依存しないというメリット（広く適用可能）があります。

5. 増加率(日感染者数/累計数)の7日間移動平均も感染拡大・収束をよく反映します。
6. ロジステック関数の最適化で得られる増加率はさらに良好な指標といえます。

SIモデルとロジスティック関数

ロジスティック関数は，コンパートメントモデルの基本となるSIモデル（SIRモデル，SIRDモデル等で，RやDへの推移を0とした場合に相当）と同じものです(Sは感染へのsusceptible，I はinfectious)。一般化ロジスティック関数であるGompertz曲線で，関数形を論ずることもあります。以下では，これまでのロジステック関数の記述（https://ysatow.blogspot.com/2020/05/covid-19_13.html）に準じて述べます。

ロジスティック関数は，

　　　N = P/{1+exp[-r(t-T)]}   　　　　　　　　　　　　　　　　　　　　1)

で表されます。

N はある時刻 t での個体数（感染者数の累計数）。

P は環境収容力とよばれ，ここでは最大となる感染個体数を表します。

r は内的自然増加率とよばれ，１個体がとりうる最大の増加率です（SIモデルのβに相当）。

T は変曲点（N が P の半分となる，またはNの増加数が最大となる点）の時刻です。

ロジスティック関数は，関数の微分の形であるロジスティック方程式

　　　dN/dｔ = r N (1-N/P) 　　　　　　　　　　　　　　　   　　　　　　2)

の微分方程式の解です。ここでは，時刻を日の単位としているので，ロジスティック方程式の値は各時刻（日）における増加数を表しています。

ロジスティック関数の最大値は漸近的に環境収容力P（私の以前の記述では記号Kを用いています）に等しくなります。また，増加数がピークとなった変曲点の時刻はT です。

感染者数の時間プロファイルは複数の感染態様の集計であることから，プロファイルを複数のロジスティック関数の和として表すと，実際の現象をより的確に表現することができます。

日毎の増加率

1つのロジスティック関数で表される感染のプロファイルについて，ある日i の日感染者数D_i は N_i - N_i-1で，(dN/dｔ)_i です。日毎の増加率 m として D_i /N_i を定めると，

　　　m = D_i /N_i = (1/N_i）(dN/dｔ)_i

です。これに2)式を代入すると，次式のようになります。

　　　m = r (1-N_i/P)　　　　　　　　　　　　　　　　　　　　　　　　　　　　3)

さらに1)式も代入すると，

　　　m = r exp[-r(t-T)]/{1+exp[-r(t-T)]} = r /{1+exp[r(t-T)]}　　　　4)

が得られます。

m の値は，3)式を見て解るように，累計数N に関しては直線です。4)式より，t << T では m は r，ｔ = T では r/2，t >> T では 0 となり，増加率の時間とともに減少する性質を有します。ｔ ≈ T の感染の累計数の変曲点の近傍では，m は時刻 t に関してほぼ直線的になります。4)式のように表してみると，m には，環境収容力 P の項が含まれていないので，感染の集団の大きさに依存しないという優れた特性があります。

ここで，感染者が他者に感染させる期間を τ とすると，τm は実効再生産数の見積もり（厳密な値ではなく，目安値）を与えることになります。新型コロナウイルスでは，τ は7-9日程度と考えられ，ここでは τ = 7として，以下を考察します。なお，ある日の前後7日間の実際のm値（D_i/N_i）の移動平均に7を乗じた値は同様に実効再生産数の見積もりを与えるものと考えられ，これをここでは m7 とよびます。

なお，実効再生産数は感染した日付を基にしており，ここでの感染者数は診断確定日を基にしており，1週間程度の時間の遅れがあります。また，感染したときの感染者数プロファイルに診断確定日のプロファイルは正確には一致しません。ご留意ください。

K値は 1-N_i-7 /N_i で定義され，特定の関数形やモデルを前提にしてはいません。ここでは，1)式を分子と分母に代入すると，やはり感染集団の大きさに依存しない式の形となります。また，(日感染者数の7日間移動平均)/(ある日の累計感染者数）なので，m7値へ近似した値と考えることができます。このことについては，次のブログ"K値を厳密に取り扱う"をご覧ください。
　

図1のグラフは，東京都の確定日別による陽性者数について，累計値（グラフでは累計obs），K値，τm値（増加率×7），m7値（平均×7）をプロットしました。同時にプロットした累計SIcalcは，私が累計値に対して最適化した２つのロジステック関数の和の値です。ここでの τm値はこの最適化した関数から4)式を用いて算出した値2つの和です。なお，累計の基準日は，感染の最初の確定日の1月22日としています。

Fig.1. K値，τm値（増加率×7），m7値（平均×7）の比較

τm値（増加率×7）は感染の主ピークの前方では，実効再生産数1に近い値から滑らかに減少し，4月7日には0.75，4月13日には0.5，4月28日には0.1まで低下しています。これら値は東京都での経過をよく反映していると考えます。

m7値（平均×7）は，3月25日からは τm値によく合致した曲線となっており，4月26日頃まではほぼ直線的です。大まかな挙動はK値のグラフに対応していますが，4月13日より前は，大きな値でかつ4日ほど挙動が先行しています（これは平均の日々の差3.5日による）。4月13日以降はK値と同様にほぼτm値に沿っています。

K値とm7値のグラフは，最適化した1変数のτm値とは異なり，細かな変化を反映しています。この変化は，累計obsと累計SIcalcの偏差に現れていますが，7倍してはいないので小さくて見難くなっていますが，図3に示した日別obsと日別calcの偏差をみると，日々の変化がより判別しやすくなります。

累計の基準日を変えてみる

図2のグラフに，累計数の基準日（以降の累計数からこの日7日前の累計数を差引く）を3月30日，そして5月23日として，K値をプロットしました。τm値については，2つのプロファイルに対応させて，基準日を最初の感染日と5月23日としています。最初の基準日の後では，図1に比べてK値はより直線的に値は減少し，次第に0に近づいています。しかし，5月23日以後の直後はK値は大きな値で，τm値との偏差も大きく，時間の経過とともにτm値に近づいています。図1の3月のm7値はK値よりも大きな値を示していて，5月23日以降ではいっそうこの傾向は大きくなっています。感染の初期のK値とm7値の解釈には注意が必要です。

Fig.2. K値とτm値（増加率×7）の比較(基準日の選択)

K値は，7×(日感染者数の7日間移動平均)/(累計感染者数）なので，τm値と7日間移動平均のm7値とは厳密には対応しませんが，τm値すなわち実効再生産数の変化を反映した指標といえます。したがって，τm値との対応が良くなる感染の収束の時期に，K値を収束の指標とすることは理にかなっています。

この考察も，たとえばK値の適用も，過去のデータに従えば，現状と将来はこのようになると述べています。将来の結果が解析と異なってきた，結果が合致しない，などに気が付くことができれば，大きな意義があります。日々の報告に接して，有意な変異に気が付くためには指標は重要です。

東京都の最近の状況

図3に，東京都の確定日別による陽性者数について，日感染者数(日別obs)，２つのロジスティック関数の和で算出したdN/dｔ を日別SI calc として，図2のτm値（増加率×7）とともに示しました。 

最初の大きな感染者数の分布を見てみましょう。その日感染者数，あるいは累計感染者数のプロファイルを，実効再生産数を反映するτm値の振舞から完全に説明はできるのですが，これら相互の関係をひと目で把握することはなかなか難しいものです。これら数量は，微分値，積分値，増幅値（密度効果を担う）となっているからです。それでも，たとえば，実効再生産数の目安値0.05となったときの，日感染者数，累計感染者数はいかなる値かは参照できるので，これらの値（場合によってはこれら値から計算できる週平均m7値やK値など）に着目してもよさそうです。

5月23日以降の感染者数の増加は，4月のピークの残り火が消えきらずに勢いを強めているようです。この増加は6月27日時点ではまだピークを迎えていないようですが，実効再生産数の目安値0.45から徐々に減少しており，関数形は奇関数なので，まだまだ現状の日感染者数が大きい状態は続きます。このまま経過するならば，実効再生産数の目安値は現時点の0.3から次第に減少するので，大幅な感染者数の増加は無さそうです。

一方，5月23日以降の感染者の増加のロジスティック関数の最適化では，環境収容力Pは3,000に近い値になっています。5月23日以降の累計感染者数は1,000名を既に超えており（6月29日現在），"このまま感染が推移する"と緩やかに経過するとしても，2,000名という大きな数の累計感染者数がこの先に発生するのではないかと危惧されます。

Fig.3. 東京都の日感染者数: データ(obs)と計算値(calc)および τm値（増加率×7）

Data source: 東京都の公開データ "確定日別による陽性者数の推移"
https://stopcovid19.metro.tokyo.lg.jp/?tab=reference

なお，τm値はこの先2か月程度までは計算してグラフにプロットできますが，変曲点の前方のデータだけから算出した値は誤差が大きいので，将来の予測に用いるときには注意が必要です。 予測性については"東京都の感染収束傾向は4月下旬には推測できた?"をご覧ください。

検討した感染の時間プロファイルは，大きな１つの集団にまとめたものです。この集団の中には多数の小さな集団（クラスター）が存在しています。東京都や日本全体などを対象とする大局的な検討は全体の傾向の把握に役立つますが，一括りにした結果では埋もれてしまって顕在化しない集団が多くあり得ます。例えば，先に示した福岡県，北海道，神奈川県などの感染のプロファイルは，日本全体のプロファイルにはほとんど顕れてきません。より小さな地域や集団についての検討が必要で，このブログでの取り組みは皆さんの理解に役立つものと考えています。

ここで示したロジスティック関数で算出したτm値（増加率×7）などは，感染のプロファイルの広い領域に渡って，感染の進行の理解に役立つ指標といえます。

大阪大学の中野貴志教授らによる論文の preprint: Novel indicator of change in COVID-19 spread status: by Takashi Nakano, Yoichi Ikeda (doi: https://doi.org/10.1101/2020.04.25.20080200)

この論文については，次のwebサイト[中野貴志教授（核物理研究センター）による論文等（K値について）］に詳しく記載されています。https://www.osaka-u.ac.jp/ja/news/info/corona/corona_info/from_members/rcnp_nakano