20/04/21 00:05:49 bowf2rRh.net
適当にSEIRモデル�
401:g張してシミュレーションすれば分かるけどかなり自粛しても感染爆発は起こるよ
402:１３２人目の素数さん
20/04/21 02:51:28 7mnZKVUh.net
URLﾘﾝｸ(www.youtube.com)
東大数学科卒の高橋洋一氏が人との接触7割減」の根拠を数式から
語ってるがどう見る？
403:１３２人目の素数さん
20/04/21 06:50:58.83 J5/+FVIQ.net
ガンマ分布：「一定期間に1回起きると期待されるランダムな事象が複数回起きるまでの時間の分布」
複数回でなくて1回だと指数分布だと思う。
何度も感染した症例を扱っているのではないのだから、ガンマ分布のパラメータを求める意味が理解できない。
404:１３２人目の素数さん
20/04/21 06:56:00.72 J5/+FVIQ.net
>>380
SEIRモデルだと感染爆発させた方が早期に収束する。死者や感染者は増えるけど。
シミュレーションでは鎖国している前提で4割が感染すればオリンピックは可能だった。
SEIRモデルはEは感染力なし、Rは再感染しないというモデルだからなぁ。
モデルを修正してR->Sの再感染が0.１%あるだけで終焉しなかったな。
しかも、外部から感染力の強いキャリアー（保菌者）が入ってくることはモデルには組み込まれていない。
405:１３２人目の素数さん
20/04/21 06:59:16.79 J5/+FVIQ.net
>>379
結局、一定期間に1回起きると期待されるランダムな事象として、その一定期間を推測しているんだろ？
406:１３２人目の素数さん
20/04/21 07:47:46 J5/+FVIQ.net
>>381
SIRモデルって鎖国モデルだから結論は信用できん。
407:１３２人目の素数さん
20/04/21 11:00:38 LfxF6Y+G.net
統計できる人尊敬するわ
MCMCとか訳分からん
408:１３２人目の素数さん
20/04/21 11:28:48 J5/+FVIQ.net
>>386
統計ってある統計量がほんにゃら分布に従うというのを黙って受容しないと次に進めないよね。
郡内分散と郡間分散の比がF分布に従うとか言われても
どうしてかは理解していない。
stanのNUTSとかfrog leapとかわからんままにMCMCさせている。
409:イナ ◆/7jUdUKiSM
20/04/21 13:30:52 0ilKIHza.net
こういうときこそコンピューターとかに頼らずに一つ一つ数をかぞえて方程式を立て、微分するべきだと思う。
∥∩∩∥ □ ∥ ＼
（-_-)） ∥______∥
（っγυ ｡∥╂─╂∥
■`(_)_)ц~ ∥╂─╂∥
＼■υυ■_∩∩､＼＼∥
＼＼＼＼⊂(_ _ )`⌒つ)
＼＼＼＼＼＼＼`υ､＼/|
＼＼＼＼＼`.,､､､＼`/ |
__＼＼＼＼彡`-`ﾐっ/ L
￣|＼_＼＼_U,~⌒ヾ /
］| ∥￣￣￣￣U~~U / /
__| ∥ □ □ ∥ |/ /
___`∥________∥／_/
￣￣￣￣￣￣￣￣￣∥ /
__________________∥/
410:１３２人目の素数さん
20/04/21 14:53:08 pHehMVs5.net
>>386
sirモデルでは罹患した患者が回復するまでの機会中常に一定の確率(感染率)で遭遇した感受受宿主に感染を広めると仮定して立式してる。
ザックリした傾向を見るならそれで十分だけど、実際のモデルではそんな事はありえない。
感染する確率は患者の体内で繁殖しているウィルス量の増加に従って増えるからその効果を勘案してΓ分布(というかt^c d^x なる形の関数)に応じて感染率が罹患した時点から変化するとするんでしょ？
もちろんコレでも大体の傾向見るにはコレで十分というモデルを臨床例から適当に選んでるだけでしょ？
実際には未来永劫感染力を持ち続けるなんて事があるはずないし、罹患初期のウィルス量の増え方はおそらく指数関数的に増えるものを採用すべきだろうし。
411:１３２人目の素数さん
20/04/21 14:54:06 pHehMVs5.net
mcmc
412:１３２人目の素数さん
20/04/22 08:56:14 0tlpvLKp.net
>>388
無理無理。
カブトガニVSシオマネキ論争はコンピュータの計算結果で決着がついたよ。
413:１３２人目の素数さん
20/04/22 13:23:53 0tlpvLKp.net
週末は休むとか週間変動の影響を除くために１週間の移動平均で線形回帰して片対数
414:グラフにすると https://i.imgur.com/7VwfswD.png 自粛の効果がでてきているな。 多分、検査自粛の効果だろうな。
415:イナ
20/04/22 19:51:00.75 iq1GZOqA.net
∥∩∩ ∥ □ ∥前>>388
（(-_-)∥ ∥______
（っ⌒⌒ﾞ ｡∥╂─╂
■`(_)_)ц~ ∥╂─╂
＼■υυ■_∩∩､＼＼＼
＼＼＼＼⊂(_ _ )`⌒つ､
＼＼＼＼＼＼＼`υ､＼＼＼＼＼＼＼＼＼＼＼＼＼＼＼＼`>>391カブトガニでもシオマネキでもシュクメルリでも微分するのがいちばん強力だと思う。
416:１３２人目の素数さん
20/04/23 04:58:29 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:１３２人目の素数さん
20/04/23 08:28:17 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:１３２人目の素数さん
20/04/23 10:37:41 3vcirHk0.net
中国には特異度１００％の検査キットがあるんだってね。
すべて陰性にでるようにセットされていると。
419:１３２人目の素数さん
20/04/23 13:51:17 3vcirHk0.net
新型コロナ患者を治療している病院に100人の職員がいる。
検体採取器具は5人分、試薬は1回分しかないとする。
無作為抽出した5人の職員から採取した検体を混合して検査したら陽性であった。
職員の陽性者数の期待値を求めよ。
また、50人以上の感染者いる確率はいくつか？
検査の陽性率はハイリスク群に検査している東京の数値2457/6654を使って計算せよ。
URLﾘﾝｸ(i.imgur.com)
420:１３２人目の素数さん
20/04/23 13:53:25 3vcirHk0.net
>>397
こんなグラフになった。
URLﾘﾝｸ(i.imgur.com)
421:１３２人目の素数さん
20/04/23 17:36:48 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:１３２人目の素数さん
20/04/23 22:02:05.89 MtjaFZpr.net
レベル低
423:１３２人目の素数さん
20/04/24 02:55:58 juJsFFfP.net
慶応大の調査で、コロナ以外で来院した人をPCR検査したところ
4/67の確率で要請だった
東京都内1500万のうち何人くらいが感染しているか推定せよ
424:１３２人目の素数さん
20/04/24 03:18:37 XPGerQAq.net
岩田健太郎・�
425:_戸大学教授『東京はすでに20万～400万人感染の可能性』 https://leia.5ch.net/test/read.cgi/poverty/1587664106/
426:１３２人目の素数さん
20/04/24 05:36:50.27 9Fe9PNfV.net
>>401
ニュー速でやったが、ベイズでこんな感じ？(感染者率)
URLﾘﾝｸ(i.imgur.com)
427:１３２人目の素数さん
20/04/24 05:52:14 9Fe9PNfV.net
>>397
混ぜて検査する方法の最適化問題ってのもあるね。
URLﾘﾝｸ(mobile.twitter.com)
(deleted an unsolicited ad)
428:１３２人目の素数さん
20/04/24 08:13:06 0onl6lJy.net
>>401
これって東京の無作為サンプリングではなさそうだよね
病院に来たって事は熱が出てたのかもしれないし
想定できる母集団ってなんになるのだろうか
429:１３２人目の素数さん
20/04/24 12:07:56 v55OWzbu.net
新型コロナ患者を治療している病院に100人の職員がいる。
検体採取器具は10人分、試薬は1回分しかないとする。
無作為抽出した10人の職員から採取した検体を混合して検査したら陰性であった。
職員の陽性者数の期待値を求めよ。
また、50人以上の感染者いる確率はいくつか？
検査の陽性率はハイリスク群に検査している東京の数値2457/6654を使って計算せよ。
430:１３２人目の素数さん
20/04/24 12:11:36 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:１３２人目の素数さん
20/04/24 17:39:04 8oiI190P.net
>>401
まじか。
まあ、体調悪いから病院に行くわけで、バイアスかかってるとはいえ...。
432:１３２人目の素数さん
20/04/24 17:42:10 v55OWzbu.net
>>404
ラテン方陣の問題？
433:１３２人目の素数さん
20/04/24 20:48:11 v55OWzbu.net
>>403
事前分布にJefferey分布を使っているな。
URLﾘﾝｸ(i.imgur.com)
破線が事前分布、実戦が事後分布
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:１３２人目の素数さん
20/04/24 20:55:16 v55OWzbu.net
>>410
青実線が事前分布を一様分布(Beta(1,1))としたとき。
URLﾘﾝｸ(i.imgur.com)
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:１３２人目の素数さん
20/04/24 21:13:10 v55OWzbu.net
>>401
URLﾘﾝｸ(georgebest1969.typepad.jp)
に準じて 東京都民の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:１３２人目の素数さん
20/04/24 21:35:58 v55OWzbu.net
>>402
感度30-70％(最頻値0.5,標準偏差0.2のβ分布),特異度(最頻値0.9 標準偏差0.05のβ分布)に設定。
有病率の事前分布は0-1の一様分布にして
MCMCしてみると
URLﾘﾝｸ(i.imgur.com)
という結果になった。
有病率の信頼区間は広すぎｗ
mean lower upper
0.22346412787 0.00000002913 0.83202346988
438:１３２人目の素数さん
20/04/24 21:50:22 v55OWzbu.net
>>413
有病率の最頻値は
> density(prev)$x[which.max(density(prev)$y)]
[1] 0.019186
439:１３２人目の素数さん
20/04/24 22:01:25 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:１３２人目の素数さん
20/04/24 22:03:05.10 v55OWzbu.net
>>415
この最頻値はオーストリアの0.3％という値に等しいな。
441:１３２人目の素数さん
20/04/24 23:51:07 v55OWzbu.net
>>415
Jeffreys が正しいスペリングみたい。 Jefferey'sかと思っていた。
442:１３２人目の素数さん
20/04/25 13:54:43 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％信頼区間
URLﾘﾝｸ(i.imgur.com)
443:１３２人目の素数さん
20/04/25 17:53:16 MKtk31sc.net
>>409
動的計画法を使うらしい。
URLﾘﾝｸ(qiita.com)
444:１３２人目の素数さん
20/04/26 21:43:36 0lgRXnyr.net
スレチですが統計に詳しい方がいると思って伺いました。
工学系の大学院生で論文を読んでて、見慣れない記号や数式が出てきたので質問させてください。
この数式中のEは何を意味するのでしょうか？
Rの期待値的なものですか？
URLﾘﾝｸ(i.imgur.com)
445:１３２人目の素数さん
20/04/27 02:54:57 cS0ISst+.net
院生でわからないってマジか
専門外の機械学習いきなりやらされたのかな？
446:１３２人目の素数さん
20/04/27 06:45:22 S7AmHM83.net
有病率が0.3％（オーストリアの無作為抽出)や6%(慶応大学の入院患者無作為抽出）のときは
感度30～70% 特異度90～100%のPCR検査キットでは
陰性的中率は
オーストリアの例では 95％[87％～99％]
慶応の例では 94%[86～99％]になる。
すべて陰性にでる中国のイカサマキットなら
陰性的中率は
オーストリアの例では99.7％
慶応の例では94%になる。
有病率が10％くらいになればイカカマキットの方が陰性的中率の成績が劣る。
447:１３２人目の素数さん
20/04/27 07:30:42 S7AmHM83.net
某医療センターでは治療後、2回連続してPCR検査陰性であれば退院させているとする。
PCR検査は
感度0.3~0.7 : β分布(2.625,2.625)相当
特異度 0.95~1.0 : β分布(26.5014,2.3422)相当
として
ｎ回連続して陰性であれば退院とするとどれくらいの感染者が野に放たれるかを
有病率（＝検査前確率）を変化させてグラフにしてみた。
URLﾘﾝｸ(i.imgur.com)
448:１３２人目の素数さん
20/04/28 10:03:02.25 /p2virDs.net
URLﾘﾝｸ(www.researchgate.net)
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:１３２人目の素数さん
20/04/28 21:14:48.22 V+QD5qnX.net
>>394
答えなくね？
451:１３２人目の素数さん
20/04/29 16:59:02 YhsQxAnp.net
>>425
正解。正直者が誰もいないようにプログラムした。
イナ師が微分で解いてくれるのを期待していたんだが。
452:１３２人目の素数さん
20/04/29 17:43:56 YhsQxAnp.net
PCR検査は
感度0.3~0.7 : β分布(2.625,2.625)相当
特異度 0.95~1.0 : β分布(26.5014,2.3422)相当
として、
URLﾘﾝｸ(statmodeling.hatenablog.com)
の設定を踏襲して、利用できるデータを更新して
推定陽性率は
> 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:１３２人目の素数さん
20/04/29 18:02:18.86 FBYarVUq.net
>>426
正直者の確定が誰もいないというより、条件を満たす正直者・嘘つきの組み合わせが存在しないよね？
URLﾘﾝｸ(i.imgur.com)
454:１３２人目の素数さん
20/04/29 18:35:59 uFfpYtab.net
rRT-PRC検査の感度・特異度ともに90%以上です。
455:１３２人目の素数さん
20/04/29 19:40:54.19 uFfpYtab.net
通常の医療でも臨床診断の5%は誤診だと推定されていたりする。
これを実際の数として見積もると膨大になる。
456:１３２人目の素数さん
20/04/29 20:37:12 YhsQxAnp.net
>>429
採取手技や採取時期で変動すると思う。
コンタミによる偽陽性もあるだろうし。
457:１３２人目の素数さん
20/04/29 20:38:43 YhsQxAnp.net
>>428
言っていることは同じ。
全部の条件を満たす組み合わせは存在しない。
458:１３２人目の素数さん
20/04/29 20:39:49 YhsQxAnp.net
>>427
このデータは再現性がないことがわかったので撤回します。
459:１３２人目の素数さん
20/04/30 10:01:57 42TM9o6B.net
ふぉー
＜新型コロナ＞抗体検査５．９％陽性 市中感染の可能性 都内の希望者２００人調査
URLﾘﾝｸ(www.tokyo-np.co.jp)
460:１３２人目の素数さん
20/04/30 12:13:53 qOF+URFa.net
>>434
"検査結果では、一般市民の百四十七人の４・８％にあたる七人が陽性、
医療従事者五十五人のうち９・１％の五人が陽性だった。
URLﾘﾝｸ(www.tokyo-np.co.jp)
"
> 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:１３２人目の素数さん
20/04/30 19:13:48 42TM9o6B.net
>>434
これの興味深いところは、陽性率の高さもさることながら、医療関係者の感染率が一般市民のそれよりも
（多分有意に）高いことだと思う。抗体検査の精度ってほとんどデータが無いに等しいけど、医療関係者の
感染者率が平均よりも多いのならば（実際そうなのだと思う）、このテストはそれを反映したものと言えて、
抗体検査の信頼性をある程度証明できるかもしれない。
そして、市中の陽性率の高さも。
この辺、ベイズで上手く検証できないかな。
462:１３２人目の素数さん
20/04/30 20:32:35.24 qOF+URFa.net
>>436
>医療関係者の感染率が一般市民のそれよりも（多分有意に）高い
χ二乗検定では、p-value = 0.41 で 標本数が少ないから有意差がでない。
エントリーが５以下のときには信頼するなと習ったな。
まあ、ポアソンでもp-value = 0.33
>435にRの出力を貼っておいた。
463:１３２人目の素数さん
20/04/30 20:35:20.83 rbz2947p.net
>>436
55人中5人と147人中7人じゃ、有意差ないだろ。
と書こうと思ったら、 >>437氏に先を越されたｗｗｗ
464:１３２人目の素数さん
20/04/30 20:41:39 rbz2947p.net
抗体検査キットの評価をしたら特異度はみな5/5だったらしいけど、サンプルが
５つだけって意味だとすると、せいぜい90％以上くらいのことしか言えないな。
特異度95%だとしたら、５％が陽性っていわれてもねぇｗ
465:１３２人目の素数さん
20/04/30 21:20:01 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)
}
結果は、
URLﾘﾝｸ(i.imgur.com)
有病率の期待値は2.3%、最頻値は0.31% 少数データなので信頼区間が広い。
有病率が平均値50%の一様分布というのは現実離れした分布だとは思う。
466:１３２人目の素数さん
20/04/30 21:22:46 qOF+URFa.net
>>439
すべてが陰性に出るイカサマキットなら特異度が100％
>18参照
467:１３２人目の素数さん
20/04/30 21:33:35 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:１３２人目の素数さん
20/04/30 21:47:02.48 qOF+URFa.net
>>440（追記)
>有病率が平均値50%の一様分布というのは現実離れした分布だとは思う。
オーストリアのデータ0.3%を参考に有病率の事前分布が0.2-0.4％
β(10.865, 3279.34)に相当として、MCMCしてみた。
感度と特異度の事前分布は前回と同じ。
URLﾘﾝｸ(i.imgur.com)
さほど、有病率が高いという結論は引き出せないな。
470:１３２人目の素数さん
20/04/30 22:23:36 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:１３２人目の素数さん
20/04/30 23:17:24 qOF+URFa.net
β分布と乱数発生で比較。
一般人と医療従事者の確率分布に従う乱数を１０００万個発生させて比率の分布をグラフにすると
９５％信頼区間が１を挟むから有意差なし。比率が１以上となる確率も87.5％で95%を越えないから有意差なし。
URLﾘﾝｸ(i.imgur.com)
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:１３２人目の素数さん
20/05/01 00:36:14 LcUyF6Tc.net
>>441
すべてが陰性に出るなら陽性は０人だろｗ
473:１３２人目の素数さん
20/05/01 07:22:53 LJVZP1pw.net
>>446
偽陽性率０％だから特異度は１００％
474:１３２人目の素数さん
20/05/01 10:05:01 LcUyF6Tc.net
>>447
だから、陽性が何人か出てんだから感度0%のインチキじゃなかろうってこと。
475:１３２人目の素数さん
20/05/04 00:17:48 Q8iLGwO2.net
最近の報道
インドは、中国企業から購入した中共ウイルス（新型コロナウイルス）の迅速スクリーニング検査キットの精度がわずか5％だとして、約50万個の注文をキャンセルした。
....
神戸市中央区にある市立医療センター中央市民病院の医師などのグループは、ことし３月末から先月７日にかけて、新型コロナウイルス以外の理由で外来を受診した患者から無作為に1000人を選び、血液中に新型コロナウイルスに感染したあとにできる「抗体」があるか調べました。
グループによりますとその結果、3.3％にあたる33人から抗体が検出されたということです。
以上を知ったある会社が
必ず陽性がでる試薬33個と必ず陰性がでる試薬967個を混ぜた1000試薬をセットにしてインドに売り込んだ。
問題
1試薬は1回しか検査できないとして
これがイカサマキット�
476:ﾅあることを証明する手段はあるか？
477:１３２人目の素数さん
20/05/04 15:22:29 jDRWX2Ph.net
3月の宿題で(1)のみ正解の数弱@shukudai_sujaku
昨年度の大学への数学(大数)での勝率は、
学コンBコースが 1/1 = 100% ，
宿題が 3/10 = 30% でした！
宿題の勝率が低すぎると思うので、
これからは一層精進していきたいです！
URLﾘﾝｸ(twitter.com)
(deleted an unsolicited ad)
478:１３２人目の素数さん
20/05/04 22:45:28 94tg6i4k.net
高橋洋一(統計数理研究所→大蔵省)
スウェーデンの感染症対策を解説
URLﾘﾝｸ(youtu.be)
479:１３２人目の素数さん
20/05/05 06:07:12 ht6rG86e.net
URLﾘﾝｸ(toyokeizai.net)
のデータを使って
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-特異度)の二項分布に従うとする。
事後確率分布は
URLﾘﾝｸ(i.imgur.com)
陽性率P=14895/153581=0.09698465での事後確率分布の密度比(Savage-Dickey density ratio)でベイズファクターを出すと
> d1/d2 # Savage-Dickey densiti ratio = BF12
[1] 1.007722
まあ、ちょっぴり、感度30-70％のモデルの方がいいかも、という結果。
陽性数/検査数の時系列データでもあればもう少し差がでるかもしれん。
東京都のデータで計算させようかと思ったが、東京都は検査人数を隠蔽しているので使いものにならない。
480:１３２人目の素数さん
20/05/05 06:26:36.20 FpoKo6a+.net
訂正
陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする
↓
陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする
481:１３２人目の素数さん
20/05/05 11:13:29 b2IqdVzK.net
3月の宿題で(1)のみ正解の数弱@shukudai_sujaku
昨年度の大学への数学(大数)での勝率は、
学コンBコースが 1/1 = 100% ，
宿題が 3/10 = 30% でした！
宿題の勝率が低すぎると思うので、
これからは一層精進していきたいです！
URLﾘﾝｸ(twitter.com)
(deleted an unsolicited ad)
482:１３２人目の素数さん
20/05/05 23:00:17 wvtmzVA1.net
>>449
抗体のはいったサンプルを２，３０個試してみりゃわかんじゃないの？
いくらなんでも感度が３％じゃつかえねぇ。
483:１３２人目の素数さん
20/05/06 15:29:26 RG00+xls.net
イカサマキットの感度特異度の事前分布を一様分布に設定して
抗体のはいったサンプルを２０個で全部陰性であったので３０個試したら全部陰性であったとすると
感度・特異度の事後分布は
URLﾘﾝｸ(i.imgur.com)
484:１３２人目の素数さん
20/05/06 15:32:34 RG00+xls.net
ところが、有病率33/1000であったときに無作為に２０人および３０人を選んで全部陰性であっても
感度・特異度の事後分布は
URLﾘﾝｸ(i.imgur.com)
485:１３２人目の素数さん
20/05/06 15:52:13 RG00+xls.net
>>454
事後分布を一様分布に設定して、このコピペ小僧の正解率の期待値・最頻値・95％信頼区間を求めよ。
その結果、
URLﾘﾝｸ(i.imgur.com)
486:１３２人目の素数さん
20/05/06 15:58:11 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:１３２人目の素数さん
20/05/06 16:25:45 RG00+xls.net
>>454
このコピペ小僧の正解率に学コンと宿題で差があるかを検定せよ。
事前分布にJefferysを用いた結果、
URLﾘﾝｸ(i.imgur.com)
488:１３２人目の素数さん
20/05/06 17:08:09 tNZVV9ZH.net
>>456
感度が5%以下じゃつかいもんにならんわなぁ。インチキで確定。
489:１３２人目の素数さん
20/05/06 23:23:13 wABsYddm.net
URLﾘﾝｸ(www.buzzfeed.com)
大阪市立大学が一昨年の血液を使って検査したところ、偽陽性は1人もいなかったそうだ。
490:１３２人目の素数さん
20/05/07 00:59:54 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:１３２人目の素数さん
20/05/07 02:09:04.11 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
URLﾘﾝｸ(ja.wikipedia.org)
492:１３２人目の素数さん
20/05/07 02:30:17 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:１３２人目の素数さん
20/05/07 17:10:06 CHL0/p02.net
>>462
５０人で０でしょ？
だったら、特異度として保証されるのはせいぜい９８％程度じゃん。
１％の陽性率に対して、それでは不十分。
494:１３２人目の素数さん
20/05/07 18:13:41.63 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:１３２人目の素数さん
20/05/07 18:52:47.08 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:１３２人目の素数さん
20/05/07 20:13:12 VnQvkZ57.net
URLﾘﾝｸ(www.buzzfeed.com)
のデータを使って、不明なものは一様分布（ベータ分布の形状母数(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)
}
その結果
URLﾘﾝｸ(i.imgur.com)
497:１３２人目の素数さん
20/05/07 21:49:37 CHL0/p02.net
>>469
なにやってるかよくわからんので、見当外れの指摘かもしれんが、
大阪市大の3/312と神戸大の33/1000ってのは特異度も感度も異なる
であろう別種のキットによる検査結果なんで、一緒くたにしちゃ
まずいんでないの？
498:１３２人目の素数さん
20/05/07 22:36:59.75 VnQvkZ57.net
>>470
同一キットじゃないから、ご指摘の通り。
しかも神戸大の方では陰性検体での確認はされていないから、神戸大方の陽性率が高いのは偽陽性を含む可能性もあるね。
大阪市大だけのデータでやってみると。
URLﾘﾝｸ(i.imgur.com)
499:１３２人目の素数さん
20/05/07 22:46:09.05 VnQvkZ57.net
>>471
使える情報が少なすぎて信頼区間が広すぎ。
有病率:0-87%
キットの感度:0-70%
ってコイントス変わらんな。
500:１３２人目の素数さん
20/05/07 22:48:00.40 VnQvkZ57.net
結局のところ、断定的な結論は出せないねということを数字で確認しているだけｗ
501:１３２人目の素数さん
20/05/08 02:04:52 Ugc87SUM.net
>>464
>>465
有難うございます。
二項分布下の母比率の信頼区間推定問題なんですね。
A２：正規分布近似法でなくWilsonの信頼区間法で解くから
上側区間≠下側区間となり平均値まわり左右非対称信頼区間と
なるのですね。
502:１３２人目の素数さん
20/05/08 09:15:37 WmDpVhCu.net
3月の宿題で(1)のみ正解の数弱@shukudai_sujaku
昨年度の大学への数学(大数)での勝率は、
学コンBコースが 1/1 = 100% ，
宿題が 3/10 = 30% でした！
宿題の勝率が低すぎると思うので、
これからは一層精進していきたいです！
URLﾘﾝｸ(twitter.com)
(deleted an unsolicited ad)
503:１３２人目の素数さん
20/05/08 11:37:38 CfFk/Uw1.net
>>474
Rのlibrary binomを使って binom::confint(3,312)で >464の出力が得られる
信頼区間をグラフにすると
URLﾘﾝｸ(i.imgur.com)
点線は3/312
正規分布近似のasymptotic以外は非対称。
どれを使うべきか？ 好きなのを使えばいい、と思う。
但し、値が負になったり１を越えたりするのは不採用の方が賢明だとは思う。
Wilson法は値が０や１に近くても信頼できるという人がいるけど、どうやって検証するのかはよくわからん。
504:１３２人目の素数さん
20/05/08 12:41:19 CfFk/Uw1.net
>>476
binom::confintはbinom.confint(3,312)の間違い
library("binom")
binom.confint(3,312)
505:１３２人目の素数さん
20/05/09 13:22:10 oDT7bFgO.net
レムデシビルで初のRCTが出たんだけど、
URLﾘﾝｸ(www.thelancet.com)(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:１３２人目の素数さん
20/05/09 13:22:52 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:１３２人目の素数さん
20/05/09 13:25:29 oDT7bFgO.net
確率分布を図示すると
URLﾘﾝｸ(i.imgur.com)
509:１３２人目の素数さん
20/05/09 14:25:05 74hNX8Dr.net
100例もやって差がでないのなら、はっきり言って効き目なしとかわらん。
むしろ副作用がこわい。
アビガンはどうなんだろうね。なぜか国内のまとまったデータがまだ出てこない。
510:１３２人目の素数さん
20/05/09 14:26:34 74hNX8Dr.net
まあ、軽症者が８割とかだと、早期投与の効き目を判定するには相当の
人数に使わないとわかんないだろうけど。
511:１３２人目の素数さん
20/05/09 15:22:40 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:１３２人目の素数さん
20/05/09 17:15:11 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:１３２人目の素数さん
20/05/09 20:41:59 74hNX8Dr.net
あ、 >>482はレムデシじゃなくてアビガンについての話だかんね。
レムデシは重症者向けらしいけど、たぶん駄目だろ。
515:１３２人目の素数さん
20/05/09 21:35:34 Fh34v+Vz.net
【政府が隠す３．１１】 地震波形が、核実験と同じ
ｽﾚﾘﾝｸ(earth板)
sssp://o.5ch.net/1nmad.png
516:１３２人目の素数さん
20/05/11 00:31:07 mGn48zmS.net
5/10夜『NHKスペシャル』「新型コロナウイルス 出口戦略は」で
厚労省新型コロナ対策班西浦博北大教授と氏が分析中のノートPC画面が
映されていた。氏の使用解析ソフトウェア名は何というのでしょうか？
517:１３２人目の素数さん
20/05/11 01:36:00 P7wrhagU.net
西浦先生はRじゃないかな。
今週ぐらいに、実効再生産数を算出するRのコードを公開できるかもと言ってた。
火曜にはニコ生があるから、そこでたっぷりと聞けるはず。
URLﾘﾝｸ(sp.live.nicovideo.jp)
518:１３２人目の素数さん
20/05/11 03:07:28 6TI8KNr5.net
>>488
>11に数式があるけど
519:１３２人目の素数さん
20/05/11 03:58:21.10 v4rWY6J/.net
>>489
いや、与えられた何日分かの新規患者数から再生産数とか回復率とかを推定するための数式じゃないの？
520:１３２人目の素数さん
20/05/11 06:53:45 ruo6tAnf.net
>>490
SEIRより単純なReed-Frostモデルだよな
URLﾘﾝｸ(en.m.wikipedia.org)
521:１３２人目の素数さん
20/05/11 09:14:46.53 P7wrhagU.net
>>489
自分で計算してみた人ならわかるけど、日々の再生算数を出すにはA「感染してから他人にうつすまでの日数の分布」が必要。西浦先生はこれをまず求めていて、その結果は割合と世界中で使われている。
加えてB「感染してから発表の数字に載るまでの日数の分布」も必要で、発表の数字をBで逆畳み込み積分してC「感染推定日毎の数字」に変換し、CをAで畳み込み積分したものでCを割れば日々のRtが求められる。
とは言うものの、そのままやると誤差で数字が発散するし、逆畳み込みは一意ではないのでどう推定するかは色々と考えなくてはならない。多分日本独自の事情もコードに入れなきゃいけない。大変だろうなあとは思う。
522:１３２人目の素数さん
20/05/11 11:02:15.57 +2HA8Yn9.net
>>492
レスありがとうございます。
確かにSEIRモデルやSIRIモデルだと
「感染してから発表の数字に載るまでの日数の分布」って考慮されてないね。
523:１３２人目の素数さん
20/05/12 20:50:54 8HqG++pl.net
>>488
公開されたぞー。StanとRだ。
URLﾘﾝｸ(nbviewer.jupyter.org)
524:１３２人目の素数さん
20/05/12 22:32:45.27 9+i4T0w3.net
ウイルス感染のモデルも必要だと思うけど
社
525:会経済活動の影響や 自粛解除した後の死者増加や 自粛継続した場合の法人の倒産や 治療薬やワクチンができた後の社会の経済や文化面の回復 そこまで考慮してどういう政策を選択したら 損失を最小化できるかという問題は解ける？ 制約付きの最小化問題だと思う 数値化するのが難しい考慮要素もあるからその辺は仮の値とかで
526:１３２人目の素数さん
20/05/13 17:23:13 PBQsFM+W.net
>>495
モデルを精密複雑化すればするほど、必要なパラメータが増えて
それを弱情報事前分布にすると事後分布の信頼区間が広がるんだなぁと思う。
527:１３２人目の素数さん
20/05/14 08:00:23 kec+XbRE.net
>>494
これにある、
JapaneseDataCOVID19
ってどこかからダウンロードできるんだろうか？
528:１３２人目の素数さん
20/05/14 09:54:38 IN0uokww.net
>>497
URLﾘﾝｸ(github.com)
529:１３２人目の素数さん
20/05/14 11:47:52 kec+XbRE.net
>>498
ありがとうございます。
530:１３２人目の素数さん
20/05/14 11:51:11 kec+XbRE.net
親ディレクトリ（フォルダ）探せばよかったんだね。
URLﾘﾝｸ(nbviewer.jupyter.org)
531:１３２人目の素数さん
20/05/14 13:08:10.58 gLQVRit5.net
>>495 >>496
パラメータを増やせば増やすほど、どんな仮説も否定できなくなる。
実験系の学界でよく使う手口だよね。
おや、誰か来たようだ…
532:１３２人目の素数さん
20/05/14 13:58:09 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:１３２人目の素数さん
20/05/14 14:10:15 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:１３２人目の素数さん
20/05/14 16:35:35 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の分布はこんな感じ、
URLﾘﾝｸ(i.imgur.com)
URLﾘﾝｸ(i.imgur.com)
535:１３２人目の素数さん
20/05/14 17:00:26 1Kd4DHQt.net
>>501
否定できないけど実証もできないことになるのでは？
単純なモデルだと現実の問題を解決できない事もあるし
536:１３２人目の素数さん
20/05/14 17:23:23.34 kec+XbRE.net
>>494
モデルで再生産数の事前分布は 平均2.5 標準偏差2.0の正規分布に設定されていたので
平均と標準偏差を変化させて、再生産数の事後分布を描出してみた。
かなり、事前分布の影響を受けるみたい。
URLﾘﾝｸ(i.imgur.com)
537:１３２人目の素数さん
20/05/14 17:26:03.64 kec+XbRE.net
どうも、こういう境地だなぁ。
断定的な結論は出せないということを数字で確認しているだけｗ
538:１３２人目の素数さん
20/05/14 17:30:34 50GljeDM.net
基本再生算数はウイルスの特徴で決まるもので
実行再生算数は更に環境や人の行動の影響で変化するものだと理解している
実行再生算数を減らすようにするには
どう行動したらいいか
かつ経済活動を出来るだけ下げずに
539:１３２人目の素数さん
20/05/14 18:25:32 IN0uokww.net
西浦モデル、一度感染日ごとの人数を最尤法で確定して、そこからベイズを回してRtを求めてるんだよなあ。
本来は、日々の感染者数も確率的に揺れがあるはず。だからRtの誤差幅は発表よりも多く見積もるべきかもしれん。
540:１３２人目の素数さん
20/05/14 18:39:49 IN0uokww.net
コードを解説してる人
URLﾘﾝｸ(mikuhatsune.hatenadiary.com)
541:１３２人目の素数さん
20/05/14 19:09:54 kec+XbRE.net
再生産数の事前分布を色々かえて事後分布を出してみた。
URLﾘﾝｸ(i.imgur.com)
542:１３２人目の素数さん
20/05/15 08:46:23 qjCTzgxb.net
>>504
切断分布だと収束しないみたいなので、
一様分布[0,10]に近そうな正規分布[5,3]
URLﾘﾝｸ(i.imgur.com)
を事前分布にして走らせてみた。
URLﾘﾝｸ(i.imgur.com)
543:１３２人目の素数さん
20/05/15 09:56:20.25 qjCTzgxb.net
【新型コロナ】 東京0.6％、東北6県0.4％陽性・・・抗体検査1000人実施 [影のたけし軍団★]
ｽﾚﾘﾝｸ(newsplus板)
複数の検査キットの性能評価と感染状況の確認が目的でしたが、東京都で献血した500人のうち3人、
東北6県で献血した500人のうち2人がいずれかの検査キットで陽性と判定されました。
満員電車など人との接触の多い東京とそうでない東北で陽性率に有意差はあるか？
544:１３２人目の素数さん
20/05/15 12:23:12 GeDvuMme.net
>>513
検査キットの特異度が99.5%以上あることが検証されました。
ってだけの話じゃね？
まぁ、東北６県と東京が同じってことは、特異度はそれより
高くはないってことだろうね。
545:１３２人目の素数さん
20/05/15 12:26:48 GeDvuMme.net
言い換えると、真の陽性率はわかりません、ってこと。
546:１３２人目の素数さん
20/05/15 16:53:27 qjCTzgxb.net
>>513
陽性率の確率分布を一様分布にすると事後分布は
URLﾘﾝｸ(i.imgur.com)
なるけど、重なりの部分の面積が差がないことの度合いを示していると考えていいかな？
547:１３２人目の素数さん
20/05/15 18:43:40 qjCTzgxb.net
これは、味気がないな。
> Epi::twoby2(x)
2 by 2 table analysis:
------------------------------------------------------
Outcome : Col 1
Comparing : Row 1 vs. Row 2
Col 1 Col 2 P(Col 1) 95% conf. interval
Row 1 3 497 0.006 0.0019 0.0184
Row 2 2 498 0.004 0.0010 0.0158
95% conf. interval
Relative Risk: 1.5000 0.2517 8.9384
Sample Odds R
548:atio: 1.5030 0.2501 9.0339 Conditional MLE Odds Ratio: 1.5024 0.1713 18.0536 Probability difference: 0.0020 -0.0092 0.0139 Exact P-value: 1.0000 Asymptotic P-value: 0.6561 ------------------------------------------------------
549:１３２人目の素数さん
20/05/16 07:59:05 /d9XIHEO.net
再生産数を計算するRのプログラムあったんだな
URLﾘﾝｸ(www.rdocumentation.org)
550:１３２人目の素数さん
20/05/16 08:16:33 28hwAVWs.net
>>509
後者の分析、RStan本の訳者が今やっているようだ。
URLﾘﾝｸ(twitter.com)
(deleted an unsolicited ad)
551:１３２人目の素数さん
20/05/16 08:17:55 28hwAVWs.net
>>519
URLﾘﾝｸ(twitter.com)
(deleted an unsolicited ad)
552:１３２人目の素数さん
20/05/16 08:36:22 /d9XIHEO.net
>>520
アヒル本の著者だね。俺も読んだ。
553:１３２人目の素数さん
20/05/16 12:36:46 cGFgpNwX.net
4/2の時点で感染者数を６０００くらいと見積もってるね。＞アヒル本の人
実数の三倍程度。いい線かもしれない。
554:１３２人目の素数さん
20/05/16 13:39:32 XSllT5Kn.net
URLﾘﾝｸ(github.com)
555:１３２人目の素数さん
20/05/16 14:29:29 BxGcLzV+.net
int<lower = K> upper_bound;
↓
int upper_bound;
にしてもエラーがでる。
stan_dataで
upper_bound = 147 にすると動くけど、何をやってんのか自分でもよくわからん。
556:１３２人目の素数さん
20/05/16 16:47:14 UBjBTXVe.net
アヒル本ってなんですか？
557:１３２人目の素数さん
20/05/16 17:55:00 BxGcLzV+.net
>>525
名著として名高いStanの入門書
StanとRでベイズ統計モデリング (Wonderful R) (日本語) 単行本 ? 2016/10/25
表紙の色がアヒルのくちばしの色に似ているかららしい。
558:１３２人目の素数さん
20/05/16 17:56:04 BxGcLzV+.net
>>524
自己解決
select
↓
dplyr::select
でちゃんと動作した
559:１３２人目の素数さん
20/05/16 19:16:41 5/7mBSAW.net
>>526
thx
560:１３２人目の素数さん
20/05/16 19:28:52 BxGcLzV+.net
R に EpiEstimというパッケージがあって、再生産数を算出する関数が搭載されている。
結局、infecterとinfecteeが発症するまでの期間serial intervalの分布をどう設定するかで結果が変わるみたいだなぁ。
Rのヘルプファイルを解読中。
Rのヘルプファイルは不親切設計で有名（理解できている人の備忘録みたいな性格だから）。
561:１３２人目の素数さん
20/05/16 20:49:50 5T7nYzW3.net
ここに居る人達って何者なの
学部の知識超えてるよね
統計でご飯食べてる人たち？
562:１３２人目の素数さん
20/05/16 21:03:16 YPW1s8+S.net
しかし、なんでも揃ってるなRって
563:１３２人目の素数さん
20/05/16 22:57:43 BxGcLzV+.net
>>529
Stanでの西浦モデルではinfecterとinfecteeが発症するまでの期間serial intervalの分布に
## Serial interval [Nishiura et al 2020 - only certain cases]
param1_SI = 2.305,
param2_SI = 5.452,
// serial interval
vector[K] gt = pweibull(param1_SI, param2_SI, K);
として使われているので、平均値などを出してみた。
乱数発生と理論値
> x=rweibull(1e5,param1_SI,param2_SI)
> mean(x) ; param2_SI*gamma(1+1/param1_SI)
[1] 4.829273
[1] 4.830129
> var(x) ; param2_SI^2*(gamma(1+2/param1_SI)-(gamma(1+1/param1_SI))^2)
[1] 4.907755
[1] 4.940682
> median(x) ; param2_SI*(log(2)^(1/param1_SI))
[1] 4.655777
[1] 4.6505
> density(x)$x[which.
564:max(density(x)$y)] ; param2_SI*(1-1/param1_SI)^(1/param1_SI) [1] 4.116837 [1] 4.259624 > optimise(function(x) dweibull(x,param1_SI,param2_SI),c(0,10),maximum = T)$max [1] 4.259623 グラフにしてみた。 https://i.imgur.com/9vvCJuZ.png 正規分布で近似してもよさそうな感じだな。
565:１３２人目の素数さん
20/05/17 00:04:44 u5AEq3c8.net
>>532
ワイブル分布の平均 4.830129と標準偏差2.222765をそのまま正規分布のパラメータに使って、グラフを重ねてみる。
URLﾘﾝｸ(i.imgur.com)
ワイブル分布で発生させた乱数をワイブルでフィットさせてAICを出してみた
Goodness-of-fit criteria
1-mle-weibull
Akaike's Information Criterion 438377.2
Bayesian Information Criterion 438396.2
ワイブル分布で発生させた乱数を正規分布でフィットさせてAICを出してみた。
Goodness-of-fit criteria
1-mle-norm
Akaike's Information Criterion 444280.9
Bayesian Information Criterion 444299.9
まぁ、許容範囲。
これで、
library(EpiEstim)の例にある、 mean_si std_siが求まった
## Estimate R with assumptions on serial interval
res <- estimate_R(incid, method = "parametric_si",
config = make_config(list(
mean_si = 4.83, std_si = 2.22)))
domestic , imported, unobserved の分類がよくわからんが、全部足してグラフを描いてみた。
URLﾘﾝｸ(i.imgur.com)
566:１３２人目の素数さん
20/05/17 00:18:59 u5AEq3c8.net
別の論文だと対数正規分布がフィットすると西浦氏は記載している。
serila interval : infector と infecteeの発症間隔
URLﾘﾝｸ(www.ijidonline.com)(20)30119-3/pdf
その分布が平均4.7 標準偏差2.9の対数正規分布が最もフィットするのはいいんだが、
その分布を与えるパラメータの記述がほしい。
最小二乗法で求めてみた。
$par
[1] 1.3862713 0.5679836
ワイブル分布にも似るとか書いてあるがパラメータ記載なし
この対数正規分布をワイブル分布で近似してみる。
Fitting of the distribution ' weibull ' by maximum likelihood
Parameters:
estimate Std. Error
shape 1.757488 0.00392072
scale 5.316986 0.01014664
URLﾘﾝｸ(i.imgur.com)
点線が2項分布で実線がワイブル分布
567:１３２人目の素数さん
20/05/17 07:09:05 /ApxcPC4.net
“頭脳王”東大生・河野玄斗
基本的な数学でコロナウイルス検査を全員にしても意味がないことを証明してみた
URLﾘﾝｸ(www.youtube.com)
568:１３２人目の素数さん
20/05/17 07:28:25 mpH4uuPx.net
実効再生産数を感染者数の推移から推定する数理的原理をキチンと解説した本、または論文誰か知りませんか？
勉強してみたい。
569:１３２人目の素数さん
20/05/17 07:33:35 u5AEq3c8.net
>>536
RのEpiEstimの著書の論文なんかどうでしょう？
A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics
URLﾘﾝｸ(www.ncbi.nlm.nih.gov)
570:１３２人目の素数さん
20/05/17 07:55:41 u5AEq3c8.net
Rで引数なしで関数を実行させようとすると、ソースが表示される。
その関数が呼び出した関数のソースを次々に辿っていけば内部計算がわかるので、そこから原理がつかめるかも（俺には無理なのでしないけど）
内部関数のときは:が3つとか、パッケージ名:::関数、パッケージ名:::関数.default（例、ｔ検定のソース表示は stats:::t.test.default ）で表示される。
EpiEstim::estimate_R
EpiEstim:::process_si_data
EpiEstim:::process_config_si_from_data
coarseDataTools::dic.fit.mcmc
571:１３２人目の素数さん
20/05/17 08:10:41 u5AEq3c8.net
>>529
感染者に０が続くと再生産数の信頼区間幅がどんどん広くなってくる。
まあ、疫病用のソフトウェアと理解しておこう。
URLﾘﾝｸ(i.imgur.com)
infected=c(0,1,1,1,0,0,0,2,0,3,2,3,1,1,1,1,3,0,1,0,0,0,0,0,1,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
## Estimate R with assumptions on serial interval
res <- estimate_R(infected, method = "parametric_si",
config = make_config(list(
mean_si = 4.83, std_si = 2.22)))
572:１３２人目の素数さん
20/05/17 09:00:46 u5AEq3c8.net
>>539
ちゃんと、記載されていた ：P
the precision of these estimates
depends directly on the number of incident cases in the time
window [t ? τ + 1; t]. This allows us to control the precision
by adjusting the window size.
URLﾘﾝｸ(www.ncbi.nlm.nih.gov)
573:１３２人目の素数さん
20/05/17 09:34:29 u5AEq3c8.net
基本コンセプトはこれだろうな。
Therefore, in practice, we apply our method to data consisting of daily counts of onset of symptoms where the infectivity profile ws is approximated by the distribution of the serial interval.
公表された感染者数で計算するために発症から診断/公表までの時間も考慮した点が西浦モデルの優れた点ではないかと感じている。
574:１３２人目の素数さん
20/05/17 12:07:27 S4iLXC97.net
>>537
ありがとう。
読んでみます。
でも"a new frame work"っていうならもっと大本のこの道の研究者なら「Rtの推定はコレ」みたいな標準的なものがあるんですかね？
どなたかご存知ないですか？
575:１３２人目の素数さん
20/05/17 12:11:01 S4iLXC97.net
>>540
おぉ、こっちがRで採用されてる推定法なんですね。
読んでみます。
576:１３２人目の素数さん
20/05/17 12:16:27 S4iLXC97.net
と思ったら>>537と>>540は同じかorz
577:１３２人目の素数さん
20/05/17 15:31:13 u5AEq3c8.net
結局、これが核心部分
Rt =I t (number of new infections generated at time t) / Σ[s=1,t] I t-s * Ws ( = total infectiousness of infected individuals at time t)
Ws : an infectivity profile given by a probability distribution ws, dependent on time since infection of the case, s, but independent of calendar time, t.
E[I t] = Rt Σ[s=1,t] I t-s * Ws
Σ[s=1,t] I t-s * Wsの部分が畳み込み積分で
Ws ∝ serial interval ガンマ分布で近似するのが定石らしい。
578:１３２人目の素数さん
20/05/17 20:41:31 /XIaXEJI.net
>>545
それはそのモデルでのRtの定義ですね、
ではなく例えば4/1～4/30までのデータx1,x2,‥,x30までが与えられた時、各日付のCIをどう定義してるのかがわからないんです。
統計量Xが一つあってその観測値xが一つある時そのレベルλのレベルλのCIがIであるとは
P(X<x|θ)>(1-λ)/2 & P(X>x|θ)>(1-λ)/2 ⇔ θ∈I
(θがIに入らないときはxが小さすぎてそんな観測値が得られる確率が(1-λ)/2以下になるか、xが大きすぎてそんな事が起こる確率が(1-λ)/2以下になってほとんど起こり得ない)
となりますが、統計量が複数になるとこの大きすぎて、小さすぎてと二つのハズレ領域を考えるだけでは済まなくなります。
ハズレ領域の設定の仕方は色々考えられるけど統計の問題なので自分で俺様流の領域を設定するわけにもいかないのでpublicにはどう処理してるのだろうと。
しかし疫学の教科書わざわざ買うほどには興味もないしw
知っ
579:てる人いたら教えてもらおうと。 まだ論文読んでないので書いてあるかもですけど。
580:１３２人目の素数さん
20/05/18 05:56:32 hVD/4gbI.net
>>546
>523氏が挙げてくれた
実効再生産数とその周辺 に 記述があったが、publicと呼べるのかどうかは門外漢なので知らない。
URLﾘﾝｸ(github.com)
こんな記述があった
>>
Several non mathematical definitions.
A:
あるカレンダー時刻 t で起こっている 2 次感染数の 1 人あたり平均値
B:
あるカレンダー時刻 t で感染した 1 人がその後の経過で生み出す 2 次感染者数の平均値
C:
罹患率 有病割合比などから推定される 1 人あたりの 2 次感染者数（ actual reproduction number とか）
D:
予防接種など流行対策下での 1 人当たりが生み出す 2次感染者数
<<
と定義は一義的ではないみたいだが、
西浦モデルでのR1(t) R2(t)は >540の論文では Rt Rct （ｃはcohortの頭文字)として言及されているから、まあ、理論疫学者の間ではcommonな定義なんだろうと思う。
581:１３２人目の素数さん
20/05/18 06:16:59 hVD/4gbI.net
>>546
>4/1～4/30までのデータx1,x2,‥,x30
x1,...,x30が感染者数なら非負整数で与えられるから、CIを考える必要があるかな？
集計ミスで19人が117人とかあったらしいから、信頼区間を考える必要があるのかもしれないとは思うけど。
確かに、発症日に関してはファジーな部分があるとは思う。
いつから熱がありましたか？味がわからなくなりましたか？と問われても１～３日位の幅はでるだろう。
RのEpiEstim::estimate_Rをヘルプファイルの実例を使って走らせてみた。
serial intervalの分布をデータから推定させるのに
infecterの発症日の下限・上限
infecteeの発症日の下限・上限
を設定する項目が(EL,ER,SL,SR)があった。
このデータに合致する分布関数(ガンマ、ワイブル、対数正規分布が指定可能）のパラメータを算出して計算させているみたい。
西浦モデルではワイブル分布で固定。
潜伏期間にも変動があるから、誰がinfecterで誰がinfecteeかを決定するのも難しいだろうとは思う。
後から発症した人間が感染源というのもありうるし。
582:１３２人目の素数さん
20/05/18 11:15:46 iMUOaQs/.net
>>548
統計データから真の罹患日を推定する方法もあるようですが
そこではないんです。
しりたいのはCIのハズレ領域の設定です。
1変数の場合、母数θに対して分布Fθが定まっている場合、レベルλに対して[0,1]の部分集合J(λ)が決まって、観測値xに対する信頼区間I(λ,x)は
θ∈I(λ,x)⇔Fθ(x)∈J(λ)
を満たす区間として定まります。
上の方でI(λ)が上下対称に取らないのはなぜという話題がありましたが、コレがその理由です。
J(λ)の方は上下の(1-λ)/2を削って((1-λ)/2,1-(1-λ)/2)をとり、上下対称に“ハズレ領域”をとりますが、それをもとに計算されるI(λ,x)は対称とはならないからです。
問題は観測値が2変数以上ある場合“ハズレ領域”をどう設定するものかわからないのです。
私が大学で勉強した時はそこまでやらなかったので。
普通に考えればI^3の中の立方体の体積がλになるように真ん中にとるんだろうなぁと思うんですけど。
583:１３２人目の素数さん
20/05/18 13:21:08 1P7V5xJn.net
職場でも最初に発症した人が感染源のように扱われるけど
潜伏期間の分布を考えたら断定はできない。
COVID19の潜伏期間の論文
URLﾘﾝｸ(www.nejm.org)
から、
潜伏期間は対数正規分布で近似できてそのパラメータは
#--- incubation period ---
# from Li et al NEJM 2020
# log
584:normal mean = 5.2 ln.par1 = 1.434065 ln.par2 = 0.6612 という。 ある人物Xが新型コロナ肺炎に罹患したとする。 行動調査によって発症前にキャバクラに行っており接客したキャバ嬢がX発症の２日後に発症していたことがわかった。 Xがキャバ嬢から移された確率を求めよ。
585:１３２人目の素数さん
20/05/18 13:25:20 1P7V5xJn.net
>>549
Highest Density Probability Intervalを求めればいいんじゃないの？
586:１３２人目の素数さん
20/05/18 14:04:29 jxOnYm2h.net
>>551
何ですかそれは？
587:１３２人目の素数さん
20/05/18 14:08:09 1P7V5xJn.net
>>550
正確にはキャバ嬢がXより先に感染していた確率だな。
588:１３２人目の素数さん
20/05/18 14:08:59 1P7V5xJn.net
>>552
確率分布が対称じゃないときの信頼区間
589:１３２人目の素数さん
20/05/18 14:13:34 1P7V5xJn.net
>>554
こんな感じ。URLﾘﾝｸ(i.imgur.com)
590:１３２人目の素数さん
20/05/18 14:26:48 jxOnYm2h.net
>>554
ただ、それだとそもそもモードに近いとこをやってます。
信頼区間は密度関数を横に切るのではなく両裾を縦に切ってハズレ部分が1-λになるようにするので少しイメージが違うしきがします。
モードなのかメジアンなのかの違いです。
いずれにせよ、こうやればいいという拡張のための俺様ルールを設定するのはいくらでもできますが、統計の話なのでそんな俺様ルールについて語っても意味ありません。
591:１３２人目の素数さん
20/05/18 15:34:14 1P7V5xJn.net
>>556
単峰性の場合、信頼区間幅が最小になるのがHighest Density Interval
>550なら
HDI
> HDInterval::hdi(x)[1:2]
lower upper
0.5822687 12.5635525
分位数
> quantile(x,c(0.025,0.975))
2.5% 97.5%
1.148711 15.334698
HDIの方が幅が小さい。
592:１３２人目の素数さん
20/05/18 15:36:39 jxOnYm2h.net
>>557
？？？
593:１３２人目の素数さん
20/05/18 15:44:16.03 jxOnYm2h.net
ああ、わかった。HDIやCIの意味を誤解してませんか？
HDIでググって調べたらコレ↓ですよ。
URLﾘﾝｸ(rindalog.blogspot.com)
594:１３２人目の素数さん
20/05/18 19:32:28.41 4BbElJo8.net
>>179
単にお勉強ができただけ。
頭が良くないのさ。
自分の頭で物事を考えるってことができない。
595:１３２人目の素数さん
20/05/18 19:34:11.31 4BbElJo8.net
>>181
その通り。
具体的には理学部の数学科と物理学科。
工学部にも時々もの凄いのがいる。
596:１３２人目の素数さん
20/05/18 19:49:54 1P7V5xJn.net
>>559
いや、Rのパッケージ HDIntervalで計算しているから誤解していないと思う。
内部の処理コードをみると信頼区間幅が最小になるのを最尤法で出しているね。
597:１３２人目の素数さん
20/05/18 19:51:25 1P7V5xJn.net
>>561
ほんとうに頭のいい人は医学部でなく理学部か工学部に行く。
ほんとうに頭の悪い奴は底辺シリツ医大に行く。
シリツ医大には手先の器用な人もいるが、頭が器用な奴をみたことがない。
598:１３２人目の素数さん
20/05/18 20:05:42 1P7V5xJn.net
>>550(自答)
#
# 人物ｄが発症してdelay日後に濃厚接触したキャバ嬢ｃが発症
# ｃの感染がｄより先行していた確率は？
rm(list=ls())
stancode=
"
data{
real onset_delay;
real ln_par1;
real ln_par2;
}
parameters{
real <lower=0> d_incubation;
real <lower=0> c_incubation;
}
transformed parameters{
real infection_delay = onset_delay + d_incubation - c_incubation;
}
model{
d_incubation ~ lognormal(ln_par1,ln_par2);
c_incubation ~ lognormal(ln_par1,ln_par2);
}
"
model=stan_model(model_code = stancode)
#--- incubation period ---
# from Li et al NEJM 2020
# lognormal mean = 5.2
ln_par1 = 1.434065
ln_par2 = 0.6612
fn.stan <- function(delay, print=FALSE, ...){
dataList=list(onset_delay=delay,ln_par1=ln_par1,ln_par2=ln_par2)
fit=sampling(model,data=dataList, ...)
ms=rstan::extract(fit)
if(print) BEST::plotPost(ms$infection_delay,compVal=0,xlab='infection delay')
mean(ms$infection_delay < 0)
}
fn.stan(2,print=T,iter=5000,warmup=1000)
onset_delays=0:20
y=sapply(onset_delays,fn.stan)
plot(onset_delays,y, ylab='Pr[ Infected Later ])')
２日後の発症だと
> fn.stan(2,print=T,iter=5000,warmup=1000)
[1] 0.2945
３割くらいはあとから発症した方が先に感染していた可能性がある。
599:１３２人目の素数さん
20/05/18 20:48:18.33 1P7V5xJn.net
Temporal dynamics in viral shedding and transmissibility of COVID-19
URLﾘﾝｸ(www.nature.com)
のRのコード
URLﾘﾝｸ(github.com)
と
西浦モデルのコード
URLﾘﾝｸ(nbviewer.jupyter.org)
から発症間時間(serial interval)の分布を重ねてみた。
URLﾘﾝｸ(i.imgur.com)
600:１３２人目の素数さん
20/05/18 21:01:41.77 4BbElJo8.net
>>563
一方的で申し訳ないが私立大医学部は金持ちのバカ息子が行くイメージ。
601:１３２人目の素数さん
20/05/18 21:02:26.04 4BbElJo8.net
西浦さんはさんざん適当なことを言って世論を煽ってどう責任を取るのかな？
602:１３２人目の素数さん
20/05/18 21:18:52.42 AArRB0Ix.net
>>567
少なくとも西浦氏は算出コードを公開しているだけでも好感が持てる。
603:１３２人目の素数さん
20/05/18 21:20:19.41 AArRB0Ix.net
>>566
それは正しい認識。
凄いのはド底辺シリツ医の馬鹿さ加減だよ。
裏口バカと呼ばれるのがよくわかる。
URLﾘﾝｸ(imagizer.imageshack.com)
URLﾘﾝｸ(i.imgur.com)
馬鹿だという自覚がないので救いようがない。
ICU Bookの最終章の冒頭で著者がこう書いている。
In clinical matters, ignorance can be dangerous,
but ignorance of ignorance can be fatal.
「叱られないと勉強しない」の対偶を「勉強すると叱られる」
と答えるのはignorance can be dangerousの範疇だが、
ドヤ顔で
>対偶をとれば意味が逆になる例文。
というのは、まさに
ignorance of ignorance can be fatal.
604:１３２人目の素数さん
20/05/21 11:39:47 hAZkHjNF.net
西浦モデルの再生算数の事前分布を変化させてグラフにしてみた。
西浦モデルでのデフォルト
URLﾘﾝｸ(i.imgur.com)
事前分布を一様分布(0,10)近似の正規分布で近似させた場合
URLﾘﾝｸ(i.imgur.com)
再生算数の平均０、標準偏差１の場合
URLﾘﾝｸ(i.imgur.com)
印象としては、西浦モデルは頑強性robustnessのあるモデルとは言えない。
事前分布に大きく影響されるモデルだと思う。
605:１３２人目の素数さん
20/05/21 11:42:59 hAZkHjNF.net
（url訂正)
西浦モデルの再生算数の事前分布を変化させてグラフにしてみた。
西浦モデルでのデフォルト
URLﾘﾝｸ(i.imgur.com)
事前分布を一様分布(0,10)近似の正規分布で近似させた場合
URLﾘﾝｸ(i.imgur.com)
再生算数の平均０、標準偏差１の場合
URLﾘﾝｸ(i.imgur.com)
印象としては、西浦モデルは頑強性robustnessのあるモデルとは言えない。
事前分布に大きく影響されるモデルだと思う。
606:１３２人目の素数さん
20/05/21 11:59:00 hAZkHjNF.net
>>571
誤解されるのは不本意なの追記するけど、ソースを公開する西浦先生の姿勢は高く評価している。
隠蔽・改竄・破棄する安倍とは大違い。
607:１３２人目の素数さん
20/05/21 12:04:41 RwIzsMl5.net
だからベイズ推定が統計学の世界でメジャーにならんのだろ？
論理的根拠のない“事前分布”なるもので答えがひょいひょい変わるのでは社会的な影響が大きい防疫政策の決定には使えない。
普通の統計学の検定なら理論的に根拠のある数字、推論しか使わない。
計算は大変だけど。
608:１３２人目の素数さん
20/05/21 12:36:51 hAZkHjNF.net
>>573
成人の平均身長を１～２ｍに事前分布にするのは納得できるし、
生まれる子供が男子である
609:確率は０．４～０．６というのも俺は納得できる。 PCRの感度が３０～７０％として計算するのも納得できるからその設定で階層ベイズモデルを組むことには異論はないな。
610:１３２人目の素数さん
20/05/21 12:43:16 hAZkHjNF.net
こういう問題
あるタクシー会社のタクシーには１から通し番号がふられている。
タクシー会社の規模から保有タクシー台数は１００台以下とわかっている（弱情報事前分布）。
この会社のタクシーを５台みかけた。最大の番号が６０であった。
この会社の保有するタクシー台数の期待値と９５％信用区間(信頼区間)を求めよ。
をベイズで解くときは、
６０台～１００台である確率を一様分布として処理しているから
これに異論があるのは理解できるけど
日本人成人の平均身長を推定に１～２ｍを事前分布に想定するのには俺は異論はないね。
一様分布ではなくてガンマ分布にすべきだというのは議論になるとは思うけど。
611:１３２人目の素数さん
20/05/21 13:00:03 hAZkHjNF.net
>>573
ベイズ信奉者から、ベイズ論者を採用したGAFA(Google等）が成長したと教わった。
迷惑メールのフィルタリングとか雑音の除去とか日常的に役立っているというぞ。
612:１３２人目の素数さん
20/05/21 13:04:49 hAZkHjNF.net
>普通の統計学の検定なら理論的に根拠のある数字、推論しか使わない。
>計算は大変だけど。
普通の統計学こそ、無理やり既知の分布に当てはめようとするんだよね。
MCMCの方の方が応用が広い。
613:１３２人目の素数さん
20/05/21 13:12:59.01 hAZkHjNF.net
プロ野球選手の打率は？と問われたら選手次第で異なる、と誰でもわかるのに
確定不能の平均値が存在していると妄想して計算を始めるのが古典主義統計。
つまり、値は存在するけど確定できないという信仰の世界。
昨今の新コロナでいえば、PCRの検査の感度・特異度が一定と考えるのが古典主義。
プロ野球選手の打率と同じでそんなのは場面場面で変化するよ、と考えて計算するのがベイズ。
614:１３２人目の素数さん
20/05/21 13:16:04.83 EdpRqUGW.net
>>576
もちろん推定の有力な方法であるにせよ、元の仮定に何の根拠もないわけだからそれから得られる結論には論理的根拠はない、ないが、数学的に伝統的な手法で与えるた結論と大差ない事がなんらかの保証があるなら、有力になる。
それが論理的に“大差ない結論が得られる”事が示されてるなら単なる計算手法に過ぎないし、示されていなくても経験的に“よい結果ぎ得られる事が多い”ならそのジャンルではそこそこ信頼するに値するんだろう。
しかしなんらかのモデルでは答えが一意に定まらず、事前分布の選び方により大きく答えが違ってしまう場合があっても不思議はないし、そのような場合ではやはり、“ではどう計算するのが正しいかのか”の論証を待たなければ信頼するのは危険になる。
615:１３２人目の素数さん
20/05/21 13:19:12.39 hAZkHjNF.net
古典的（頻度主義）統計信者って、この計算はどうやるんだ？
俺は乱数発生させて計算できるけど、
そればベイズなのかどうかは知らんが、条件付き確率なのでベイズなんだろうな。
（開業医スレに投稿したけど、回答できるやつは０）
COVID19の潜伏期間の論文
URLﾘﾝｸ(www.nejm.org)
結論は
#--- incubation period ---
# from Li et al NEJM 2020
# lognormal mean = 5.2
ln.par1 = 1.434065
ln.par2 = 0.6612
ある開業医が新型コロナ肺炎に罹患したとする。
行動調査によって発症前にキャバクラに行っており接客したキャバ嬢が開業医発症の２日後に発症していたことがわかった。
キャバ嬢は開業医から移されたと主張して１億円の賠償を求めている。
潜伏期間には幅がありキャバ嬢から移された可能性もあると主張してその確率を計算して賠償金を値切りたい。
いくら値切れるか計算せよ。
616:１３２人目の素数さん
20/05/21 13:31:57 3Y6HuNza.net
>>580
潜伏期間なんて数学的に決定できるハズないやん？
数学がやるのは例えば感染日と実際
617:発症した日が確定できるような症例がある程度以上あって、それが従うと病理学的に信頼できる分布の族があって、その中からデータに最も“沿う”分布を選び出すことしかできない。 それは何件も実際の患者のウイルス量の統計をとったり、ウイルスが体内でどのように増えていくのかの病理学的研究データがあって初めて可能になる。
618:１３２人目の素数さん
20/05/21 17:20:41 Du+rgnHD.net
>>571
よく理解できていないので質問ですけど
事前分布とは具体的に何の分布ですか？
基本再生算数の推定値の分布？
実行再生算数の推定値の分布？
実行再生算数の事前分布は基本再生算数の分布としたらいいのかなと思いますけど
619:１３２人目の素数さん
20/05/21 17:26:10 BjgAHZWs.net
まぁこのスレは用語がめちゃくちゃだからなぁ。
620:１３２人目の素数さん
20/05/21 17:58:56 hAZkHjNF.net
モデル前提での計算できないアホ発見！
621:１３２人目の素数さん
20/05/21 18:22:49 hAZkHjNF.net
潜伏期は病原体と罹患者のパワーバランスで決まるだろうから
定数でなくてばらつきはあると思うね。
622:１３２人目の素数さん
20/05/21 18:25:39 8+K/mYXQ.net
>>572
>隠蔽・改竄・破棄する安倍とは大違い。
安倍ちゃんがやってるわけじゃないだろ。
指示もされてねーのに、官僚がやってんだよ。
そんなこともわからんようでは、あんたの解析もあてにならんな。
623:１３２人目の素数さん
20/05/21 18:26:26 hAZkHjNF.net
>>582
西浦のソースだと Rt ~ normal(2.4 ,2)
624:１３２人目の素数さん
20/05/21 18:28:21.45 hAZkHjNF.net
>>586
安倍じゃなければ官僚もまともだったんじゃねえの？
安倍らの集団に組み込まれると東大卒も馬鹿になるようだぞ。
625:１３２人目の素数さん
20/05/21 18:31:24.22 8+K/mYXQ.net
>>581
その通り。
結局、そういう本質的なデータや理論抜きでは、ベイズ推定やったって
限界があるし、結果の説得力もない。
まあ、適用限界ってものがあるのは何やっても同じだけど。
626:１３２人目の素数さん
20/05/21 18:32:55.80 8+K/mYXQ.net
>>588
>安倍らの集団に組み込まれると東大卒も馬鹿になるようだぞ。
統計的に確認してみれば？ベイズ推定でもっともらしい結果が出るんじゃね？ｗ
627:１３２人目の素数さん
20/05/21 18:33:59.08 hAZkHjNF.net
>>590
今は2/2かなｗ
加藤と西村
628:１３２人目の素数さん
20/05/21 18:40:02 hAZkHjNF.net
>>590
馬鹿になる事後確率分布はこんな結果になりました。
URLﾘﾝｸ(i.imgur.com)
629:１３２人目の素数さん
20/05/21 18:41:26 hAZkHjNF.net
95% CI は 0.5005265 1
630:１３２人目の素数さん
20/05/21 18:46:18 hAZkHjNF.net
>>592
ちなみに事前分布は一様分布に設定。
631:１３２人目の素数さん
20/05/21 18:53:16.25 hAZkHjNF.net
結局、統計ってこういう道楽なんだよなぁ。
632:１３２人目の素数さん
20/05/21 19:01:50.90 hAZkHjNF.net
>>581
いやウイルスの振る舞いも素粒子の振る舞いもばらつきがあるんじゃないの？
存在も位置も確率的にしかわからない、という議論になると思う。
633:１３２人目の素数さん
20/05/21 21:14:17 hAZkHjNF.net
>>571
最初の２つでは頑強と言えそう。
平均身長のたとえだと
事前分布に1～２ｍを選んでも１～１０ｍを選んでも同様の結果だが、０～１ｍを選ぶと現実離れした結果が返ってくるという感じだな。
634:１３２人目の素数さん
20/05/21 21:27:14 hAZkHjNF.net
次のおもちゃ
しばらく、これで遊べそう。
臨床所見からロジスティック回帰でCOVID19の確率を出すペーパーがでるだろうなと思っていた。
Real-time tracking of self-reported symptoms to predict potential COVID-19
URLﾘﾝｸ(www.nature.com)
635:１３２人目の素数さん
20/05/22 01:26:24 MSJJjK3u.net
>>591
馬鹿の定義は？単なるお前の主観じゃねーかｗ
ベイズ統計にのめり込むと馬鹿になる事後確率分布でも求めてろ。
636:１３２人目の素数さん
20/05/22 04:53:19 C0tPqEF8.net
こういうのも興味ある人多い？感染してからの日数とPCR陰性になる確率の関係。
URLﾘﾝｸ(twitter.com)
(deleted an unsolicited ad)
637:１３２人目の素数さん
20/05/22 05:39:46 s9txG4+B.net
各国のロックダウンの度合いを数値化してるところ。色々と分析に使えるかも。
URLﾘﾝｸ(ourworldindata.org)
638:１３２人目の素数さん
20/05/22 08:35:38 DQesskhT.net
>>599
あれ3/3になったのか。
それだと、95％CrIは0.5559329 1
事前分布にはJeffereysで計算
639:１３２人目の素数さん
20/05/22 09:06:15 DQesskhT.net
>>600
ありがとうございます。
おねーちゃんと濃厚接触したあとは、何日目に検査したらいいかの参考になりました。
640:１３２人目の素数さん
20/05/22 09:49:33 DQesskhT.net
>>600
論文を検索してみた、この論文からのグラフのよう。
Variation in False-Negative Rate of Reverse Transcriptase PolymeraseChain Reaction?Based SARS-CoV-2 Tests by Time Since Exposure
URLﾘﾝｸ(www.acpjournals.org)
641:１３２人目の素数さん
20/05/22 10:18:18 DQesskhT.net
Measurements: A Bayesian hierarchical model was fitted to estimate
the false-negative rate by day since exposure and symptom
onset.
とのことで、
こんな、ありふれたベイズ階層モデル
x[j,t] ~ Binomial(n[j,t],p[j,t])
logit(p[j,t]) = βj + β*1log(t) + β2*log(t)^2 + β3*(t)^3
βj ~ Normal(β0,σ2)
642:１３２人目の素数さん
20/05/22 10:36:46 DQesskhT.net
Rとstanのコードが公表されているので、遊べそう。
URLﾘﾝｸ(github.com)
643:１３２人目の素数さん
20/05/22 17:38:12 MSJJjK3u.net
>>602
1/1だよ。あんただけｗ
644:１３２人目の素数さん
20/05/22 19:01:13 DQesskhT.net
>>607
信頼区間出せない馬鹿発見！
645:１３２人目の素数さん
20/05/22 19:09:59 xdbfxKAO.net
このスレででてる信頼区間なんかほとんど統計の検定試験では出てこない俺様信頼区間だけどな。
信頼区間の定義ちゃんと言える奴の方が少ないだろ。
646:１３２人目の素数さん
20/05/22 19:11:40 DQesskhT.net
1/1でも信頼区間を出せちゃうのがベイズ
647:１３２人目の素数さん
20/05/22 19:13:01 DQesskhT.net
>>609
ベイズでは信用区間と呼んで信頼区間と区別する人もいるね。
CIでなくてCrIという記載も目にする。
648:１３２人目の素数さん
20/05/22 19:20:12 DQesskhT.net
ベイズではconfidence interval でなくてcredibility interval
649:１３２人目の素数さん
20/05/22 23:34:07 MSJJjK3u.net
>>612
あんたの独占スレになったな。1/1達成、おめでとうｗ
650:１３２人目の素数さん
20/05/23 06:44:49 COZ69QMb.net
bootstrap使うと大抵の数値の信頼区間は出せるけど1/1には無理だな。
651:１３２人目の素数さん
20/05/23 15:08:38 COZ69QMb.net
JAGSを使って1/1のときの95％ Credibility Intervalを求めるスクリプト
事前分布はJeffereys
library(rjags)
data=list(x=1,n=1)
cat('
model{
x ~ dbin(p,n) # binomial distribution
p ~ dbeta(0.5,0.5) # Jeffereys prior
}',file='tmp.txt')
jagsModel=jags.model('tmp.txt',data=data,n.chains=4)
update(jagsModel)
codaSamples=coda.samples(jagsModel,n.iter=1e5,var='p')
gelman.plot(codaSamples) # 収束を確認
js=as.data.frame(as.matrix(codaSamples))
HDInterval::hdi(js$p) # 95% Highest Densty Interval
# 既存のパッケージで確認。
HDInterval::hdi(qbeta,shape1=0.5+1,shape2=0.5)
binom::binom.bayes(1,1)
652:１３２人目の素数さん
20/05/23 17:34:11 COZ69QMb.net
エクセルのツールをCDCが公開しているね
653:。 https://www.cdc.gov/coronavirus/2019-ncov/hcp/COVIDSurge.html COVID-19Surge is a spreadsheet-based tool that hospital administrators and public health officials can use to estimate the surge in demand for hospital-based services during the COVID-19 pandemic. A user of COVID-19Surge can produce estimates of the number of COVID-19 patients that need to be hospitalized, the number requiring ICU care, and the number requiring ventilator support. The user can then compare those estimates with hospital capacity, using either existing capacity or estimates of expanded capacity. With COVID-19Surge, users define the population in the hospital “catchment area” or local jurisdiction. The user also enters the number of cases to date, and the available hospital resources (non-ICU beds, ICU beds, and mechanical ventilators). Users can assess up to three community mitigation strategies simultaneously and compare the impacts on hospital resources versus a “no intervention” scenario.