StatModeling Memorandum

StatModeling Memorandum

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

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の略です)。有症状者の濃厚接触者など、数が限られており、低い値だと考えられます。

統計モデル

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

f:id:StatModeling:20201106155552p:plain

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

  •  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人ぐらいは存在しそうという推定結果になりました。

f:id:StatModeling:20201106155545p:plain

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

f:id:StatModeling:20201106155539p:plain

コード

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

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

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

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

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