StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

西浦先生らによる実効再生産数の統計モデルを解説&拡張する試み

先日の西浦先生のニコ生の発表を聞いていない人はぜひ聞いてください。

モデルとデータを以下のリポジトリでオープンにしていただいたので、モデルについて僕が分かる範囲内で少し解説を加えたいと思います。

github.com

実効再生産数を推定するコードが2種類ありまして、最尤推定(Maximum Likelihood Estimation, MLE)を使ったMLE版(Sungmok Jungさん作成)と 、ベイズ推定版(Andrei Akhmetzhanovさん作成)があります。どちらもコンセプトはほぼ同じで、実装が若干異なります。この記事では、ベイズ推定版(以降、元コードと呼びます)の流れを簡単に説明し、その後でその拡張を試みます。

ベイズ推定版の流れ

大きく分けて「データの集計」「back projection」「実効再生産数の推定」の3つの部分からなります。

データの集計

まずは日付ごとの人数のデータを作ります。

  • 発病日(onset)があるデータ(12526人)については発病日を使います。日付ごとに人数の合計を算出します。
  • その他のデータ(2872人)については、診断日または報告日(report)を報告日と呼んで、日付ごとに人数の合計を算出します。

なお、発病日と報告日を両方持っている人は(悲しいことに)存在しません。

back projection

再生産方程式は、「ある時刻において感染者は何倍に増えたか?」を記述する式なので、「感染者の人数」と「感染から感染までの時間(生成時間, generation time)」をベースに考えます。そこで、上記の「発病日の人数」や「報告日の人数」を「感染日の人数」に変換します。この変換をback projectionと呼んでいます。

具体的には以下の式に従うと考えて、ある日 tにおける感染者数  infect_{t} を推定します。

ここで、 onset_{t} はある日 tにおける発病者数です。 inc_{t} は潜伏期間(incubation period, 感染から発病までの期間)の分布です。発病者数と潜伏期間の分布はデータとして与えます。

例えば、 t = 10 の発病者数  onset_{10} の平均は、1日目における感染者数  infect_{1} にその感染者たちが9日後に発病する確率をかけたもの、2日目における感染者数  infect_{2} にその感染者たちが8日後に発病する確率をかけたもの、…、9日目における感染者数  infect_{9} にその感染者たちが1日後に発病する確率をかけたもの、の総和で定まるだろうということを表しています。この(合計が10になるように)掛け算して総和をとる部分は畳み込み(convolution)と呼ばれ、本記事でたびたび登場します。

元コードではback projectionを実効再生産数の推定の前に行います。そこで使用している{surveillance}パッケージとその推定方法については、Med_KUさんの以下のブログ記事が詳しいのでぜひそちらを参照してください。

mikuhatsune.hatenadiary.com

英語になりますが、こちらにも詳しい説明があります。

元コードでは海外帰国と国内の感染者に分けて、それぞれに対してback projectionを行っています。まず2872人の報告日のデータについて、発病日→報告日の分布を与えて、発病日を計算するback projectionを行っています。すると、12526+2872人分の発病日のデータとなります。さらに発病日のデータについて、潜伏期間の分布を与えて、感染日を計算するback projectionを行っています。

back projectionではスムージングもあわせて実行されます。スムージングは「本来の感染者数は日付についてそこまでギザギザになっておらず、ある程度なめらかでしょう」という仮定を反映しています。スムージングの強さはハイパーパラメータとして与えます。

また、このback projectionの結果、求まる感染者数は小数の値になります。さらに変換前後で合計人数を変えないように、感染者数の合計=変換前の人数の合計 が成り立つように正規化しています。

実効再生産数の推定

元コードではback projectionを最尤推定で求めた感染者数を使って、実効再生産数Rtを推定しています。

まずは簡単のため、「感染しているけど、発病日あるいは報告日のデータとしてまだ上がってきていない」という「報告遅れ」を考慮しないモデルを考えます。具体的には、以下の再生産のモデル式に従うと考えて、実効再生産数を推定します。

ここで、 domestic~infect_{t} はある日  t における国内の感染者数です。 cases~infect_{t} はある日  t における海外帰国の感染者数と国内の感染者数の合計です。 p~gt_{t} は生成時間の分布です。感染者数と生成時間の分布はデータとして与えます。

ポアソン分布の平均が、実効再生産数と「感染者数と生成時間の分布の畳み込み」の積になっている点に注目です。ある日  t の実効再生産数が1程度であることは、「その時点より過去の感染者数」と「生成時間の分布」から計算される「現在の感染者数」が同程度で、感染が広がっていない(狭まってもいない)ことを表します。

左辺が小数の値をとるので厳密にはポアソン分布ではないのですが、階乗をガンマ関数で表現することでポアソン分布と同じ形の尤度を考えることができます。*1

西浦先生が説明されていましたが、無症状者の割合を時刻によらず一定と仮定すれば、無症状者まで含めた感染者数に対して、上記の再生産のモデル式がそのまま成り立ちます。よって、推定した実効再生産数をもとにした議論はそのまま適用可能です。

Stanによる実装

元コードにおいてモデルはStanで実装されています。元コードでは以下の部分が少し難解です。

  • 畳み込みに使う分布をStanコード内で定義している。
  • 連続値を許すポアソン分布の代わりに、rateパラメータが1のガンマ分布を用いている。
  • 尤度がガンマ分布の二度漬けになっている。
  • 報告遅れを平均ベースで考慮している。

そこで、ここでは元コードの意図と推定結果をほとんど変えないままで、少しシンプルにしたStanコードを紹介します。

  • 2~10行目: すべての日(ただし初日は除く)における、xyの畳み込みを求める関数です。計算をラクにするため、あらかじめyを逆順に並べたものを引数として渡します。
  • 12~16行目: 連続値を許すポアソン分布の定義です。lamがゼロを含む場合にはもう少し定義を頑張る必要がありますが、ここではこれで大丈夫です。正規化定数はclosed formで表現できないのですが、1に近いのでこれで近似する意味です(参考URL1, 参考URL2)。
  • 23行目: p_gt_revは生成時間の分布の値(逆順)です。Stan内で計算するよりRで計算して渡す方がラクなのでそうしています。
  • 35行目: (1.2)式に相当する部分です。.*vector型どうしの要素ごとの積です。

推定時間はiter=10000で、私のX1 carbonで1chainあたりおよそ20秒でした。

back projectionして求めた感染者数は元コードと同じものを使い、上記のStanコードで推定した結果の図が以下になります。実線と灰帯はそれぞれ実効再生産数の中央値と95%信用区間(左軸)、棒グラフはback projectionで求めた感染者数(右軸)です。

以下の図はMLE版の推定結果です。インターネット上でよく議論されていた図です。この図と大体同じ結果が得られていると分かります(上の図では横軸を05/03まで伸ばしている点は異なりますが)。

モデルの拡張

西浦先生の話と質疑応答にもありましたが、以下のような拡張は考えうると思います。

  • 検査のバイアス・陽性率の変動は考慮していない → 今のところ難しい。他のデータが必要そう。
  • 生成時間の分布をあらかじめ与えている → 分布の算出に使用したデータが少ないので検討の余地あり

この記事ではひとまず上記の拡張はさておき、元コードの意図を極力保持しながら、仮定を緩和するような拡張方法を考えます。具体的には、「back projectionのモデル化」と「報告遅れのモデル化」を行いましたのでそれらを説明します。

back projectionのモデル化

元コードにおける上記の流れの懸念点として、back projectionで確定的に求めた感染者数を、実効再生産数の推定においてデータとして与えてしまっていることが挙げられます。これをやると、back projection自体の不確定性を考慮していないことに相当し、信頼区間を過小評価することにつながります。

そこで、この節ではback projectionも実効再生産数のモデルに組み込みます。上記で説明したback projectionの(1.1)式をモデルに組み込んだだけで、残りの部分は元コードに忠実です。Stanコードは以下になります。

  • 26行目: p_inc_revは潜伏期間の分布の値(逆順)です。Rコード参照。
  • 27行目: p_rep_revは感染から報告までの時間の分布の値(逆順)です。Rコード参照。
  • 32行目: 感染者数に対して1階差分のスムージングをかけています。s_smoothはそのばらつきの大きさでデータから推定されます。
  • 33~36, 40~43行目: back projectionの変換前後で合計人数が変わらないようにするために、要素の合計が変換前の合計人数と一致するようなパラメータを定義しています。simplexにスケーリングを掛けることで実現するテクニックを使っています。国内domesticと海外帰国imported、発病onsetと報告report、の合計4種類のback projectionがあることに注意。
  • 51~55行目: simplexの空間でスムージングしています。スムージングのばらつきs_smoothを共有しています。
  • 57~61行目: back projectionの(1.1)式に相当する部分です。

上記をコンパイルするときに、「~の左辺がパラメータの非線形変換を含んでいたらヤコビアンを調整してね」という旨のお知らせが表示されますが、ここでは変換は定数をかけたり足したりする線形変換なので、ヤコビアンの調整は不要と思います。

推定を実行するRコードも載せておきます。推定時間はiter=6000で、私のX1 carbonで1chainあたりおよそ12分でした。

  • 10~34行目: データの集計部分です。
  • 37~59行目: 生成時間、潜伏期間、感染から報告までの時間のそれぞれ分布を定義しています。本来は連続値をとる確率分布を、日付ごとの離散的な分布に直すため、累積分布関数から計算している点に注意してください。特に最後の、感染から報告までの時間の分布は、潜伏期間+発病から報告までの時間 という二つの確率変数の和の分布になるので、ここでも分布の畳み込みが出てきます。互いに独立な確率変数の和の累積分布関数については、例えばこの講義ノートを参照してください。

結果は以下の図です。棒グラフは推定された感染者数の中央値、エラーバーは95%信用区間です。先ほどと比べて、実効再生産数の信頼区間が全体的に少し広くなっていることが分かります。

報告遅れのモデル化

元コードでは、「back projectionで求めた感染者数」を「感染から報告までの時間の累積分布関数の値」で割ることで、まだどのデータにも上がってきていない感染者数の分を増やして実効再生産数の推定に用いていました。その方法では、人数の多寡にかかわらず常に確定的に水増しされることになり、少人数の場合の不確定性をやはり考慮できていません。結果として、まだどのデータにも上がってきていない人数が含まれるであろう、ここ2週間ぐらいの実効再生産数の信頼区間が過小評価されることにつながります。

そこで、この節では報告の遅れも実効再生産数のモデルに組み込みます。ただし、元コードと同じ方法ではなく、以下の式に従うと仮定しました。

ここで、 infect_{t} はback projectionで求めた、ある日  t における感染者数です。 infect~latent_{t} は報告が遅れている人まで含めた、 t における感染者数です。 prob~report_{t} t においてどれぐらいの人数がデータとして上がってきているかを表す確率です。

Stanコードは以下になります。

  • 28行目: 感染から報告までの時間の累積分布関数の値です。
  • 38~41, 53~55行目: 報告が遅れている人も含めた感染者数を定義しています。
  • 71~73行目: (1.3)式に相当する部分です。

Rコードの説明は省略します。推定時間はiter=6000で、私のX1 carbonで1chainあたりおよそ30分でした。

結果は以下の図です。棒グラフは報告が遅れている人まで含めた感染者数を表します。先ほどと比べて、ここ2週間ぐらいの信頼区間が少し広くなっていることが分かります。

ここ2週間の不確定性を含めた推定結果が重要だと思うので、個人的には最後のモデルを推したいです。

enjoy!!

*1:WinBUGSではデフォルトでポアソン分布がそのように拡張されていました。

COVID-19 日本国内の潜在的な陽性者数を推定する試み

日本国内の潜在的な陽性者数を推定することは有益ですが、簡単ではありません。PCR検査がランダムになっていないことが推定を難しくしています。有症状者が検査されやすいというselection biasがあるからです。この記事ではいくつか仮定を置いて潜在的な陽性者数を推定したいと思います。

仮定

全国民のうち潜在的に陽性になっている割合

この割合は年代によらず一定と仮定します。ここでは  p_{pos} と書きます(posはpositiveの略)。例えば0.0001なら日本人約1億2千万人中、おおよそ12000人が潜在的に陽性になっている計算です。

なお、国民の年代別人口の値はこのページ令和2年3月報 (令和元年10月確定値,令和2年3月概算値) (PDF:301KB) の「2019年10月1日現在(確定値)」の総人口 男女計の値を使用しました。

陽性者中の有症状者の割合

若年層で無症状が多いなど、年代で異なることが知られています。そこで、全員が検査されて診断されているダイアモンドプリンセスのデータをモデルに組み込みます。ここでは  q[a] と書きます。 aは年代の添字(1: 0~9歳、 2: 10~19歳、…、 10: 90~99歳)を表します。

有症状者が厚生省の報道発表資料上に観測される割合

この割合は年代によらず一定と仮定します。ここでは p_{obs.sym} と書きます(obsはobservation, symはsymptomaticの略です)。そこそこ高い値であると考えられます。*1

なお、厚生省の報道発表資料はこのページの一番下の方の新型コロナウイルス感染症の国内発生動向(2020年4月3日掲載分)の2ページ目のスライドの右のグラフの値を使用しました。4/2 18:00時点の値になります。

無症状者が厚生省の報道発表資料上に観測される割合

この割合も年代によらず一定と仮定します。ここでは p_{obs.asy} と書きます(obsはobservation, asyはasymptomaticの略です)。有症状者の濃厚接触者など、数が限られており、低い値だと考えられます。

統計モデル

使用した統計モデルは以下になります。

ここで頭文字が大文字の変数はデータが与えられている変数です。各変数の説明は以下です。

  •  N[a] は年代別の国民の人数
  •  n_{pos}[a] は年代別の陽性者数(推定したいもの)
  •  n_{sym}[a] は年代別の有症状者数
  •  n_{asy}[a] は年代別の無症状者数
  •  Y_{sym}[a] は厚生省の報道発表資料の年代別の有症状者数
  •  Y_{asy}[a] は厚生省の報道発表資料の年代別の無症状者数
  •  Y_{sym.DP}[a] は年代別のダイアモンドプリンセスの有症状者数
  •  N_{pos.DP}[a] は年代別のダイアモンドプリンセスの陽性者数

各事前分布の説明は以下です。

  • (1.8)式は p_{pos}が低い値であることへの事前分布
  • (1.9)式は p_{obs.sym}が高い値であることへの事前分布
  • (1.10)式は p_{obs.asy}が低い値であることへの事前分布
  • (1.11)式は q[a]が年代に対してある程度滑らかであることへの事前分布

(1.8)式は楽観的、すなわち推定される陽性者数を減らす方向に働いています。色々変えて推定してみると良いと思います。

推定結果と考察

推定総陽性者数の分布は以下になります。分布の中央値でも7000人ぐらいは存在しそうという推定結果になりました。

年代別の推定総陽性者数は以下になります。黒点は中央値、灰帯は95%ベイズ信頼区間です。厚生省の報道発表資料と比べると、無症状かつ陽性者がかなり存在していると予想されます。特に10代以下の潜在的な陽性者はもっと多いと推測されています。

コード

Stanは整数値を推定するのが少し手間であるため、今回はJAGSを用いました。統計モデルを表すファイルは以下です。

JAGSコードでは、二項分布(dbinom)の引数は確率を先に書く点、正規分布dnorm)の引数は標準偏差ではなく精度であることに注意です。あとは統計モデルのままなので説明を省きます。

推定を実行して図を描くRコードは以下です。ここでは悲観的に「症状確認中」を「有症状」に入れてますが、「無症状」に入れて実行してみるのも良いと思います。「無症状」としてカウントすると、推定総陽性者数の分布の中央値は6800人ぐらいでそこまで変わりませんでした。

その他の推定結果は以下です。

*1:ちなみにこの割合が年代で異なるモデルでも結果に大きな差はありませんでした。

『わけがわかる機械学習』中谷秀洋(著)の書評

僕が中谷さんと初めて会ったのはみどりぼんの読書会で、初めて話したのは岩波DSの打ち合わせだったと思います。今でもそんなに親しくはないと思います。しかし、中谷さんのブログは10年ぐらい前から読んでいました。自然言語処理を中心とする機械学習に関連する理論(の解釈)・論文レビュー・数値実験の記事が多く、他のブログでは見られない独特かつ理解の深い内容で、毎日勉強させてもらっていました。今でも何度も読むべきブログです。その中谷さんが機械学習についてまるごと一冊書いたものが本書になります。もともと買うつもりでしたが、献本いただいたので簡単にご紹介いたします。

わけがわかる機械学習 ── 現実の問題を解くために、しくみを理解する

わけがわかる機械学習 ── 現実の問題を解くために、しくみを理解する

目次は以下になります。

  • 0章: はじめに
  • 1章: 機械学習ことはじめ
  • 2章: 確率
  • 3章: 連続確率と正規分布
  • 4章: 線形回帰
  • 5章: ベイズ確率
  • 6章: ベイズ線形回帰
  • 7章: 分類問題
  • 8章: 最適化
  • 9章: モデル選択
  • 10章: おわりに
  • 付録A: 本書で用いる数学

特長

一読して最初の感想は前半(2章~3章, 5章)が「確率および分布の初中級者向けの完ぺきな入門」、後半(4章~7章)が「わけがわかるPRML」あるいは「声に出して読みたいPRML」です。

残りの部分は、深層学習を含む機械学習全般(ところどころ深層学習にも言及しています)で必須の内容(モデルとは、交差検証、など)を扱っています。

理論の本にありがちな「天下りで定義や式が与えられて、わけわからんけど式変形していく」ではなく、「経験から確率は〇〇という条件を満たしてほしいよね。仮に〇〇として議論を進めるとどうだろう?」という感じに、ボトムアップで話が展開されていくので納得しながら読みやすいです。

5章は4章を拡張するために導入され、8章は7章の問題を解くために導入されます。章のタイトルだけ見ると一見関係性が不明ですが、流れるように議論をすすめるためにこの配置になっています。

数式に関しては0.1節に

本書は「数式がなくてもわかる本」ではありません。

と書いてあるように大学学部レベルの数式が出てきます。しかし、遠慮なく式が出てくるのは4章以降ですし、式変形の内容や考え方が非常に丁寧に書いてあるので、じっくり読めば理解できて力になることは間違いありません。

以降では部分的に補足します。

2~3章

0, 1, 9, 10章は合計して25ページぐらいなのに、2章と3章は合計して65ページもあって力の入れようが分かります。一見普通の章タイトルですが、内容はめちゃくちゃいいです。特に例題をまじえて「同時分布・周辺分布・条件付き確率分布」を説明しているあたりがいいです。

僕のブログの以前の記事

色々読んでみましたが、現在決定版と言えるものは存在しないように思えました。個人的には、シグマと積分の復習、場合の数・数え上げの方法、確率、確率変数、確率密度、度数分布とヒストグラム、代表値・平均・分散、確率分布、同時分布、周辺分布、確率変数の変数変換、検定、散布図と箱ひげ図、回帰、相関あたりをRやPythonなどを使いながらシンプルに説明していく本があるといいと思うのですが、なかなかバランスのとれたいい本がありません。初歩の初歩しか説明してない、グラフが少ない、検定にページを割きすぎ、分厚い、ちょっと難しいなどの不満点があります。立ち読みして自分にあった本を選ぶのがいいと思います。

と書きましたが、「確率、確率変数、確率密度、代表値・平均・分散、確率分布、同時分布、周辺分布、確率変数の変数変換」に関して僕が想定する内容と比べて200%ぐらいの完成度です。これからは「そのあたりがわかんなかったら、この本の2,3章読んでおいてね」と言うことにします。

4~7章

4.1節 最小二乗法の最後の段落(p.88)で

ここまで、自然な流れに従って「正解」を出したように感じるかもしれません。しかしこの短い話の中で大きな仮定を4つも使っています。仮定1:関数の形が1次式。 仮定2:二乗和誤差が最小=良い答え。仮定3:誤差独立の仮定。仮定4…(略)

というように前提条件をきっちり説明したうえで、「どこが拡張しやすいか」「拡張したらどうなるか」を続けて論じていきます。一貫してこのような理屈の説明が出てくるのが本書の特徴と思います。

7章は中谷さんの専門の自然言語処理の例題で話がすすむためか、モデルの考え方や苦労などもちらほらあって、勉強になるのはもちろんとても楽しく読める章になっています。

4~7章は、PRMLの3~4章を読もうとしたけど挫折した人、読んだけどイマイチ分かってない人に特にすすめたいです。

Enjoy!