読んだ。
Ascertainment rate of novel coronavirus disease (COVID-19) in Japan | medRxiv
ascertainment rate という、感染者数（PCR陽性ベース）がどれくらいか、つまり、1だと実際の報告数が潜在的な患者数と同一で、>1だと過剰に報告されている、<1だと過小報告されている、と考えるならば、0.44（95%CI 0.37-0.50）で、軽症患者については実際の患者数は2倍くらい（実際の0.44分しか報告されていない）だろう、ということらしい。
2020年2月28日までの疫学データから、RT-PCRで陽性確定となった患者がほぼ毎日厚生労働省のHPで見れるので、それを都道府県、10歳きざみの年齢、そして重症度でカウントする。重症度の定義は、酸素療法を要して肺炎もしくは挿管されている患者か、ICUに入室した患者、となっている。
厚生労働省のページをいちいち見に行ってもよいが、広島県のHPが時系列で厚生労働省の発表をまとめているので、そこから毎日データを見に行って確認した。2020年2月28日までは症例198まで番号がついている。論文ではNは言及がないが、図を見ると足すと200人前後っぽいのでたぶんいいのかもしれない。
患者の発生状況等 - 広島市公式ホームページ
都道府県ベースで10歳きざみのデータが報告されている。そして各自治体レベルで、重症かそうかの報告がされている。ただし、「胸部X線およびCTで肺炎像」というのがはたしてすべて重症肺炎かというとそうでもない。抗生剤さえ開始すれば酸素療法はいらない肺炎は多数存在する。そもそもコロナウイルスは抗生剤が効くようなものではないので、細菌性肺炎が合併するならまだしも、抗生剤はいらない。
重症度については各自治体で報告様式が異なっている。例えば東京都は、「症例○は重篤である」と書いてあるし、北海道では「症例○は非侵襲性呼吸補助療法を要した」「挿管された」などと書いてある。
図からはsevere な症例について21症例がみてとれるが、今回自分でデータを漁ったところ、12症例しかわからなかった。ただし、ただの肺炎（画像診断で肺炎、となっている症例）は除外したのでこれが過小カウントになってしまったかもしれない。
都道府県ごとの人口は平成26年度版の5歳階級データを拝借した。
https://www.e-stat.go.jp/dbview?sid=0003104197
モデルとしては、各都道府県の各年齢階級
の人口
について、非重症患者数
と重症患者数
（これらは各都道府県の各年齢階級のインデックスがつく）について、確率
のポアソン過程で患者が発生する、としている。
ここで、非重症患者は重症患者より の割合で患者報告数が多い（これは年齢階級のみのインデックス）、というパラメータをいれている。ここのパラメータは中国の初期段階での患者報告データから流用している。おそらくtable 1である。
Clinical Characteristics of Coronavirus Disease 2019 in China. - PubMed - NCBI
さて、肝心のascertainment rate は として、対数尤度関数を
と定義するが、これはポアソン分布の確率密度関数が
であり、非重症患者では、重症患者では
とすればよい。
ここまで出来たのでrstanでやってみる。
結論から言うとascertainment rate の推定はそれなりによかった。
2.5% 50% 97.5% 0.3640505 0.4345467 0.5170327
しかし、肝心の推定患者数はグダグダだった。非重症患者はもとより、重症患者が異常に多く推定されてしまった。
グダグダだった理由として、2020年2月28日までの都道府県別年齢階級別のデータがほんとうに論文で解析したデータとあっているのか不明だし、 の値も不明だし、重症患者の定義が
severe dyspnea that required oxygen support plus pneumonia or intubation
と書いてあるが、肺炎であることが必ずしも重症ではないし、重症患者のカウントの仕方が不明だった。
こんな短い論文ですら再現出来ないのだから自分の技量は（ここで記事が途絶えている
# 都道府県別の年齢階級別の人口 pop <- read.table(text=" 10歳未満 10代 20代 30代 40代 50代 60代 70代 80代 北海道 397 462 511 643 738 693 853 641 461 青森県 95 123 109 147 172 182 210 163 120 岩手県 98 119 108 144 160 172 196 160 128 宮城県 191 217 256 305 316 296 329 239 180 秋田県 70 88 73 110 123 144 172 139 118 山形県 88 105 92 127 136 151 174 135 123 福島県 149 190 171 221 243 265 293 218 185 茨城県 239 279 281 363 411 366 444 324 213 栃木県 165 186 189 253 280 254 301 207 145 群馬県 163 194 179 239 283 242 296 221 159 埼玉県 603 673 778 953 1137 858 1031 808 398 千葉県 505 560 636 798 953 738 907 715 385 東京都 1020 1041 1662 2072 2220 1583 1611 1339 841 神奈川県 760 824 999 1246 1497 1080 1186 957 547 新潟県 179 213 201 268 300 295 356 272 228 富山県 84 99 88 125 148 128 167 130 101 石川県 97 112 115 139 161 138 171 127 97 福井県 68 78 70 92 104 100 116 89 75 山梨県 66 83 79 94 118 108 122 95 76 長野県 176 204 171 245 289 258 306 250 211 岐阜県 174 203 190 240 283 246 303 238 165 静岡県 316 349 331 452 528 460 548 431 292 愛知県 682 727 835 1012 1147 858 976 773 446 三重県 154 177 170 219 257 224 263 211 149 滋賀県 135 145 154 185 202 167 191 139 97 京都府 208 238 301 321 371 294 374 303 200 大阪府 723 824 950 1131 1366 999 1231 1047 566 兵庫県 472 532 542 674 815 664 798 634 411 奈良県 110 134 133 156 193 166 210 168 105 和歌山県 75 91 81 102 130 123 150 125 94 鳥取県 48 55 49 67 71 74 88 65 60 島根県 57 63 54 76 82 87 110 85 82 岡山県 165 184 194 230 255 223 279 220 173 広島県 247 265 276 346 396 332 417 319 234 山口県 111 128 118 155 178 167 229 182 143 徳島県 58 68 65 87 96 97 122 92 78 香川県 82 92 84 116 130 118 154 114 93 愛媛県 112 128 117 159 180 175 219 167 138 高知県 54 66 57 81 92 92 119 92 84 福岡県 455 472 557 655 688 609 736 536 381 佐賀県 77 86 76 97 102 107 124 90 77 長崎県 117 132 117 150 171 185 216 163 136 熊本県 160 171 166 207 217 230 265 202 177 大分県 98 106 105 135 144 145 183 140 116 宮崎県 100 108 91 127 133 146 173 129 107 鹿児島県 148 157 143 188 195 224 250 190 172 沖縄県 166 162 157 188 195 182 168 116 85 ", check.name=FALSE ) * 1000 # 非重症患者 infec <- read.table(text=" 10歳未満 10代 20代 30代 40代 50代 60代 70代 80代 北海道 4 2 5 5 8 9 13 10 7 青森県 0 0 0 0 0 0 0 0 0 岩手県 0 0 0 0 0 0 0 0 0 宮城県 0 0 0 0 0 0 0 0 0 秋田県 0 0 0 0 0 0 0 0 0 山形県 0 0 0 0 0 0 0 0 0 福島県 0 0 0 0 0 0 0 0 0 茨城県 0 0 0 0 0 0 0 0 0 栃木県 0 0 0 0 0 0 1 0 0 群馬県 0 0 0 0 0 0 0 0 0 埼玉県 0 0 0 1 0 0 0 0 0 千葉県 0 0 2 0 1 2 4 2 0 東京都 0 0 1 2 4 7 2 8 2 神奈川県 0 0 3 2 2 6 2 1 4 新潟県 0 0 0 0 0 0 1 0 0 富山県 0 0 0 0 0 0 0 0 0 石川県 0 0 0 0 0 2 1 0 0 福井県 0 0 0 0 0 0 0 0 0 山梨県 0 0 0 0 0 0 0 0 0 長野県 0 0 0 0 0 0 1 0 0 岐阜県 0 0 0 0 0 1 0 0 0 静岡県 0 0 0 0 0 0 1 0 0 愛知県 0 0 1 0 2 2 14 8 1 三重県 0 0 0 0 0 1 0 0 0 滋賀県 0 0 0 0 0 1 0 0 0 京都府 0 0 2 0 0 0 0 0 0 大阪府 0 0 0 0 5 1 0 0 0 兵庫県 0 0 0 0 0 0 0 0 0 奈良県 0 0 0 0 0 0 1 0 0 和歌山県 0 0 0 1 1 5 2 1 1 鳥取県 0 0 0 0 0 0 0 0 0 島根県 0 0 0 0 0 0 0 0 0 岡山県 0 0 0 0 0 0 0 0 0 広島県 0 0 0 0 0 0 0 0 0 山口県 0 0 0 0 0 0 0 0 0 徳島県 0 0 0 0 0 0 0 0 0 香川県 0 0 0 0 0 0 0 0 0 愛媛県 0 0 0 0 0 0 0 0 0 高知県 0 0 0 1 0 0 0 0 0 福岡県 0 0 0 0 0 0 2 0 0 佐賀県 0 0 0 0 0 0 0 0 0 長崎県 0 0 0 0 0 0 0 0 0 熊本県 0 0 1 0 0 2 2 0 0 大分県 0 0 0 0 0 0 0 0 0 宮崎県 0 0 0 0 0 0 0 0 0 鹿児島県 0 0 0 0 0 0 0 0 0 沖縄県 0 0 0 0 0 0 2 0 1 ", check.name=FALSE ) # 重症患者 infec_severe <- read.table(text=" 10歳未満 10代 20代 30代 40代 50代 60代 70代 80代 北海道 0 0 1 0 0 1 0 0 2 青森県 0 0 0 0 0 0 0 0 0 岩手県 0 0 0 0 0 0 0 0 0 宮城県 0 0 0 0 0 0 0 0 0 秋田県 0 0 0 0 0 0 0 0 0 山形県 0 0 0 0 0 0 0 0 0 福島県 0 0 0 0 0 0 0 0 0 茨城県 0 0 0 0 0 0 0 0 0 栃木県 0 0 0 0 0 0 0 0 0 群馬県 0 0 0 0 0 0 0 0 0 埼玉県 0 0 0 0 0 0 0 0 0 千葉県 0 0 0 0 0 0 0 0 0 東京都 0 0 0 0 0 2 1 1 1 神奈川県 0 0 0 0 0 0 0 0 1 新潟県 0 0 0 0 0 0 0 0 0 富山県 0 0 0 0 0 0 0 0 0 石川県 0 0 0 0 0 0 0 0 0 福井県 0 0 0 0 0 0 0 0 0 山梨県 0 0 0 0 0 0 0 0 0 長野県 0 0 0 0 0 0 0 0 0 岐阜県 0 0 0 0 0 0 0 0 0 静岡県 0 0 0 0 0 0 0 0 0 愛知県 0 0 0 0 0 0 0 0 0 三重県 0 0 0 0 0 0 0 0 0 滋賀県 0 0 0 0 0 0 0 0 0 京都府 0 0 0 0 0 0 0 0 0 大阪府 0 0 0 0 0 0 0 0 0 兵庫県 0 0 0 0 0 0 0 0 0 奈良県 0 0 0 0 0 0 0 0 0 和歌山県 0 0 0 0 0 0 0 1 0 鳥取県 0 0 0 0 0 0 0 0 0 島根県 0 0 0 0 0 0 0 0 0 岡山県 0 0 0 0 0 0 0 0 0 広島県 0 0 0 0 0 0 0 0 0 山口県 0 0 0 0 0 0 0 0 0 徳島県 0 0 0 0 0 0 0 0 0 香川県 0 0 0 0 0 0 0 0 0 愛媛県 0 0 0 0 0 0 0 0 0 高知県 0 0 0 0 0 0 0 0 0 福岡県 0 0 0 0 0 0 0 0 0 佐賀県 0 0 0 0 0 0 0 0 0 長崎県 0 0 0 0 0 0 0 0 0 熊本県 0 0 0 0 0 0 0 0 0 大分県 0 0 0 0 0 0 0 0 0 宮崎県 0 0 0 0 0 0 0 0 0 鹿児島県 0 0 0 0 0 0 0 0 0 沖縄県 0 0 0 0 0 0 1 0 0 ", check.name=FALSE ) # nejm のfa f1 <- rep(c(8, 490, 241, 109), c(2, 3, 2, 2))/848 f2 <- rep(c(1, 67, 51, 44), c(2, 3, 2, 2))/163 fa <- f1/f2 # http://www.ourphn.org.au/wp-content/uploads/20200225-Article-COVID-19.pdf # の論文にある fa # これを使ってもいい再現にはならない f1 <- c(0.01,1,7,18,38,130,309,213,208) # death f2 <- c(416,549,3619,7600,8571,10008,8583,3918,1408) # confirmed case fa <- f2/f1 library(vioplot) library(rstan) rstan_options(auto_write=TRUE) options(mc.cores=parallel::detectCores()) code <- " data{ int<lower=0> A; // age group int<lower=0> X; // prefecture matrix<lower=0>[X, A] N; // population int<lower=0> Dn[X, A]; // non-severe int<lower=0> Ds[X, A]; // severe matrix<lower=0>[X, A] f; // ratio of non-severe to severe } parameters{ matrix<lower=0, upper=0.3>[X, A] p; real<lower=0> k; } transformed parameters{ matrix<lower=0>[X, A] lambda1; matrix<lower=0>[X, A] lambda2; lambda1 = (k * N .* p) .* f; lambda2 = (N .* p); } model{ for(a in 1:A){ for(x in 1:X){ Dn[x, a] ~ poisson(lambda1[x, a]); Ds[x, a] ~ poisson(lambda2[x, a]); } } } " m0 <- stan_model(model_code=code) standata <- list(X=nrow(pop), A=ncol(pop), N=pop, Dn=infec, Ds=infec_severe, f=t(replicate(nrow(pop), fa))) fit <- sampling(m0, standata, iter=10000, warmup=3000, chain=4) ex <- extract(fit, pars=head(fit@model_pars, -1)) alpha <- 0.05 cia <- c(alpha/2, 0.5, 1-alpha/2) quantile(ex$k, cia) # k の推定区間 est <- list("non-severe"=t(mapply(function(z) colSums(ex$lambda1[z,,]), 1:dim(ex$lambda1)[1])), severe=t(mapply(function(z) colSums(ex$lambda2[z,,]), 1:dim(ex$lambda2)[1]))) Ninfec <- lapply(list(infec, infec_severe), colSums) xl <- c(0.5, ncol(pop)+0.5) yl <- c(0, max(unlist(est), unlist(Ninfec))) par(mfrow=c(2, 1), las=1) for(i in seq(Ninfec)){ plot(Ninfec[[i]], type="n", col="red", xaxt="n", xlim=xl, ylim=yl, xlab="Age group", ylab="count") vioplot(as.data.frame(est[[i]]), ylim=yl, add=TRUE) points(Ninfec[[i]], pch=15, col="red") axis(1, at=seq(ncol(pop)), labels=colnames(pop)) legend("topleft", legend=c("Estimate", "Data"), pch=15, col=c(grey(0.3), "red")) title(names(est)[i]) }