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

1 名前：１３２人目の素数さん [2020/02/29(土) 02:18:41 ID:twdO677Q.net] 東大数学科卒の元官僚はこう分析してるが、お前らはどうなると思う？ www.zakzak.co.jp/soc/news/200220/dom2002200003-n2.html 中国国外感染者の中国国内との比率をみると、 １月２０日の数字公表以降は、０・８〜２・６％で比較的安定している。 これは、新型肺炎の感染者のほとんどは中国国内、それも湖北省に集中しているからだ。 ちなみに中国国外での感染者数は、中国国内の１・１％だ（２月１６日現在）。 本コラムで紹介したが、現時点では、最終的な中国国内の感染者数は２０万人超と筆者は推計している。 となると、中国国外の感染者は数千人程度になるだろう。 中国国外のうち日本の比率は１割弱なので、日本の感染者数は数百人程度であろう。 その場合、死者も数人から１０人程度になるだろう。 こうした推計をすると、今の感染者は氷山の一角だと思われるが、今後の増加ペースはどうなるだろうか。 新型コロナウイルスの検査は簡単に行えるので、今後、日本での感染者数は増えていくだろう。 ある時点ではそれがネズミ算的に増えるかのように思える局面もあるだろうが、 筆者の推計が正しければ、現時点ではせいぜい数百人が一つのメドだ。 492 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 02:30:17 ID:VnQvkZ57.net] Wilsonの式での95％信頼区間幅を1％以内にしたいなら サンプルサイズは38411が必要という計算になった。 Rでの算出プログラム library(binom) x=3 n=312 binom.wilson(x,n) fn <- function(n){ x=0:n l=binom.wilson(x,n)[,5] u=binom.wilson(x,n)[,6] max(u-l) } fn=Vectorize(fn) n=seq(1000,50000,by=1000) plot(n,fn(n)) abline(h=0.01,lty=3) uniroot(function(x,u0=0.01) fn(x)-u0, c(10000,50000)) 493 名前：１３２人目の素数さん [2020/05/07(木) 17:10:06 ID:CHL0/p02.net] >>462 ５０人で０でしょ？ だったら、特異度として保証されるのはせいぜい９８％程度じゃん。 １％の陽性率に対して、それでは不十分。 494 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 18:13:41.63 ID:VnQvkZ57.net] >>466 95％信頼区間の下限境界は 事前分布をJeffreysで > qbeta(0.95,0.5+50,0.5,lower=F) [1] 0.9624989 事前分布を一様分布で > qbeta(0.95,1+50,1,lower=F) [1] 0.942952 495 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 18:52:47.08 ID:VnQvkZ57.net] 特異度の95%の信頼区間の下限値を0.99にするのに必要なサンプルサイズは 事前分布を一様分布で > uniroot(function(x) fn(x)-0.99, c(100,500))$root [1] 297.0728 Jeffreysで > uniroot(function(x) fn(x,0.5,0.5)-0.99, c(100,500))$root [1] 190.8606 fn <- function(x,shape1=1,shape2=1){ qbeta(0.95,shape1 + x, shape2, lower=F) } n=50:500 plot(n,fn(n),type='l', ylab='95%CI.lower') abline(h=0.99,lty=3) uniroot(function(x) fn(x)-0.99, c(100,500))$root uniroot(function(x) fn(x,0.5,0.5)-0.99, c(100,500))$root 496 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 20:13:12 ID:VnQvkZ57.net] https://www.buzzfeed.com/jp/naokoiwanaga/covid-19-antibody-test のデータを使って、不明なものは一様分布（ベータ分布の形状母数(1,1)）に事前分布を設定してMCMCしてみる。 x=c(3,33) n=c(312,1000) m=50 N=length(n) shape1=1 shape2=1 model{ for(i in 1:N){ x[i] ~ dbin(p,n[i]) # 二項分布 } p <- prev*sen+(1-prev)*(1-spc) # 陽性=真陽性+偽陽性 sen ~ dbeta(shape1,shape2) spc ~ dbeta(shape1+m,shape2) prev ~ dbeta(shape1,shape1) } その結果 https://i.imgur.com/B7p825B.png 497 名前：１３２人目の素数さん [2020/05/07(木) 21:49:37 ID:CHL0/p02.net] >>469 なにやってるかよくわからんので、見当外れの指摘かもしれんが、 大阪市大の3/312と神戸大の33/1000ってのは特異度も感度も異なる であろう別種のキットによる検査結果なんで、一緒くたにしちゃ まずいんでないの？ 498 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 22:36:59.75 ID:VnQvkZ57.net] >>470 同一キットじゃないから、ご指摘の通り。 しかも神戸大の方では陰性検体での確認はされていないから、神戸大方の陽性率が高いのは偽陽性を含む可能性もあるね。 大阪市大だけのデータでやってみると。 https://i.imgur.com/njVQtRZ.png 499 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 22:46:09.05 ID:VnQvkZ57.net] >>471 使える情報が少なすぎて信頼区間が広すぎ。 有病率:0-87% キットの感度:0-70% ってコイントス変わらんな。 500 名前：１３２人目の素数さん mailto:sage [2020/05/07(木) 22:48:00.40 ID:VnQvkZ57.net] 結局のところ、断定的な結論は出せないねということを数字で確認しているだけｗ 501 名前：１３２人目の素数さん mailto:sage [2020/05/08(金) 02:04:52 ID:Ugc87SUM.net] >>464 >>465 有難うございます。 二項分布下の母比率の信頼区間推定問題なんですね。 A２：正規分布近似法でなくWilsonの信頼区間法で解くから 　　上側区間≠下側区間となり平均値まわり左右非対称信頼区間と 　　なるのですね。 502 名前：１３２人目の素数さん [2020/05/08(金) 09:15:37 ID:WmDpVhCu.net] 3月の宿題で(1)のみ正解の数弱@shukudai_sujaku 昨年度の大学への数学(大数)での勝率は、 学コンBコースが 1/1 = 100% ， 宿題が 3/10 = 30% でした！ 宿題の勝率が低すぎると思うので、 これからは一層精進していきたいです！ https://twitter.com/shukudai_sujaku (deleted an unsolicited ad) 503 名前：１３２人目の素数さん mailto:sage [2020/05/08(金) 11:37:38 ID:CfFk/Uw1.net] >>474 Rのlibrary binomを使って　binom::confint(3,312)で　>464の出力が得られる 信頼区間をグラフにすると https://i.imgur.com/Mlvv7bB.png 点線は3/312 正規分布近似のasymptotic以外は非対称。 どれを使うべきか？　好きなのを使えばいい、と思う。 但し、値が負になったり１を越えたりするのは不採用の方が賢明だとは思う。 Wilson法は値が０や１に近くても信頼できるという人がいるけど、どうやって検証するのかはよくわからん。 504 名前：１３２人目の素数さん mailto:sage [2020/05/08(金) 12:41:19 ID:CfFk/Uw1.net] >>476 binom::confintはbinom.confint(3,312)の間違い library("binom") binom.confint(3,312) 505 名前：１３２人目の素数さん mailto:sage [2020/05/09(土) 13:22:10 ID:oDT7bFgO.net] レムデシビルで初のRCTが出たんだけど、 https://www.thelancet.com/pdfs/journals/lancet/PIIS0140-6736(20)31022-9.pdf レムデシビルRCTで有意差だせず、症例数不足で検出力不足とか。 確かに、効力がわずからなら検出力不足 > n1=158 > n2=78 > pwr.2p2n.test(h=c(0.8,0.5,0.2),n1=n1,n2=n2,sig.level = 0.05,power=NULL, + alternative = "two.sided") difference of proportion power calculation for binomial distribution (arcsine transformation) h = 0.8, 0.5, 0.2 n1 = 158 n2 = 78 sig.level = 0.05 power = 0.9999336, 0.9508568, 0.3037150 alternative = two.sided NOTE: different sample sizes " Remdesivir group (n=158)Placebo group (n=78)Difference Clinical improvement rates Day 7 4 (3%) 2 (3%) 0.0% (-4.3 to 4.2) Day 14 42 (27%) 18 (23%) 3.5% (-8.1 to 15.1) Day 28 103 (65%) 45 (58%) 7.5% (-5.7 to 20.7) " どの週においても両群で症状改善率は不変 症状改善率はどちらも一様分布を事前分布に仮定して 症状改善率の差δ=:０というモデルとδ≠０というモデルでのベイズファクター（周辺尤度比＝エビデンス比）を求めると > (BF01=d.post/d.prio) [1] 1.01472 なので、どちらのモデルが得られたデータを説明するのに適しているかは判断できず 506 名前：、という結果になった。 検出力不足か、差がないのかには決着つかず。 [] [ここ壊れてます] 507 名前：１３２人目の素数さん mailto:sage [2020/05/09(土) 13:22:52 ID:oDT7bFgO.net] library(pwr) library(rjags) library(polspline) par(bty='l') n1=158 n2=78 pwr.2p2n.test(h=c(0.8,0.5,0.2),n1=n1,n2=n2,sig.level = 0.05,power=NULL, alternative = "two.sided") pwr.2p2n.test(h=NULL,n1=n1,n2=n2,sig.level = 0.05,power=0.80, alternative = "two.sided") s1=c(4,42,103) # improvement with Remdesivir s2=c(2,18,45) # improvement with placebo # s1=103 ; s2=45 # Day 28 N=length(s1) # 3 shape1=1 # prior parameter in dbeta shape2=1 data=list(n1=n1,n2=n2,s1=s2,N=N,shape1=shape1,shape2=shape2) modelString=' model{ # data for(i in 1:N){ s1[i] ~ dbin(p1,n1) s2[i] ~ dbin(p2,n2) } # parameter delta <- p1 - p2 # priors p1 ~ dbeta(shape1, shape2) p2 ~ dbeta(shape1, shape2) # sampling from priors pri.p1 ~ dbeta(shape1, shape2) pri.p2 ~ dbeta(shape1, shape2) pri.delta <- pri.p1 - pri.p2 } ' writeLines(modelString,'tmp.txt') n.iter=1e5 ; thin=1 jagsModel=jags.model('tmp.txt',data,inits = NULL,n.chains=4,n.adapt=n.iter/5) update(jagsModel) codaSamples=coda.samples(jagsModel,c('delta','pri.delta'),n.iter,thin) # plot(codaSamples) coda2summary(codaSamples) js=as.data.frame(as.matrix(codaSamples)) BEST::plotPost(js$delta,compVal=0,xlab=bquote(delta)) fit.post=logspline(js$delta) fit.prio=logspline(js$pri.delta) curve(dlogspline(x, fit.post), -2,2, lty=2, xlab=bquote(delta),ylab='Density') curve(dlogspline(x, fit.prio), add=T) (d.post=dlogspline(0,fit.post)) (d.prio=dlogspline(0,fit.prio)) legend('topright', bty='n', legend=c('δ≠0','δ=0'), lty=1:2) abline(v=0,col=8) points(c(0,0),c(d.post,d.prio),pch=c(1,19)) (BF01=d.post/d.prio) 508 名前：１３２人目の素数さん mailto:sage [2020/05/09(土) 13:25:29 ID:oDT7bFgO.net] 確率分布を図示すると https://i.imgur.com/hfH7gZO.png 509 名前：１３２人目の素数さん [2020/05/09(土) 14:25:05 ID:74hNX8Dr.net] 100例もやって差がでないのなら、はっきり言って効き目なしとかわらん。 むしろ副作用がこわい。 アビガンはどうなんだろうね。なぜか国内のまとまったデータがまだ出てこない。 510 名前：１３２人目の素数さん [2020/05/09(土) 14:26:34 ID:74hNX8Dr.net] まあ、軽症者が８割とかだと、早期投与の効き目を判定するには相当の 人数に使わないとわかんないだろうけど。 511 名前：１３２人目の素数さん mailto:sage [2020/05/09(土) 15:22:40 ID:oDT7bFgO.net] 死亡率に関してはベイズファクターで推論できた。 Remdesivir group (n=158) Placebo group (n=78) Difference Day 28 mortality 22 (14%) 10 (13%) 1.1(-8.1% to 10.3%)" 事前確率分布をどうするかだが、なんの情報もないので 実薬群も対照群も一様分布(もしくはJeffreys分布)とする。 死亡率の差δの事前分布と事後分布の確率分布曲線を描いてδ=0での確率密度比(Savage-Dickey density ratio)が ベイズファクター（周辺尤度の比)になる。　 これをJAGSを用いたプログラムで計算すると、一様分布のとき8.34 Jeffrey分布のとき6.40となる。、 10を超えないいのでまあ、中程度の確信で検出力不足によるのではなく、もともと差がないのであろうと推論できる。 サンプルサイズを増やしてもレムデシビルは軽症・早期投与でも致死率を改善しないであろうと予想。 512 名前：１３２人目の素数さん mailto:sage [2020/05/09(土) 17:15:11 ID:oDT7bFgO.net] >どの週においても両群で症状改善率は不変 この前提はおかしいので週ごとに症状改善率が変わるようにモデルを変更 model{ # data for(i in 1:N){ s1[i] ~ dbin(p1[i],n1) s2[i] ~ dbin(p2[i],n2) } # parameter alpha <- delta*sigma for(i in 1:N){ p1[i] <- phi(x1[i]) # probit:pnorm(x) p2[i] <- phi(x2[i]) # p1[i] <- ilogit(x1[i]) # logit:exp(x)/(1+exp(x)) # p2[i] <- ilogit(x2[i]) } # model for(i in 1:N){ x1[i] ~ dnorm(mu + alpha/2, pow(sigma,-2)) # alpha:difference of mean x2[i] ~ dnorm(mu - alpha/2, pow(sigma,-2)) } # priors delta ~ dnorm(0,1) mu ~ dnorm(0,1) sigma ~ dt(0,1,1)T(0,) } ベイズファクターは > (BF01=d.post/d.prio) [1] 0.997627 513 名前：4 症状改善率に差がないというモデルも差があるというモデルも周辺尤度（エビデンス）はほぼ同じという結果。 [] [ここ壊れてます] 514 名前：１３２人目の素数さん [2020/05/09(土) 20:41:59 ID:74hNX8Dr.net] あ、 >>482はレムデシじゃなくてアビガンについての話だかんね。 レムデシは重症者向けらしいけど、たぶん駄目だろ。 515 名前：１３２人目の素数さん [2020/05/09(土) 21:35:34 ID:Fh34v+Vz.net] 【政府が隠す３．１１】　地震波形が、核実験と同じ rio2016.5ch.net/test/read.cgi/earth/1575170033/l50 sssp://o.5ch.net/1nmad.png 516 名前：１３２人目の素数さん [2020/05/11(月) 00:31:07 ID:mGn48zmS.net] 5/10夜『NHKスペシャル』「新型コロナウイルス　出口戦略は」で 厚労省新型コロナ対策班西浦博北大教授と氏が分析中のノートPC画面が 映されていた。氏の使用解析ソフトウェア名は何というのでしょうか？ 517 名前：１３２人目の素数さん mailto:sage [2020/05/11(月) 01:36:00 ID:P7wrhagU.net] 西浦先生はRじゃないかな。 今週ぐらいに、実効再生産数を算出するRのコードを公開できるかもと言ってた。 火曜にはニコ生があるから、そこでたっぷりと聞けるはず。 https://sp.live.nicovideo.jp/watch/lv325833316 518 名前：１３２人目の素数さん mailto:sage [2020/05/11(月) 03:07:28 ID:6TI8KNr5.net] >>488 >11に数式があるけど 519 名前：１３２人目の素数さん mailto:sage [2020/05/11(月) 03:58:21.10 ID:v4rWY6J/.net] >>489 いや、与えられた何日分かの新規患者数から再生産数とか回復率とかを推定するための数式じゃないの？ 520 名前：１３２人目の素数さん mailto:sage [2020/05/11(月) 06:53:45 ID:ruo6tAnf.net] >>490 SEIRより単純なReed-Frostモデルだよな https://en.m.wikipedia.org/wiki/Reed%E2%80%93Frost_model 521 名前：１３２人目の素数さん mailto:sage [2020/05/11(月) 09:14:46.53 ID:P7wrhagU.net] >>489 自分で計算してみた人ならわかるけど、日々の再生算数を出すにはA「感染してから他人にうつすまでの日数の分布」が必要。西浦先生はこれをまず求めていて、その結果は割合と世界中で使われている。 加えてB「感染してから発表の数字に載るまでの日数の分布」も必要で、発表の数字をBで逆畳み込み積分してC「感染推定日毎の数字」に変換し、CをAで畳み込み積分したものでCを割れば日々のRtが求められる。 とは言うものの、そのままやると誤差で数字が発散するし、逆畳み込みは一意ではないのでどう推定するかは色々と考えなくてはならない。多分日本独自の事情もコードに入れなきゃいけない。大変だろうなあとは思う。 522 名前：１３２人目の素数さん mailto:sage [2020/05/11(月) 11:02:15.57 ID:+2HA8Yn9.net] >>492 レスありがとうございます。 確かにSEIRモデルやSIRIモデルだと 「感染してから発表の数字に載るまでの日数の分布」って考慮されてないね。 523 名前：１３２人目の素数さん mailto:sage [2020/05/12(火) 20:50:54 ID:8HqG++pl.net] >>488 公開されたぞー。StanとRだ。 https://nbviewer.jupyter.org/github/contactmodel/COVID19-Japan-Reff/blob/master/scripts/C.%20Calculating%20the%20Rt%20in%20Stan.ipynb 524 名前：１３２人目の素数さん mailto:sage [2020/05/12(火) 22:32:45.27 ID:9+i4T0w3.net] ウイルス感染のモデルも必要だと思うけど 社 525 名前：会経済活動の影響や 自粛解除した後の死者増加や 自粛継続した場合の法人の倒産や 治療薬やワクチンができた後の社会の経済や文化面の回復 そこまで考慮してどういう政策を選択したら 損失を最小化できるかという問題は解ける？ 制約付きの最小化問題だと思う 数値化するのが難しい考慮要素もあるからその辺は仮の値とかで [] [ここ壊れてます] 526 名前：１３２人目の素数さん mailto:sage [2020/05/13(水) 17:23:13 ID:PBQsFM+W.net] >>495 モデルを精密複雑化すればするほど、必要なパラメータが増えて それを弱情報事前分布にすると事後分布の信頼区間が広がるんだなぁと思う。 527 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 08:00:23 ID:kec+XbRE.net] >>494 これにある、 JapaneseDataCOVID19 ってどこかからダウンロードできるんだろうか？ 528 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 09:54:38 ID:IN0uokww.net] >>497 https://github.com/contactmodel/COVID19-Japan-Reff/tree/master/data 529 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 11:47:52 ID:kec+XbRE.net] >>498 ありがとうございます。 530 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 11:51:11 ID:kec+XbRE.net] 親ディレクトリ（フォルダ）探せばよかったんだね。 https://nbviewer.jupyter.org/github/contactmodel/COVID19-Japan-Reff/tree/master/ 531 名前：１３２人目の素数さん [2020/05/14(木) 13:08:10.58 ID:gLQVRit5.net] >>495 >>496 パラメータを増やせば増やすほど、どんな仮説も否定できなくなる。 実験系の学界でよく使う手口だよね。 おや、誰か来たようだ… 532 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 13:58:09 ID:kec+XbRE.net] >>494 これを走らせてみたい人いますか？ > datestar = as.Date("2020-05-10") > datemin = as.Date("2019-12-25") # particular choice > (tstar = as.numeric(datestar - datemin)) [1] 137 > (K = nrow(df_cases)) # 147 [1] 147 なので、 data { int<lower = 1> K; //number of days vector<lower = 0>[K] imported_backproj; vector<lower = 0>[K] domestic_backproj; int<lower = K> upper_bound; の int<lower = K> upper_bound; だ Error in new_CppObject_xp(fields$.module, fields$.pointer, ...) : Exception: modelf2b01ab215_fit_infection_namespace::modelf2b01ab215_fit_infection: upper_bound is 137, but must be greater than or equal to 147 (in 'modelf2b01ab215_fit_infection' at line 53) upper_bound 137の下限を K 147にしているからみたい。 533 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 14:10:15 ID:kec+XbRE.net] upper_boundの制限を外して data { // int<lower = K> upper_bound; int upper_bound; 再生数の平均値を以下の出すブロックを加えて走らせてみた。 transformed parameters{ real mean_Rt; real mean_Rt_adj; mean_Rt = mean(Rt); mean_Rt_adj = mean(Rt_adj); } その結果、 Inference for Stan model: fit_infection. 2 chains, each with iter=10000; warmup=2000; thin=5; post-warmup draws per chain=1600, total post-warmup draws=3200. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat mean_Rt 1.79 0 0.07 1.66 1.74 1.79 1.83 1.92 3122 1 mean_Rt_adj 1.79 0 0.07 1.66 1.74 1.79 1.84 1.93 3123 1 534 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 16:35:35 ID:kec+XbRE.net] 再生算数を０〜１０人の一様分布にすると、収束しない。 > print(fit_u) Inference for Stan model: fit_infection_u. 4 chains, each with iter=10000; warmup=5000; thin=5; post-warmup draws per chain=1000, total post-warmup draws=4000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat mean_Rt 2.15 0.05 0.10 1.98 2.06 2.16 2.23 2.30 3 1.99 mean_Rt_adj 2.13 0.07 0.11 1.93 2.06 2.13 2.21 2.36 3 2.27 lp__ -789.19 7.18 14.43 -820.19 -797.69 -787.44 -779.42 -762.53 4 1.89 traceplotやchainの分布はこんな感じ、 https://i.imgur.com/z0RL1KW.png https://i.imgur.com/VRrrsKw.png 535 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 17:00:26 ID:1Kd4DHQt.net] >>501 否定できないけど実証もできないことになるのでは？ 単純なモデルだと現実の問題を解決できない事もあるし 536 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 17:23:23.34 ID:kec+XbRE.net] >>494 モデルで再生産数の事前分布は　平均2.5　標準偏差2.0の正規分布に設定されていたので 平均と標準偏差を変化させて、再生産数の事後分布を描出してみた。 かなり、事前分布の影響を受けるみたい。 https://i.imgur.com/OwqsFC1.png 537 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 17:26:03.64 ID:kec+XbRE.net] どうも、こういう境地だなぁ。 断定的な結論は出せないということを数字で確認しているだけｗ 538 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 17:30:34 ID:50GljeDM.net] 基本再生算数はウイルスの特徴で決まるもので 実行再生算数は更に環境や人の行動の影響で変化するものだと理解している 実行再生算数を減らすようにするには どう行動したらいいか かつ経済活動を出来るだけ下げずに 539 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 18:25:32 ID:IN0uokww.net] 西浦モデル、一度感染日ごとの人数を最尤法で確定して、そこからベイズを回してRtを求めてるんだよなあ。 本来は、日々の感染者数も確率的に揺れがあるはず。だからRtの誤差幅は発表よりも多く見積もるべきかもしれん。 540 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 18:39:49 ID:IN0uokww.net] コードを解説してる人 mikuhatsune.hatenadiary.com/entry/2020/05/13/205824 541 名前：１３２人目の素数さん mailto:sage [2020/05/14(木) 19:09:54 ID:kec+XbRE.net] 再生産数の事前分布を色々かえて事後分布を出してみた。 https://i.imgur.com/CT2TRbg.png 542 名前：１３２人目の素数さん mailto:sage [2020/05/15(金) 08:46:23 ID:qjCTzgxb.net] >>504 切断分布だと収束しないみたいなので、 一様分布[0,10]に近そうな正規分布[5,3] https://i.imgur.com/h8vMZUM.png を事前分布にして走らせてみた。 https://i.imgur.com/O5s0Y8a.png 543 名前：１３２人目の素数さん mailto:sage [2020/05/15(金) 09:56:20.25 ID:qjCTzgxb.net] 【新型コロナ】 東京0.6％、東北6県0.4％陽性・・・抗体検査1000人実施 [影のたけし軍団★] https://asahi.5ch.net/test/read.cgi/newsplus/1589502801/ 複数の検査キットの性能評価と感染状況の確認が目的でしたが、東京都で献血した500人のうち3人、 東北6県で献血した500人のうち2人がいずれかの検査キットで陽性と判定されました。 満員電車など人との接触の多い東京とそうでない東北で陽性率に有意差はあるか？ 544 名前：１３２人目の素数さん [2020/05/15(金) 12:23:12 ID:GeDvuMme.net] >>513 検査キットの特異度が99.5%以上あることが検証されました。 ってだけの話じゃね？ まぁ、東北６県と東京が同じってことは、特異度はそれより 高くはないってことだろうね。 545 名前：１３２人目の素数さん [2020/05/15(金) 12:26:48 ID:GeDvuMme.net] 言い換えると、真の陽性率はわかりません、ってこと。 546 名前：１３２人目の素数さん mailto:sage [2020/05/15(金) 16:53:27 ID:qjCTzgxb.net] >>513 陽性率の確率分布を一様分布にすると事後分布は https://i.imgur.com/YF6m869.png なるけど、重なりの部分の面積が差がないことの度合いを示していると考えていいかな？ 547 名前：１３２人目の素数さん mailto:sage [2020/05/15(金) 18:43:40 ID: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 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 07:59:05 ID:/d9XIHEO.net] 再生産数を計算するRのプログラムあったんだな https://www.rdocumentation.org/packages/EpiEstim/versions/2.2-1/topics/estimate_R 550 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 08:16:33 ID:28hwAVWs.net] >>509 後者の分析、RStan本の訳者が今やっているようだ。 https://twitter.com/hankagosa/status/1261430169283125248 (deleted an unsolicited ad) 551 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 08:17:55 ID:28hwAVWs.net] >>519 https://twitter.com/hankagosa/status/1261426113374416897 (deleted an unsolicited ad) 552 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 08:36:22 ID:/d9XIHEO.net] >>520 アヒル本の著者だね。俺も読んだ。 553 名前：１３２人目の素数さん [2020/05/16(土) 12:36:46 ID:cGFgpNwX.net] 4/2の時点で感染者数を６０００くらいと見積もってるね。＞アヒル本の人 実数の三倍程度。いい線かもしれない。 554 名前：１３２人目の素数さん [2020/05/16(土) 13:39:32 ID:XSllT5Kn.net] https://github.com/contactmodel/COVID19-Japan-Reff/blob/master/nishiura_Rt%E4%BC%9A%E8%AD%B0_12May2020.pdf 555 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 14:29:29 ID:BxGcLzV+.net] int<lower = K> upper_bound; ↓ int　upper_bound; にしてもエラーがでる。 stan_dataで upper_bound = 147 にすると動くけど、何をやってんのか自分でもよくわからん。 556 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 16:47:14 ID:UBjBTXVe.net] アヒル本ってなんですか？ 557 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 17:55:00 ID:BxGcLzV+.net] >>525 名著として名高いStanの入門書 StanとRでベイズ統計モデリング (Wonderful R) (日本語) 単行本 ? 2016/10/25 表紙の色がアヒルのくちばしの色に似ているかららしい。 558 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 17:56:04 ID:BxGcLzV+.net] >>524 自己解決 select ↓ dplyr::select でちゃんと動作した 559 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 19:16:41 ID:5/7mBSAW.net] >>526 thx 560 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 19:28:52 ID:BxGcLzV+.net] R に　EpiEstimというパッケージがあって、再生産数を算出する関数が搭載されている。 結局、infecterとinfecteeが発症するまでの期間serial intervalの分布をどう設定するかで結果が変わるみたいだなぁ。 Rのヘルプファイルを解読中。 Rのヘルプファイルは不親切設計で有名（理解できている人の備忘録みたいな性格だから）。 561 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 20:49:50 ID:5T7nYzW3.net] ここに居る人達って何者なの 学部の知識超えてるよね 統計でご飯食べてる人たち？ 562 名前：１３２人目の素数さん [2020/05/16(土) 21:03:16 ID:YPW1s8+S.net] しかし、なんでも揃ってるなRって 563 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 22:57:43 ID: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 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 00:04:44 ID:u5AEq3c8.net] >>532 ワイブル分布の平均 4.830129と標準偏差2.222765をそのまま正規分布のパラメータに使って、グラフを重ねてみる。 https://i.imgur.com/TnzGwWx.png ワイブル分布で発生させた乱数をワイブルでフィットさせて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 の分類がよくわからんが、全部足してグラフを描いてみた。 https://i.imgur.com/rKBeWgq.png 566 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 00:18:59 ID:u5AEq3c8.net] 別の論文だと対数正規分布がフィットすると西浦氏は記載している。 serila interval : infector と infecteeの発症間隔 https://www.ijidonline.com/article/S1201-9712(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 https://i.imgur.com/Uzg6u84.png 点線が2項分布で実線がワイブル分布 567 名前：１３２人目の素数さん [2020/05/17(日) 07:09:05 ID:/ApxcPC4.net] “頭脳王”東大生・河野玄斗 基本的な数学でコロナウイルス検査を全員にしても意味がないことを証明してみた https://www.youtube.com/watch?v=jMIScCb04qs 568 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 07:28:25 ID:mpH4uuPx.net] 実効再生産数を感染者数の推移から推定する数理的原理をキチンと解説した本、または論文誰か知りませんか？ 勉強してみたい。 569 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 07:33:35 ID:u5AEq3c8.net] >>536 RのEpiEstimの著書の論文なんかどうでしょう？ A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3816335/ 570 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 07:55:41 ID: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 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 08:10:41 ID:u5AEq3c8.net] >>529 感染者に０が続くと再生産数の信頼区間幅がどんどん広くなってくる。 まあ、疫病用のソフトウェアと理解しておこう。 https://i.imgur.com/QbwNydN.png 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 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 09:00:46 ID: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. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3816335/pdf/kwt133.pdf 573 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 09:34:29 ID: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 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 12:07:27 ID:S4iLXC97.net] >>537 ありがとう。 読んでみます。 でも"a new frame work"っていうならもっと大本のこの道の研究者なら「Rtの推定はコレ」みたいな標準的なものがあるんですかね？ どなたかご存知ないですか？ 575 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 12:11:01 ID:S4iLXC97.net] >>540 おぉ、こっちがRで採用されてる推定法なんですね。 読んでみます。 576 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 12:16:27 ID:S4iLXC97.net] と思ったら>>537と>>540は同じかorz 577 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 15:31:13 ID: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 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 20:41:31 ID:/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 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 05:56:32 ID:hVD/4gbI.net] >>546 >523氏が挙げてくれた 実効再生産数とその周辺　に　記述があったが、publicと呼べるのかどうかは門外漢なので知らない。 https://github.com/contactmodel/COVID19-Japan-Reff/blob/master/nishiura_Rt%E4%BC%9A%E8%AD%B0_12May2020.pdf こんな記述があった >> 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 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 06:16:59 ID: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 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 11:15:46 ID: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 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 13:21:08 ID:1P7V5xJn.net] 職場でも最初に発症した人が感染源のように扱われるけど 潜伏期間の分布を考えたら断定はできない。 COVID19の潜伏期間の論文 https://www.nejm.org/doi/full/10.1056/NEJMoa2001316 から、 潜伏期間は対数正規分布で近似できてそのパラメータは #--- 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 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 13:25:20 ID:1P7V5xJn.net] >>549 Highest Density Probability Intervalを求めればいいんじゃないの？ 586 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 14:04:29 ID:jxOnYm2h.net] >>551 何ですかそれは？ 587 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 14:08:09 ID:1P7V5xJn.net] >>550 正確にはキャバ嬢がXより先に感染していた確率だな。 588 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 14:08:59 ID:1P7V5xJn.net] >>552 確率分布が対称じゃないときの信頼区間 589 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 14:13:34 ID:1P7V5xJn.net] >>554 こんな感じ。https://i.imgur.com/C8jOPlx.jpg 590 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 14:26:48 ID:jxOnYm2h.net] >>554 ただ、それだとそもそもモードに近いとこをやってます。 信頼区間は密度関数を横に切るのではなく両裾を縦に切ってハズレ部分が1-λになるようにするので少しイメージが違うしきがします。 モードなのかメジアンなのかの違いです。 いずれにせよ、こうやればいいという拡張のための俺様ルールを設定するのはいくらでもできますが、統計の話なのでそんな俺様ルールについて語っても意味ありません。 591 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 15:34:14 ID: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 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 15:36:39 ID:jxOnYm2h.net] >>557 ？？？

---

---

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

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