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

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

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

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

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

もし、筆者の伝えたいことが、「陽性判定は、即、感染者ということにはならない」　つまり、
「間違えることもある」という点にあるのであれば、あなたの言い分は通るかもしれないが、筆者の力点は、

>>正解は、約26％だ。判定が「陽性」でも、本当の感染者は10人中3人以下ということだ。

を見て判る通り、「陽性的中率が低い」というというところにある。
そのために、感度や特異度、真の感染率などに、具体的な数字を与えているし、問題設定の中では、わざわざ

>> （検査に至った経緯や、発熱・咳など他の所見は無視する）。
等と断り、終始定量的な説明が加えられている。

ならばこそ、なおさら、「クラスター」という言葉は使うべきでは無かった。
クラスター周辺での調査と、市中での無作為検査では、陽性的中率が変わってしまうのは、
全体を通して、筆者が言いたいことであっただろうに、にもかかわらず、問題の解説で
前提を崩してしまう「クラスター」という言葉を使ったのだから。

そもそも、ニュースでは、「PCRで陽性が○○人」等という使い方をしていただろうか？
単に「感染者○○人」だと思う。もちろん、この感染者の中には、偽陽性も含まれているだろうが、
検査陽性者数と感染者数の違いに注目を与えかねない「陽性者○○人」のような報道は記憶に無い。
そう考えると、「誤解の種」を自ら蒔いて、刈り取るかのような記事に見えてきた。

>>367
クラスターの有病率の事前確率を1%として計算してんじゃないの？
でなければ
>判定が「陽性」でも、本当の感染者は10人中3人以下ということだ。
ということにはならない。

クラスターで１０人が陽性として検査した人数と陽性的中率PPVとの関係をグラフにしてみた。
灰色実線は95%信頼区間境界、灰色点線はPPV＝0.26の線

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

全人口の有病率をクラスター内の有病率にすり替えて、１０人陽性でも感染しているのは３人以下という間違った結論を出している。

わかっていて書いているのか、馬鹿なのか、どちらかは不明。

>>368
罹患の定義による。
他人のゲノムでコンタミネーションが起こったりしていなければ、
あるゲノムのシークアンスが検出されたら罹患というなら、罹患率は１００％
ウイルスとして増殖能力を有しているかは不明、死骸の一部を検出しているだけかもしれない。

>>369
東京都の発表では　感染者数とせず、陽性患者数と表現している。

正確を期すなら陽性者数とすべきだろうな。

https://stopcovid19.metro.tokyo.lg.jp/

Natureのこの論文はエクセルとRのコードがついていて自分で再現できるので入力の手間が省ける。
https://www.nature.com/articles/s41591-020-0869-5
感染させる確率分布をガンマ分布を平行移動させた分布として想定して、データから最尤法でパラメータ算出している。

>>374
そのプログラムをみていくと
感染させる確率分布を
#--- infectiousness, gamma distribution ---
# gpar[1:2]: hyper-parameters (gamma)
# x : infection time of infectee w.r.t onset time of infector
f.Xc = function(x, gpar) { dgamma(x, gpar[1], gpar[2]) }
ガンマ分布にしているけど
一人が一人を感染させるモデルだから
形状パラメータgpar[1]は1で固定、つまり、指数分布だと思うんだがどうだろ？待ち時間の分布と同じ考え。

w.r.t = with reference to らしい

>>375
補足

ひとりが別の一人に一回感染させるというモデルなので待ち時間の分布の指数分布でいいと思う。
ガンマ分布の形状パラメータ＝１とおけば指数分布になる。
プログラムを書き直してグラフにすると
https://i.imgur.com/paErvfa.png
発症前に感染させている確率は47％と原著とあまりかわらないが、そのピークは発症1.6日前という結果になった。

>>369

>>正解は、約26％だ。判定が「陽性」でも、本当の感染者は10人中3人以下ということだ。

ああ、確かにそこは問題だな。前段の「クラスターで陽性が１０人」とリンクすると
思われても仕方がない。著者もうかつだとは思うが、有病率により的中率が変わると
いう話の内容は間違いではない。

>>371
元記事のグラフの有病率0.1のところをみれば的中率が80％越えになってるわけで、
単純に「クラスター」って言葉を軽んじてただけなんだろうね。
https://webronza.asahi.com/photo/photo.html?photo=/S2010/upload/2020032600006_3.jpg

>>375
たにんへの感染力の分布なんか数学的にわかるわけないやん？罹患して気道部にウィルスを放出できるくらいのウィルス量が繁殖するのは罹患していきなりなわけがない。

適当にSEIRモデル拡張してシミュレーションすれば分かるけどかなり自粛しても感染爆発は起こるよ

https://www.youtube.com/watch?v=4B5LTAdM0NE

東大数学科卒の高橋洋一氏が人との接触7割減」の根拠を数式から
語ってるがどう見る？

ガンマ分布：「一定期間に1回起きると期待されるランダムな事象が複数回起きるまでの時間の分布」
複数回でなくて1回だと指数分布だと思う。
何度も感染した症例を扱っているのではないのだから、ガンマ分布のパラメータを求める意味が理解できない。

>>380
SEIRモデルだと感染爆発させた方が早期に収束する。死者や感染者は増えるけど。
シミュレーションでは鎖国している前提で4割が感染すればオリンピックは可能だった。
SEIRモデルはEは感染力なし、Rは再感染しないというモデルだからなぁ。
モデルを修正してR->Sの再感染が0.１%あるだけで終焉しなかったな。
しかも、外部から感染力の強いキャリアー（保菌者）が入ってくることはモデルには組み込まれていない。

>>379
結局、一定期間に1回起きると期待されるランダムな事象として、その一定期間を推測しているんだろ？

>>381
SIRモデルって鎖国モデルだから結論は信用できん。

統計できる人尊敬するわ
MCMCとか訳分からん

>>386
統計ってある統計量がほんにゃら分布に従うというのを黙って受容しないと次に進めないよね。
郡内分散と郡間分散の比がF分布に従うとか言われても
どうしてかは理解していない。
stanのNUTSとかfrog leapとかわからんままにMCMCさせている。

>>386
sirモデルでは罹患した患者が回復するまでの機会中常に一定の確率(感染率)で遭遇した感受受宿主に感染を広めると仮定して立式してる。
ザックリした傾向を見るならそれで十分だけど、実際のモデルではそんな事はありえない。
感染する確率は患者の体内で繁殖しているウィルス量の増加に従って増えるからその効果を勘案してΓ分布(というかt^c d^x なる形の関数)に応じて感染率が罹患した時点から変化するとするんでしょ？
もちろんコレでも大体の傾向見るにはコレで十分というモデルを臨床例から適当に選んでるだけでしょ？
実際には未来永劫感染力を持ち続けるなんて事があるはずないし、罹患初期のウィルス量の増え方はおそらく指数関数的に増えるものを採用すべきだろうし。

mcmc

>>388
無理無理。
カブトガニVSシオマネキ論争はコンピュータの計算結果で決着がついたよ。

週末は休むとか週間変動の影響を除くために１週間の移動平均で線形回帰して片対数グラフにすると

https://i.imgur.com/7VwfswD.png

自粛の効果がでてきているな。　多分、検査自粛の効果だろうな。

>>393
これを微分で解いたら、ネ申か狂人だろな。

AからHの8人はそれぞれ正直者か嘘つきであり、誰が正直者か嘘つきかはお互いに知っている。
A,B,C,D,Eは嘘つきなら必ず嘘をつくが、F,G,Hは嘘つきでも正しいことを言う場合がある。
次の証言から確実に正直者と断定できるのは誰か？

A「嘘つきの方が正直者より多い」
B「Hは嘘つきである」
C「Bは嘘つきである」
D「CもFも嘘つきである」
E「8人の中に、少なくとも1人嘘つきがいる」
F「8人の中に、少なくとも2人嘘つきがいる」
G「Eは嘘つきである」
H「AもFも正直者である」

CTとPCRの一致係数（κ値）を信頼区間つきでMCMCしようかと思ってスクリプトを書いたはいいが
肝心なデータがない　：（
事前分布を一様分布にするのには異論があるかもしれん。

library(rjags)
kappa.model='
model{
A ~ dbeta(1,1) # A:Pr[CT+], 1-A:Pr[CT-]
B ~ dbeta(1,1) # B:Pr[PCR+|CT+]
C ~ dbeta(1,1) # C:Pr[PCR-|CT-]
p[1]=A*B # CT+PCR-
p[2]=A*(1-B) # CT+PCR-
p[3]=(1-A)*(1-C) # CT-PCR+
p[4]=(1-A)*C # CT-PCR-
y[1:4] ~ dmulti(p[],n) # multinominal distribution
po=(p[1]+p[4])/n # observed agreement
pe=(p[1]+p[2])/n*(p[1]+p[3])/n + (p[3]+p[4])/n*(p[2]+p[4])/n # coincidence
kappa=(po-pe)/(1-pe)
PABAK=2*po-1 # Prevalence Adjusted Bias Adjusted Kappa
}
'
writeLines(kappa.model,'kappaj.txt')

中国には特異度１００％の検査キットがあるんだってね。
すべて陰性にでるようにセットされていると。

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

>>397
こんなグラフになった。

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

>>397
これであってるかな？

> # 期待値
> integrate(function(x) x*pdf(x),0,100)$value
[1] 37.13
> # 50人以上の確率
> integrate(pdf,50,100)$value
[1] 0.0041903
> c(HPDI.lower=lwr,HPDI.upper=upr) # HPDI
HPDI.lower HPDI.upper
27.778 46.558

レベル低

慶応大の調査で、コロナ以外で来院した人をPCR検査したところ
4/67の確率で要請だった
東京都内1500万のうち何人くらいが感染しているか推定せよ

岩田健太郎・神戸大学教授『東京はすでに20万〜400万人感染の可能性』
https://leia.5ch.net/test/read.cgi/poverty/1587664106/

>>401
ニュー速でやったが、ベイズでこんな感じ？(感染者率)
https://i.imgur.com/Nfs1vCm.jpg

>>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