数学 統計に詳しい人が語るコロナウイルス

東大数学科卒の元官僚はこう分析してるが、お前らはどうなると思う？

http://www.zakzak.co.jp/soc/news/200220/dom2002200003-n2.html

中国国外感染者の中国国内との比率をみると、
１月２０日の数字公表以降は、０・８〜２・６％で比較的安定している。
これは、新型肺炎の感染者のほとんどは中国国内、それも湖北省に集中しているからだ。
ちなみに中国国外での感染者数は、中国国内の１・１％だ（２月１６日現在）。
本コラムで紹介したが、現時点では、最終的な中国国内の感染者数は２０万人超と筆者は推計している。
となると、中国国外の感染者は数千人程度になるだろう。
中国国外のうち日本の比率は１割弱なので、日本の感染者数は数百人程度であろう。
その場合、死者も数人から１０人程度になるだろう。

こうした推計をすると、今の感染者は氷山の一角だと思われるが、今後の増加ペースはどうなるだろうか。
新型コロナウイルスの検査は簡単に行えるので、今後、日本での感染者数は増えていくだろう。
ある時点ではそれがネズミ算的に増えるかのように思える局面もあるだろうが、
筆者の推計が正しければ、現時点ではせいぜい数百人が一つのメドだ。

>>397
混ぜて検査する方法の最適化問題ってのもあるね。

https://mobile.twitter.com/p_gotcha/status/1243501000943702017
https://twitter.com/5chan_nel (5ch newer account)

>>401
これって東京の無作為サンプリングではなさそうだよね
病院に来たって事は熱が出てたのかもしれないし
想定できる母集団ってなんになるのだろうか

新型コロナ患者を治療している病院に100人の職員がいる。
検体採取器具は10人分、試薬は1回分しかないとする。
無作為抽出した10人の職員から採取した検体を混合して検査したら陰性であった。
職員の陽性者数の期待値を求めよ。
また、50人以上の感染者いる確率はいくつか？
検査の陽性率はハイリスク群に検査している東京の数値2457/6654を使って計算せよ。

>>401
95%CIで
> binom::binom.confint(4,67)
method x n mean lower upper
1 agresti-coull 4 67 0.05970149 0.019131154 0.1480232
2 asymptotic 4 67 0.05970149 0.002968439 0.1164345
3 bayes 4 67 0.06617647 0.015026904 0.1253507
4 cloglog 4 67 0.05970149 0.019283398 0.1337560
5 exact 4 67 0.05970149 0.016504404 0.1458632
6 logit 4 67 0.05970149 0.022588780 0.1485238
7 probit 4 67 0.05970149 0.020905075 0.1402573
8 profile 4 67 0.05970149 0.018970462 0.1332788
9 lrt 4 67 0.05970149 0.018929939 0.1332756
10 prop.test 4 67 0.05970149 0.019297952 0.1534709
11 wilson 4 67 0.05970149 0.023459351 0.1436950

>>401
まじか。
まあ、体調悪いから病院に行くわけで、バイアスかかってるとはいえ...。

>>404
ラテン方陣の問題？

>>403
事前分布にJefferey分布を使っているな。

https://i.imgur.com/m7kdY8Z.png

破線が事前分布、実戦が事後分布

curve(dbeta(x,0.5+4,0.5+67-4),bty='l',xlab='probability',ylab='density')
curve(dbeta(x,0.5,0.5),add=T,lty=2)

>>410
青実線が事前分布を一様分布(Beta(1,1))としたとき。

https://i.imgur.com/qUF1Fue.png

Jeffereyの方が95%CI幅が小さいな。

> binom::binom.bayes(4,67,prior.shape1 = 0.5,prior.shape2 = 0.5)
method x n shape1 shape2 mean lower upper sig
1 bayes 4 67 4.5 63.5 0.06617647 0.0150269 0.1253507 0.04999999
> binom::binom.bayes(4,67,prior.shape1 = 1, prior.shape2 = 1)
method x n shape1 shape2 mean lower upper sig
1 bayes 4 67 5 64 0.07246377 0.01876916 0.1338218 0.04999999

>>401

https://georgebest1969.typepad.jp/blog/2020/04/%E6%85%B6%E5%BF%9C%E3%81%AEpcr6%E3%81%AE%E6%84%8F%E5%91%B3.html

に準じて　東京都民の1395万人に当てはめると

> data.frame(method=ci[,1],round(ci[,4:6]*pop))
method mean lower upper
1 agresti-coull 832836 266880 2064924
2 asymptotic 832836 41410 1624262
3 bayes 923162 209625 1748642
4 cloglog 832836 269003 1865896
5 exact 832836 230236 2034792
6 logit 832836 315113 2071907
7 probit 832836 291626 1956590
8 profile 832836 264638 1859240
9 lrt 832836 264073 1859195
10 prop.test 832836 269206 2140919
11 wilson 832836 327258 2004545

岩田の計算は
5 exact 832836 230236 2034792

>>402
感度30-70％(最頻値0.5,標準偏差0.2のβ分布),特異度(最頻値0.9 標準偏差0.05のβ分布)に設定。
有病率の事前分布は0-1の一様分布にして　

MCMCしてみると

https://i.imgur.com/VfDTj51.png

という結果になった。

有病率の信頼区間は広すぎｗ

mean lower upper
0.22346412787 0.00000002913 0.83202346988

>>413
有病率の最頻値は
> density(prev)$x[which.max(density(prev)$y)]
[1] 0.019186

>>413
事前分布をJeffereyにしたら、

> js=PCRj4(67,4,SEN=0.5,SD1=0.2,SPC=0.9,SD2=0.1,N.ITER=1e6)$js
mean lower upper
1.4972e-01 6.3120e-14 8.5733e-01

> summary(prev)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.0111 0.0472 0.1497 0.1442 1.0000

> density(prev)$x[which.max(density(prev)$y)] # mode
[1] 0.0032278

>>415
この最頻値はオーストリアの0.3％という値に等しいな。

>>415
Jeffreys が正しいスペリングみたい。　Jefferey'sかと思っていた。

某医療センターでは治療後、2回連続してPCR検査陰性であれば退院させているとする。
PCR検査は
感度0.3~0.7 : β分布(2.625,2.625)相当
特異度 0.95~1.0 : β分布(26.5014,2.3422)相当
有病率は一様分布として
2回連続してPCR検査陰性の患者の有病率の期待値を求めよ。

２回連続検査陰性の患者が感染者である確率と有病率との関係をグラフにしてみた。
灰色点線は95％信頼区間

https://i.imgur.com/2U5HBmb.png

>>409
動的計画法を使うらしい。

https://qiita.com/p_gotcha/items/e2741b12b4f606b96136

スレチですが統計に詳しい方がいると思って伺いました。
工学系の大学院生で論文を読んでて、見慣れない記号や数式が出てきたので質問させてください。
この数式中のEは何を意味するのでしょうか？
Rの期待値的なものですか？

https://i.imgur.com/oVwva8r.jpg

院生でわからないってマジか
専門外の機械学習いきなりやらされたのかな？

有病率が0.3％（オーストリアの無作為抽出)や6%(慶応大学の入院患者無作為抽出）のときは

感度30〜70%　特異度90〜100%のPCR検査キットでは

陰性的中率は
オーストリアの例では　95％[87％〜99％]
慶応の例では 94%[86〜99％]になる。

すべて陰性にでる中国のイカサマキットなら
陰性的中率は
オーストリアの例では99.7％
慶応の例では94%になる。

有病率が10％くらいになればイカカマキットの方が陰性的中率の成績が劣る。

某医療センターでは治療後、2回連続してPCR検査陰性であれば退院させているとする。
PCR検査は
感度0.3~0.7 : β分布(2.625,2.625)相当
特異度 0.95~1.0 : β分布(26.5014,2.3422)相当
として

ｎ回連続して陰性であれば退院とするとどれくらいの感染者が野に放たれるかを
有病率（＝検査前確率）を変化させてグラフにしてみた。


https://i.imgur.com/fvNvXhX.png

https://www.researchgate.net/publication/340869568_Saliva_Sample_as_a_Non-Invasive_Specimen_for_the_Diagnosis_of_Coronavirus_Disease-2019_COVID-19_a_Cross-Sectional_Study/link/5ea19006458515ec3aff8f36/download
https://i.imgur.com/3GWmZz3.png
に有病率の低い群で行った唾液と鼻咽頭スワッブ200例のcross-section tableがあったので
これでKappa係数とPABAK（prevalence ad-justed bias adjusted kappa）とその信頼区間をstanのMCMCで出してみた。
swabと唾液でまあ、結果が合致している。

> print(fit.kappa,pars=c('kappa','PABAK'))
Inference for Stan model: kappa.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
kappa 0.81 0 0.07 0.65 0.77 0.81 0.86 0.92 3540 1
PABAK 0.93 0 0.02 0.88 0.92 0.93 0.95 0.97 3342 1

JAGSでMCMCしてもほぼ同じ結果。
https://i.imgur.com/XKT53Fy.png

>>394
答えなくね？

>>425
正解。正直者が誰もいないようにプログラムした。
イナ師が微分で解いてくれるのを期待していたんだが。

PCR検査は
感度0.3~0.7 : β分布(2.625,2.625)相当
特異度 0.95~1.0 : β分布(26.5014,2.3422)相当
として、
http://statmodeling.hatenablog.com/entry/covid19-estimate-total-number-of-positives-in-japan
の設定を踏襲して、利用できるデータを更新して

推定陽性率は
> summary(p.pos) ; HDInterval::hdi(p.pos)[1:2]
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000031 0.000067 0.004681 0.005229 0.009325 0.016658
lower upper
0.00003082 0.01408833

推定有病率は
> summary(prevalence) ; HDInterval::hdi(prevalence)[1:2]
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00000 0.00166 0.00529 0.01040 0.01129 0.16078
lower upper
0.0000005693 0.0305640211

>>426
正直者の確定が誰もいないというより、条件を満たす正直者・嘘つきの組み合わせが存在しないよね？
https://i.imgur.com/NJw3GJN.jpg

rRT-PRC検査の感度・特異度ともに90%以上です。

通常の医療でも臨床診断の5%は誤診だと推定されていたりする。
これを実際の数として見積もると膨大になる。

>>429
採取手技や採取時期で変動すると思う。
コンタミによる偽陽性もあるだろうし。

>>428
言っていることは同じ。
全部の条件を満たす組み合わせは存在しない。

>>427
このデータは再現性がないことがわかったので撤回します。

ふぉー

＜新型コロナ＞抗体検査５．９％陽性　市中感染の可能性　都内の希望者２００人調査
https://www.tokyo-np.co.jp/s/article/2020043090070748.html

>>434

"検査結果では、一般市民の百四十七人の４・８％にあたる七人が陽性、
医療従事者五十五人のうち９・１％の五人が陽性だった。
https://www.tokyo-np.co.jp/s/article/2020043090070748.html

"

> r1=5;r2=7;n1=55;n2=147
> prop.test(c(r1,r2),c(n1,n2))

2-sample test for equality of proportions with continuity
correction

data: c(r1, r2) out of c(n1, n2)
X-squared = 0.679, df = 1, p-value = 0.41
alternative hypothesis: two.sided
95 percent confidence interval:
-0.052613 0.139194
sample estimates:
prop 1 prop 2
0.090909 0.047619

Warning message:
In prop.test(c(r1, r2), c(n1, n2)) :
Chi-squared approximation may be incorrect
> poisson.test(c(r1,r2),c(n1,n2))

Comparison of Poisson rates

data: c(r1, r2) time base: c(n1, n2)
count1 = 5, expected count1 = 3.27, p-value = 0.33
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
0.47778 6.98763
sample estimates:
rate ratio
1.9091

>>434
これの興味深いところは、陽性率の高さもさることながら、医療関係者の感染率が一般市民のそれよりも
（多分有意に）高いことだと思う。抗体検査の精度ってほとんどデータが無いに等しいけど、医療関係者の
感染者率が平均よりも多いのならば（実際そうなのだと思う）、このテストはそれを反映したものと言えて、
抗体検査の信頼性をある程度証明できるかもしれない。
そして、市中の陽性率の高さも。

この辺、ベイズで上手く検証できないかな。

>>436
>医療関係者の感染率が一般市民のそれよりも（多分有意に）高い

χ二乗検定では、p-value = 0.41 で　標本数が少ないから有意差がでない。
エントリーが５以下のときには信頼するなと習ったな。

まあ、ポアソンでもp-value = 0.33

>435にRの出力を貼っておいた。

>>436
55人中5人と147人中7人じゃ、有意差ないだろ。

と書こうと思ったら、 >>437氏に先を越されたｗｗｗ

抗体検査キットの評価をしたら特異度はみな5/5だったらしいけど、サンプルが
５つだけって意味だとすると、せいぜい90％以上くらいのことしか言えないな。
特異度95%だとしたら、５％が陽性っていわれてもねぇｗ

>>436
事前分布にかなり影響をうけるが、

感度特異度とも50-70％(最頻値60%標準偏差10％のβ分布),
有病率は一様分布、
検査陽性数は陽性確率が　有病率*感度＋（１-有病率)*(1-特異度)の二項分布に従う

というモデルでプログラムを組むと

model
{
for(i in 1:N) {
x[i] ~ dbin(p,n[i])
}
p = prev*sen + (1-prev)*(1-spc)
PPV=sen*prev/(sen*prev+(1-prev)*(1-spc))
NPV=(1-prev)*spc/((1-prev)*spc+prev*(1-sen))
precision=(prev*sen+(1-prev)*spc)/
((prev*sen+(1-prev)*spc + (1-prev)*(1-spc)+(prev*(1-sen))))
pLR=sen/(1-spc)
nLR=(1-sen)/spc
DOR=pLR/nLR
sen ~ dbeta(sn[1],sn[2])
spc ~ dbeta(sp[1],sp[2])
prev ~ dbeta(shape1,shape2)
}

結果は、
https://i.imgur.com/eKXLUXZ.png

有病率の期待値は2.3%、最頻値は0.31%　少数データなので信頼区間が広い。
有病率が平均値50%の一様分布というのは現実離れした分布だとは思う。

>>439
すべてが陰性に出るイカサマキットなら特異度が100％
>18参照

>55人中5人と147人中7人じゃ、有意差ないだろ。

その比率のまま、約4倍(3.87倍)になればχ二乗検定で有意差がでるね。
r1=5;r2=7;n1=55;n2=147
mat=matrix(c(r1,r2,n1,n2),2,b=T)
fn <- function(x) chisq.test(mat*x)$p.value
fn=Vectorize(fn)
uniroot(function(x) fn(x)-0.05,c(1,10))$root
fn(4)

> fn(4)
[1] 0.045678

>>440（追記)
>有病率が平均値50%の一様分布というのは現実離れした分布だとは思う。
オーストリアのデータ0.3%を参考に有病率の事前分布が0.2-0.4％
β(10.865, 3279.34)に相当として、MCMCしてみた。
感度と特異度の事前分布は前回と同じ。

https://i.imgur.com/dHcOxsR.png

さほど、有病率が高いという結論は引き出せないな。

自己流な検定

A群55人中5人とB群147人中7人に有意差はあるか？

モンテカルロ法で検証

A群もB群も同じ陽性率で、0.0594と仮定する
∵ 陽性率は、 12/(55+147) = 12/202 = 0.0594

モンテカルロ法100万回したら
B群の陽性者が7名となったのは、125238回
その内 A群の陽性が5名以上は、28293回

∴ P(55人中5人以上｜147人中7人) = 28293/125238 = 0.226
∴ 有意差なし

EXCEL VBAソースコード概略

Dim I As Long
Dim J As Long
Dim K As Long

K = 1
For I = 1 To 1000000
A = 0
For J = 1 To 55
If Rnd(1) < 0.0594 Then '陽性
A = A + 1
End If
Next

B = 0
For J = 1 To 147
If Rnd(1) < 0.0594 Then '陽性
B = B + 1
End If
Next

If B = 7 Then 'B群の陽性者が7名の場合
Cells(K, "A") = A 'A群の陽性者の人数
K = K + 1
End If
B = 0

Next

β分布と乱数発生で比較。

一般人と医療従事者の確率分布に従う乱数を１０００万個発生させて比率の分布をグラフにすると
９５％信頼区間が１を挟むから有意差なし。比率が１以上となる確率も87.5％で95%を越えないから有意差なし。
https://i.imgur.com/8zFKODI.png

a=0.5 ; b=0.5
r1=5 ; r2=7 ; n1=55 ; n2=147
layout(1)
layout(matrix(1:2,2))
curve(dbeta(x,a+r2,b+n2-r2),0,0.3,bty='l',ann=F,lwd=2)
curve(dbeta(x,a+r1,b+n1-r1),col=2,add=T,lwd=2)
legend('center',bty='n',legend=c('general','medical'),lwd=2,col=1:2)
k=1e7
general=rbeta(k,a+r2,b+n2-r2)
medical=rbeta(k,a+r1,b+n1-r1)
BEST::plotPost(medical/general,compVal = 1)

>>441
すべてが陰性に出るなら陽性は０人だろｗ

>>446
偽陽性率０％だから特異度は１００％

>>447
だから、陽性が何人か出てんだから感度0%のインチキじゃなかろうってこと。

最近の報道

インドは、中国企業から購入した中共ウイルス（新型コロナウイルス）の迅速スクリーニング検査キットの精度がわずか5％だとして、約50万個の注文をキャンセルした。

....
神戸市中央区にある市立医療センター中央市民病院の医師などのグループは、ことし３月末から先月７日にかけて、新型コロナウイルス以外の理由で外来を受診した患者から無作為に1000人を選び、血液中に新型コロナウイルスに感染したあとにできる「抗体」があるか調べました。
グループによりますとその結果、3.3％にあたる33人から抗体が検出されたということです。

以上を知ったある会社が
必ず陽性がでる試薬33個と必ず陰性がでる試薬967個を混ぜた1000試薬をセットにしてインドに売り込んだ。

問題

1試薬は1回しか検査できないとして
これがイカサマキットであることを証明する手段はあるか？

3月の宿題で(1)のみ正解の数弱@shukudai_sujaku

昨年度の大学への数学(大数)での勝率は、

学コンBコースが 1/1 = 100% ，

宿題が 3/10 = 30% でした！

宿題の勝率が低すぎると思うので、

これからは一層精進していきたいです！

https://twitter.com/shukudai_sujaku
https://twitter.com/5chan_nel (5ch newer account)

高橋洋一(統計数理研究所→大蔵省)
スウェーデンの感染症対策を解説
https://youtu.be/2bsYRuorwMI?t=192

https://toyokeizai.net/sp/visual/tko/covid19/
のデータを使って
P=14895/153581 　#　2020/05/04
PCR検査の感度30-70％のモデルM1と感度70-90％のモデルM2のどちらが信憑性があるか、ベイズファクターで計算してみる。
M1は感度が最頻値60%標準偏差10%、M2は最頻値80%標準偏差10%のベータ分布に設定
特異度はいずれも最頻値95%標準偏差2.5%に設定し、有病率は一様分布を仮定
陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする。
事後確率分布は
https://i.imgur.com/81bK4KE.png
陽性率P=14895/153581=0.09698465での事後確率分布の密度比(Savage-Dickey density ratio)でベイズファクターを出すと
> d1/d2 # Savage-Dickey densiti ratio = BF12
[1] 1.007722

まあ、ちょっぴり、感度30-70％のモデルの方がいいかも、という結果。
陽性数/検査数の時系列データでもあればもう少し差がでるかもしれん。
東京都のデータで計算させようかと思ったが、東京都は検査人数を隠蔽しているので使いものにならない。

訂正

陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする
↓
陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする

3月の宿題で(1)のみ正解の数弱@shukudai_sujaku

昨年度の大学への数学(大数)での勝率は、

学コンBコースが 1/1 = 100% ，

宿題が 3/10 = 30% でした！

宿題の勝率が低すぎると思うので、

これからは一層精進していきたいです！

https://twitter.com/shukudai_sujaku
https://twitter.com/5chan_nel (5ch newer account)

>>449
抗体のはいったサンプルを２，３０個試してみりゃわかんじゃないの？
いくらなんでも感度が３％じゃつかえねぇ。

イカサマキットの感度特異度の事前分布を一様分布に設定して
抗体のはいったサンプルを２０個で全部陰性であったので３０個試したら全部陰性であったとすると
感度・特異度の事後分布は

https://i.imgur.com/6ZYn34k.png

ところが、有病率33/1000であったときに無作為に２０人および３０人を選んで全部陰性であっても
感度・特異度の事後分布は
https://i.imgur.com/jWg4RgH.png

>>454
事後分布を一様分布に設定して、このコピペ小僧の正解率の期待値・最頻値・95％信頼区間を求めよ。
その結果、
https://i.imgur.com/M9SdCSL.png

>>458
ベータ分布の理論値

> HDInterval::hdi(qbeta,shape1=5,shape2=8)[1:2]
lower upper
0.1406542 0.6377277
> 5/(5+8) # mean
[1] 0.3846154
> (5-1)/(5-1+8-1) # mode
[1] 0.3636364

>>454
このコピペ小僧の正解率に学コンと宿題で差があるかを検定せよ。
事前分布にJefferysを用いた結果、
https://i.imgur.com/2yYagtX.png

>>456
感度が5%以下じゃつかいもんにならんわなぁ。インチキで確定。

https://www.buzzfeed.com/jp/naokoiwanaga/covid-19-antibody-test
大阪市立大学が一昨年の血液を使って検査したところ、偽陽性は1人もいなかったそうだ。

>>463
Wilsonの式を使ったんだろうね。

method x n mean lower upper
1 agresti-coull 3 312 0.0096 0.0019 0.029
2 asymptotic 3 312 0.0096 -0.0012 0.020
3 bayes 3 312 0.0112 0.0016 0.023
4 cloglog 3 312 0.0096 0.0027 0.026
5 exact 3 312 0.0096 0.0020 0.028
6 logit 3 312 0.0096 0.0031 0.029
7 probit 3 312 0.0096 0.0029 0.027
8 profile 3 312 0.0096 0.0024 0.025
9 lrt 3 312 0.0096 0.0024 0.025
10 prop.test 3 312 0.0096 0.0025 0.030
11 wilson 3 312 0.0096 0.0033 0.028

https://ja.wikipedia.org/wiki/%E3%82%A6%E3%82%A3%E3%83%AB%E3%82%BD%E3%83%B3%E3%81%AE%E4%BF%A1%E9%A0%BC%E5%8C%BA%E9%96%93

Wilsonの式での95％信頼区間幅を1％以内にしたいなら
サンプルサイズは38411が必要という計算になった。

Rでの算出プログラム
library(binom)
x=3
n=312
binom.wilson(x,n)
fn <- function(n){
x=0:n
l=binom.wilson(x,n)[,5]
u=binom.wilson(x,n)[,6]
max(u-l)
}
fn=Vectorize(fn)
n=seq(1000,50000,by=1000)
plot(n,fn(n))
abline(h=0.01,lty=3)
uniroot(function(x,u0=0.01) fn(x)-u0, c(10000,50000))

>>462
５０人で０でしょ？
だったら、特異度として保証されるのはせいぜい９８％程度じゃん。
１％の陽性率に対して、それでは不十分。

>>466
95％信頼区間の下限境界は

事前分布をJeffreysで
> qbeta(0.95,0.5+50,0.5,lower=F)
[1] 0.9624989
事前分布を一様分布で
> qbeta(0.95,1+50,1,lower=F)
[1] 0.942952

特異度の95%の信頼区間の下限値を0.99にするのに必要なサンプルサイズは
事前分布を一様分布で
> uniroot(function(x) fn(x)-0.99, c(100,500))$root
[1] 297.0728
Jeffreysで
> uniroot(function(x) fn(x,0.5,0.5)-0.99, c(100,500))$root
[1] 190.8606

fn <- function(x,shape1=1,shape2=1){
qbeta(0.95,shape1 + x, shape2, lower=F)
}
n=50:500
plot(n,fn(n),type='l', ylab='95%CI.lower')
abline(h=0.99,lty=3)
uniroot(function(x) fn(x)-0.99, c(100,500))$root
uniroot(function(x) fn(x,0.5,0.5)-0.99, c(100,500))$root

https://www.buzzfeed.com/jp/naokoiwanaga/covid-19-antibody-test
のデータを使って、不明なものは一様分布（ベータ分布の形状母数(1,1)）に事前分布を設定してMCMCしてみる。

x=c(3,33)
n=c(312,1000)
m=50
N=length(n)
shape1=1
shape2=1

model{
for(i in 1:N){
x[i] ~ dbin(p,n[i]) # 二項分布
}
p <- prev*sen+(1-prev)*(1-spc) # 陽性=真陽性+偽陽性
sen ~ dbeta(shape1,shape2)
spc ~ dbeta(shape1+m,shape2)
prev ~ dbeta(shape1,shape1)
}

その結果
https://i.imgur.com/B7p825B.png

>>469
なにやってるかよくわからんので、見当外れの指摘かもしれんが、
大阪市大の3/312と神戸大の33/1000ってのは特異度も感度も異なる
であろう別種のキットによる検査結果なんで、一緒くたにしちゃ
まずいんでないの？

>>470
同一キットじゃないから、ご指摘の通り。
しかも神戸大の方では陰性検体での確認はされていないから、神戸大方の陽性率が高いのは偽陽性を含む可能性もあるね。

大阪市大だけのデータでやってみると。

https://i.imgur.com/njVQtRZ.png

>>471
使える情報が少なすぎて信頼区間が広すぎ。
有病率:0-87%
キットの感度:0-70%
ってコイントス変わらんな。

結局のところ、断定的な結論は出せないねということを数字で確認しているだけｗ

3月の宿題で(1)のみ正解の数弱@shukudai_sujaku

昨年度の大学への数学(大数)での勝率は、

学コンBコースが 1/1 = 100% ，

宿題が 3/10 = 30% でした！

宿題の勝率が低すぎると思うので、

これからは一層精進していきたいです！

https://twitter.com/shukudai_sujaku
https://twitter.com/5chan_nel (5ch newer account)

>>474
Rのlibrary binomを使って　binom::confint(3,312)で　>464の出力が得られる
信頼区間をグラフにすると
https://i.imgur.com/Mlvv7bB.png
点線は3/312
正規分布近似のasymptotic以外は非対称。

どれを使うべきか？　好きなのを使えばいい、と思う。
但し、値が負になったり１を越えたりするのは不採用の方が賢明だとは思う。

Wilson法は値が０や１に近くても信頼できるという人がいるけど、どうやって検証するのかはよくわからん。

>>476
binom::confintはbinom.confint(3,312)の間違い

library("binom")
binom.confint(3,312)

レムデシビルで初のRCTが出たんだけど、
https://www.thelancet.com/pdfs/journals/lancet/PIIS0140-6736(20)31022-9.pdf
レムデシビルRCTで有意差だせず、症例数不足で検出力不足とか。

確かに、効力がわずからなら検出力不足
> n1=158
> n2=78
> pwr.2p2n.test(h=c(0.8,0.5,0.2),n1=n1,n2=n2,sig.level = 0.05,power=NULL,
+ alternative = "two.sided")

difference of proportion power calculation for binomial distribution (arcsine transformation)

h = 0.8, 0.5, 0.2
n1 = 158
n2 = 78
sig.level = 0.05
power = 0.9999336, 0.9508568, 0.3037150
alternative = two.sided

NOTE: different sample sizes



"
Remdesivir group (n=158)Placebo group (n=78)Difference
Clinical improvement rates
Day 7 4 (3%) 2 (3%) 0.0% (-4.3 to 4.2)
Day 14 42 (27%) 18 (23%) 3.5% (-8.1 to 15.1)
Day 28 103 (65%) 45 (58%) 7.5% (-5.7 to 20.7)
"

どの週においても両群で症状改善率は不変
症状改善率はどちらも一様分布を事前分布に仮定して

症状改善率の差δ=:０というモデルとδ≠０というモデルでのベイズファクター（周辺尤度比＝エビデンス比）を求めると
> (BF01=d.post/d.prio)
[1] 1.01472
なので、どちらのモデルが得られたデータを説明するのに適しているかは判断できず、という結果になった。
検出力不足か、差がないのかには決着つかず。

library(pwr)
library(rjags)
library(polspline)
par(bty='l')

n1=158
n2=78
pwr.2p2n.test(h=c(0.8,0.5,0.2),n1=n1,n2=n2,sig.level = 0.05,power=NULL,
alternative = "two.sided")
pwr.2p2n.test(h=NULL,n1=n1,n2=n2,sig.level = 0.05,power=0.80,
alternative = "two.sided")

s1=c(4,42,103) # improvement with Remdesivir
s2=c(2,18,45) # improvement with placebo
# s1=103 ; s2=45 # Day 28
N=length(s1) # 3
shape1=1 # prior parameter in dbeta
shape2=1

data=list(n1=n1,n2=n2,s1=s2,N=N,shape1=shape1,shape2=shape2)

modelString='
model{
# data
for(i in 1:N){
s1[i] ~ dbin(p1,n1)
s2[i] ~ dbin(p2,n2)
}
# parameter
delta <- p1 - p2
# priors
p1 ~ dbeta(shape1, shape2)
p2 ~ dbeta(shape1, shape2)
# sampling from priors
pri.p1 ~ dbeta(shape1, shape2)
pri.p2 ~ dbeta(shape1, shape2)
pri.delta <- pri.p1 - pri.p2
}
'
writeLines(modelString,'tmp.txt')
n.iter=1e5 ; thin=1
jagsModel=jags.model('tmp.txt',data,inits = NULL,n.chains=4,n.adapt=n.iter/5)
update(jagsModel)
codaSamples=coda.samples(jagsModel,c('delta','pri.delta'),n.iter,thin)
# plot(codaSamples)
coda2summary(codaSamples)
js=as.data.frame(as.matrix(codaSamples))
BEST::plotPost(js$delta,compVal=0,xlab=bquote(delta))
fit.post=logspline(js$delta)
fit.prio=logspline(js$pri.delta)
curve(dlogspline(x, fit.post), -2,2, lty=2, xlab=bquote(delta),ylab='Density')
curve(dlogspline(x, fit.prio), add=T)
(d.post=dlogspline(0,fit.post))
(d.prio=dlogspline(0,fit.prio))
legend('topright', bty='n', legend=c('δ≠0','δ=0'), lty=1:2)
abline(v=0,col=8)
points(c(0,0),c(d.post,d.prio),pch=c(1,19))
(BF01=d.post/d.prio)

確率分布を図示すると

https://i.imgur.com/hfH7gZO.png

100例もやって差がでないのなら、はっきり言って効き目なしとかわらん。
むしろ副作用がこわい。

アビガンはどうなんだろうね。なぜか国内のまとまったデータがまだ出てこない。

まあ、軽症者が８割とかだと、早期投与の効き目を判定するには相当の
人数に使わないとわかんないだろうけど。

死亡率に関してはベイズファクターで推論できた。

Remdesivir group (n=158) Placebo group (n=78) Difference
Day 28 mortality 22 (14%) 10 (13%) 1.1(-8.1% to 10.3%)"

事前確率分布をどうするかだが、なんの情報もないので
実薬群も対照群も一様分布(もしくはJeffreys分布)とする。

死亡率の差δの事前分布と事後分布の確率分布曲線を描いてδ=0での確率密度比(Savage-Dickey density ratio)が
ベイズファクター（周辺尤度の比)になる。　

これをJAGSを用いたプログラムで計算すると、一様分布のとき8.34 Jeffrey分布のとき6.40となる。、
10を超えないいのでまあ、中程度の確信で検出力不足によるのではなく、もともと差がないのであろうと推論できる。

サンプルサイズを増やしてもレムデシビルは軽症・早期投与でも致死率を改善しないであろうと予想。

>どの週においても両群で症状改善率は不変
この前提はおかしいので週ごとに症状改善率が変わるようにモデルを変更

model{
# data
for(i in 1:N){
s1[i] ~ dbin(p1[i],n1)
s2[i] ~ dbin(p2[i],n2)
}

# parameter
alpha <- delta*sigma
for(i in 1:N){
p1[i] <- phi(x1[i]) # probit:pnorm(x)
p2[i] <- phi(x2[i])
# p1[i] <- ilogit(x1[i]) # logit:exp(x)/(1+exp(x))
# p2[i] <- ilogit(x2[i])
}
# model
for(i in 1:N){
x1[i] ~ dnorm(mu + alpha/2, pow(sigma,-2)) # alpha:difference of mean
x2[i] ~ dnorm(mu - alpha/2, pow(sigma,-2))
}

# priors
delta ~ dnorm(0,1)
mu ~ dnorm(0,1)
sigma ~ dt(0,1,1)T(0,)
}

ベイズファクターは
> (BF01=d.post/d.prio)
[1] 0.9976274

症状改善率に差がないというモデルも差があるというモデルも周辺尤度（エビデンス）はほぼ同じという結果。

あ、 >>482はレムデシじゃなくてアビガンについての話だかんね。
レムデシは重症者向けらしいけど、たぶん駄目だろ。

【政府が隠す３．１１】　地震波形が、核実験と同じ
http://rio2016.5ch.net/test/read.cgi/earth/1575170033/l50

5/10夜『NHKスペシャル』「新型コロナウイルス　出口戦略は」で
厚労省新型コロナ対策班西浦博北大教授と氏が分析中のノートPC画面が
映されていた。氏の使用解析ソフトウェア名は何というのでしょうか？

西浦先生はRじゃないかな。
今週ぐらいに、実効再生産数を算出するRのコードを公開できるかもと言ってた。

火曜にはニコ生があるから、そこでたっぷりと聞けるはず。
https://sp.live.nicovideo.jp/watch/lv325833316

>>488
>11に数式があるけど

>>489
いや、与えられた何日分かの新規患者数から再生産数とか回復率とかを推定するための数式じゃないの？

>>490
SEIRより単純なReed-Frostモデルだよな
https://en.m.wikipedia.org/wiki/Reed%E2%80%93Frost_model

>>489
自分で計算してみた人ならわかるけど、日々の再生算数を出すにはA「感染してから他人にうつすまでの日数の分布」が必要。西浦先生はこれをまず求めていて、その結果は割合と世界中で使われている。
加えてB「感染してから発表の数字に載るまでの日数の分布」も必要で、発表の数字をBで逆畳み込み積分してC「感染推定日毎の数字」に変換し、CをAで畳み込み積分したものでCを割れば日々のRtが求められる。
とは言うものの、そのままやると誤差で数字が発散するし、逆畳み込みは一意ではないのでどう推定するかは色々と考えなくてはならない。多分日本独自の事情もコードに入れなきゃいけない。大変だろうなあとは思う。

>>492
レスありがとうございます。
確かにSEIRモデルやSIRIモデルだと
「感染してから発表の数字に載るまでの日数の分布」って考慮されてないね。

>>488
公開されたぞー。StanとRだ。

https://nbviewer.jupyter.org/github/contactmodel/COVID19-Japan-Reff/blob/master/scripts/C.%20Calculating%20the%20Rt%20in%20Stan.ipynb

ウイルス感染のモデルも必要だと思うけど
社会経済活動の影響や
自粛解除した後の死者増加や
自粛継続した場合の法人の倒産や
治療薬やワクチンができた後の社会の経済や文化面の回復
そこまで考慮してどういう政策を選択したら
損失を最小化できるかという問題は解ける？

制約付きの最小化問題だと思う
数値化するのが難しい考慮要素もあるからその辺は仮の値とかで

>>495
モデルを精密複雑化すればするほど、必要なパラメータが増えて
それを弱情報事前分布にすると事後分布の信頼区間が広がるんだなぁと思う。

>>494
これにある、
JapaneseDataCOVID19
ってどこかからダウンロードできるんだろうか？

>>497

https://github.com/contactmodel/COVID19-Japan-Reff/tree/master/data

>>498
ありがとうございます。

親ディレクトリ（フォルダ）探せばよかったんだね。

https://nbviewer.jupyter.org/github/contactmodel/COVID19-Japan-Reff/tree/master/

>>495 >>496
パラメータを増やせば増やすほど、どんな仮説も否定できなくなる。
実験系の学界でよく使う手口だよね。

おや、誰か来たようだ…

>>494
これを走らせてみたい人いますか？

> datestar = as.Date("2020-05-10")
> datemin = as.Date("2019-12-25") # particular choice
> (tstar = as.numeric(datestar - datemin))
[1] 137

> (K = nrow(df_cases)) # 147
[1] 147

なので、
data {
int<lower = 1> K; //number of days
vector<lower = 0>[K] imported_backproj;
vector<lower = 0>[K] domestic_backproj;
int<lower = K> upper_bound;
の

int<lower = K> upper_bound;
だ
Error in new_CppObject_xp(fields$.module, fields$.pointer, ...) :
Exception: modelf2b01ab215_fit_infection_namespace::modelf2b01ab215_fit_infection: upper_bound is 137, but must be greater than or equal to 147 (in 'modelf2b01ab215_fit_infection' at line 53)

upper_bound 137の下限を K 147にしているからみたい。

upper_boundの制限を外して

data {
// int<lower = K> upper_bound;
int upper_bound;

再生数の平均値を以下の出すブロックを加えて走らせてみた。
transformed parameters{
real mean_Rt;
real mean_Rt_adj;
mean_Rt = mean(Rt);
mean_Rt_adj = mean(Rt_adj);
}

その結果、
Inference for Stan model: fit_infection.
2 chains, each with iter=10000; warmup=2000; thin=5;
post-warmup draws per chain=1600, total post-warmup draws=3200.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mean_Rt 1.79 0 0.07 1.66 1.74 1.79 1.83 1.92 3122 1
mean_Rt_adj 1.79 0 0.07 1.66 1.74 1.79 1.84 1.93 3123 1

再生算数を０〜１０人の一様分布にすると、収束しない。

> print(fit_u)
Inference for Stan model: fit_infection_u.
4 chains, each with iter=10000; warmup=5000; thin=5;
post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mean_Rt 2.15 0.05 0.10 1.98 2.06 2.16 2.23 2.30 3 1.99
mean_Rt_adj 2.13 0.07 0.11 1.93 2.06 2.13 2.21 2.36 3 2.27
lp__ -789.19 7.18 14.43 -820.19 -797.69 -787.44 -779.42 -762.53 4 1.89

traceplotやchainの分布はこんな感じ、
https://i.imgur.com/z0RL1KW.png
https://i.imgur.com/VRrrsKw.png