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

1 名前：１３２人目の素数さん [2020/02/29(土) 02:18:41 ID:twdO677Q.net] 東大数学科卒の元官僚はこう分析してるが、お前らはどうなると思う？ www.zakzak.co.jp/soc/news/200220/dom2002200003-n2.html 中国国外感染者の中国国内との比率をみると、 １月２０日の数字公表以降は、０・８〜２・６％で比較的安定している。 これは、新型肺炎の感染者のほとんどは中国国内、それも湖北省に集中しているからだ。 ちなみに中国国外での感染者数は、中国国内の１・１％だ（２月１６日現在）。 本コラムで紹介したが、現時点では、最終的な中国国内の感染者数は２０万人超と筆者は推計している。 となると、中国国外の感染者は数千人程度になるだろう。 中国国外のうち日本の比率は１割弱なので、日本の感染者数は数百人程度であろう。 その場合、死者も数人から１０人程度になるだろう。 こうした推計をすると、今の感染者は氷山の一角だと思われるが、今後の増加ペースはどうなるだろうか。 新型コロナウイルスの検査は簡単に行えるので、今後、日本での感染者数は増えていくだろう。 ある時点ではそれがネズミ算的に増えるかのように思える局面もあるだろうが、 筆者の推計が正しければ、現時点ではせいぜい数百人が一つのメドだ。 377 名前：１３２人目の素数さん mailto:sage [2020/04/18(土) 17:17:45 ID:PWtsyGc9.net] >>ニュースなどで「クラスターで陽性が10人」＝「10人が感染」と素直に受け取ってしまう向きも多いだろう。 >>１万人を検査すると、70人+198人が陽性と判定される。本当の陽性は70人だから、結果が正しい確率は70÷268＝26％ >>正解は、約26％だ。判定が「陽性」でも、本当の感染者は10人中3人以下ということだ。 明らかに誤っている部分がある。 市中感染率を1%と仮定し、その中から無作為に調査を行って、10人の陽性と判断される者を見つけた場合、 確かに、本当の感染者は2.6人程度。 しかし、クラスターが発見され、その周辺から、感染者を探す場合は、事前確率は1%ではない。 どの範囲まで、調査対象を広げるかにもよるが、50%とかが期待される。 仮に50%だとすると、87.5%、陽性と判断された結果は正しい。 378 名前：１３２人目の素数さん mailto:sage [2020/04/18(土) 17:17:48 ID:vWIoYYH+.net] 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 計算はあってる。面倒だから解説は読まない。 379 名前：１３２人目の素数さん mailto:sage [2020/04/18(土) 17:37:48 ID:vWIoYYH+.net] クラスターから何人検査して何人陽性であったのかによって結果が違ってくるね。 検査の感度を最頻値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 380 名前：１３２人目の素数さん mailto:sage [2020/04/18(土) 17:41:14.84 ID:vWIoYYH+.net] クラスターから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 381 名前：１３２人目の素数さん mailto:sage [2020/04/18(土) 17:55:40.37 ID:vWIoYYH+.net] >358の指摘の通り、陽性的中率はその集団の有病率の影響されるから、クラスターの有病率を１％にするのは間違い。 情報がないから、クラスター内の有病率の確率分布を一様分布として計算するとか、母集団の有病率の10倍以内とか設定すればクラスター内の有病率の確率分布が出せる。 クラスタ100人検査で10人陽性、クラスター内の有病率は母集団の10倍以内と設定すると 陽性的中率は、0.684　95％CI（0.497-0.857)と計算された。 382 名前：１３２人目の素数さん mailto:sage [2020/04/18(土) 18:10:34.06 ID:vWIoYYH+.net] 結局、誤った解説を読む羽目になったｗ 383 名前：マイヤービートリス [2020/04/18(土) 18:18:36 ID:Joe2nOPR.net] >>358 有難うございます。 >>しかし、クラスターが発見され、その周辺から、感染者を探す場合は、事前確率は1%ではない。 この見積もりがキーですね。 >>359,360,361,362 計算有難うございます。 384 名前：１３２人目の素数さん [2020/04/18(土) 20:40:48.15 ID:iquzLVJz.net] >>358 おいおい、なにも間違ったことは書いてないでしょ。 君の引用の仕方に問題があるんだよ。 >>ニュースなどで「クラスターで陽性が10人」＝「10人が感染」と素直に受け取ってしまう向きも多いだろう。 で一旦段落は終わってるんだから、１万人を検査すると云々という、設定された 「問題」の答えとは直接リンクしてない。君のように勘違いする人はいるかも しれないが、きちんと読めば間違ったことは書いてないことがわかるはず。 どこにもクラスターの有病率が１％だなどとは書いてない。 で、仮に有病率が50%でも、 >「クラスターで陽性が10人」＝「10人が感染」 ではないわけで、引用元の著者の主張はまっとうなものだと思える。 >>364 だから、有病率が高いと見込まれる集団での検査であればいいわけで、 むやみに検査せずにある程度絞り込むべきだ（たとえば、クラスター 感染が疑われる集団を対象にせよ）とまで記事の中で述べられていれば よかったんだろうけど、そこまで書かなかったのは片手落ちだったかもね。 385 名前：１３２人目の素数さん mailto:sage [2020/04/18(土) 21:33:23 ID:PWtsyGc9.net] もし、記事が、 「ニュースなどで「陽性判定が10人」＝「10人が感染」と素直に受け取ってしまう向きも多いだろう。」 と書かれていたなら、あなたの指摘は正しいだろう。しかし、実際は >> ニュースなどで「クラスターで陽性が10人」＝「10人が感染」と素直に受け取ってしまう向きも多いだろう。 となっている。 クラスターとは、集団の存在、つまり、濃度の濃い部分があることを示唆したものであり、 母集団の感染率をそのまま用いるべきではない点を指摘した。 だから、非クラスターなら、正解だか、クラスターなら正解とは言えないという主旨で書いている。 換言すれば、問題の設定では、 「真の感染率＝1％とする。（検査に至った経緯や、発熱・咳など他の所見は無視する）。」 と書いているのに対し、問題の解説では、「クラスター」を持ち出して解説しているのは、明らかに不適当。 386 名前：１３２人目の素数さん [2020/04/19(日) 00:22:33.17 ID:eKu2VjGM.net] >>366 クラスターですら、「陽性判定＝感染」ってことではないんだから、その 一文で言いたいことに間違いはない。 でもって、著者は、クラスターの有病率が１％なんてことは一言も言って ないし、そもそも、クラスターについての具体的な解説もしていない。単に ニュースで扱われる文言の例としてそこで１回だけクラスターという言葉が 使われてるだけ。クラスターの検査を否定してる内容でもない。 有病率があがると的中率も大きくあがることについては、記事中のグラフで 示してあるから、情報が示されてないわけでもない。記事の主旨とは直接関係 ないので文章中では触れなかっただけでしょうね。 387 名前：１３２人目の素数さん mailto:sage [2020/04/19(日) 01:26:51.78 ID:R/sv/z3n.net] 統計全然わからないんだけど、今行われてるPCR検査で陽性と判断されて、本当に罹患してる確率ってどの程度なの？ 388 名前：１３２人目の素数さん mailto:sage [2020/04/19(日) 01:56:19.99 ID:gn6lsHTI.net] もし、筆者の伝えたいことが、「陽性判定は、即、感染者ということにはならない」　つまり、 「間違えることもある」という点にあるのであれば、あなたの言い分は通るかもしれないが、筆者の力点は、 >>正解は、約26％だ。判定が「陽性」でも、本当の感染者は10人中3人以下ということだ。 を見て判る通り、「陽性的中率が低い」というというところにある。 そのために、感度や特異度、真の感染率などに、具体的な数字を与えているし、問題設定の中では、わざわざ >> （検査に至った経緯や、発熱・咳など他の所見は無視する）。 等と断り、終始定量的な説明が加えられている。 ならばこそ、なおさら、「クラスター」という言葉は使うべきでは無かった。 クラスター周辺での調査と、市中での無作為検査では、陽性的中率が変わってしまうのは、 全体を通して、筆者が言いたいことであっただろうに、にもかかわらず、問題の解説で 前提を崩してしまう「クラスター」という言葉を使ったのだから。 そもそも、ニュースでは、「PCRで陽性が○○人」等という使い方をしていただろうか？ 単に「感染者○○人」だと思う。もちろん、この感染者の中には、偽陽性も含まれているだろうが、 検査陽性者数と感染者数の違いに注目を与えかねない「陽性者○○人」のような報道は記憶に無い。 そう考えると、「誤解の種」を自ら蒔いて、刈り取るかのような記事に見えてきた。 389 名前：１３２人目の素数さん mailto:sage [2020/04/19(日) 02:31 ] [ここ壊れてます] 390 名前：:27.50 ID:Czp86qrf.net mailto: >>367 クラスターの有病率の事前確率を1%として計算してんじゃないの？ でなければ >判定が「陽性」でも、本当の感染者は10人中3人以下ということだ。 ということにはならない。 [] [ここ壊れてます] 391 名前：１３２人目の素数さん mailto:sage [2020/04/19(日) 07:50:33.21 ID:Czp86qrf.net] クラスターで１０人が陽性として検査した人数と陽性的中率PPVとの関係をグラフにしてみた。 灰色実線は95%信頼区間境界、灰色点線はPPV＝0.26の線 https://i.imgur.com/T6TDh2r.png 全人口の有病率をクラスター内の有病率にすり替えて、１０人陽性でも感染しているのは３人以下という間違った結論を出している。 わかっていて書いているのか、馬鹿なのか、どちらかは不明。 392 名前：１３２人目の素数さん mailto:sage [2020/04/19(日) 07:56:38.54 ID:Czp86qrf.net] >>368 罹患の定義による。 他人のゲノムでコンタミネーションが起こったりしていなければ、 あるゲノムのシークアンスが検出されたら罹患というなら、罹患率は１００％ ウイルスとして増殖能力を有しているかは不明、死骸の一部を検出しているだけかもしれない。 393 名前：１３２人目の素数さん mailto:sage [2020/04/19(日) 08:01:23.85 ID:Czp86qrf.net] >>369 東京都の発表では　感染者数とせず、陽性患者数と表現している。 正確を期すなら陽性者数とすべきだろうな。 https://stopcovid19.metro.tokyo.lg.jp/ 394 名前：１３２人目の素数さん mailto:sage [2020/04/19(日) 20:39:24 ID:Czp86qrf.net] Natureのこの論文はエクセルとRのコードがついていて自分で再現できるので入力の手間が省ける。 https://www.nature.com/articles/s41591-020-0869-5 感染させる確率分布をガンマ分布を平行移動させた分布として想定して、データから最尤法でパラメータ算出している。 395 名前：１３２人目の素数さん mailto:sage [2020/04/20(月) 17:39:51 ID:Db3kUO+J.net] >>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 らしい 396 名前：１３２人目の素数さん mailto:sage [2020/04/20(月) 19:39:48 ID:Db3kUO+J.net] >>375 補足 ひとりが別の一人に一回感染させるというモデルなので待ち時間の分布の指数分布でいいと思う。 ガンマ分布の形状パラメータ＝１とおけば指数分布になる。 プログラムを書き直してグラフにすると https://i.imgur.com/paErvfa.png 発症前に感染させている確率は47％と原著とあまりかわらないが、そのピークは発症1.6日前という結果になった。 397 名前：１３２人目の素数さん [2020/04/20(月) 21:40:46 ID:LmPkmRXS.net] >>369 >>正解は、約26％だ。判定が「陽性」でも、本当の感染者は10人中3人以下ということだ。 ああ、確かにそこは問題だな。前段の「クラスターで陽性が１０人」とリンクすると 思われても仕方がない。著者もうかつだとは思うが、有病率により的中率が変わると いう話の内容は間違いではない。 398 名前：１３２人目の素数さん [2020/04/20(月) 21:52:20 ID:LmPkmRXS.net] >>371 元記事のグラフの有病率0.1のところをみれば的中率が80％越えになってるわけで、 単純に「クラスター」って言葉を軽んじてただけなんだろうね。 https://webronza.asahi.com/photo/photo.html?photo=/S2010/upload/2020032600006_3.jpg 399 名前：１３２人目の素数さん mailto:sage [2020/04/20(月) 23:44:00 ID:IwCIr2qd.net] >>375 たにんへの感染力の分布なんか数学的にわかるわけないやん？罹患して気道部にウィルスを放出できるくらいのウィルス量が繁殖するのは罹患していきなりなわけがない。 400 名前：１３２人目の素数さん [2020/04/21(火) 00:05:49 ID:bowf2rRh.net] 適当にSEIRモデル� 401 名前：g張してシミュレーションすれば分かるけどかなり自粛しても感染爆発は起こるよ [] [ここ壊れてます] 402 名前：１３２人目の素数さん [2020/04/21(火) 02:51:28 ID:7mnZKVUh.net] https://www.youtube.com/watch?v=4B5LTAdM0NE 東大数学科卒の高橋洋一氏が人との接触7割減」の根拠を数式から 語ってるがどう見る？ 403 名前：１３２人目の素数さん mailto:sage [2020/04/21(火) 06:50:58.83 ID:J5/+FVIQ.net] ガンマ分布：「一定期間に1回起きると期待されるランダムな事象が複数回起きるまでの時間の分布」 複数回でなくて1回だと指数分布だと思う。 何度も感染した症例を扱っているのではないのだから、ガンマ分布のパラメータを求める意味が理解できない。 404 名前：１３２人目の素数さん mailto:sage [2020/04/21(火) 06:56:00.72 ID:J5/+FVIQ.net] >>380 SEIRモデルだと感染爆発させた方が早期に収束する。死者や感染者は増えるけど。 シミュレーションでは鎖国している前提で4割が感染すればオリンピックは可能だった。 SEIRモデルはEは感染力なし、Rは再感染しないというモデルだからなぁ。 モデルを修正してR->Sの再感染が0.１%あるだけで終焉しなかったな。 しかも、外部から感染力の強いキャリアー（保菌者）が入ってくることはモデルには組み込まれていない。 405 名前：１３２人目の素数さん mailto:sage [2020/04/21(火) 06:59:16.79 ID:J5/+FVIQ.net] >>379 結局、一定期間に1回起きると期待されるランダムな事象として、その一定期間を推測しているんだろ？ 406 名前：１３２人目の素数さん mailto:sage [2020/04/21(火) 07:47:46 ID:J5/+FVIQ.net] >>381 SIRモデルって鎖国モデルだから結論は信用できん。 407 名前：１３２人目の素数さん mailto:sage [2020/04/21(火) 11:00:38 ID:LfxF6Y+G.net] 統計できる人尊敬するわ MCMCとか訳分からん 408 名前：１３２人目の素数さん mailto:sage [2020/04/21(火) 11:28:48 ID:J5/+FVIQ.net] >>386 統計ってある統計量がほんにゃら分布に従うというのを黙って受容しないと次に進めないよね。 郡内分散と郡間分散の比がF分布に従うとか言われても どうしてかは理解していない。 stanのNUTSとかfrog leapとかわからんままにMCMCさせている。 409 名前：イナ ◆/7jUdUKiSM mailto:sage [2020/04/21(火) 13:30:52 ID:0ilKIHza.net] こういうときこそコンピューターとかに頼らずに一つ一つ数をかぞえて方程式を立て、微分するべきだと思う。 ‖∩∩‖ □ ‖　　 ＼ （-_-)）　　‖______‖ （っγυ　 ｡‖╂─╂‖ ■`(_)_)ц~ ‖╂─╂‖ ＼■υυ■_∩∩､＼＼‖ ＼＼＼＼⊂(_ _ )`⌒つ) ＼＼＼＼＼＼＼`υ､＼/| ＼＼＼＼＼`.,､､､＼`/ | __＼＼＼＼彡`-`ﾐっ/　L ￣|＼_＼＼_U,~⌒ヾ　/ ］| ‖￣￣￣￣U~~U / / __| ‖ □　□ ‖ |/ / ___`‖________‖／_/ ￣￣￣￣￣￣￣￣￣‖ / __________________‖/ 410 名前：１３２人目の素数さん mailto:sage [2020/04/21(火) 14:53:08 ID:pHehMVs5.net] >>386 sirモデルでは罹患した患者が回復するまでの機会中常に一定の確率(感染率)で遭遇した感受受宿主に感染を広めると仮定して立式してる。 ザックリした傾向を見るならそれで十分だけど、実際のモデルではそんな事はありえない。 感染する確率は患者の体内で繁殖しているウィルス量の増加に従って増えるからその効果を勘案してΓ分布(というかt^c d^x なる形の関数)に応じて感染率が罹患した時点から変化するとするんでしょ？ もちろんコレでも大体の傾向見るにはコレで十分というモデルを臨床例から適当に選んでるだけでしょ？ 実際には未来永劫感染力を持ち続けるなんて事があるはずないし、罹患初期のウィルス量の増え方はおそらく指数関数的に増えるものを採用すべきだろうし。 411 名前：１３２人目の素数さん mailto:sage [2020/04/21(火) 14:54:06 ID:pHehMVs5.net] mcmc 412 名前：１３２人目の素数さん mailto:sage [2020/04/22(水) 08:56:14 ID:0tlpvLKp.net] >>388 無理無理。 カブトガニVSシオマネキ論争はコンピュータの計算結果で決着がついたよ。 413 名前：１３２人目の素数さん mailto:sage [2020/04/22(水) 13:23:53 ID:0tlpvLKp.net] 週末は休むとか週間変動の影響を除くために１週間の移動平均で線形回帰して片対数 414 名前：グラフにすると https://i.imgur.com/7VwfswD.png 自粛の効果がでてきているな。　多分、検査自粛の効果だろうな。 [] [ここ壊れてます] 415 名前：イナ mailto:sage [2020/04/22(水) 19:51:00.75 ID:iq1GZOqA.net] ‖∩∩ ‖ □ ‖前>>388 （(-_-)‖　　‖______ （っ⌒⌒ﾞ　 ｡‖╂─╂ ■`(_)_)ц~　‖╂─╂ ＼■υυ■_∩∩､＼＼＼ ＼＼＼＼⊂(_ _ )`⌒つ､ ＼＼＼＼＼＼＼`υ､＼＼＼＼＼＼＼＼＼＼＼＼＼＼＼＼`>>391カブトガニでもシオマネキでもシュクメルリでも微分するのがいちばん強力だと思う。 416 名前：１３２人目の素数さん mailto:sage [2020/04/23(木) 04:58:29 ID:3vcirHk0.net] >>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も正直者である」 417 名前：１３２人目の素数さん mailto:sage [2020/04/23(木) 08:28:17 ID:3vcirHk0.net] 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') 418 名前：１３２人目の素数さん mailto:sage [2020/04/23(木) 10:37:41 ID:3vcirHk0.net] 中国には特異度１００％の検査キットがあるんだってね。 すべて陰性にでるようにセットされていると。 419 名前：１３２人目の素数さん mailto:sage [2020/04/23(木) 13:51:17 ID:3vcirHk0.net] 新型コロナ患者を治療している病院に100人の職員がいる。 検体採取器具は5人分、試薬は1回分しかないとする。 無作為抽出した5人の職員から採取した検体を混合して検査したら陽性であった。 職員の陽性者数の期待値を求めよ。 また、50人以上の感染者いる確率はいくつか？ 検査の陽性率はハイリスク群に検査している東京の数値2457/6654を使って計算せよ。 https://i.imgur.com/zYK75Lo.jpg 420 名前：１３２人目の素数さん mailto:sage [2020/04/23(木) 13:53:25 ID:3vcirHk0.net] >>397 こんなグラフになった。 https://i.imgur.com/gortYDb.png 421 名前：１３２人目の素数さん mailto:sage [2020/04/23(木) 17:36:48 ID:3vcirHk0.net] >>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 422 名前：１３２人目の素数さん mailto:sage [2020/04/23(木) 22:02:05.89 ID:MtjaFZpr.net] レベル低 423 名前：１３２人目の素数さん [2020/04/24(金) 02:55:58 ID:juJsFFfP.net] 慶応大の調査で、コロナ以外で来院した人をPCR検査したところ 4/67の確率で要請だった 東京都内1500万のうち何人くらいが感染しているか推定せよ 424 名前：１３２人目の素数さん [2020/04/24(金) 03:18:37 ID:XPGerQAq.net] 岩田健太郎・� 425 名前：_戸大学教授『東京はすでに20万〜400万人感染の可能性』 https://leia.5ch.net/test/read.cgi/poverty/1587664106/ [] [ここ壊れてます] 426 名前：１３２人目の素数さん mailto:sage [2020/04/24(金) 05:36:50.27 ID:9Fe9PNfV.net] >>401 ニュー速でやったが、ベイズでこんな感じ？(感染者率) https://i.imgur.com/Nfs1vCm.jpg 427 名前：１３２人目の素数さん mailto:sage [2020/04/24(金) 05:52:14 ID:9Fe9PNfV.net] >>397 混ぜて検査する方法の最適化問題ってのもあるね。 https://mobile.twitter.com/p_gotcha/status/1243501000943702017 (deleted an unsolicited ad) 428 名前：１３２人目の素数さん mailto:sage [2020/04/24(金) 08:13:06 ID:0onl6lJy.net] >>401 これって東京の無作為サンプリングではなさそうだよね 病院に来たって事は熱が出てたのかもしれないし 想定できる母集団ってなんになるのだろうか 429 名前：１３２人目の素数さん mailto:sage [2020/04/24(金) 12:07:56 ID:v55OWzbu.net] 新型コロナ患者を治療している病院に100人の職員がいる。 検体採取器具は10人分、試薬は1回分しかないとする。 無作為抽出した10人の職員から採取した検体を混合して検査したら陰性であった。 職員の陽性者数の期待値を求めよ。 また、50人以上の感染者いる確率はいくつか？ 検査の陽性率はハイリスク群に検査している東京の数値2457/6654を使って計算せよ。 430 名前：１３２人目の素数さん mailto:sage [2020/04/24(金) 12:11:36 ID:v55OWzbu.net] >>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 431 名前：１３２人目の素数さん [2020/04/24(金) 17:39:04 ID:8oiI190P.net] >>401 まじか。 まあ、体調悪いから病院に行くわけで、バイアスかかってるとはいえ...。 432 名前：１３２人目の素数さん mailto:sage [2020/04/24(金) 17:42:10 ID:v55OWzbu.net] >>404 ラテン方陣の問題？ 433 名前：１３２人目の素数さん mailto:sage [2020/04/24(金) 20:48:11 ID:v55OWzbu.net] >>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) 434 名前：１３２人目の素数さん mailto:sage [2020/04/24(金) 20:55:16 ID:v55OWzbu.net] >>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 435 名前：１３２人目の素数さん mailto:sage [2020/04/24(金) 21:13:10 ID:v55OWzbu.net] >>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 436 名前： 10 prop.test 832836 269206 2140919 11 wilson 832836 327258 2004545 岩田の計算は 5 exact 832836 230236 2034792 [] [ここ壊れてます] 437 名前：１３２人目の素数さん mailto:sage [2020/04/24(金) 21:35:58 ID:v55OWzbu.net] >>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 438 名前：１３２人目の素数さん mailto:sage [2020/04/24(金) 21:50:22 ID:v55OWzbu.net] >>413 有病率の最頻値は > density(prev)$x[which.max(density(prev)$y)] [1] 0.019186 439 名前：１３２人目の素数さん mailto:sage [2020/04/24(金) 22:01:25 ID:v55OWzbu.net] >>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 440 名前：１３２人目の素数さん mailto:sage [2020/04/24(金) 22:03:05.10 ID:v55OWzbu.net] >>415 この最頻値はオーストリアの0.3％という値に等しいな。 441 名前：１３２人目の素数さん mailto:sage [2020/04/24(金) 23:51:07 ID:v55OWzbu.net] >>415 Jeffreys が正しいスペリングみたい。　Jefferey'sかと思っていた。 442 名前：１３２人目の素数さん mailto:sage [2020/04/25(土) 13:54:43 ID:R/UD6QQG.net] 某医療センターでは治療後、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 443 名前：１３２人目の素数さん mailto:sage [2020/04/25(土) 17:53:16 ID:MKtk31sc.net] >>409 動的計画法を使うらしい。 https://qiita.com/p_gotcha/items/e2741b12b4f606b96136 444 名前：１３２人目の素数さん mailto:sage [2020/04/26(日) 21:43:36 ID:0lgRXnyr.net] スレチですが統計に詳しい方がいると思って伺いました。 工学系の大学院生で論文を読んでて、見慣れない記号や数式が出てきたので質問させてください。 この数式中のEは何を意味するのでしょうか？ Rの期待値的なものですか？ https://i.imgur.com/oVwva8r.jpg 445 名前：１３２人目の素数さん mailto:sage [2020/04/27(月) 02:54:57 ID:cS0ISst+.net] 院生でわからないってマジか 専門外の機械学習いきなりやらされたのかな？ 446 名前：１３２人目の素数さん mailto:sage [2020/04/27(月) 06:45:22 ID:S7AmHM83.net] 有病率が0.3％（オーストリアの無作為抽出)や6%(慶応大学の入院患者無作為抽出）のときは 感度30〜70%　特異度90〜100%のPCR検査キットでは 陰性的中率は オーストリアの例では　95％[87％〜99％] 慶応の例では 94%[86〜99％]になる。 すべて陰性にでる中国のイカサマキットなら 陰性的中率は オーストリアの例では99.7％ 慶応の例では94%になる。 有病率が10％くらいになればイカカマキットの方が陰性的中率の成績が劣る。 447 名前：１３２人目の素数さん mailto:sage [2020/04/27(月) 07:30:42 ID:S7AmHM83.net] 某医療センターでは治療後、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 448 名前：１３２人目の素数さん mailto:sage [2020/04/28(火) 10:03:02.25 ID:/p2virDs.net] 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 449 名前：_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 [] [ここ壊れてます] 450 名前：１３２人目の素数さん mailto:sage [2020/04/28(火) 21:14:48.22 ID:V+QD5qnX.net] >>394 答えなくね？ 451 名前：１３２人目の素数さん mailto:sage [2020/04/29(水) 16:59:02 ID:YhsQxAnp.net] >>425 正解。正直者が誰もいないようにプログラムした。 イナ師が微分で解いてくれるのを期待していたんだが。 452 名前：１３２人目の素数さん mailto:sage [2020/04/29(水) 17:43:56 ID:YhsQxAnp.net] PCR検査は 感度0.3~0.7 : β分布(2.625,2.625)相当 特異度 0.95~1.0 : β分布(26.5014,2.3422)相当 として、 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 453 名前：１３２人目の素数さん mailto:sage [2020/04/29(水) 18:02:18.86 ID:FBYarVUq.net] >>426 正直者の確定が誰もいないというより、条件を満たす正直者・嘘つきの組み合わせが存在しないよね？ https://i.imgur.com/NJw3GJN.jpg 454 名前：１３２人目の素数さん [2020/04/29(水) 18:35:59 ID:uFfpYtab.net] rRT-PRC検査の感度・特異度ともに90%以上です。 455 名前：１３２人目の素数さん [2020/04/29(水) 19:40:54.19 ID:uFfpYtab.net] 通常の医療でも臨床診断の5%は誤診だと推定されていたりする。 これを実際の数として見積もると膨大になる。 456 名前：１３２人目の素数さん mailto:sage [2020/04/29(水) 20:37:12 ID:YhsQxAnp.net] >>429 採取手技や採取時期で変動すると思う。 コンタミによる偽陽性もあるだろうし。 457 名前：１３２人目の素数さん mailto:sage [2020/04/29(水) 20:38:43 ID:YhsQxAnp.net] >>428 言っていることは同じ。 全部の条件を満たす組み合わせは存在しない。 458 名前：１３２人目の素数さん mailto:sage [2020/04/29(水) 20:39:49 ID:YhsQxAnp.net] >>427 このデータは再現性がないことがわかったので撤回します。 459 名前：１３２人目の素数さん mailto:sage [2020/04/30(木) 10:01:57 ID:42TM9o6B.net] ふぉー ＜新型コロナ＞抗体検査５．９％陽性　市中感染の可能性　都内の希望者２００人調査 https://www.tokyo-np.co.jp/s/article/2020043090070748.html 460 名前：１３２人目の素数さん mailto:sage [2020/04/30(木) 12:13:53 ID:qOF+URFa.net] >>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 461 名前：１３２人目の素数さん mailto:sage [2020/04/30(木) 19:13:48 ID:42TM9o6B.net] >>434 これの興味深いところは、陽性率の高さもさることながら、医療関係者の感染率が一般市民のそれよりも （多分有意に）高いことだと思う。抗体検査の精度ってほとんどデータが無いに等しいけど、医療関係者の 感染者率が平均よりも多いのならば（実際そうなのだと思う）、このテストはそれを反映したものと言えて、 抗体検査の信頼性をある程度証明できるかもしれない。 そして、市中の陽性率の高さも。 この辺、ベイズで上手く検証できないかな。 462 名前：１３２人目の素数さん mailto:sage [2020/04/30(木) 20:32:35.24 ID:qOF+URFa.net] >>436 >医療関係者の感染率が一般市民のそれよりも（多分有意に）高い χ二乗検定では、p-value = 0.41 で　標本数が少ないから有意差がでない。 エントリーが５以下のときには信頼するなと習ったな。 まあ、ポアソンでもp-value = 0.33 >435にRの出力を貼っておいた。 463 名前：１３２人目の素数さん [2020/04/30(木) 20:35:20.83 ID:rbz2947p.net] >>436 55人中5人と147人中7人じゃ、有意差ないだろ。 と書こうと思ったら、 >>437氏に先を越されたｗｗｗ 464 名前：１３２人目の素数さん [2020/04/30(木) 20:41:39 ID:rbz2947p.net] 抗体検査キットの評価をしたら特異度はみな5/5だったらしいけど、サンプルが ５つだけって意味だとすると、せいぜい90％以上くらいのことしか言えないな。 特異度95%だとしたら、５％が陽性っていわれてもねぇｗ 465 名前：１３２人目の素数さん mailto:sage [2020/04/30(木) 21:20:01 ID:qOF+URFa.net] >>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%の一様分布というのは現実離れした分布だとは思う。 466 名前：１３２人目の素数さん mailto:sage [2020/04/30(木) 21:22:46 ID:qOF+URFa.net] >>439 すべてが陰性に出るイカサマキットなら特異度が100％ >18参照 467 名前：１３２人目の素数さん mailto:sage [2020/04/30(木) 21:33:35 ID:qOF+URFa.net] >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) uniro 468 名前：ot(function(x) fn(x)-0.05,c(1,10))$root fn(4) > fn(4) [1] 0.045678 [] [ここ壊れてます] 469 名前：１３２人目の素数さん mailto:sage [2020/04/30(木) 21:47:02.48 ID:qOF+URFa.net] >>440（追記) >有病率が平均値50%の一様分布というのは現実離れした分布だとは思う。 オーストリアのデータ0.3%を参考に有病率の事前分布が0.2-0.4％ β(10.865, 3279.34)に相当として、MCMCしてみた。 感度と特異度の事前分布は前回と同じ。 https://i.imgur.com/dHcOxsR.png さほど、有病率が高いという結論は引き出せないな。 470 名前：１３２人目の素数さん [2020/04/30(木) 22:23:36 ID:gdjT5b3G.net] 自己流な検定 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 471 名前：１３２人目の素数さん mailto:sage [2020/04/30(木) 23:17:24 ID:qOF+URFa.net] β分布と乱数発生で比較。 一般人と医療従事者の確率分布に従う乱数を１０００万個発生させて比率の分布をグラフにすると ９５％信頼区間が１を挟むから有意差なし。比率が１以上となる確率も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) 472 名前：１３２人目の素数さん [2020/05/01(金) 00:36:14 ID:LcUyF6Tc.net] >>441 すべてが陰性に出るなら陽性は０人だろｗ 473 名前：１３２人目の素数さん mailto:sage [2020/05/01(金) 07:22:53 ID:LJVZP1pw.net] >>446 偽陽性率０％だから特異度は１００％ 474 名前：１３２人目の素数さん [2020/05/01(金) 10:05:01 ID:LcUyF6Tc.net] >>447 だから、陽性が何人か出てんだから感度0%のインチキじゃなかろうってこと。 475 名前：１３２人目の素数さん mailto:sage [2020/05/04(月) 00:17:48 ID:Q8iLGwO2.net] 最近の報道 インドは、中国企業から購入した中共ウイルス（新型コロナウイルス）の迅速スクリーニング検査キットの精度がわずか5％だとして、約50万個の注文をキャンセルした。 .... 神戸市中央区にある市立医療センター中央市民病院の医師などのグループは、ことし３月末から先月７日にかけて、新型コロナウイルス以外の理由で外来を受診した患者から無作為に1000人を選び、血液中に新型コロナウイルスに感染したあとにできる「抗体」があるか調べました。 グループによりますとその結果、3.3％にあたる33人から抗体が検出されたということです。 以上を知ったある会社が 必ず陽性がでる試薬33個と必ず陰性がでる試薬967個を混ぜた1000試薬をセットにしてインドに売り込んだ。 問題 1試薬は1回しか検査できないとして これがイカサマキット� 476 名前：ﾅあることを証明する手段はあるか？ [] [ここ壊れてます] 477 名前：１３２人目の素数さん [2020/05/04(月) 15:22:29 ID:jDRWX2Ph.net] 3月の宿題で(1)のみ正解の数弱@shukudai_sujaku 昨年度の大学への数学(大数)での勝率は、 学コンBコースが 1/1 = 100% ， 宿題が 3/10 = 30% でした！ 宿題の勝率が低すぎると思うので、 これからは一層精進していきたいです！ https://twitter.com/shukudai_sujaku (deleted an unsolicited ad)

---

---

[ 続きを読む ] / [ 携帯版 ]

read.cgi ver5.27 [feat.BBS2 +1.6] / e.0.2 (02/09/03) / eucaly.net products.
担当:undef