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

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

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

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

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

>>346
ロスマンの疫学を読んだことのある人なら答えられるし、
感度特異度の意味を理解していれば自分で答えが出せる。

サパイラ 身体診察のアートとサイエンスにも載っているが、こっちは100%の所見ではないな。

>>347
君が書いたらいい

サパイラ 身体診察のアートとサイエンスには
あらゆる疾患に感度の高い検査（所見）として　10　finger test が挙げられている。
すなわち、指が１０本あれば疾患ありと判断する。

この検査の感度は１００％に近い。新型コロナでも高感度。

身体診察のアートとサイエンス　の　翻訳者のひとりには岩田健太郎がいる。
俺は旧版を原著で読んだ。　

引用すると、

A test with high sensitivity is not necessarily a useful test. The sign “10 fingers” would be
extremely sensitive for almost any disease because most patients with the disease will have ten
fingers. Very few patients with the disease will have a different number of fingers. Thus, the ratio of
true positives (number of patients with the disease who have ten fingers) to the sum of true positives
plus false negatives (where false negatives are people who have the disease and do not have ten
fingers) will usually be greater than 0.99 (except in a sanitarium for Hansen disease). Yet, common
sense tells us that the possession of ten fingers, however sensitive on paper, is not of great use to the
diagnostician. Why not? The reason is that most of the people in the world have ten fingers but do not
have the disease.

>>334
各国の検査数と陽性数のデータがあれば、
β分布のパラメータを以下のように設定してMCMCすれば出せる。

beta(13.6991,9.4661)でmode 0.6, sd = 0.1
beta(36.172,4.908)でmode 0.9 sd = 0.05

東京都で検査を受けたハイリスクグループでの陽性率を出そうしたんだが、
東京都は行政検査数しか公表しないので出す術がない。

>>353
＞コロナが苦手なもの

知性だな。
台湾の知性が制御している。

>>343
オープンなリソース管理はコロナの弱点だろうな
公務員が数ちょろまかしてるような国では強くなれる

ベイズを議論している人に教えてほしいのだが
https://webronza.asahi.com/science/articles/2020032600006.html?iref=wrp_rnavi_new
でこの人の計算や解釈は正しいの？

>>ニュースなどで「クラスターで陽性が10人」＝「10人が感染」と素直に受け取ってしまう向きも多いだろう。
>>１万人を検査すると、70人+198人が陽性と判定される。本当の陽性は70人だから、結果が正しい確率は70÷268＝26％
>>正解は、約26％だ。判定が「陽性」でも、本当の感染者は10人中3人以下ということだ。

明らかに誤っている部分がある。

市中感染率を1%と仮定し、その中から無作為に調査を行って、10人の陽性と判断される者を見つけた場合、
確かに、本当の感染者は2.6人程度。

しかし、クラスターが発見され、その周辺から、感染者を探す場合は、事前確率は1%ではない。
どの範囲まで、調査対象を広げるかにもよるが、50%とかが期待される。
仮に50%だとすると、87.5%、陽性と判断された結果は正しい。

pr2pv <- function( # prevalence to predicative value
pr ,# prevalence
sn=0.7, # sensitibity=TP/(TP+FN)
sp=0.9) # specificity=TN/(TN+FP)
{
N=1 # polutaion million, billion,or any proper unit
si=pr*N # sick population
he=(1-pr)*N # healthy population
TP=si*sn
FN=si*(1-sn)
TN=he*sp
FP=he*(1-sp)
PPV=TP/(TP+FP)
NPV=TN/(TN+FN)
PV=c(PPV=PPV,NPV=NPV)
return(PV)
}



> pr2pv(0.01,0.7,0.98)
PPV NPV
0.2611940 0.9969174
計算はあってる。面倒だから解説は読まない。

クラスターから何人検査して何人陽性であったのかによって結果が違ってくるね。


検査の感度を最頻値0.7標準偏差0.1
特異度を最頻値0.98 標準偏差0.01
有病率は(0,1)の一様分布
を弱情報事前分布(情弱事前分布w)として

クラスターから10人検査したら10人陽性であったとき　stanでMCMCした結果

https://i.imgur.com/AaOB2LX.png
https://i.imgur.com/yeIgnCX.png

クラスターから100人検査して10人陽性だった場合


mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
prev 0.13475 0.00044 0.05797 0.04287 0.12754 0.27095 17334 1.00010
sen 0.66409 0.00074 0.10423 0.44686 0.67054 0.84855 20036 1.00010
spc 0.97580 0.00007 0.01006 0.95256 0.97714 0.99146 20272 1.00012
p 0.10797 0.00019 0.03044 0.05633 0.10530 0.17455 24621 1.00001
PPV 0.76071 0.00072 0.09797 0.53774 0.77243 0.91536 18317 1.00011
NPV 0.95986 0.00013 0.01749 0.91811 0.96260 0.98585 17185 1.00010
precision 0.94212 0.00013 0.01766 0.90229 0.94405 0.97089 19310 1.00010
pLR 33.28155 0.14736 18.25227 12.71591 28.84717 80.12622 15341 1.00013
nLR 0.34429 0.00075 0.10692 0.15527 0.33782 0.56712 20099 1.00010
DOR 114.74459 0.78412 97.10981 25.95828 89.02272 356.05482 15338 0.99996
lp__ -75.90415 0.01168 1.27866 -79.22234 -75.58340 -74.42549 11979 1.00024

PPVは0.76

>358の指摘の通り、陽性的中率はその集団の有病率の影響されるから、クラスターの有病率を１％にするのは間違い。
情報がないから、クラスター内の有病率の確率分布を一様分布として計算するとか、母集団の有病率の10倍以内とか設定すればクラスター内の有病率の確率分布が出せる。

クラスタ100人検査で10人陽性、クラスター内の有病率は母集団の10倍以内と設定すると
陽性的中率は、0.684　95％CI（0.497-0.857)と計算された。

結局、誤った解説を読む羽目になったｗ

>>358
有難うございます。
>>しかし、クラスターが発見され、その周辺から、感染者を探す場合は、事前確率は1%ではない。
この見積もりがキーですね。
>>359,360,361,362
計算有難うございます。

>>358
おいおい、なにも間違ったことは書いてないでしょ。
君の引用の仕方に問題があるんだよ。

>>ニュースなどで「クラスターで陽性が10人」＝「10人が感染」と素直に受け取ってしまう向きも多いだろう。

で一旦段落は終わってるんだから、１万人を検査すると云々という、設定された
「問題」の答えとは直接リンクしてない。君のように勘違いする人はいるかも
しれないが、きちんと読めば間違ったことは書いてないことがわかるはず。
どこにもクラスターの有病率が１％だなどとは書いてない。

で、仮に有病率が50%でも、
>「クラスターで陽性が10人」＝「10人が感染」
ではないわけで、引用元の著者の主張はまっとうなものだと思える。

>>364
だから、有病率が高いと見込まれる集団での検査であればいいわけで、
むやみに検査せずにある程度絞り込むべきだ（たとえば、クラスター
感染が疑われる集団を対象にせよ）とまで記事の中で述べられていれば
よかったんだろうけど、そこまで書かなかったのは片手落ちだったかもね。

もし、記事が、

「ニュースなどで「陽性判定が10人」＝「10人が感染」と素直に受け取ってしまう向きも多いだろう。」

と書かれていたなら、あなたの指摘は正しいだろう。しかし、実際は

>> ニュースなどで「クラスターで陽性が10人」＝「10人が感染」と素直に受け取ってしまう向きも多いだろう。

となっている。
クラスターとは、集団の存在、つまり、濃度の濃い部分があることを示唆したものであり、
母集団の感染率をそのまま用いるべきではない点を指摘した。
だから、非クラスターなら、正解だか、クラスターなら正解とは言えないという主旨で書いている。

換言すれば、問題の設定では、
「真の感染率＝1％とする。（検査に至った経緯や、発熱・咳など他の所見は無視する）。」
と書いているのに対し、問題の解説では、「クラスター」を持ち出して解説しているのは、明らかに不適当。

>>366
クラスターですら、「陽性判定＝感染」ってことではないんだから、その
一文で言いたいことに間違いはない。

でもって、著者は、クラスターの有病率が１％なんてことは一言も言って
ないし、そもそも、クラスターについての具体的な解説もしていない。単に
ニュースで扱われる文言の例としてそこで１回だけクラスターという言葉が
使われてるだけ。クラスターの検査を否定してる内容でもない。

有病率があがると的中率も大きくあがることについては、記事中のグラフで
示してあるから、情報が示されてないわけでもない。記事の主旨とは直接関係
ないので文章中では触れなかっただけでしょうね。

統計全然わからないんだけど、今行われてるPCR検査で陽性と判断されて、本当に罹患してる確率ってどの程度なの？

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

>>正解は、約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
すべてが陰性に出るなら陽性は０人だろｗ