X
トップページ⇒数学
1002コメント563KB

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

■ このスレッドは過去ログ倉庫に格納されています
0001１３２人目の素数さん
垢版 |
2020/02/29(土) 02:18:41.53ID:twdO677Q
東大数学科卒の元官僚はこう分析してるが、お前らはどうなると思う？

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

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

こうした推計をすると、今の感染者は氷山の一角だと思われるが、今後の増加ペースはどうなるだろうか。
新型コロナウイルスの検査は簡単に行えるので、今後、日本での感染者数は増えていくだろう。
ある時点ではそれがネズミ算的に増えるかのように思える局面もあるだろうが、
筆者の推計が正しければ、現時点ではせいぜい数百人が一つのメドだ。
0358１３２人目の素数さん
垢版 |
2020/04/18(土) 17:17:45.60ID:PWtsyGc9
>>ニュースなどで「クラスターで陽性が10人」＝「10人が感染」と素直に受け取ってしまう向きも多いだろう。
>>１万人を検査すると、70人+198人が陽性と判定される。本当の陽性は70人だから、結果が正しい確率は70÷268＝26％
>>正解は、約26％だ。判定が「陽性」でも、本当の感染者は10人中3人以下ということだ。

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

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

しかし、クラスターが発見され、その周辺から、感染者を探す場合は、事前確率は1%ではない。
どの範囲まで、調査対象を広げるかにもよるが、50%とかが期待される。
仮に50%だとすると、87.5%、陽性と判断された結果は正しい。
0359１３２人目の素数さん
垢版 |
2020/04/18(土) 17:17:48.47ID:vWIoYYH+
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
計算はあってる。面倒だから解説は読まない。
0360１３２人目の素数さん
垢版 |
2020/04/18(土) 17:37:48.25ID:vWIoYYH+
クラスターから何人検査して何人陽性であったのかによって結果が違ってくるね。


検査の感度を最頻値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
0361１３２人目の素数さん
垢版 |
2020/04/18(土) 17:41:14.84ID:vWIoYYH+
クラスターから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
0362１３２人目の素数さん
垢版 |
2020/04/18(土) 17:55:40.37ID:vWIoYYH+
>358の指摘の通り、陽性的中率はその集団の有病率の影響されるから、クラスターの有病率を１％にするのは間違い。
情報がないから、クラスター内の有病率の確率分布を一様分布として計算するとか、母集団の有病率の10倍以内とか設定すればクラスター内の有病率の確率分布が出せる。

クラスタ100人検査で10人陽性、クラスター内の有病率は母集団の10倍以内と設定すると
陽性的中率は、0.684　95％CI（0.497-0.857)と計算された。
0364マイヤービートリス
垢版 |
2020/04/18(土) 18:18:36.32ID:Joe2nOPR
>>358
有難うございます。
>>しかし、クラスターが発見され、その周辺から、感染者を探す場合は、事前確率は1%ではない。
この見積もりがキーですね。
>>359,360,361,362
計算有難うございます。
0365１３２人目の素数さん
垢版 |
2020/04/18(土) 20:40:48.15ID:iquzLVJz
>>358
おいおい、なにも間違ったことは書いてないでしょ。
君の引用の仕方に問題があるんだよ。

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

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

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

>>364
だから、有病率が高いと見込まれる集団での検査であればいいわけで、
むやみに検査せずにある程度絞り込むべきだ（たとえば、クラスター
感染が疑われる集団を対象にせよ）とまで記事の中で述べられていれば
よかったんだろうけど、そこまで書かなかったのは片手落ちだったかもね。
0366１３２人目の素数さん
垢版 |
2020/04/18(土) 21:33:23.51ID:PWtsyGc9
もし、記事が、

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

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

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

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

換言すれば、問題の設定では、
「真の感染率＝1％とする。（検査に至った経緯や、発熱・咳など他の所見は無視する）。」
と書いているのに対し、問題の解説では、「クラスター」を持ち出して解説しているのは、明らかに不適当。
0367１３２人目の素数さん
垢版 |
2020/04/19(日) 00:22:33.17ID:eKu2VjGM
>>366
クラスターですら、「陽性判定＝感染」ってことではないんだから、その
一文で言いたいことに間違いはない。

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

有病率があがると的中率も大きくあがることについては、記事中のグラフで
示してあるから、情報が示されてないわけでもない。記事の主旨とは直接関係
ないので文章中では触れなかっただけでしょうね。
0368１３２人目の素数さん
垢版 |
2020/04/19(日) 01:26:51.78ID:R/sv/z3n
統計全然わからないんだけど、今行われてるPCR検査で陽性と判断されて、本当に罹患してる確率ってどの程度なの？
0369１３２人目の素数さん
垢版 |
2020/04/19(日) 01:56:19.99ID:gn6lsHTI
もし、筆者の伝えたいことが、「陽性判定は、即、感染者ということにはならない」　つまり、
「間違えることもある」という点にあるのであれば、あなたの言い分は通るかもしれないが、筆者の力点は、

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

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

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

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

そもそも、ニュースでは、「PCRで陽性が○○人」等という使い方をしていただろうか？
単に「感染者○○人」だと思う。もちろん、この感染者の中には、偽陽性も含まれているだろうが、
検査陽性者数と感染者数の違いに注目を与えかねない「陽性者○○人」のような報道は記憶に無い。
そう考えると、「誤解の種」を自ら蒔いて、刈り取るかのような記事に見えてきた。
0370１３２人目の素数さん
垢版 |
2020/04/19(日) 02:31:27.50ID:Czp86qrf
>>367
クラスターの有病率の事前確率を1%として計算してんじゃないの？
でなければ
>判定が「陽性」でも、本当の感染者は10人中3人以下ということだ。
ということにはならない。
0371１３２人目の素数さん
垢版 |
2020/04/19(日) 07:50:33.21ID:Czp86qrf
クラスターで１０人が陽性として検査した人数と陽性的中率PPVとの関係をグラフにしてみた。
灰色実線は95%信頼区間境界、灰色点線はPPV＝0.26の線

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

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

わかっていて書いているのか、馬鹿なのか、どちらかは不明。
0372１３２人目の素数さん
垢版 |
2020/04/19(日) 07:56:38.54ID:Czp86qrf
>>368
罹患の定義による。
他人のゲノムでコンタミネーションが起こったりしていなければ、
あるゲノムのシークアンスが検出されたら罹患というなら、罹患率は１００％
ウイルスとして増殖能力を有しているかは不明、死骸の一部を検出しているだけかもしれない。
0374１３２人目の素数さん
垢版 |
2020/04/19(日) 20:39:24.23ID:Czp86qrf
Natureのこの論文はエクセルとRのコードがついていて自分で再現できるので入力の手間が省ける。
https://www.nature.com/articles/s41591-020-0869-5
感染させる確率分布をガンマ分布を平行移動させた分布として想定して、データから最尤法でパラメータ算出している。
0375１３２人目の素数さん
垢版 |
2020/04/20(月) 17:39:51.08ID:Db3kUO+J
>>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 らしい
0376１３２人目の素数さん
垢版 |
2020/04/20(月) 19:39:48.89ID:Db3kUO+J
>>375
補足

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

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

ああ、確かにそこは問題だな。前段の「クラスターで陽性が１０人」とリンクすると
思われても仕方がない。著者もうかつだとは思うが、有病率により的中率が変わると
いう話の内容は間違いではない。
0379１３２人目の素数さん
垢版 |
2020/04/20(月) 23:44:00.74ID:IwCIr2qd
>>375
たにんへの感染力の分布なんか数学的にわかるわけないやん？罹患して気道部にウィルスを放出できるくらいのウィルス量が繁殖するのは罹患していきなりなわけがない。
0380１３２人目の素数さん
垢版 |
2020/04/21(火) 00:05:49.98ID:bowf2rRh
適当にSEIRモデル拡張してシミュレーションすれば分かるけどかなり自粛しても感染爆発は起こるよ
0381１３２人目の素数さん
垢版 |
2020/04/21(火) 02:51:28.64ID:7mnZKVUh
https://www.youtube.com/watch?v=4B5LTAdM0NE

東大数学科卒の高橋洋一氏が人との接触7割減」の根拠を数式から
語ってるがどう見る？
0382１３２人目の素数さん
垢版 |
2020/04/21(火) 06:50:58.83ID:J5/+FVIQ
ガンマ分布：「一定期間に1回起きると期待されるランダムな事象が複数回起きるまでの時間の分布」
複数回でなくて1回だと指数分布だと思う。
何度も感染した症例を扱っているのではないのだから、ガンマ分布のパラメータを求める意味が理解できない。
0383１３２人目の素数さん
垢版 |
2020/04/21(火) 06:56:00.72ID:J5/+FVIQ
>>380
SEIRモデルだと感染爆発させた方が早期に収束する。死者や感染者は増えるけど。
シミュレーションでは鎖国している前提で4割が感染すればオリンピックは可能だった。
SEIRモデルはEは感染力なし、Rは再感染しないというモデルだからなぁ。
モデルを修正してR->Sの再感染が0.１%あるだけで終焉しなかったな。
しかも、外部から感染力の強いキャリアー（保菌者）が入ってくることはモデルには組み込まれていない。
0384１３２人目の素数さん
垢版 |
2020/04/21(火) 06:59:16.79ID:J5/+FVIQ
>>379
結局、一定期間に1回起きると期待されるランダムな事象として、その一定期間を推測しているんだろ？
0387１３２人目の素数さん
垢版 |
2020/04/21(火) 11:28:48.72ID:J5/+FVIQ
>>386
統計ってある統計量がほんにゃら分布に従うというのを黙って受容しないと次に進めないよね。
郡内分散と郡間分散の比がF分布に従うとか言われても
どうしてかは理解していない。
stanのNUTSとかfrog leapとかわからんままにMCMCさせている。
0388イナ ◆/7jUdUKiSM
垢版 |
2020/04/21(火) 13:30:52.20ID:0ilKIHza
こういうときこそコンピューターとかに頼らずに一つ一つ数をかぞえて方程式を立て、微分するべきだと思う。
‖∩∩‖ □ ‖　　 ＼
（-_-)）　　‖______‖
（っγυ　 ｡‖╂─╂‖
■`(_)_)ц~ ‖╂─╂‖
＼■υυ■_∩∩､＼＼‖
＼＼＼＼⊂(_ _ )`⌒つ)
＼＼＼＼＼＼＼`υ､＼/|
＼＼＼＼＼`.,､､､＼`/ |
__＼＼＼＼彡`-`ﾐっ/　L
￣|＼_＼＼_U,~⌒ヾ　/
］| ‖￣￣￣￣U~~U / /
__| ‖ □　□ ‖ |/ /
___`‖________‖／_/
￣￣￣￣￣￣￣￣￣‖ /
__________________‖/
0389１３２人目の素数さん
垢版 |
2020/04/21(火) 14:53:08.35ID:pHehMVs5
>>386
sirモデルでは罹患した患者が回復するまでの機会中常に一定の確率(感染率)で遭遇した感受受宿主に感染を広めると仮定して立式してる。
ザックリした傾向を見るならそれで十分だけど、実際のモデルではそんな事はありえない。
感染する確率は患者の体内で繁殖しているウィルス量の増加に従って増えるからその効果を勘案してΓ分布(というかt^c d^x なる形の関数)に応じて感染率が罹患した時点から変化するとするんでしょ？
もちろんコレでも大体の傾向見るにはコレで十分というモデルを臨床例から適当に選んでるだけでしょ？
実際には未来永劫感染力を持ち続けるなんて事があるはずないし、罹患初期のウィルス量の増え方はおそらく指数関数的に増えるものを採用すべきだろうし。
0391１３２人目の素数さん
垢版 |
2020/04/22(水) 08:56:14.63ID:0tlpvLKp
>>388
無理無理。
カブトガニVSシオマネキ論争はコンピュータの計算結果で決着がついたよ。
0392１３２人目の素数さん
垢版 |
2020/04/22(水) 13:23:53.12ID:0tlpvLKp
週末は休むとか週間変動の影響を除くために１週間の移動平均で線形回帰して片対数グラフにすると

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

自粛の効果がでてきているな。　多分、検査自粛の効果だろうな。
0393イナ ◆/7jUdUKiSM
垢版 |
2020/04/22(水) 19:51:00.75ID:iq1GZOqA
‖∩∩ ‖ □ ‖前>>388
（(-_-)‖　　‖______
（っ⌒⌒ﾞ　 ｡‖╂─╂
■`(_)_)ц~　‖╂─╂
＼■υυ■_∩∩､＼＼＼
＼＼＼＼⊂(_ _ )`⌒つ､
＼＼＼＼＼＼＼`υ､＼＼＼＼＼＼＼＼＼＼＼＼＼＼＼＼`>>391カブトガニでもシオマネキでもシュクメルリでも微分するのがいちばん強力だと思う。
0394１３２人目の素数さん
垢版 |
2020/04/23(木) 04:58:29.35ID:3vcirHk0
>>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も正直者である」
0395１３２人目の素数さん
垢版 |
2020/04/23(木) 08:28:17.25ID:3vcirHk0
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')
0396１３２人目の素数さん
垢版 |
2020/04/23(木) 10:37:41.66ID:3vcirHk0
中国には特異度１００％の検査キットがあるんだってね。
すべて陰性にでるようにセットされていると。
0397１３２人目の素数さん
垢版 |
2020/04/23(木) 13:51:17.19ID:3vcirHk0
新型コロナ患者を治療している病院に100人の職員がいる。
検体採取器具は5人分、試薬は1回分しかないとする。
無作為抽出した5人の職員から採取した検体を混合して検査したら陽性であった。
職員の陽性者数の期待値を求めよ。
また、50人以上の感染者いる確率はいくつか？
検査の陽性率はハイリスク群に検査している東京の数値2457/6654を使って計算せよ。
https://i.imgur.com/zYK75Lo.jpg
0399１３２人目の素数さん
垢版 |
2020/04/23(木) 17:36:48.57ID:3vcirHk0
>>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
0401１３２人目の素数さん
垢版 |
2020/04/24(金) 02:55:58.61ID:juJsFFfP
慶応大の調査で、コロナ以外で来院した人をPCR検査したところ
4/67の確率で要請だった
東京都内1500万のうち何人くらいが感染しているか推定せよ
0405１３２人目の素数さん
垢版 |
2020/04/24(金) 08:13:06.02ID:0onl6lJy
>>401
これって東京の無作為サンプリングではなさそうだよね
病院に来たって事は熱が出てたのかもしれないし
想定できる母集団ってなんになるのだろうか
0406１３２人目の素数さん
垢版 |
2020/04/24(金) 12:07:56.92ID:v55OWzbu
新型コロナ患者を治療している病院に100人の職員がいる。
検体採取器具は10人分、試薬は1回分しかないとする。
無作為抽出した10人の職員から採取した検体を混合して検査したら陰性であった。
職員の陽性者数の期待値を求めよ。
また、50人以上の感染者いる確率はいくつか？
検査の陽性率はハイリスク群に検査している東京の数値2457/6654を使って計算せよ。
0407１３２人目の素数さん
垢版 |
2020/04/24(金) 12:11:36.28ID:v55OWzbu
>>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
0408１３２人目の素数さん
垢版 |
2020/04/24(金) 17:39:04.40ID:8oiI190P
>>401
まじか。
まあ、体調悪いから病院に行くわけで、バイアスかかってるとはいえ...。
0411１３２人目の素数さん
垢版 |
2020/04/24(金) 20:55:16.38ID:v55OWzbu
>>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
0412１３２人目の素数さん
垢版 |
2020/04/24(金) 21:13:10.41ID:v55OWzbu
>>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
0413１３２人目の素数さん
垢版 |
2020/04/24(金) 21:35:58.21ID:v55OWzbu
>>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
0415１３２人目の素数さん
垢版 |
2020/04/24(金) 22:01:25.57ID:v55OWzbu
>>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
0418１３２人目の素数さん
垢版 |
2020/04/25(土) 13:54:43.50ID:R/UD6QQG
某医療センターでは治療後、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
0420１３２人目の素数さん
垢版 |
2020/04/26(日) 21:43:36.22ID:0lgRXnyr
スレチですが統計に詳しい方がいると思って伺いました。
工学系の大学院生で論文を読んでて、見慣れない記号や数式が出てきたので質問させてください。
この数式中のEは何を意味するのでしょうか？
Rの期待値的なものですか？

https://i.imgur.com/oVwva8r.jpg
0421１３２人目の素数さん
垢版 |
2020/04/27(月) 02:54:57.87ID:cS0ISst+
院生でわからないってマジか
専門外の機械学習いきなりやらされたのかな？
0422１３２人目の素数さん
垢版 |
2020/04/27(月) 06:45:22.65ID:S7AmHM83
有病率が0.3％（オーストリアの無作為抽出)や6%(慶応大学の入院患者無作為抽出）のときは

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

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

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

有病率が10％くらいになればイカカマキットの方が陰性的中率の成績が劣る。
0423１３２人目の素数さん
垢版 |
2020/04/27(月) 07:30:42.68ID:S7AmHM83
某医療センターでは治療後、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
0424１３２人目の素数さん
垢版 |
2020/04/28(火) 10:03:02.25ID:/p2virDs
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
0426１３２人目の素数さん
垢版 |
2020/04/29(水) 16:59:02.34ID:YhsQxAnp
>>425
正解。正直者が誰もいないようにプログラムした。
イナ師が微分で解いてくれるのを期待していたんだが。
0427１３２人目の素数さん
垢版 |
2020/04/29(水) 17:43:56.54ID:YhsQxAnp
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
0429１３２人目の素数さん
垢版 |
2020/04/29(水) 18:35:59.17ID:uFfpYtab
rRT-PRC検査の感度・特異度ともに90%以上です。
0430１３２人目の素数さん
垢版 |
2020/04/29(水) 19:40:54.19ID:uFfpYtab
通常の医療でも臨床診断の5%は誤診だと推定されていたりする。
これを実際の数として見積もると膨大になる。
0435１３２人目の素数さん
垢版 |
2020/04/30(木) 12:13:53.58ID:qOF+URFa
>>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
0436１３２人目の素数さん
垢版 |
2020/04/30(木) 19:13:48.86ID:42TM9o6B
>>434
これの興味深いところは、陽性率の高さもさることながら、医療関係者の感染率が一般市民のそれよりも
（多分有意に）高いことだと思う。抗体検査の精度ってほとんどデータが無いに等しいけど、医療関係者の
感染者率が平均よりも多いのならば（実際そうなのだと思う）、このテストはそれを反映したものと言えて、
抗体検査の信頼性をある程度証明できるかもしれない。
そして、市中の陽性率の高さも。

この辺、ベイズで上手く検証できないかな。
0437１３２人目の素数さん
垢版 |
2020/04/30(木) 20:32:35.24ID:qOF+URFa
>>436
>医療関係者の感染率が一般市民のそれよりも（多分有意に）高い

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

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

>435にRの出力を貼っておいた。
0438１３２人目の素数さん
垢版 |
2020/04/30(木) 20:35:20.83ID:rbz2947p
>>436
55人中5人と147人中7人じゃ、有意差ないだろ。

と書こうと思ったら、 >>437氏に先を越されたｗｗｗ
0439１３２人目の素数さん
垢版 |
2020/04/30(木) 20:41:39.46ID:rbz2947p
抗体検査キットの評価をしたら特異度はみな5/5だったらしいけど、サンプルが
５つだけって意味だとすると、せいぜい90％以上くらいのことしか言えないな。
特異度95%だとしたら、５％が陽性っていわれてもねぇｗ
0440１３２人目の素数さん
垢版 |
2020/04/30(木) 21:20:01.97ID:qOF+URFa
>>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%の一様分布というのは現実離れした分布だとは思う。
0442１３２人目の素数さん
垢版 |
2020/04/30(木) 21:33:35.19ID:qOF+URFa
>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
0443１３２人目の素数さん
垢版 |
2020/04/30(木) 21:47:02.48ID:qOF+URFa
>>440（追記)
>有病率が平均値50%の一様分布というのは現実離れした分布だとは思う。
オーストリアのデータ0.3%を参考に有病率の事前分布が0.2-0.4％
β(10.865, 3279.34)に相当として、MCMCしてみた。
感度と特異度の事前分布は前回と同じ。

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

さほど、有病率が高いという結論は引き出せないな。
0444１３２人目の素数さん
垢版 |
2020/04/30(木) 22:23:36.86ID:gdjT5b3G
自己流な検定

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
0445１３２人目の素数さん
垢版 |
2020/04/30(木) 23:17:24.52ID:qOF+URFa
β分布と乱数発生で比較。

一般人と医療従事者の確率分布に従う乱数を１０００万個発生させて比率の分布をグラフにすると
９５％信頼区間が１を挟むから有意差なし。比率が１以上となる確率も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)
0446１３２人目の素数さん
垢版 |
2020/05/01(金) 00:36:14.80ID:LcUyF6Tc
>>441
すべてが陰性に出るなら陽性は０人だろｗ
0448１３２人目の素数さん
垢版 |
2020/05/01(金) 10:05:01.58ID:LcUyF6Tc
>>447
だから、陽性が何人か出てんだから感度0%のインチキじゃなかろうってこと。
0449１３２人目の素数さん
垢版 |
2020/05/04(月) 00:17:48.34ID:Q8iLGwO2
最近の報道

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

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

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

問題

1試薬は1回しか検査できないとして
これがイカサマキットであることを証明する手段はあるか？
0450１３２人目の素数さん
垢版 |
2020/05/04(月) 15:22:29.91ID:jDRWX2Ph
3月の宿題で(1)のみ正解の数弱@shukudai_sujaku

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

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

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

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

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

https://twitter.com/shukudai_sujaku
https://twitter.com/5chan_nel (5ch newer account)
0451１３２人目の素数さん
垢版 |
2020/05/04(月) 22:45:28.89ID:94tg6i4k
高橋洋一(統計数理研究所→大蔵省)
スウェーデンの感染症対策を解説
https://youtu.be/2bsYRuorwMI?t=192
0452１３２人目の素数さん
垢版 |
2020/05/05(火) 06:07:12.42ID:ht6rG86e
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％のモデルの方がいいかも、という結果。
陽性数/検査数の時系列データでもあればもう少し差がでるかもしれん。
東京都のデータで計算させようかと思ったが、東京都は検査人数を隠蔽しているので使いものにならない。
0453１３２人目の素数さん
垢版 |
2020/05/05(火) 06:26:36.20ID:FpoKo6a+
訂正

陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする
↓
陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする
0454１３２人目の素数さん
垢版 |
2020/05/05(火) 11:13:29.85ID:b2IqdVzK
3月の宿題で(1)のみ正解の数弱@shukudai_sujaku

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

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

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

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

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

https://twitter.com/shukudai_sujaku
https://twitter.com/5chan_nel (5ch newer account)
0455１３２人目の素数さん
垢版 |
2020/05/05(火) 23:00:17.52ID:wvtmzVA1
>>449
抗体のはいったサンプルを２，３０個試してみりゃわかんじゃないの？
いくらなんでも感度が３％じゃつかえねぇ。
0456１３２人目の素数さん
垢版 |
2020/05/06(水) 15:29:26.96ID:RG00+xls
イカサマキットの感度特異度の事前分布を一様分布に設定して
抗体のはいったサンプルを２０個で全部陰性であったので３０個試したら全部陰性であったとすると
感度・特異度の事後分布は

https://i.imgur.com/6ZYn34k.png
■ このスレッドは過去ログ倉庫に格納されています
ニューススポーツなんでも実況