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

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

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

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

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

こうした推計をすると、今の感染者は氷山の一角だと思われるが、今後の増加ペースはどうなるだろうか。
新型コロナウイルスの検査は簡単に行えるので、今後、日本での感染者数は増えていくだろう。
ある時点ではそれがネズミ算的に増えるかのように思える局面もあるだろうが、
筆者の推計が正しければ、現時点ではせいぜい数百人が一つのメドだ。
0399１３２人目の素数さん
垢版 |
2020/04/23(木) 17:36:48.57ID:3vcirHk0
>>397
これであってるかな？

> # 期待値
> integrate(function(x) x*pdf(x),0,100)$value
[1] 37.13
> # 50人以上の確率
> integrate(pdf,50,100)$value
[1] 0.0041903
> c(HPDI.lower=lwr,HPDI.upper=upr) # HPDI
HPDI.lower HPDI.upper
27.778 46.558
0401１３２人目の素数さん
垢版 |
2020/04/24(金) 02:55:58.61ID:juJsFFfP
慶応大の調査で、コロナ以外で来院した人をPCR検査したところ
4/67の確率で要請だった
東京都内1500万のうち何人くらいが感染しているか推定せよ
0405１３２人目の素数さん
垢版 |
2020/04/24(金) 08:13:06.02ID:0onl6lJy
>>401
これって東京の無作為サンプリングではなさそうだよね
病院に来たって事は熱が出てたのかもしれないし
想定できる母集団ってなんになるのだろうか
0406１３２人目の素数さん
垢版 |
2020/04/24(金) 12:07:56.92ID:v55OWzbu
新型コロナ患者を治療している病院に100人の職員がいる。
検体採取器具は10人分、試薬は1回分しかないとする。
無作為抽出した10人の職員から採取した検体を混合して検査したら陰性であった。
職員の陽性者数の期待値を求めよ。
また、50人以上の感染者いる確率はいくつか？
検査の陽性率はハイリスク群に検査している東京の数値2457/6654を使って計算せよ。
0407１３２人目の素数さん
垢版 |
2020/04/24(金) 12:11:36.28ID:v55OWzbu
>>401
95%CIで
> binom::binom.confint(4,67)
method x n mean lower upper
1 agresti-coull 4 67 0.05970149 0.019131154 0.1480232
2 asymptotic 4 67 0.05970149 0.002968439 0.1164345
3 bayes 4 67 0.06617647 0.015026904 0.1253507
4 cloglog 4 67 0.05970149 0.019283398 0.1337560
5 exact 4 67 0.05970149 0.016504404 0.1458632
6 logit 4 67 0.05970149 0.022588780 0.1485238
7 probit 4 67 0.05970149 0.020905075 0.1402573
8 profile 4 67 0.05970149 0.018970462 0.1332788
9 lrt 4 67 0.05970149 0.018929939 0.1332756
10 prop.test 4 67 0.05970149 0.019297952 0.1534709
11 wilson 4 67 0.05970149 0.023459351 0.1436950
0408１３２人目の素数さん
垢版 |
2020/04/24(金) 17:39:04.40ID:8oiI190P
>>401
まじか。
まあ、体調悪いから病院に行くわけで、バイアスかかってるとはいえ...。
0411１３２人目の素数さん
垢版 |
2020/04/24(金) 20:55:16.38ID:v55OWzbu
>>410
青実線が事前分布を一様分布(Beta(1,1))としたとき。

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

Jeffereyの方が95%CI幅が小さいな。

> binom::binom.bayes(4,67,prior.shape1 = 0.5,prior.shape2 = 0.5)
method x n shape1 shape2 mean lower upper sig
1 bayes 4 67 4.5 63.5 0.06617647 0.0150269 0.1253507 0.04999999
> binom::binom.bayes(4,67,prior.shape1 = 1, prior.shape2 = 1)
method x n shape1 shape2 mean lower upper sig
1 bayes 4 67 5 64 0.07246377 0.01876916 0.1338218 0.04999999
0412１３２人目の素数さん
垢版 |
2020/04/24(金) 21:13:10.41ID:v55OWzbu
>>401

https://georgebest1969.typepad.jp/blog/2020/04/%E6%85%B6%E5%BF%9C%E3%81%AEpcr6%E3%81%AE%E6%84%8F%E5%91%B3.html

に準じて　東京都民の1395万人に当てはめると

> data.frame(method=ci[,1],round(ci[,4:6]*pop))
method mean lower upper
1 agresti-coull 832836 266880 2064924
2 asymptotic 832836 41410 1624262
3 bayes 923162 209625 1748642
4 cloglog 832836 269003 1865896
5 exact 832836 230236 2034792
6 logit 832836 315113 2071907
7 probit 832836 291626 1956590
8 profile 832836 264638 1859240
9 lrt 832836 264073 1859195
10 prop.test 832836 269206 2140919
11 wilson 832836 327258 2004545

岩田の計算は
5 exact 832836 230236 2034792
0413１３２人目の素数さん
垢版 |
2020/04/24(金) 21:35:58.21ID:v55OWzbu
>>402
感度30-70％(最頻値0.5,標準偏差0.2のβ分布),特異度(最頻値0.9 標準偏差0.05のβ分布)に設定。
有病率の事前分布は0-1の一様分布にして　

MCMCしてみると

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

という結果になった。

有病率の信頼区間は広すぎｗ

mean lower upper
0.22346412787 0.00000002913 0.83202346988
0415１３２人目の素数さん
垢版 |
2020/04/24(金) 22:01:25.57ID:v55OWzbu
>>413
事前分布をJeffereyにしたら、

> js=PCRj4(67,4,SEN=0.5,SD1=0.2,SPC=0.9,SD2=0.1,N.ITER=1e6)$js
mean lower upper
1.4972e-01 6.3120e-14 8.5733e-01

> summary(prev)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.0111 0.0472 0.1497 0.1442 1.0000

> density(prev)$x[which.max(density(prev)$y)] # mode
[1] 0.0032278
0418１３２人目の素数さん
垢版 |
2020/04/25(土) 13:54:43.50ID:R/UD6QQG
某医療センターでは治療後、2回連続してPCR検査陰性であれば退院させているとする。
PCR検査は
感度0.3~0.7 : β分布(2.625,2.625)相当
特異度 0.95~1.0 : β分布(26.5014,2.3422)相当
有病率は一様分布として
2回連続してPCR検査陰性の患者の有病率の期待値を求めよ。

２回連続検査陰性の患者が感染者である確率と有病率との関係をグラフにしてみた。
灰色点線は95％信頼区間

https://i.imgur.com/2U5HBmb.png
0420１３２人目の素数さん
垢版 |
2020/04/26(日) 21:43:36.22ID:0lgRXnyr
スレチですが統計に詳しい方がいると思って伺いました。
工学系の大学院生で論文を読んでて、見慣れない記号や数式が出てきたので質問させてください。
この数式中のEは何を意味するのでしょうか？
Rの期待値的なものですか？

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

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

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

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

有病率が10％くらいになればイカカマキットの方が陰性的中率の成績が劣る。
0423１３２人目の素数さん
垢版 |
2020/04/27(月) 07:30:42.68ID:S7AmHM83
某医療センターでは治療後、2回連続してPCR検査陰性であれば退院させているとする。
PCR検査は
感度0.3~0.7 : β分布(2.625,2.625)相当
特異度 0.95~1.0 : β分布(26.5014,2.3422)相当
として

ｎ回連続して陰性であれば退院とするとどれくらいの感染者が野に放たれるかを
有病率（＝検査前確率）を変化させてグラフにしてみた。


https://i.imgur.com/fvNvXhX.png
0424１３２人目の素数さん
垢版 |
2020/04/28(火) 10:03:02.25ID:/p2virDs
https://www.researchgate.net/publication/340869568_Saliva_Sample_as_a_Non-Invasive_Specimen_for_the_Diagnosis_of_Coronavirus_Disease-2019_COVID-19_a_Cross-Sectional_Study/link/5ea19006458515ec3aff8f36/download
https://i.imgur.com/3GWmZz3.png
に有病率の低い群で行った唾液と鼻咽頭スワッブ200例のcross-section tableがあったので
これでKappa係数とPABAK（prevalence ad-justed bias adjusted kappa）とその信頼区間をstanのMCMCで出してみた。
swabと唾液でまあ、結果が合致している。

> print(fit.kappa,pars=c('kappa','PABAK'))
Inference for Stan model: kappa.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
kappa 0.81 0 0.07 0.65 0.77 0.81 0.86 0.92 3540 1
PABAK 0.93 0 0.02 0.88 0.92 0.93 0.95 0.97 3342 1

JAGSでMCMCしてもほぼ同じ結果。
https://i.imgur.com/XKT53Fy.png
0426１３２人目の素数さん
垢版 |
2020/04/29(水) 16:59:02.34ID:YhsQxAnp
>>425
正解。正直者が誰もいないようにプログラムした。
イナ師が微分で解いてくれるのを期待していたんだが。
0427１３２人目の素数さん
垢版 |
2020/04/29(水) 17:43:56.54ID:YhsQxAnp
PCR検査は
感度0.3~0.7 : β分布(2.625,2.625)相当
特異度 0.95~1.0 : β分布(26.5014,2.3422)相当
として、
http://statmodeling.hatenablog.com/entry/covid19-estimate-total-number-of-positives-in-japan
の設定を踏襲して、利用できるデータを更新して

推定陽性率は
> summary(p.pos) ; HDInterval::hdi(p.pos)[1:2]
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000031 0.000067 0.004681 0.005229 0.009325 0.016658
lower upper
0.00003082 0.01408833

推定有病率は
> summary(prevalence) ; HDInterval::hdi(prevalence)[1:2]
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00000 0.00166 0.00529 0.01040 0.01129 0.16078
lower upper
0.0000005693 0.0305640211
0429１３２人目の素数さん
垢版 |
2020/04/29(水) 18:35:59.17ID:uFfpYtab
rRT-PRC検査の感度・特異度ともに90%以上です。
0430１３２人目の素数さん
垢版 |
2020/04/29(水) 19:40:54.19ID:uFfpYtab
通常の医療でも臨床診断の5%は誤診だと推定されていたりする。
これを実際の数として見積もると膨大になる。
0435１３２人目の素数さん
垢版 |
2020/04/30(木) 12:13:53.58ID:qOF+URFa
>>434

"検査結果では、一般市民の百四十七人の４・８％にあたる七人が陽性、
医療従事者五十五人のうち９・１％の五人が陽性だった。
https://www.tokyo-np.co.jp/s/article/2020043090070748.html

"

> r1=5;r2=7;n1=55;n2=147
> prop.test(c(r1,r2),c(n1,n2))

2-sample test for equality of proportions with continuity
correction

data: c(r1, r2) out of c(n1, n2)
X-squared = 0.679, df = 1, p-value = 0.41
alternative hypothesis: two.sided
95 percent confidence interval:
-0.052613 0.139194
sample estimates:
prop 1 prop 2
0.090909 0.047619

Warning message:
In prop.test(c(r1, r2), c(n1, n2)) :
Chi-squared approximation may be incorrect
> poisson.test(c(r1,r2),c(n1,n2))

Comparison of Poisson rates

data: c(r1, r2) time base: c(n1, n2)
count1 = 5, expected count1 = 3.27, p-value = 0.33
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
0.47778 6.98763
sample estimates:
rate ratio
1.9091
0436１３２人目の素数さん
垢版 |
2020/04/30(木) 19:13:48.86ID:42TM9o6B
>>434
これの興味深いところは、陽性率の高さもさることながら、医療関係者の感染率が一般市民のそれよりも
（多分有意に）高いことだと思う。抗体検査の精度ってほとんどデータが無いに等しいけど、医療関係者の
感染者率が平均よりも多いのならば（実際そうなのだと思う）、このテストはそれを反映したものと言えて、
抗体検査の信頼性をある程度証明できるかもしれない。
そして、市中の陽性率の高さも。

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

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

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

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

と書こうと思ったら、 >>437氏に先を越されたｗｗｗ
0439１３２人目の素数さん
垢版 |
2020/04/30(木) 20:41:39.46ID:rbz2947p
抗体検査キットの評価をしたら特異度はみな5/5だったらしいけど、サンプルが
５つだけって意味だとすると、せいぜい90％以上くらいのことしか言えないな。
特異度95%だとしたら、５％が陽性っていわれてもねぇｗ
0440１３２人目の素数さん
垢版 |
2020/04/30(木) 21:20:01.97ID:qOF+URFa
>>436
事前分布にかなり影響をうけるが、

感度特異度とも50-70％(最頻値60%標準偏差10％のβ分布),
有病率は一様分布、
検査陽性数は陽性確率が　有病率*感度＋（１-有病率)*(1-特異度)の二項分布に従う

というモデルでプログラムを組むと

model
{
for(i in 1:N) {
x[i] ~ dbin(p,n[i])
}
p = prev*sen + (1-prev)*(1-spc)
PPV=sen*prev/(sen*prev+(1-prev)*(1-spc))
NPV=(1-prev)*spc/((1-prev)*spc+prev*(1-sen))
precision=(prev*sen+(1-prev)*spc)/
((prev*sen+(1-prev)*spc + (1-prev)*(1-spc)+(prev*(1-sen))))
pLR=sen/(1-spc)
nLR=(1-sen)/spc
DOR=pLR/nLR
sen ~ dbeta(sn[1],sn[2])
spc ~ dbeta(sp[1],sp[2])
prev ~ dbeta(shape1,shape2)
}

結果は、
https://i.imgur.com/eKXLUXZ.png

有病率の期待値は2.3%、最頻値は0.31%　少数データなので信頼区間が広い。
有病率が平均値50%の一様分布というのは現実離れした分布だとは思う。
0442１３２人目の素数さん
垢版 |
2020/04/30(木) 21:33:35.19ID:qOF+URFa
>55人中5人と147人中7人じゃ、有意差ないだろ。

その比率のまま、約4倍(3.87倍)になればχ二乗検定で有意差がでるね。
r1=5;r2=7;n1=55;n2=147
mat=matrix(c(r1,r2,n1,n2),2,b=T)
fn <- function(x) chisq.test(mat*x)$p.value
fn=Vectorize(fn)
uniroot(function(x) fn(x)-0.05,c(1,10))$root
fn(4)

> fn(4)
[1] 0.045678
0443１３２人目の素数さん
垢版 |
2020/04/30(木) 21:47:02.48ID:qOF+URFa
>>440（追記)
>有病率が平均値50%の一様分布というのは現実離れした分布だとは思う。
オーストリアのデータ0.3%を参考に有病率の事前分布が0.2-0.4％
β(10.865, 3279.34)に相当として、MCMCしてみた。
感度と特異度の事前分布は前回と同じ。

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

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

A群55人中5人とB群147人中7人に有意差はあるか？

モンテカルロ法で検証

A群もB群も同じ陽性率で、0.0594と仮定する
∵ 陽性率は、 12/(55+147) = 12/202 = 0.0594

モンテカルロ法100万回したら
B群の陽性者が7名となったのは、125238回
その内 A群の陽性が5名以上は、28293回

∴ P(55人中5人以上｜147人中7人) = 28293/125238 = 0.226
∴ 有意差なし

EXCEL VBAソースコード概略

Dim I As Long
Dim J As Long
Dim K As Long

K = 1
For I = 1 To 1000000
A = 0
For J = 1 To 55
If Rnd(1) < 0.0594 Then '陽性
A = A + 1
End If
Next

B = 0
For J = 1 To 147
If Rnd(1) < 0.0594 Then '陽性
B = B + 1
End If
Next

If B = 7 Then 'B群の陽性者が7名の場合
Cells(K, "A") = A 'A群の陽性者の人数
K = K + 1
End If
B = 0

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

一般人と医療従事者の確率分布に従う乱数を１０００万個発生させて比率の分布をグラフにすると
９５％信頼区間が１を挟むから有意差なし。比率が１以上となる確率も87.5％で95%を越えないから有意差なし。
https://i.imgur.com/8zFKODI.png

a=0.5 ; b=0.5
r1=5 ; r2=7 ; n1=55 ; n2=147
layout(1)
layout(matrix(1:2,2))
curve(dbeta(x,a+r2,b+n2-r2),0,0.3,bty='l',ann=F,lwd=2)
curve(dbeta(x,a+r1,b+n1-r1),col=2,add=T,lwd=2)
legend('center',bty='n',legend=c('general','medical'),lwd=2,col=1:2)
k=1e7
general=rbeta(k,a+r2,b+n2-r2)
medical=rbeta(k,a+r1,b+n1-r1)
BEST::plotPost(medical/general,compVal = 1)
0446１３２人目の素数さん
垢版 |
2020/05/01(金) 00:36:14.80ID:LcUyF6Tc
>>441
すべてが陰性に出るなら陽性は０人だろｗ
0448１３２人目の素数さん
垢版 |
2020/05/01(金) 10:05:01.58ID:LcUyF6Tc
>>447
だから、陽性が何人か出てんだから感度0%のインチキじゃなかろうってこと。
0449１３２人目の素数さん
垢版 |
2020/05/04(月) 00:17:48.34ID:Q8iLGwO2
最近の報道

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

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

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

問題

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

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

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

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

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

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

https://twitter.com/shukudai_sujaku
https://twitter.com/5chan_nel (5ch newer account)
0451１３２人目の素数さん
垢版 |
2020/05/04(月) 22:45:28.89ID:94tg6i4k
高橋洋一(統計数理研究所→大蔵省)
スウェーデンの感染症対策を解説
https://youtu.be/2bsYRuorwMI?t=192
0452１３２人目の素数さん
垢版 |
2020/05/05(火) 06:07:12.42ID:ht6rG86e
https://toyokeizai.net/sp/visual/tko/covid19/
のデータを使って
P=14895/153581 　#　2020/05/04
PCR検査の感度30-70％のモデルM1と感度70-90％のモデルM2のどちらが信憑性があるか、ベイズファクターで計算してみる。
M1は感度が最頻値60%標準偏差10%、M2は最頻値80%標準偏差10%のベータ分布に設定
特異度はいずれも最頻値95%標準偏差2.5%に設定し、有病率は一様分布を仮定
陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする。
事後確率分布は
https://i.imgur.com/81bK4KE.png
陽性率P=14895/153581=0.09698465での事後確率分布の密度比(Savage-Dickey density ratio)でベイズファクターを出すと
> d1/d2 # Savage-Dickey densiti ratio = BF12
[1] 1.007722

まあ、ちょっぴり、感度30-70％のモデルの方がいいかも、という結果。
陽性数/検査数の時系列データでもあればもう少し差がでるかもしれん。
東京都のデータで計算させようかと思ったが、東京都は検査人数を隠蔽しているので使いものにならない。
0453１３２人目の素数さん
垢版 |
2020/05/05(火) 06:26:36.20ID:FpoKo6a+
訂正

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

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

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

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

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

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

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

https://i.imgur.com/6ZYn34k.png
0459１３２人目の素数さん
垢版 |
2020/05/06(水) 15:58:11.98ID:RG00+xls
>>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
0461１３２人目の素数さん
垢版 |
2020/05/06(水) 17:08:09.44ID:tNZVV9ZH
>>456
感度が5%以下じゃつかいもんにならんわなぁ。インチキで確定。
0463１３２人目の素数さん
垢版 |
2020/05/07(木) 00:59:54.31ID:0vbdeA0/
>>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は過少では？過少に伴う推定誤差加算が必要では？
0464１３２人目の素数さん
垢版 |
2020/05/07(木) 02:09:04.11ID:VnQvkZ57
>>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
0465１３２人目の素数さん
垢版 |
2020/05/07(木) 02:30:17.69ID:VnQvkZ57
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))
0466１３２人目の素数さん
垢版 |
2020/05/07(木) 17:10:06.19ID:CHL0/p02
>>462
５０人で０でしょ？
だったら、特異度として保証されるのはせいぜい９８％程度じゃん。
１％の陽性率に対して、それでは不十分。
0467１３２人目の素数さん
垢版 |
2020/05/07(木) 18:13:41.63ID:VnQvkZ57
>>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
0468１３２人目の素数さん
垢版 |
2020/05/07(木) 18:52:47.08ID:VnQvkZ57
特異度の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
0469１３２人目の素数さん
垢版 |
2020/05/07(木) 20:13:12.78ID:VnQvkZ57
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
0470１３２人目の素数さん
垢版 |
2020/05/07(木) 21:49:37.22ID:CHL0/p02
>>469
なにやってるかよくわからんので、見当外れの指摘かもしれんが、
大阪市大の3/312と神戸大の33/1000ってのは特異度も感度も異なる
であろう別種のキットによる検査結果なんで、一緒くたにしちゃ
まずいんでないの？
0471１３２人目の素数さん
垢版 |
2020/05/07(木) 22:36:59.75ID:VnQvkZ57
>>470
同一キットじゃないから、ご指摘の通り。
しかも神戸大の方では陰性検体での確認はされていないから、神戸大方の陽性率が高いのは偽陽性を含む可能性もあるね。

大阪市大だけのデータでやってみると。

https://i.imgur.com/njVQtRZ.png
0472１３２人目の素数さん
垢版 |
2020/05/07(木) 22:46:09.05ID:VnQvkZ57
>>471
使える情報が少なすぎて信頼区間が広すぎ。
有病率:0-87%
キットの感度:0-70%
ってコイントス変わらんな。
0473１３２人目の素数さん
垢版 |
2020/05/07(木) 22:48:00.40ID:VnQvkZ57
結局のところ、断定的な結論は出せないねということを数字で確認しているだけｗ
0474１３２人目の素数さん
垢版 |
2020/05/08(金) 02:04:52.03ID:Ugc87SUM
>>464
>>465
有難うございます。
二項分布下の母比率の信頼区間推定問題なんですね。
A２：正規分布近似法でなくWilsonの信頼区間法で解くから
　　上側区間≠下側区間となり平均値まわり左右非対称信頼区間と
　　なるのですね。
0475１３２人目の素数さん
垢版 |
2020/05/08(金) 09:15:37.41ID:WmDpVhCu
3月の宿題で(1)のみ正解の数弱@shukudai_sujaku

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

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

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

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

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

https://twitter.com/shukudai_sujaku
https://twitter.com/5chan_nel (5ch newer account)
0476１３２人目の素数さん
垢版 |
2020/05/08(金) 11:37:38.79ID:CfFk/Uw1
>>474
Rのlibrary binomを使って　binom::confint(3,312)で　>464の出力が得られる
信頼区間をグラフにすると
https://i.imgur.com/Mlvv7bB.png
点線は3/312
正規分布近似のasymptotic以外は非対称。

どれを使うべきか？　好きなのを使えばいい、と思う。
但し、値が負になったり１を越えたりするのは不採用の方が賢明だとは思う。

Wilson法は値が０や１に近くても信頼できるという人がいるけど、どうやって検証するのかはよくわからん。
0478１３２人目の素数さん
垢版 |
2020/05/09(土) 13:22:10.71ID:oDT7bFgO
レムデシビルで初の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
なので、どちらのモデルが得られたデータを説明するのに適しているかは判断できず、という結果になった。
検出力不足か、差がないのかには決着つかず。
0479１３２人目の素数さん
垢版 |
2020/05/09(土) 13:22:52.37ID:oDT7bFgO
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)
0481１３２人目の素数さん
垢版 |
2020/05/09(土) 14:25:05.91ID:74hNX8Dr
100例もやって差がでないのなら、はっきり言って効き目なしとかわらん。
むしろ副作用がこわい。

アビガンはどうなんだろうね。なぜか国内のまとまったデータがまだ出てこない。
0482１３２人目の素数さん
垢版 |
2020/05/09(土) 14:26:34.05ID:74hNX8Dr
まあ、軽症者が８割とかだと、早期投与の効き目を判定するには相当の
人数に使わないとわかんないだろうけど。
0483１３２人目の素数さん
垢版 |
2020/05/09(土) 15:22:40.72ID:oDT7bFgO
死亡率に関してはベイズファクターで推論できた。

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を超えないいのでまあ、中程度の確信で検出力不足によるのではなく、もともと差がないのであろうと推論できる。

サンプルサイズを増やしてもレムデシビルは軽症・早期投与でも致死率を改善しないであろうと予想。
0484１３２人目の素数さん
垢版 |
2020/05/09(土) 17:15:11.02ID:oDT7bFgO
>どの週においても両群で症状改善率は不変
この前提はおかしいので週ごとに症状改善率が変わるようにモデルを変更

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.9976274

症状改善率に差がないというモデルも差があるというモデルも周辺尤度（エビデンス）はほぼ同じという結果。
0485１３２人目の素数さん
垢版 |
2020/05/09(土) 20:41:59.18ID:74hNX8Dr
あ、 >>482はレムデシじゃなくてアビガンについての話だかんね。
レムデシは重症者向けらしいけど、たぶん駄目だろ。
0487１３２人目の素数さん
垢版 |
2020/05/11(月) 00:31:07.85ID:mGn48zmS
5/10夜『NHKスペシャル』「新型コロナウイルス　出口戦略は」で
厚労省新型コロナ対策班西浦博北大教授と氏が分析中のノートPC画面が
映されていた。氏の使用解析ソフトウェア名は何というのでしょうか？
0490１３２人目の素数さん
垢版 |
2020/05/11(月) 03:58:21.10ID:v4rWY6J/
>>489
いや、与えられた何日分かの新規患者数から再生産数とか回復率とかを推定するための数式じゃないの？
0492１３２人目の素数さん
垢版 |
2020/05/11(月) 09:14:46.53ID:P7wrhagU
>>489
自分で計算してみた人ならわかるけど、日々の再生算数を出すにはA「感染してから他人にうつすまでの日数の分布」が必要。西浦先生はこれをまず求めていて、その結果は割合と世界中で使われている。
加えてB「感染してから発表の数字に載るまでの日数の分布」も必要で、発表の数字をBで逆畳み込み積分してC「感染推定日毎の数字」に変換し、CをAで畳み込み積分したものでCを割れば日々のRtが求められる。
とは言うものの、そのままやると誤差で数字が発散するし、逆畳み込みは一意ではないのでどう推定するかは色々と考えなくてはならない。多分日本独自の事情もコードに入れなきゃいけない。大変だろうなあとは思う。
0493１３２人目の素数さん
垢版 |
2020/05/11(月) 11:02:15.57ID:+2HA8Yn9
>>492
レスありがとうございます。
確かにSEIRモデルやSIRIモデルだと
「感染してから発表の数字に載るまでの日数の分布」って考慮されてないね。
0495１３２人目の素数さん
垢版 |
2020/05/12(火) 22:32:45.27ID:9+i4T0w3
ウイルス感染のモデルも必要だと思うけど
社会経済活動の影響や
自粛解除した後の死者増加や
自粛継続した場合の法人の倒産や
治療薬やワクチンができた後の社会の経済や文化面の回復
そこまで考慮してどういう政策を選択したら
損失を最小化できるかという問題は解ける？

制約付きの最小化問題だと思う
数値化するのが難しい考慮要素もあるからその辺は仮の値とかで
0496１３２人目の素数さん
垢版 |
2020/05/13(水) 17:23:13.01ID:PBQsFM+W
>>495
モデルを精密複雑化すればするほど、必要なパラメータが増えて
それを弱情報事前分布にすると事後分布の信頼区間が広がるんだなぁと思う。
■ このスレッドは過去ログ倉庫に格納されています
ニューススポーツなんでも実況