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.
654:１３２人目の素数さん
20/05/24 08:44:40 6T3JvGbL.net
@kazu_fujisawa
西浦先生が一番の功労者だというのは僕も同意見だな。あと、細かいこと言うと、
あの限られたデータしかない中で、８割削減なら1ヶ月で収束、接触6-7割削減なら2ヶ月で収束、という予測精度もお見事としか言いようがない。
655:１３２人目の素数さん
20/05/24 10:47:51 uCyuaPog.net
アレは単に再生産数の推定値から出しただけやろ。
656:１３２人目の素数さん
20/05/26 15:13:08 SaUF2oVL.net
わかった後からだと
「簡単じゃん！」とか言えるけど
先んじて発表した西浦先生は
すごくまともな研究者だと思います。
657:１３２人目の素数さん
20/05/26 16:32:28 +yxrXOhy.net
中国・武漢全域で行われた新型コロナウイルスのPCR検査で、189人が無症状の感染者と確認されました。
武漢では35日ぶりに新たな感染者が確認されたため、14日から「10日間大戦争」と銘打ち、延べ900万人にPCR検査を行いました。中国メディアによりますと、そのうち657万人の結果が出て、189人が無症状の感染者でした。
10万人のうち2.87人が無症状感染者の計算になるとしています。中国のSNSには、検査人数のあまりの多さに看護師が泣き叫ぶ動画が投稿されていて、
一日の検査能力である最大10万件を大幅に超えた検査の精度を疑問視する声もあります。10日間ですべての検査が終わらなかったので、検査は26日午後も続いています。
[2020/05/26 13:22]
URLﾘﾝｸ(news.tv-asahi.co.jp)
URLﾘﾝｸ(news.tv-asahi.co.jp)
658:１３２人目の素数さん
20/05/26 16:47:55.61 +yxrXOhy.net
189/6570000は
有病率、感度、特異度の事前分布はどうするかな？
有病率：一様分布
感度：最頻値0.6標準偏差0.1のベータ分布
特異度：最頻値0.9標準偏差0.05のベータ分布
として検査陽性数は有病率＊感度＋（１－有病率）＊（１－特異度）の確率に従う二項分布
でいいか?
659:１３２人目の素数さん
20/05/26 16:50:21.89 +yxrXOhy.net
>>618
再生産数から収束に至るまでの期間の計算法って公開されてたっけ？
660:１３２人目の素数さん
20/05/26 18:11:34 5HIYV6Gm.net
>>622 それは簡単に計算できるだろ？逆は難しくて、それが西浦先生の功績だと思うけど。
662:１３２人目の素数さん
20/05/26 18:34:15 aYF++qy3.net
>>622
いっぱんに再生産数Rの伝染病が収束を始める免疫の比率は1-1/R。
R=2.5なら免疫獲得者が1-1/2.5=60%になれば拡大はしない。
なので60%になれば平衡状態になり、最低でもそれ以上減らさないとダメなのは明らか。
集団免疫の話読んだ事ある人間なら常識。
オレみにわかではあるが。
感染率が70%減、80%減になったとき、どれくらいの速度で落ちていくかもsirでだいたいわかる。
663:１３２人目の素数さん
20/05/26 21:08:44.19 hVGZ6ofy.net
え？ >>617って皮肉じゃないの？ｗ
西浦さんの８割削減って言葉がメディアに乗るころには拡大再生産数は
1を切ってたんだよね（4/1頃）。拡大再生産数が2を越えて高止まりして
たのは3/15-3/25頃で、そこから１週間ほどで急落して1を切ってんだよね。
東京、大阪で自粛が叫ばれ、危機感がつのってたころだ。
664:１３２人目の素数さん
20/05/26 21:10:56.84 hVGZ6ofy.net
×拡大再生産数
○実効再生産数
665:１３２人目の素数さん
20/05/26 21:46:56 aYF++qy3.net
>>625
それは今でも良くわかってないだろ。
何せ４月の初め頃って検査リソースの限界で東京とか陽性率が８割に近かった、もうほとんど陽性確定の人しか検査してなかった頃。
この“認知されたケース(confirmed case”となったケースのピークと実際の感染者数のピークは当然ずれてる。
実際の東京の罹患者は8万人を超えてるらしいので、相当数が検査されずに見過ごされてる。
おそらく見過ごされた多くのケースは陽性者数がピークを迎えた(陽性率が８割ぐらいあった)4/15前後に発症したケースだろうけど、それは明らかにconfirmed casesのピークより遅い。
666:１３２人目の素数さん
20/05/27 00:41:21 SJ9C7hro.net
>>620
5～10人の検体を混ぜて検査して、陽性が出たロットだけ全部検査してるらしい。
URLﾘﾝｸ(www.sankei.com)
陽性はほとんどいないから、1日に50～100万人分検査できる。
667:１３２人目の素数さん
20/05/27 00:42:18.53 VNGePIH1.net
>>627
４月の初め頃の検査対象はその１週間前に感染してんだから、
まさに専門家委員会が示したconfirmed casesの感染者数の
ピークの頃だよ。
668:１３２人目の素数さん
20/05/27 01:15:53.98 sG8aLkL9.net
>>629
confermed casesのピークはもっと遅い４/15前後だよ。
その前と後に二つの最大のピークがある。
いずれにせよ、検査拒否の陽性率の最大がそのあたりにあって、どのくらいの見逃しがどのくらいあったのか科学的に推し量る方法などない。
その当時どんなコールセンターへの相談があり、どの程度怪しいケースを弾いていたのかなどコールセンターの資料をツブサに調べていくしかない。
そんなもんここで議論するような話ではない永遠の水かけ論にしかならない。
669:１３２人目の素数さん
20/05/27 04:32:30.99 DONHjyrT.net
>>624
収束する時期は再生産数だけじゃでは計算できないのではという疑問。
670:１３２人目の素数さん
20/05/27 07:13:32.08 DONHjyrT.net
>>628
有病率を１０万分の３として、何人の検体を混ぜるのが一番効率がいいんだろ？
１０万を３でわった３万３千人くらいかな？
671:１３２人目の素数さん
20/05/27 08:58:32.23 olFI7/IU.net
>>632
n/k+n/k(1-(1-p)^k)kを最小にするkだから、k=2
672:Productlog(-0.5√-log(1-p))/log(1-p)だね。 p=3/100000を代入してk=183。 ちなみにk=10となるのはp=0.0111の時。そんなに陽性率が高いと見積もったわけでは無いだろうから、10ぐらいが混ぜて検査の限界なのだろう。 もっと究極的に検査回数を少なくすると>>419 の方法に行き着くと思う。
673:１３２人目の素数さん
20/05/27 09:20:13.89 DONHjyrT.net
>>633
早速のレスありがとうございました。
674:１３２人目の素数さん
20/05/27 12:14:38 VNGePIH1.net
>>630
だ、か、ら、検査結果の判明した日付では4/15頃がピークだけど、検査の
待ち時間や発症までの潜伏期間を含めると、感染したのはその２週間以上
前になるっていうのが、専門家委員会の見積もり。
URLﾘﾝｸ(www.mhlw.go.jp)
675:１３２人目の素数さん
20/05/27 12:23:58 cpzgnD1W.net
>>635
confirmed cases の罹患→発症の平均時間と
all cases の罹患→発症の平均時間にはズレが出る。
confirmed cases のそれが2週間以上でそのピークが４月より前でも
all cases のそれはもっと小さいのでそのピークはもっと後ろになる。
676:１３２人目の素数さん
20/05/27 12:48:12.43 0C2sk63u.net
具体例で実験してみればいい。
山がある曲線かく。
山をy軸方向に0.7倍したもの(軽症段階で相談した人)と0.3倍したもの(重篤化して認知されたもの)を作り、後者はx軸反転する。
前者は右に1週間、後者は左手に一週間ずらす。
この二つの曲線の差がその日に検査要請がかかった件数。
しかし検査リソースには限界があるのでx軸より下の重症者は全部検査されたとして、そこから一定幅より上の部分は検査拒否されたと考える、つまりy軸より下の部分を検査リソースの限界分だけ上に平行移動したものより下が実際に検査された軽症者の人数。
もとまった検査された軽症者の曲線を改めて一週間前にずらし、重症者のカステラの曲線と足してみる、コレがconfirmed casesに限った推定発症日の曲線。
元の山とできた山のピーク比べてみる。
ズレてるから。
677:１３２人目の素数さん
20/05/27 13:01:39.14 VNGePIH1.net
>>636
だ、か、ら、>>627で検査が滞って隠れ感染者が相対的に増えてたと
君がいう４月上旬の検査対象は３月下旬の感染者なんだから、全感染
者数のピーク時期はconfirmed casesの4/1頃のピークの前にこそなれ、
後ろにはならんでしょ？
678:１３２人目の素数さん
20/05/27 13:04:06.21 Rn3481n9.net
>>637
なるよ。
>>635の実験やってみた？
必ずconfirmed casesのピークはall casesのピークより前にズレる。
つまりall casesのピークはconfirmed casesのピークより後ろにズレる。
679:１３２人目の素数さん
20/05/27 14:10:34 VNGePIH1.net
>>637
前提に問題がありすぎ。検査数は逐次増えてて、感染者のピーク前後では
ほぼ検査数に比例して陽性者数も増えてる（陽性率は十数％前後で大きな
変化はない）。キャパ一定で推移してるわけではない。
680:１３２人目の素数さん
20/05/27 14:36:45.45 m+sqB0p7.net
>>640
前にも書いたけどその話をしだすとココでは絶対答えが出ない水掛け論にしかならない。
結局症状が出てても検査を断られる４月頭時点での検査体制では、実際どれくらいのズレが発生したのかはハッキリとは言えない。
それは後年の調査である程度はわかるかもしれないけど。
何せ１日数万件もあった相談と検査拒否の実態を掴むのは容易ではない。
いずれにせよ、実際の感染のピークが
confirmed casesのピーク日
-conformed cases内の(確定日-推定罹患日)の平均
より後ろにズレるのは確実。
何故ならどの程度の検査拒否があったにせよ、この相談日-罹患日の値が小さいグループ(症状の軽いグループ)の方が検査拒否された割合が有意に高いので。
まぁ実際confirmed casesのピーク日-(xxx)
のxxxに入れるべき数値は2週間よりはずっと小さい。
681:１３２人目の素数さん
20/05/27 18:01:10.48 DONHjyrT.net
職場で最初に発症した人が感染源と決めつけてはならないことを
数学的には確認。
あるクリニックで院長が新コロナを発症。
翌日２人の看護師が発症、その次の日に３人の看護師が発症。
各人が最初に感染していた確率を計算せよ。
こういう計算ができると、誰が移したとか責めたりせずにすむ。
（春節ウェルカムした安倍が悪い、で職員の心が一つになるｗ
計算に必要な数値は
#--- incubation period ---
# URLﾘﾝｸ(www.nejm.org)
# Li et al NEJM 2020
# 潜伏期間は対数正規分布でパラメータは
ln_par1 = 1.434065
ln_par2 = 0.6612
682:１３２人目の素数さん
20/05/27 18:03:06.39 DONHjyrT.net
RとStanでやってみた。
> data.frame(p=round(sapply(1:n,fi),2))
p
1 0.26
2 0.18
3 0.18
4 0.13
5 0.13
6 0.13
コードはここ
ｽﾚﾘﾝｸ(hosp板:199番)
683:１３２人目の素数さん
20/05/27 18:05:18.29 1E963DQR.net
>>632
ウイルスを検出できる最低限度の濃度とかあるんじゃないかな
制約付きで何人分を混ぜるかを最大化する問題になると思う
684:１３２人目の素数さん
20/05/27 19:07:41.79 VNGePIH1.net
>>641
しつこいな。どこからそんな式が出てくるんだよｗ
現状で分かってる事実から推測して全感染者数の
ピークが後ろにくる理由はないよ。むしろ前にくるほうが自然。
685:１３２人目の素数さん
20/05/27 19:13:02 VNGePIH1.net
>>642
発症前から感染力があり、潜伏期間にも幅があるんだから、
その情報だけでだれが感染源かなんて断定できるわけない
のは当たり前だろ。それでも最初の発症者が一番怪しい
ってことに変わりはない。
常識で考えろ。
686:１３２人目の素数さん
20/05/27 19:21:49 v/kXgPR4.net
>>645
>>637の実験やってみた？
687:１３２人目の素数さん
20/05/27 19:34:58.61 DONHjyrT.net
>>646
院長と３日めの発症した看護師とどれくらい確率が違うかは
常識では出せないだろ。
こういう計算には数値化が必要。
ある開業医が新型コロナ肺炎に罹患したとする。
行動調査によって発症前にキャバクラに行っており接客したキャバ嬢が開業医発症の２日後に発症していたことがわかった。
キャバ嬢は開業医から移されたと主張して１億円の賠償を求めている。
潜伏期間には幅がありキャバ嬢から移された可能性もあると主張してその確率を計算して賠償金を値切りたい。
いくら値切れるか計算せよ。
688:１３２人目の素数さん
20/05/27 20:06:50 olFI7/IU.net
>>642
Mathematicaだと
cond = Distributed[{a, b, c, d, e, f},ProductDistribution[{LogNormalDistribution[1.434065, 0.6612], 6}]];
NProbability[a > Max[b - 1, c - 1, d - 2, e - 2, f - 2], cond]
NProbability[b - 1 > Max[a, c - 1, d - 2, e - 2, f - 2], cond]
NProbability[d - 2 > Max[a, b - 1, c - 1, e - 2, f - 2], cond]
>0.257392
>0.178947
>0.128264
689:１３２人目の素数さん
20/05/27 20:50:23 DONHjyrT.net
>>649
レスありがとうございます。
有効数字２桁で同じ結果なので自分のコードの検証になりました。
690:１３２人目の素数さん
20/05/28 18:53:42 9QoKXLHk.net
中国湖北省武漢市（人口約１１００万人）が今月中旬から全市民を対象に実施していたＰＣＲ検査（遺伝子検査）がおおむね終了した。
地元紙「湖北日報」（電子版）によると、直近に検査を済ませた人や乳児を除く９００万人以上から検体を採取。
武漢では１日に１００万人以上の検査が可能で、２４日までに約６５０万人の検査が終わり、無症状感染者が２１８人見つかった。
事前分布に有病率は一様分布、感度は40-60%、特異度は90－100%のベータ分布を設定してStanでMCMCしてみた。
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
spc 0.999977 0.000000 0.000008 0.999965 0.999976 0.999994 12593 1.00001
sen 0.452612 0.002116 0.220226 0.073638 0.444200 0.871391 10834 1.00000
prev 0.000041 0.000001 0.000118 0.000001 0.000023 0.000183 8206 1.00004
p 0.000034 0.000000 0.000002 0.000030 0.000034 0.000039 36412 0.99998
spc:特異度、sen:感度 prev:有病率 p:陽性率
900万人での陽性者数の予想
最頻値 306.15
95%信頼区間
lower upper
266.46 346.70
1100万人のうちの有病者数予想
> summary(ms$prev*11e6) ; MODE(ms$prev*11e6)[1] ; hdi(ms$prev*11e6)[1:2]
Min. 1st Qu. Median Mean 3rd Qu. Max.
0 108 251 449 483 137796
x
179.24
lower upper
0.020213 1360.523395
691:１３２人目の素数さん
20/05/29 17:18:00.15 OOwyRgpM.net
岩田健太郎 Kentaro Iwata
@georgebest1969
日本が第一波をかなりうまく乗り切ろうとしているのだけど、最大の功労者の一人は西浦博先生だよ。それは絶対に間違いない。
午後5:44 2020年5月20日
692:１３２人目の素数さん
20/05/29 17:34:00 2wb0CO5Y.net
>>652
俺は志村けんとフィリピンパブ親父だと思う。
693:１３２人目の素数さん
20/05/29 19:29:59 7LO8JYfk.net
政府：「接触８割減を目指して下さい。」
マスコミ：「今朝の○○駅の映像をご覧下さい。人出８割減にはほど遠いですね．．．」
政府広報：
「接触８割減には、(1-x)^2=1-0.8 → x=0.552786 だから 人出55%減でいいんだけど、
馬○なマスコミが、勝手に目標を上げてくれてる。あえて指摘する必要ないよね．．．」
694:１３２人目の素数さん
20/05/29 19:53:22.15 iNWKoruo.net
>>653
藤浪も仲間に入れてやって。
岡江久美子も５月１０日前後の感染急減には貢献したかも。
しかし、岩田健太郎に擁護されてもむしろ逆効果でお気の毒＞８割おじさん
695:１３２人目の素数さん
20/05/29 19:56:02.55 iNWKoruo.net
>>648
賠償金ゼロまで値切れる。
俺じゃない可能性が十分あるのならそれでいい。
期待値で賠償金を計算されちゃかなわん。
696:１３２人目の素数さん
20/05/29 21:32:01 2wb0CO5Y.net
>>656
期待値を求めよ、では食いつきが悪いと思って金の問題にアレンジして開業医スレに書いてみたんだが、
答られる開業医は０だったよ。
697:１３２人目の素数さん
20/05/29 21:33:11 2wb0CO5Y.net
>>655
石田と志村の順が逆だったら、自粛は捗らなかったのではと思うな。
698:１３２人目の素数さん
20/05/29 21:44:18 2wb0CO5Y.net
潜伏期間が従うとされる対数正規分布（パラメータは、ln_par1 = 1.434065 ln_par2 = 0.6612)で乱数発生させて
その差が２以上になる割合を求めれば数値解は出せるけど、これって解析解は可能だのだろうか？
パラメータが同じ対数正規分布の差の分布は数式で表せるのだろうか？引き算したら０になってしまう気がする。
699:１３２人目の素数さん
20/05/29 21:51:25.78 2wb0CO5Y.net
>>654
感染の起こりうる濃厚接触するときの人数って期待値はどれくらいで、どんな分布に従うんだろうね？
700:１３２人目の素数さん
20/05/30 10:10:43.57 oTxR9Gdg.net
>>659(自己解決)
mu = 1.434065
sg = 0.6612
pdf1 <- function(x) dlnorm(x,mu,sg)
pdf2 <- function(y) dlnorm(y,mu,sg)
f <- function(x,y) pdf1(x+y)*pdf2(y)
vf=Vectorize(f,vectorize.args = 'y')
pdf <- function(x) integrate(function(y) vf(x,y),-Inf,
701:Inf)$value pdf=Vectorize(pdf) curve(pdf(x),-30,30)
702:１３２人目の素数さん
20/06/05 09:35:47.94 a3w5V/C1.net
COVID19関連での分布とそのパラメータ
Distributional fits to key COVID-19 distributions.
URLﾘﾝｸ(i.imgur.com)
URLﾘﾝｸ(www.medrxiv.org)
703:１３２人目の素数さん
20/06/05 10:32:03.26 a3w5V/C1.net
学校閉鎖の功罪
School closure and management practices during coronavirus outbreaks including COVID-19: a rapid systematic review
URLﾘﾝｸ(www.thelancet.com)(20)30095-X.pdf
704:１３２人目の素数さん
20/06/05 10:41:43.54 Vcc0zApK.net
>>654
西浦博教授が新しい試算を発表されたようです。
1影のたけし軍団 ★2020/06/03(水) 18:41:19.54ID:v4jJ+LDc9
西浦博教授らは、海外から新型コロナウイルスの感染者が入国すると、
国内で再び大規模な流行が起きるとする試算をまとめた。
1日当たり10人の場合、3カ月以内に再流行する確率は98.7%に達すると・・・・・・・
【コロナ】 西浦博教授、試算・・・海外から1日10人入国で再流行 [影のたけし軍団★]
ｽﾚﾘﾝｸ(newsplus板)
705:１３２人目の素数さん
20/06/05 11:24:15 a3w5V/C1.net
>>664
この試算コードは公表されている？
706:１３２人目の素数さん
20/06/05 14:20:57 PWYfRUOi.net
m3（医療関係者のサイト？）に割と詳しく載ってるらしい
707:１３２人目の素数さん
20/06/05 22:30:49 a3w5V/C1.net
>>666
みてきた。未知のパラメータが多すぎるのではという印象。
（3）どうやって計算したか
90日間の入国者総数をN人とし、感染率（感染している者の割合）をp、停留やPCRの効果（すり抜ける感染者の相対的減少率）をeとします。
そのとき、侵入者は
n = (1-e)pN
と計算されます。
上述の1人の感染者が侵入したときの「絶滅確率」をqとし、n人が独立して入国してきたとすると、
分岐過程に基づく大規模流行の確率Xは以下で与えられます。
X = 1 - q^n = 1 - q^((1-e)pN)
708:１３２人目の素数さん
20/06/05 22:35:27 a3w5V/C1.net
Rt のときと違ってあまりにも簡単な式でびっくりした。
709:１３２人目の素数さん
20/06/05 23:43:21.30 PWYfRUOi.net
みほちゃんというアカウントの物理学者がTwitterで批判してる 正しい気がするんだが、どうなん？
710:１３２人目の素数さん
20/06/06 09:52:55.88 dnuHAH8y.net
>>667
グラフにしてみた。
URLﾘﾝｸ(i.imgur.com)
711:１３２人目の素数さん
20/06/06 11:36:10.82 n1RssdRL.net
URLﾘﾝｸ(webcache.googleusercontent.com)
m3の中味
URLﾘﾝｸ(togetter.com)
ツイートまとめ
712:１３２人目の素数さん
20/06/09 09:32:36.39 olJVeaWc.net
ここ数日流れている、7割の人は誰にもうつさないの元論文。
URLﾘﾝｸ(www.researchsquare.com)
figure2Bが個人の再生産数の分布で、負の二項分布で近似可
713:１３２人目の素数さん
20/06/09 10:15:44.14 vlglV9Ra.net
>>672
abstract読んだだけではわからんかった。
どういうことか詳しくお願い。
ってか、そういうのも考慮にいれて実効再生産数の値が決まるんでは？
714:１３２人目の素数さん
20/06/09 10:49:09.07 olJVeaWc.net
>>673
これ、まさに>>667あたりの議論に関係してくるんだよ。実効再生産数が1を超えていても、スーパースプレッダーが拡散していてその値なら、運が良ければ大規模拡散は無い。逆に全員が均等に1人以上にうつしてしていたら、確実に大規模拡散につながる。
拡散能力の個人差を決めるのが、論文で言うkだと思う。
715:１３２人目の素数さん
20/06/09 11:03:38.13 BdGxyACB.net
>>674
拡散能力は何で決まるかわかっている？
個人の行動で決まるなら行動変容は必要だろう
マスクや咳エチケット、消毒、物理的距離や
換気など
716:１３２人目の素数さん
20/06/09 11:18:48 olJVeaWc.net
>>675
だから、それらがいい加減でRt>1でも拡散しない確率があるという話。帰国者のように少数が散発的に爆撃して行く場合には、Rt値だけでは判断できないわけだ。
東京も、日々の感染者数が１人とか２人になってくるとこの辺りの話が効いてくるだろうね。
717:１３２人目の素数さん
20/06/09 11:25:34.75 QHj7P+P7.net
そうそう。西浦先生の >>671 で面白いのは、帰国者が大規模感染を起こす云々より、むしろ再生産数が1を超えているのに感染が止まる確率があることなんだよね。
718:１３２人目の素数さん
20/06/09 12:01:48.99 vlglV9Ra.net
>>674
よくわかんないんだけど、感染者数が全体で数十人レベルなら隔離対策で
全部のスーパースプレッダーをたまたま捕縛できれば確かに終息できる
だろうけど、数百人レベルに拡がってそういうことができなくなれば、
一定数以上のスーパースプレッダーは常に残ることになって、結局感染
拡大は防げないのでは？（SS1人当たりの再生算数が４だと、3割がSSだと
して、実効再生産数は単純に0.3×4=1.2とかにならんの？）
719:１３２人目の素数さん
20/06/09 13:33:10.98 BdGxyACB.net
>>677
実行再生算数は過去のデータを反映した値で
Rt>1でも患者を全員隔離できたら収束可能ではあるから
矛盾はしないと思う
それを実現できる可能性は低いだろうけど
720:１３２人目の素数さん
20/06/09 20:27:46.64 olJVeaWc.net
>>678
そう。だから>>671の(4)の表で入国者数×感染率が数百単位の場合は大流行確率がほぼ100%になってる。あくまでも少人数しか感染していない場合の話。
>>679
隔離の話ではなく、感染力のバラツキに起因する感染消滅。サイトに載ってる考え方でシミュレーションを回して自然消滅確率q=0.9226は再現できたよ。
721:１３２人目の素数さん
20/06/09 23:46:15.61 Rz+Wm47s.net
>>672
理論的根拠があって負の二項分布を選んだのかと思って読んでみたらRであてはまりのよい分布を探したみたい。
Calculation of serial interval and empirical offspring distribution
We calculated the median serial interval as the difference between the symptom onset dates of each infector-infectee pair, excluding asymptomatic cases, and fitted normal, lognormal, gamma and Weibull distributions using the R package “fitdistrplus" and maximum-likelihoo
722:d methods, excluding eight non-positive data points for the latter three distributions. We generated the empirical offspring distribution from the observed number of secondary cases per individual infector and similarly fit negative binomial, geometric and Poisson distributions as before.
723:１３２人目の素数さん
20/06/10 05:29:56.44 uBmNcOzn.net
>>681
関連論文をいくつか読んでみたのだが、個人ごとのRtに注目する研究では、その違いに負の二項分布が使われるようだ。
バラツキを表すdispersion parameter のk(西浦が0.1と置いてるのはこれ)を成功回数、1/(1+Rt/k)を成功確率とした負の二項分布をどの論文でも使っていて、これで分布の平均はRtになる。
本来整数であるべき成功回数を実数としているのだから、負の二項分布である意味はなく本当に当てはめただけなのだろう。
724:１３２人目の素数さん
20/06/11 23:28:19.41 /K72s0+6.net
一人を感染させる時間=serial interval(SI)
Rt人感染させるまでに要する時間は
形状パラメータ=Rt
尺度パラメータ=SI
のガンマ分布gamma(Rt, SI)に従うはず、
λがこのガンマ分布に従う変数として平均λのポアソン分布pois(λ)を考えれば負の二項分布になるのだけど
これはRtの分布じゃなくてRt人感染させる時間の分布だよなぁ、と考えたりしたけど
理詰めじゃなくて、AICが最も小さくなる既知の分布を探したってことなのだろうな。
725:１３２人目の素数さん
20/06/12 15:38:09.41 SHffzP0W.net
ソフトバンクグループは、社員や医療関係者らを対象に独自に行った新型コロナウイルスの抗体検査の結果を速報値として公表しました。4万4000人余りのうち、0.43％に当たる191人が陽性だったということです。
ソフトバンクグループは、先月中旬から今月8日までに、社員や取引先、それに医療機関の関係者ら4万4066人を対象に新型コロナウイルスに感染したことがあるかを調べる抗体検査を行い、結果を公表しました。
それによりますと、191人、0.43％が陽性だったということです。
このうち、医療機関の関係者は5850人のうち105人、1.79％が陽性、社員や取引先の関係者は3万8216人のうち86人、0.23％が陽性でした。
以下全文はソース先で
2020年6月9日 23時20分
URLﾘﾝｸ(www3.nhk.or.jp)
東大の7/1000の0.7%に近い感じだな。
一方、ロシアでは
ロシア・モスクワで保健当局が市民を対象に新型コロナウイルスの抗体検査を行ったところ、
およそ６人に１人にあたる１７．４％が陽性だったことがわかりました。
発表によりますと、モスクワの保健当局は１１日までの２週間で市民１６万５千人を対象に抗体検査を実施しました。
このうち、およそ６人に１人にあたる１７．４％の人に抗体があることが判明したということです。
先月４日から２１日まで行われた前回の抗体検査では１２．５％の人から抗体が見つかっていて、
前回から今回にかけて、陽性者がおよそ５ポイント増えた形です。
ロシアの感染者は１１日、５０万人を超えましたが、モスクワでは感染が抑制されてきているとして、
９日には外出制限が解除されています。
URLﾘﾝｸ(www.news24.jp)
726:１３２人目の素数さん
20/06/12 18:24:28.46 6ksURmJ8.net
>>684
特異度がわかんないけど、少なくとも99.8%くらいはありそうだな。
したがって、医療関係者が2%足らず陽性ってのは本物っぽい。
だから、どうしたって感じだけどね。
727:１３２人目の素数さん
20/06/12 21:51:01.76 DZUOxvlk.net
>>685
まだ未感染の人の割合が多いから集団免疫ができていなくて感染爆発のリスクが比�
728:r的高い 予防対策が引き続き必要
729:１３２人目の素数さん
20/06/12 22:01:46.43 opvlT+26.net
中和抗体の存在が実証されるまでは集団免疫神話だと思っている。
730:１３２人目の素数さん
20/06/12 22:16:50 6ksURmJ8.net
>>686
感染者の多い東京に限っても抗体陽性率は0.6%以下だってことは一ヶ月前から
分かってた話なんだから、今回0.23%以下と出たからといって、本質的になんも
かわらんでしょ。しいていえば、致死率の下限が上がったってことくらいか。
感染者が0.23%以下だとすれば、現在の死者数はやく千人なので、致死率は0.4%以上
となる。
731:１３２人目の素数さん
20/06/16 09:59:45.11 RYoHf/j0.net
外来患者0.5％に抗体、山形大 付属病院の1009人を検査
山形大（山形市）は15日、医学部付属病院の外来患者1009人について新型コロナウイルスの抗体検査をした結果、0.5％に当たる5人が陽性だったと発表した。山形県内でPCR検査により感染確認されたのが69人なのに対し、統計学上は県民約107万人のうち、670～1万人が感染した経験があると推計されるという。
大学によると、1～4日に外来で受診した人の血液を、検査キットより精度の高い機器と試薬で調べた。ただ検査数が千人余りに限られており、推計される陽性率は0.063～0.937％にばらつくという。
732:１３２人目の素数さん
20/06/16 10:10:06.88 RYoHf/j0.net
>>689
この95％信頼区間は正規分布での近似で計算しているな。
> p=5/1009
> q=1-p
> n=1009
> p-qnorm(0.975)*sqrt(p*q/n)
[1] 0.0006226557
> p+qnorm(0.975)*sqrt(p*q/n)
[1] 0.009288147
733:１３２人目の素数さん
20/06/16 10:12:02.29 RYoHf/j0.net
各種の95％信頼区間
> binom.ci(5,1009)
method x n mean lower upper
1 agresti-coull 5 1009 0.004955401 0.0017596471 0.011906321
2 asymptotic 5 1009 0.004955401 0.0006226557 0.009288147
3 bayes 5 1009 0.005445545 0.0014689159 0.010037024
4 cloglog 5 1009 0.004955401 0.0019155107 0.011096185
5 exact 5 1009 0.004955401 0.0016109041 0.011526082
6 logit 5 1009 0.004955401 0.0020640649 0.011848825
7 probit 5 1009 0.004955401 0.0019822732 0.011396609
8 profile 5 1009 0.004955401 0.0017805162 0.010620264
9 lrt 5 1009 0.004955401 0.0017936701 0.010616984
10 prop.test 5 1009 0.004955401 0.0018257774 0.012233747
11 wilson 5 1009 0.004955401 0.0021184531 0.011547515
734:１３２人目の素数さん
20/06/16 10:25:31.57 IaE5USJ0.net
厚労省がやってた抗体検査も結果出た。
URLﾘﾝｸ(www.mhlw.go.jp)
735:１３２人目の素数さん
20/06/16 12:01:44.37 RYoHf/j0.net
>>692
ロッシェとアボトットという業界では一流ブランドなのに一致率が低いのに驚き
TYO=matrix(c(2,4,2,1963),2,b=T)
> OSA=matrix(c(5,5,11,2949),2,b=T)
> MYA=matrix(c(1,6,2,3000),2,b=T)
> (dat=TYO+OSA+MYA)
[,1] [,2]
[1,] 8 15
[2,] 15 7912
> library('epiR')
> epi.kappa(dat)
$prop.agree
obs exp
1 0.9962264 0.9942306
$pindex
est se lower upper
1 -0.9942138 0.0008513625 -0.9958825 -0.9925452
$bindex
est se lower upper
1 0 0.0008518883 -0.00166967 0.00166967
$pabak
est lower upper
1 0.9924528 0.9892346 0.9949051
$kappa
est se lower upper
1 0.3459338 0.01121544 0.323952 0.3679157
$z
test.statistic p.value
1 30.84442 6.655295e-209
$mcnemar
test.statistic df p.value
1 0 1 1
736:１３２人目の素数さん
20/06/16 12:
737:19:16.92 ID:RYoHf/j0.net
738:１３２人目の素数さん
20/06/16 13:07:27.70 RYoHf/j0.net
>>694
これだけ一致率が低いデータで東京大阪宮城で陽性率に有意差があるを検定することに意味があるんだろうか？
739:１３２人目の素数さん
20/06/16 13:24:20.56 RYoHf/j0.net
>>692
多重比較になるけど、補正なしでも都市間の有意差はでないね。
むしろ、有名メーカーの検査キットの不一致の方が意外だった。
> positives=c(TYO=2,OSA=5,MYA=1)
> subjects =c(TYO=1971,OSA=2970,MYA=3009)
> pairwise.prop.test(positives,subjects,p.adjust.method = 'none')
Pairwise comparisons using Pairwise comparison of proportions
data: positives out of subjects
TYO OSA
OSA 0.82 -
MYA 0.71 0.21
P value adjustment method: none
Warning messages:
1: In prop.test(x[c(i, j)], n[c(i, j)], ...) :
Chi-squared approximation may be incorrect
2: In prop.test(x[c(i, j)], n[c(i, j)], ...) :
Chi-squared approximation may be incorrect
3: In prop.test(x[c(i, j)], n[c(i, j)], ...) :
Chi-squared approximation may be incorrect
> library(fmsb)
> pairwise.fisher.test(positives,subjects,p.adjust.method = 'none')
Pairwise comparisons using Pairwise comparison of proportions (Fisher)
data: positives out of subjects
TYO OSA
OSA 0.71 -
MYA 0.57 0.12
P value adjustment method: none
740:１３２人目の素数さん
20/06/16 21:54:44.89 tyACuXlX.net
まあ、特異度99.5%でも偽陽性が15人ぐらいでるから、陽性者はほとんどそれなんだろうな。なので両方陽性の人のみを数えてるのは正しいと思う。
日本の抗体保有率が低すぎたと。
741:１３２人目の素数さん
20/06/16 22:19:18.40 5zvBa6Bd.net
統計勉強しててこのスレ見つけたんですけどみなさんRつかってるんですね
pythonや他の言語ではないということは統計界隈はRが主流なのでしょうか？
742:１３２人目の素数さん
20/06/17 06:59:24.38 ozGKItOe.net
>>697
偽陽性率 ＝ 1 - 特異度だから全部が偽陽性の可能性もあるなぁ。
その確率は計算できるのだろうか？
743:１３２人目の素数さん
20/06/17 07:14:15.69 Sjc3FLqU.net
>>699
URLﾘﾝｸ(www.magojibi.jp)
URLﾘﾝｸ(www.magojibi.jp)
ここを信じる限り、特異度は99.6%と99.8%。3000人の検査で6～12人偽陽性が出る感じで、あってるように思う。
両方同時に偽陽性になる検体が1つ以上出る確率は2.3%だから、2つ調べればほぼ確実ということだろうね。まあ、互いに独立した事象ならという仮定だが。
744:１３２人目の素数さん
20/06/17 10:20:29.77 dTTrdB49.net
2キットが両方陽性のときのみ検査陽性とすると
特異度は
1-((1-0.996)×(1-0.998))
=0.999992
でいいのかな？
745:１３２人目の素数さん
20/06/17 10:23:14.49 dTTrdB49.net
偽陽性率は100万分の8でいいのか。
すると全部偽陽性というわけでもなさそう。
746:１３２人目の素数さん
20/06/17 16:40:39.13 Sjc3FLqU.net
>>702
あくまで、偽陽性の発生が2つの試験で独立の場合。
実際には、判定メカニズムに共通部分があるだろうから、完全に独立ではなくある程度相関はあるかも。相関が低いほど最終的な特異度は高くなる。
そういう意味では、カッパ値が低いことはむしろ利点となる。
747:１３２人目の素数さん
20/06/17 18:20:14.90 ozGKItOe.net
レスありがとうございます。
有病率prevalenceが低いとカッパ値は低めにでるらしいので
これを補正した
PABAK = ２✕(陽性合致数＋陰性合致数）／ 検査総数 − 1
は>693にある通り0.992と高いので
２キットは相関しているような気もする。
748:１３２人目の素数さん
20/06/18 11:44:36.30 PfxBa8Vk.net
>>698
もともと統計畑ではRを使っていた人が多数派だったと思う
Pythonは統計では後発で、pandasはRのデータフレームによく似ている
749:１３２人目の素数さん
20/06/18 12:50:29.49 +4oV8c92.net
>>705
統計のパッケージが一番揃っているよね。
最新理論のパッケージもまず、R版がリリースされるし。
例：Bayesian Network Meta-analysisとか。
こんなのあったらいいなぁと思って探すとみつかることが多い。
むかし、ヨンクヒール・タプストラ検定のスクリプト作ったけど、パッケージの方が高速でエラー処理も優れていたなぁ。
Rで??Jonckheereと入力するとそれらしいのがヒットする。
750:１３２人目の素数さん
20/06/18 20:20:41.28 +4oV8c92.net
> (dat=TYO+OSA+MYA)
[,1] [,2]
[1,] 8 15
[2,] 15 7912
で検査総数は7950
全員が抗体不保持でキットは全部偽陽性と仮定。
２キットが陽性で検査陽性判定なら8/7950が偽陽性で
特異度 1 - 8 / 7,950 = 0.9989937
１キットでも陽性で検査陽性判定なら
特異度 7,912 / 7,950 = 0.9952201
751:１３２人目の素数さん
20/06/19 15:35:57.95 qWJYCYqQ.net
抗体保有率＝０という帰無仮説で検定するにはどうすればいいんだろう？棄却できないと気がするんだが。
どちらのキットも陽性が23人で陽性率は
23 / 7,950 = 0.0028931
公称の偽陽性率(= 1 - 特異度)程度の陽性率にとどまっている。
752:１３２人目の素数さん
20/06/20 14:39:39.05 KyfQupfS.net
６月末に全国初の治験開始へ『大阪発の新型コロナワクチン』大阪府と阪大などが進める
URLﾘﾝｸ(www.mbs.jp)
>新型コロナウイルスのワクチンの治験は全国で初めてで、安全性が確認できれば、
>今年１０月に数百人規模で治験を行った上で、年内に２０万人分のワクチンを製造し、
>２０２１年春か秋の実用化を目指すということです。
有病率が0.5%程度だから数百人規模の治験じゃ有意差出せないだろうな。
753:１３２人目の素数さん
20/06/20 15:18:10.16 S9Mc7SGy.net
患者数が減って標本数が足りないって問題があるみたい
754:１３２人目の素数さん
20/06/20 19:20:16.10 tzFG6qz7.net
公表された抗体陽性割合
8/7950 # 厚労省
7/1000 # 東大
191/44066 # ソフトバンク
一番高い0.7％を採用して
>今年１０月に数百人規模で治験
を実薬500人、偽薬500人として
偽薬の500人中4人が陽性、ワクチン投与で500人中0人にできたとすると
> prop.test(c(4,0),c(500,500))
2-sample test for equality of proportions with continuity correction
data: c(4, 0) out of c(500, 500)
X-squared = 2.259, df = 1, p-value = 0.1328
alternative hypothesis: two.sided
95 percent confidence interval:
-0.001808434 0.017808434
sample estimates:
prop 1 prop 2
0.008 0.000
Warning message:
In prop.test(c(4, 0), c(500, 500)) :
Chi-squared approximation may be incorrect
> Fisher.test(c(4,0),c(500,500))
Fisher's Exact Test for Count Data
data: cbind(hit, shot - hit)
p-value = 0.1242
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.6617083 Inf
sample estimates:
odds ratio
Inf
> poisson.test(c(4,0),c(500,500))
Comparison of Poisson rates
data: c(4, 0) time base: c(500, 500)
count1 = 4, expected count1 = 2, p-value = 0.125
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
0.660124 Inf
sample estimates:
rate ratio
Inf
どれでも、有意差がでない。
どういうデザインで治験するんだろうなぁ？
755:１３２人目の素数さん
20/06/20 20:58:51 S9Mc7SGy.net
ワクチンと治療薬だと治験の方法が違うかも
安全性と効果を確認するだろうけど
ワクチンの場合は
健康な人がワクチン接種で病気などにならない事と
ワクチン接種後に
756:できた抗体でウイルスなどを中和出来ることを試験管なとで確認したら良さそう 治療薬は患者に投与して 副作用があまりない事の確認と治癒などの効果確認をするだろう
757:１３２人目の素数さん
20/06/20 21:56:40.97 tzFG6qz7.net
>>712
なるほど。
エンドポイントを感染予防でなくて抗体獲得という代理エンドポイントにするということね。
ありがとうございました。
758:１３２人目の素数さん
20/06/21 21:48:59.46 oYQHWaj+.net
岩手は、やはり抗体も検出されず
URLﾘﾝｸ(www3.nhk.or.jp)
759:１３２人目の素数さん
20/06/24 17:52:58.44 uHXQOLLh.net
こういう問題もあるようだな。
>>
ワクチンによって正常な抗体が作られる保証もない
中途半端な抗体ができるとADEとなり、コロナに感染した時に重篤化する
サイトカインストームみたいな状態になる
この場合、ワクチンの副作用として失敗作となる
今の段階はワクチンを接種しただけで、単体の副作用ないということが確認出来ただけ
760:１３２人目の素数さん
20/06/30 07:15:26.87 vzuOiySp.net
>>697
東大の児玉龍彦氏がその計算法に異議を唱えているね。
URLﾘﾝｸ(youtu.be)
27分30秒辺りから
761:１３２人目の素数さん
20/06/30 07:18:54.51 vzuOiySp.net
ワクチン神話を疑え！SARSで17年ワクチンができないわけ【新型コロナと闘う 児玉龍彦×金子勝】
で動画検索すると出てくる。
762:１３２人目の素数さん
20/06/30 07:36:35.39 HSIVZLVa.net
金子勝ってまだ生きてたのか
763:１３２人目の素数さん
20/06/30 08:06:57.85 vzuOiySp.net
>>716
デング熱ワクチンのことを調べてみたらこんな記述にであった。
URLﾘﾝｸ(www.forth.go.jp)
予防接種
2015年12月に、流行地域の9－45歳に向けて、サノフィ・パスツールによる初めてのデング熱ワクチンDengvaxia(CYD-TDV)が、20ヶ国で承認登録されました。
2016年4月、WHOはデング熱が住民の血清抗体保有率が70％以上の高度に常在している地域でのワクチン使用の条件付き勧告を行いました。
2017年11月に、予防接種時の血清抗体保有状況を後ろ向き研究で明らかにするための追加分析の結果が発表されました。
この分析からは、最初の予防接種の時点で血清抗体陰性であることが推測される研究の参加者では、予防接種を受けていない人と比較して、より重症なデング熱と入院のリスクがより高くなるということが示されました。
764:１３２人目の素数さん
20/06/30 09:22:31.83 dful2LQF.net
抗体を効果で分類すると3つあって
役なし、善玉、悪玉
善玉が中和抗体で期待するもの
役なしは変化なし
悪玉は抗体がある方が症状が重くなる
765:１３２人目の素数さん
20/06/30 11:35:22.28 vzuOiySp.net
>>720
役なしは感染既往の指標にはなるよ。
B型肝炎でいうとHBc抗体くらいの位置づけ。
766:１３２人目の素数さん
20/06/30 11:36:33.02 ouOildtQ.net
>>717
金子勝ってだけで見る気がせんわ。アホすぎ。
767:１３２人目の素数さん
20/06/30 11:40:04.49 ouOildtQ.net
>>720
仮にワクチンで抗体ができても、その３つのどれに該当するかを判定するのは
デング熱の場合よりはるかに難しそうだな。
768:１３２人目の素数さん
20/06/30 11:46:44.85 vzuOiySp.net
>>723
とりわけ有病率が低い国で有意差出すのは至難の業だと思う。
>>722
俺も児玉氏のところしか聞いていない。
769:１３２人目の素数さん
20/06/30 13:16:14.51 vzuOiySp.net
>>716
こっちの方が再生できるようだ。
://www.youtube.com/watch
770:?v=y6W83Y85zJs&t=1634s ワクチン神話を疑え！SARSで17年ワクチンができないわけ 【新型コロナと闘う 児玉龍彦×金子勝】
771:１３２人目の素数さん
20/06/30 13:48:18 b0dQmLIf.net
アラビア人のついた嘘だ
数学掲示板群 URLﾘﾝｸ(x0000.net)
学術の巨大掲示板群 - アルファ・ラボ URLﾘﾝｸ(x0000.net)<)
微分幾何学入門
URLﾘﾝｸ(x0000.net)
772:１３２人目の素数さん
20/07/03 01:23:32.12 DYGsod0z.net
新規感染者が増えかかっているがどう見る？
773:１３２人目の素数さん
20/07/03 02:35:31.91 90y63y3Z.net
>>727
西浦氏が再評価される。
774:１３２人目の素数さん
20/07/03 10:13:27.22 hNK4yhyJ.net
ｽﾚﾘﾝｸ(sisou板:133番)-135
すまん寝＠秩序回復・財務省廃止・反財政再建 @sumannne
【新型コロナウィルス】藤井さん、そのグラフ変じゃないですか？
URLﾘﾝｸ(himorjp.blog.fc2.com)
こっちで改めて計算したところ、「自粛要請は滅茶苦茶効果あった」という結論になっちゃったんで、藤井さん、どうにかした方が良いと思うよ。多分専門家はもうツッコんですらくれないだろうし。
ちなみにデータソースは感染研に1か月遅れですがありましたので、怪しげなデータを使っているわけではありません。
むしろ藤井さんが使ってるデータソースが不明。
藤井さんの統計のいじり方で一番分からんのが、対数グラフを取るでもなく、普通に差分するでもなく、「対前日比」という謎の計算で「速度」「加速度」を名乗ることなんだよなぁ。
この程度のオーダーで対数グラフ取る意味が分からんし、実際藤井さんは対数取ってるわけじゃないので、そうであれば変化率も加速度も差分しないと正しく見えないと思うが。
対数グラフを取る、もしくは対前日比でモノを言うのは、日ごとに指数増大していることがほぼ明らかな状況でなら分かるが、少なくともデータ上はそうではないから適切な計算方法とはちょっと思えない。
藤井さんの記事の「加速度」、縦軸を見てもらえばわかるが上下が0.01の幅で動いてる。
計算は「対前日比」でやってるのでほとんど0ですよこれは。
数値計算で、厳密解が0のものが数値誤差でちょっと荒れてるのを拡大顕微鏡で見て騒いでるような感じがする。
(実際、指数増大しているわけでないものを「比」で取ったら本質的に比率1しか出ません。)
極めて不適切だと思う。
まぁ、今回はさすがに時系列データに相関分析かけたりはしなかったから、その点は成長してるんじゃないっすか。
「対前日比」がいかにおかしいかの例題。
2次多項式のオーダーで増大する場合でも、「対前日比」を取っちゃうと増減が全然分からない。
URLﾘﾝｸ(pbs.twimg.com)
URLﾘﾝｸ(pbs.twimg.com)
775:１３２人目の素数さん
20/07/03 10:13:49.21 GtC8ue70.net
西浦モデルはもう使われないんじゃないの？
とりあえず夜の飲食業は営業停止で。
776:１３２人目の素数さん
20/07/03 10:16:46 GtC8ue70.net
>>729
藤井聡太が迷惑しそうな名前のオヤジだなｗ
もう藤井恥
777:に改名したほうがよろしい。
778:１３２人目の素数さん
20/07/03 21:04:53.02 I3C9PAZr.net
数学系でかつ鉄道趣味のある人で通勤電車の中にどのくらい元気な感染者が存在しうるか
なんて計算してる人っているんだろうか？ （元気な感染者＝感染力の強い時期の人）
あるいはそんな研究あるんだろうか？ -素朴な疑問です 他意はありません
779:１３２人目の素数さん
20/07/03 21:52:40 UvmxXx4m.net
URLﾘﾝｸ(forbesjapan.com)
上の記事、新型インフルエンザで申し訳ないが、感染研が2011年に電車に1人乗ってただけでも感染爆発するというシミュレーション出してる。
780:１３２人目の素数さん
20/07/03 22:04:30.84 GtC8ue70.net
コロナの飛沫感染に関しては、マスクをしてるだけで感染力は
大幅に下がるんでねぇの？
インフルエンザとの違いはそこだね。
781:１３２人目の素数さん
20/07/03 22:05:15.79 GtC8ue70.net
したがって、電車に乗る人はマスク着用を義務付けるべきなんだよな。
782:１３２人目の素数さん
20/07/03 22:45:08.73 90y63y3Z.net
>>733
今は削除されているのでWebarchiveにある。
曰く
首都圏の鉄道に新型インフルエンザを発症した人が1人乗ったと仮定して、まったく対策を取らなかった場合、重症軽症を含めた感染者数は1週間で12万人に拡大するというシミュレーションがあります。
URLﾘﾝｸ(web.archive.org)
783:１３２人目の素数さん
20/07/03 22:47:39.63 90y63y3Z.net
コロナが始まって日本人はマスクを買い増したが、アメリカ人が買い増したのはライフル銃の弾。
マスクに関してはこれが面白かった。
A modelling framework to assess the likely effectiveness of facemasks in combination with ‘lock-down’ in managing the COVID-19 pandemic
URLﾘﾝｸ(royalsocietypublishing.org)
Identifying airborne transmission as the dominant route for the spread of COVID-19
URLﾘﾝｸ(www.pnas.org)
URLﾘﾝｸ(www.pnas.org)
飛沫の可視化
URLﾘﾝｸ(youtu.be)
784:１３２人目の素数さん
20/07/03 22:56:23.10 90y63y3Z.net
>>729
ここのK値というのも、意味が理解できなかった。
ｽﾚﾘﾝｸ(newsplus板)
「（感染増の）波のとらえ方に、K値を採用しようと思っています。中野先生のおっしゃるK値モデルは、
感染の収束速度を計算して感染者数を想定するうえでは、かなり正確だと思う」
K値について詳しい説明はこちらで確認
URLﾘﾝｸ(www.dailyshincho.jp)
URLﾘﾝｸ(www.dailyshincho.jp)
785:１３２人目の素数さん
20/07/03 23:50:48 UvmxXx4m.net
>>738
K値ってのは簡単に言うと、適当に数値の比をとったもの。疫学的な意味はない。
786:１３２人目の素数さん
20/07/04 07:13:01 kFD5yPDP.net
>>739
直線でに線形回帰みたいなもの？
787:１３２人目の素数さん
20/07/04 07:55:47.68 2sU7d1SU.net
>>740
適当に線引いたらそれっぽい数字が出たという数値
788:１３２人目の素数さん
20/07/04 08:21:47.85 kFD5yPDP.net
こういうパラメータも既知の分布によく当て嵌まるかで決定なのだろうな (>681参照)
Distributional fits to key COVID-19 distributions.
URLﾘﾝｸ(i.imgur.com)
URLﾘﾝｸ(www.medrxiv.org)
789:１３２人目の素数さん
20/07/04 13:22:16.78 Y1vg+l8y.net
K値ってのは１週間の感染者数を総感染者数で割っただけなんだが、そうすると既に大流行が終わったところに2波目が来ても不当に小さい値になってしまう。
そこで適当なところで「リセット」出来ることになっているのだが、それが恣意的な時刻で可能なのが最大の問題。
発明者は「ほら、2波目を予測できてるでしょ」と言うのだが、いやいやそれ都合の良いとこでリセットかけてるだけやんとしか。
790:１３２人目の素数さん
20/07/04 20:49:26.54 kFD5yPDP.net
>>716
URLﾘﾝｸ(www.youtube.com)
15分辺りで２キットの検査の違いを解説している。
791:１３２人目の素数さん
20/07/04 21:33:00.37 YBorZJOe.net
>>733
ありがとうございます
コロナの現状（マスクあり）でも誰か
やっててもよさそうなもんですね
不安を煽るとかで規制あるのかな
792:１３２人目の素数さん
20/07/05 08:07:43.49 4xu+FlwC.net
>>743
解説ありがとうございます。
都合よくあう点を選んで接線を引いて予測どおりというみたいなものという理解でいいですか？