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

1 名前：１３２人目の素数さん [2020/02/29(土) 02:18:41 ID:twdO677Q.net] 東大数学科卒の元官僚はこう分析してるが、お前らはどうなると思う？ www.zakzak.co.jp/soc/news/200220/dom2002200003-n2.html 中国国外感染者の中国国内との比率をみると、 １月２０日の数字公表以降は、０・８〜２・６％で比較的安定している。 これは、新型肺炎の感染者のほとんどは中国国内、それも湖北省に集中しているからだ。 ちなみに中国国外での感染者数は、中国国内の１・１％だ（２月１６日現在）。 本コラムで紹介したが、現時点では、最終的な中国国内の感染者数は２０万人超と筆者は推計している。 となると、中国国外の感染者は数千人程度になるだろう。 中国国外のうち日本の比率は１割弱なので、日本の感染者数は数百人程度であろう。 その場合、死者も数人から１０人程度になるだろう。 こうした推計をすると、今の感染者は氷山の一角だと思われるが、今後の増加ペースはどうなるだろうか。 新型コロナウイルスの検査は簡単に行えるので、今後、日本での感染者数は増えていくだろう。 ある時点ではそれがネズミ算的に増えるかのように思える局面もあるだろうが、 筆者の推計が正しければ、現時点ではせいぜい数百人が一つのメドだ。 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) 478 名前：１３２人目の素数さん [2020/05/04(月) 22:45:28 ID:94tg6i4k.net] 高橋洋一(統計数理研究所→大蔵省) スウェーデンの感染症対策を解説 https://youtu.be/2bsYRuorwMI?t=192 479 名前：１３２人目の素数さん mailto:sage [2020/05/05(火) 06:07:12 ID:ht6rG86e.net] 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％のモデルの方がいいかも、という結果。 陽性数/検査数の時系列データでもあればもう少し差がでるかもしれん。 東京都のデータで計算させようかと思ったが、東京都は検査人数を隠蔽しているので使いものにならない。 480 名前：１３２人目の素数さん mailto:sage [2020/05/05(火) 06:26:36.20 ID:FpoKo6a+.net] 訂正 陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする ↓ 陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする 481 名前：１３２人目の素数さん [2020/05/05(火) 11:13:29 ID:b2IqdVzK.net] 3月の宿題で(1)のみ正解の数弱@shukudai_sujaku 昨年度の大学への数学(大数)での勝率は、 学コンBコースが 1/1 = 100% ， 宿題が 3/10 = 30% でした！ 宿題の勝率が低すぎると思うので、 これからは一層精進していきたいです！ https://twitter.com/shukudai_sujaku (deleted an unsolicited ad) 482 名前：１３２人目の素数さん [2020/05/05(火) 23:00:17 ID:wvtmzVA1.net] >>449 抗体のはいったサンプルを２，３０個試してみりゃわかんじゃないの？ いくらなんでも感度が３％じゃつかえねぇ。 483 名前：１３２人目の素数さん mailto:sage [2020/05/06(水) 15:29:26 ID:RG00+xls.net] イカサマキットの感度特異度の事前分布を一様分布に設定して 抗体のはいったサンプルを２０個で全部陰性であったので３０個試したら全部陰性であったとすると 感度・特異度の事後分布は https://i.imgur.com/6ZYn34k.png 484 名前：１３２人目の素数さん mailto:sage [2020/05/06(水) 15:32:34 ID:RG00+xls.net] ところが、有病率33/1000であったときに無作為に２０人および３０人を選んで全部陰性であっても 感度・特異度の事後分布は https://i.imgur.com/jWg4RgH.png 485 名前：１３２人目の素数さん mailto:sage [2020/05/06(水) 15:52:13 ID:RG00+xls.net] >>454 事後分布を一様分布に設定して、このコピペ小僧の正解率の期待値・最頻値・95％信頼区間を求めよ。 その結果、 https://i.imgur.com/M9SdCSL.png 486 名前：１３２人目の素数さん mailto:sage [2020/05/06(水) 15:58:11 ID:RG00+xls.net] >>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 487 名前：１３２人目の素数さん mailto:sage [2020/05/06(水) 16:25:45 ID:RG00+xls.net] >>454 このコピペ小僧の正解率に学コンと宿題で差があるかを検定せよ。 事前分布にJefferysを用いた結果、 https://i.imgur.com/2yYagtX.png 488 名前：１３２人目の素数さん [2020/05/06(水) 17:08:09 ID:tNZVV9ZH.net] >>456 感度が5%以下じゃつかいもんにならんわなぁ。インチキで確定。 489 名前：１３２人目の素数さん mailto:sage [2020/05/06(水) 23:23:13 ID:wABsYddm.net] https://www.buzzfeed.com/jp/naokoiwanaga/covid-19-antibody-test 大阪市立大学が一昨年の血液を使って検査したところ、偽陽性は1人もいなかったそうだ。 490 名前：１３２人目の素数さん mailto:age [2020/05/07(木) 00:59:54 ID:0vbdeA0/.net] >>462 >2020年4月中の2日間に同大学の付属病院を、新型コロナウイルス感染症の診療以外で受診した >患者を対象に、そこから無作為に312人を抽出・・・・312 人（年齢中央値 66.5 歳、 >男性：女性=154：158）のうち、3人が陽性であることがわかった。約1%の陽性率だ。 >統計的な誤差を考慮すると、95%の確率で0.33〜2.8%の間に入る・・・・・・・・・・・・ 教えてエロい人 Ｑ１統計的推定を述べるなら「95%の確率で」でなく「信頼率95%で」と書くのが適切ではないか？ Ｑ２無作為抽出陽性率3/312=0.96%の母集団陽性率上側/下側信頼区間=0.96-0.33/2.80-0.96と 　　上側区間≠下側区間で左右非対称分布想定は何故？ Ｑ３サンプルサイズ312は過少では？過少に伴う推定誤差加算が必要では？ 491 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 02:09:04.11 ID:VnQvkZ57.net] >>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 492 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 02:30:17 ID:VnQvkZ57.net] 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)) 493 名前：１３２人目の素数さん [2020/05/07(木) 17:10:06 ID:CHL0/p02.net] >>462 ５０人で０でしょ？ だったら、特異度として保証されるのはせいぜい９８％程度じゃん。 １％の陽性率に対して、それでは不十分。 494 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 18:13:41.63 ID:VnQvkZ57.net] >>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 495 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 18:52:47.08 ID:VnQvkZ57.net] 特異度の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 496 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 20:13:12 ID:VnQvkZ57.net] 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 497 名前：１３２人目の素数さん [2020/05/07(木) 21:49:37 ID:CHL0/p02.net] >>469 なにやってるかよくわからんので、見当外れの指摘かもしれんが、 大阪市大の3/312と神戸大の33/1000ってのは特異度も感度も異なる であろう別種のキットによる検査結果なんで、一緒くたにしちゃ まずいんでないの？ 498 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 22:36:59.75 ID:VnQvkZ57.net] >>470 同一キットじゃないから、ご指摘の通り。 しかも神戸大の方では陰性検体での確認はされていないから、神戸大方の陽性率が高いのは偽陽性を含む可能性もあるね。 大阪市大だけのデータでやってみると。 https://i.imgur.com/njVQtRZ.png 499 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 22:46:09.05 ID:VnQvkZ57.net] >>471 使える情報が少なすぎて信頼区間が広すぎ。 有病率:0-87% キットの感度:0-70% ってコイントス変わらんな。 500 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 22:48:00.40 ID:VnQvkZ57.net] 結局のところ、断定的な結論は出せないねということを数字で確認しているだけｗ 501 名前：１３２人目の素数さん mailto:sage [2020/05/08(金) 02:04:52 ID:Ugc87SUM.net] >>464 >>465 有難うございます。 二項分布下の母比率の信頼区間推定問題なんですね。 A２：正規分布近似法でなくWilsonの信頼区間法で解くから 　　上側区間≠下側区間となり平均値まわり左右非対称信頼区間と 　　なるのですね。 502 名前：１３２人目の素数さん [2020/05/08(金) 09:15:37 ID:WmDpVhCu.net] 3月の宿題で(1)のみ正解の数弱@shukudai_sujaku 昨年度の大学への数学(大数)での勝率は、 学コンBコースが 1/1 = 100% ， 宿題が 3/10 = 30% でした！ 宿題の勝率が低すぎると思うので、 これからは一層精進していきたいです！ https://twitter.com/shukudai_sujaku (deleted an unsolicited ad) 503 名前：１３２人目の素数さん mailto:sage [2020/05/08(金) 11:37:38 ID:CfFk/Uw1.net] >>474 Rのlibrary binomを使って　binom::confint(3,312)で　>464の出力が得られる 信頼区間をグラフにすると https://i.imgur.com/Mlvv7bB.png 点線は3/312 正規分布近似のasymptotic以外は非対称。 どれを使うべきか？　好きなのを使えばいい、と思う。 但し、値が負になったり１を越えたりするのは不採用の方が賢明だとは思う。 Wilson法は値が０や１に近くても信頼できるという人がいるけど、どうやって検証するのかはよくわからん。 504 名前：１３２人目の素数さん mailto:sage [2020/05/08(金) 12:41:19 ID:CfFk/Uw1.net] >>476 binom::confintはbinom.confint(3,312)の間違い library("binom") binom.confint(3,312) 505 名前：１３２人目の素数さん mailto:sage [2020/05/09(土) 13:22:10 ID:oDT7bFgO.net] レムデシビルで初の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 なので、どちらのモデルが得られたデータを説明するのに適しているかは判断できず 506 名前：、という結果になった。 検出力不足か、差がないのかには決着つかず。 [] [ここ壊れてます] 507 名前：１３２人目の素数さん mailto:sage [2020/05/09(土) 13:22:52 ID:oDT7bFgO.net] 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) 508 名前：１３２人目の素数さん mailto:sage [2020/05/09(土) 13:25:29 ID:oDT7bFgO.net] 確率分布を図示すると https://i.imgur.com/hfH7gZO.png 509 名前：１３２人目の素数さん [2020/05/09(土) 14:25:05 ID:74hNX8Dr.net] 100例もやって差がでないのなら、はっきり言って効き目なしとかわらん。 むしろ副作用がこわい。 アビガンはどうなんだろうね。なぜか国内のまとまったデータがまだ出てこない。 510 名前：１３２人目の素数さん [2020/05/09(土) 14:26:34 ID:74hNX8Dr.net] まあ、軽症者が８割とかだと、早期投与の効き目を判定するには相当の 人数に使わないとわかんないだろうけど。 511 名前：１３２人目の素数さん mailto:sage [2020/05/09(土) 15:22:40 ID:oDT7bFgO.net] 死亡率に関してはベイズファクターで推論できた。 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を超えないいのでまあ、中程度の確信で検出力不足によるのではなく、もともと差がないのであろうと推論できる。 サンプルサイズを増やしてもレムデシビルは軽症・早期投与でも致死率を改善しないであろうと予想。 512 名前：１３２人目の素数さん mailto:sage [2020/05/09(土) 17:15:11 ID:oDT7bFgO.net] >どの週においても両群で症状改善率は不変 この前提はおかしいので週ごとに症状改善率が変わるようにモデルを変更 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.997627 513 名前：4 症状改善率に差がないというモデルも差があるというモデルも周辺尤度（エビデンス）はほぼ同じという結果。 [] [ここ壊れてます] 514 名前：１３２人目の素数さん [2020/05/09(土) 20:41:59 ID:74hNX8Dr.net] あ、 >>482はレムデシじゃなくてアビガンについての話だかんね。 レムデシは重症者向けらしいけど、たぶん駄目だろ。 515 名前：１３２人目の素数さん [2020/05/09(土) 21:35:34 ID:Fh34v+Vz.net] 【政府が隠す３．１１】　地震波形が、核実験と同じ rio2016.5ch.net/test/read.cgi/earth/1575170033/l50 sssp://o.5ch.net/1nmad.png 516 名前：１３２人目の素数さん [2020/05/11(月) 00:31:07 ID:mGn48zmS.net] 5/10夜『NHKスペシャル』「新型コロナウイルス　出口戦略は」で 厚労省新型コロナ対策班西浦博北大教授と氏が分析中のノートPC画面が 映されていた。氏の使用解析ソフトウェア名は何というのでしょうか？ 517 名前：１３２人目の素数さん mailto:sage [2020/05/11(月) 01:36:00 ID:P7wrhagU.net] 西浦先生はRじゃないかな。 今週ぐらいに、実効再生産数を算出するRのコードを公開できるかもと言ってた。 火曜にはニコ生があるから、そこでたっぷりと聞けるはず。 https://sp.live.nicovideo.jp/watch/lv325833316 518 名前：１３２人目の素数さん mailto:sage [2020/05/11(月) 03:07:28 ID:6TI8KNr5.net] >>488 >11に数式があるけど 519 名前：１３２人目の素数さん mailto:sage [2020/05/11(月) 03:58:21.10 ID:v4rWY6J/.net] >>489 いや、与えられた何日分かの新規患者数から再生産数とか回復率とかを推定するための数式じゃないの？ 520 名前：１３２人目の素数さん mailto:sage [2020/05/11(月) 06:53:45 ID:ruo6tAnf.net] >>490 SEIRより単純なReed-Frostモデルだよな https://en.m.wikipedia.org/wiki/Reed%E2%80%93Frost_model 521 名前：１３２人目の素数さん mailto:sage [2020/05/11(月) 09:14:46.53 ID:P7wrhagU.net] >>489 自分で計算してみた人ならわかるけど、日々の再生算数を出すにはA「感染してから他人にうつすまでの日数の分布」が必要。西浦先生はこれをまず求めていて、その結果は割合と世界中で使われている。 加えてB「感染してから発表の数字に載るまでの日数の分布」も必要で、発表の数字をBで逆畳み込み積分してC「感染推定日毎の数字」に変換し、CをAで畳み込み積分したものでCを割れば日々のRtが求められる。 とは言うものの、そのままやると誤差で数字が発散するし、逆畳み込みは一意ではないのでどう推定するかは色々と考えなくてはならない。多分日本独自の事情もコードに入れなきゃいけない。大変だろうなあとは思う。 522 名前：１３２人目の素数さん mailto:sage [2020/05/11(月) 11:02:15.57 ID:+2HA8Yn9.net] >>492 レスありがとうございます。 確かにSEIRモデルやSIRIモデルだと 「感染してから発表の数字に載るまでの日数の分布」って考慮されてないね。 523 名前：１３２人目の素数さん mailto:sage [2020/05/12(火) 20:50:54 ID:8HqG++pl.net] >>488 公開されたぞー。StanとRだ。 https://nbviewer.jupyter.org/github/contactmodel/COVID19-Japan-Reff/blob/master/scripts/C.%20Calculating%20the%20Rt%20in%20Stan.ipynb 524 名前：１３２人目の素数さん mailto:sage [2020/05/12(火) 22:32:45.27 ID:9+i4T0w3.net] ウイルス感染のモデルも必要だと思うけど 社 525 名前：会経済活動の影響や 自粛解除した後の死者増加や 自粛継続した場合の法人の倒産や 治療薬やワクチンができた後の社会の経済や文化面の回復 そこまで考慮してどういう政策を選択したら 損失を最小化できるかという問題は解ける？ 制約付きの最小化問題だと思う 数値化するのが難しい考慮要素もあるからその辺は仮の値とかで [] [ここ壊れてます] 526 名前：１３２人目の素数さん mailto:sage [2020/05/13(水) 17:23:13 ID:PBQsFM+W.net] >>495 モデルを精密複雑化すればするほど、必要なパラメータが増えて それを弱情報事前分布にすると事後分布の信頼区間が広がるんだなぁと思う。 527 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 08:00:23 ID:kec+XbRE.net] >>494 これにある、 JapaneseDataCOVID19 ってどこかからダウンロードできるんだろうか？ 528 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 09:54:38 ID:IN0uokww.net] >>497 https://github.com/contactmodel/COVID19-Japan-Reff/tree/master/data 529 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 11:47:52 ID:kec+XbRE.net] >>498 ありがとうございます。 530 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 11:51:11 ID:kec+XbRE.net] 親ディレクトリ（フォルダ）探せばよかったんだね。 https://nbviewer.jupyter.org/github/contactmodel/COVID19-Japan-Reff/tree/master/ 531 名前：１３２人目の素数さん [2020/05/14(木) 13:08:10.58 ID:gLQVRit5.net] >>495 >>496 パラメータを増やせば増やすほど、どんな仮説も否定できなくなる。 実験系の学界でよく使う手口だよね。 おや、誰か来たようだ… 532 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 13:58:09 ID:kec+XbRE.net] >>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にしているからみたい。 533 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 14:10:15 ID:kec+XbRE.net] 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 534 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 16:35:35 ID:kec+XbRE.net] 再生算数を０〜１０人の一様分布にすると、収束しない。 > 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 535 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 17:00:26 ID:1Kd4DHQt.net] >>501 否定できないけど実証もできないことになるのでは？ 単純なモデルだと現実の問題を解決できない事もあるし 536 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 17:23:23.34 ID:kec+XbRE.net] >>494 モデルで再生産数の事前分布は　平均2.5　標準偏差2.0の正規分布に設定されていたので 平均と標準偏差を変化させて、再生産数の事後分布を描出してみた。 かなり、事前分布の影響を受けるみたい。 https://i.imgur.com/OwqsFC1.png 537 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 17:26:03.64 ID:kec+XbRE.net] どうも、こういう境地だなぁ。 断定的な結論は出せないということを数字で確認しているだけｗ 538 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 17:30:34 ID:50GljeDM.net] 基本再生算数はウイルスの特徴で決まるもので 実行再生算数は更に環境や人の行動の影響で変化するものだと理解している 実行再生算数を減らすようにするには どう行動したらいいか かつ経済活動を出来るだけ下げずに 539 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 18:25:32 ID:IN0uokww.net] 西浦モデル、一度感染日ごとの人数を最尤法で確定して、そこからベイズを回してRtを求めてるんだよなあ。 本来は、日々の感染者数も確率的に揺れがあるはず。だからRtの誤差幅は発表よりも多く見積もるべきかもしれん。

---

---

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

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