数学 統計に詳しい人が語るコロナウイルス
■ このスレッドは過去ログ倉庫に格納されています
東大数学科卒の元官僚はこう分析してるが、お前らはどうなると思う? http://www.zakzak.co.jp/soc/news/200220/dom2002200003-n2.html 中国国外感染者の中国国内との比率をみると、 1月20日の数字公表以降は、0・8〜2・6%で比較的安定している。 これは、新型肺炎の感染者のほとんどは中国国内、それも湖北省に集中しているからだ。 ちなみに中国国外での感染者数は、中国国内の1・1%だ(2月16日現在)。 本コラムで紹介したが、現時点では、最終的な中国国内の感染者数は20万人超と筆者は推計している。 となると、中国国外の感染者は数千人程度になるだろう。 中国国外のうち日本の比率は1割弱なので、日本の感染者数は数百人程度であろう。 その場合、死者も数人から10人程度になるだろう。 こうした推計をすると、今の感染者は氷山の一角だと思われるが、今後の増加ペースはどうなるだろうか。 新型コロナウイルスの検査は簡単に行えるので、今後、日本での感染者数は増えていくだろう。 ある時点ではそれがネズミ算的に増えるかのように思える局面もあるだろうが、 筆者の推計が正しければ、現時点ではせいぜい数百人が一つのメドだ。 >>454 事後分布を一様分布に設定して、このコピペ小僧の正解率の期待値・最頻値・95%信頼区間を求めよ。 その結果、 https://i.imgur.com/M9SdCSL.png >>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 >>454 このコピペ小僧の正解率に学コンと宿題で差があるかを検定せよ。 事前分布にJefferysを用いた結果、 https://i.imgur.com/2yYagtX.png >>456 感度が5%以下じゃつかいもんにならんわなぁ。インチキで確定。 >>462 >2020年4月中の2日間に同大学の付属病院を、新型コロナウイルス感染症の診療以外で受診した >患者を対象に、そこから無作為に312人を抽出・・・・312 人(年齢中央値 66.5 歳、 >男性:女性=154:158)のうち、3人が陽性であることがわかった。約1%の陽性率だ。 >統計的な誤差を考慮すると、95%の確率で0.33〜2.8%の間に入る・・・・・・・・・・・・ 教えてエロい人 Q1統計的推定を述べるなら「95%の確率で」でなく「信頼率95%で」と書くのが適切ではないか? Q2無作為抽出陽性率3/312=0.96%の母集団陽性率上側/下側信頼区間=0.96-0.33/2.80-0.96と 上側区間≠下側区間で左右非対称分布想定は何故? Q3サンプルサイズ312は過少では?過少に伴う推定誤差加算が必要では? >>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 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)) >>462 50人で0でしょ? だったら、特異度として保証されるのはせいぜい98%程度じゃん。 1%の陽性率に対して、それでは不十分。 >>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 特異度の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 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 >>469 なにやってるかよくわからんので、見当外れの指摘かもしれんが、 大阪市大の3/312と神戸大の33/1000ってのは特異度も感度も異なる であろう別種のキットによる検査結果なんで、一緒くたにしちゃ まずいんでないの? >>470 同一キットじゃないから、ご指摘の通り。 しかも神戸大の方では陰性検体での確認はされていないから、神戸大方の陽性率が高いのは偽陽性を含む可能性もあるね。 大阪市大だけのデータでやってみると。 https://i.imgur.com/njVQtRZ.png >>471 使える情報が少なすぎて信頼区間が広すぎ。 有病率:0-87% キットの感度:0-70% ってコイントス変わらんな。 結局のところ、断定的な結論は出せないねということを数字で確認しているだけw >>464 >>465 有難うございます。 二項分布下の母比率の信頼区間推定問題なんですね。 A2:正規分布近似法でなくWilsonの信頼区間法で解くから 上側区間≠下側区間となり平均値まわり左右非対称信頼区間と なるのですね。 3月の宿題で(1)のみ正解の数弱@shukudai_sujaku 昨年度の大学への数学(大数)での勝率は、 学コンBコースが 1/1 = 100% , 宿題が 3/10 = 30% でした! 宿題の勝率が低すぎると思うので、 これからは一層精進していきたいです! https://twitter.com/shukudai_sujaku https://twitter.com/5chan_nel (5ch newer account) >>474 Rのlibrary binomを使って binom::confint(3,312)で >464の出力が得られる 信頼区間をグラフにすると https://i.imgur.com/Mlvv7bB.png 点線は3/312 正規分布近似のasymptotic以外は非対称。 どれを使うべきか? 好きなのを使えばいい、と思う。 但し、値が負になったり1を越えたりするのは不採用の方が賢明だとは思う。 Wilson法は値が0や1に近くても信頼できるという人がいるけど、どうやって検証するのかはよくわからん。 >>476 binom::confintはbinom.confint(3,312)の間違い library("binom") binom.confint(3,312) レムデシビルで初の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) " どの週においても両群で症状改善率は不変 症状改善率はどちらも一様分布を事前分布に仮定して 症状改善率の差δ=:0というモデルとδ≠0というモデルでのベイズファクター(周辺尤度比=エビデンス比)を求めると > (BF01=d.post/d.prio) [1] 1.01472 なので、どちらのモデルが得られたデータを説明するのに適しているかは判断できず、という結果になった。 検出力不足か、差がないのかには決着つかず。 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) 100例もやって差がでないのなら、はっきり言って効き目なしとかわらん。 むしろ副作用がこわい。 アビガンはどうなんだろうね。なぜか国内のまとまったデータがまだ出てこない。 まあ、軽症者が8割とかだと、早期投与の効き目を判定するには相当の 人数に使わないとわかんないだろうけど。 死亡率に関してはベイズファクターで推論できた。 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を超えないいのでまあ、中程度の確信で検出力不足によるのではなく、もともと差がないのであろうと推論できる。 サンプルサイズを増やしてもレムデシビルは軽症・早期投与でも致死率を改善しないであろうと予想。 >どの週においても両群で症状改善率は不変 この前提はおかしいので週ごとに症状改善率が変わるようにモデルを変更 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 症状改善率に差がないというモデルも差があるというモデルも周辺尤度(エビデンス)はほぼ同じという結果。 あ、 >>482 はレムデシじゃなくてアビガンについての話だかんね。 レムデシは重症者向けらしいけど、たぶん駄目だろ。 5/10夜『NHKスペシャル』「新型コロナウイルス 出口戦略は」で 厚労省新型コロナ対策班西浦博北大教授と氏が分析中のノートPC画面が 映されていた。氏の使用解析ソフトウェア名は何というのでしょうか? 西浦先生はRじゃないかな。 今週ぐらいに、実効再生産数を算出するRのコードを公開できるかもと言ってた。 火曜にはニコ生があるから、そこでたっぷりと聞けるはず。 https://sp.live.nicovideo.jp/watch/lv325833316 >>489 いや、与えられた何日分かの新規患者数から再生産数とか回復率とかを推定するための数式じゃないの? >>489 自分で計算してみた人ならわかるけど、日々の再生算数を出すにはA「感染してから他人にうつすまでの日数の分布」が必要。西浦先生はこれをまず求めていて、その結果は割合と世界中で使われている。 加えてB「感染してから発表の数字に載るまでの日数の分布」も必要で、発表の数字をBで逆畳み込み積分してC「感染推定日毎の数字」に変換し、CをAで畳み込み積分したものでCを割れば日々のRtが求められる。 とは言うものの、そのままやると誤差で数字が発散するし、逆畳み込みは一意ではないのでどう推定するかは色々と考えなくてはならない。多分日本独自の事情もコードに入れなきゃいけない。大変だろうなあとは思う。 >>492 レスありがとうございます。 確かにSEIRモデルやSIRIモデルだと 「感染してから発表の数字に載るまでの日数の分布」って考慮されてないね。 ウイルス感染のモデルも必要だと思うけど 社会経済活動の影響や 自粛解除した後の死者増加や 自粛継続した場合の法人の倒産や 治療薬やワクチンができた後の社会の経済や文化面の回復 そこまで考慮してどういう政策を選択したら 損失を最小化できるかという問題は解ける? 制約付きの最小化問題だと思う 数値化するのが難しい考慮要素もあるからその辺は仮の値とかで >>495 モデルを精密複雑化すればするほど、必要なパラメータが増えて それを弱情報事前分布にすると事後分布の信頼区間が広がるんだなぁと思う。 >>494 これにある、 JapaneseDataCOVID19 ってどこかからダウンロードできるんだろうか? >>495 >>496 パラメータを増やせば増やすほど、どんな仮説も否定できなくなる。 実験系の学界でよく使う手口だよね。 おや、誰か来たようだ… >>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にしているからみたい。 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 再生算数を0〜10人の一様分布にすると、収束しない。 > 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 >>501 否定できないけど実証もできないことになるのでは? 単純なモデルだと現実の問題を解決できない事もあるし >>494 モデルで再生産数の事前分布は 平均2.5 標準偏差2.0の正規分布に設定されていたので 平均と標準偏差を変化させて、再生産数の事後分布を描出してみた。 かなり、事前分布の影響を受けるみたい。 https://i.imgur.com/OwqsFC1.png どうも、こういう境地だなぁ。 断定的な結論は出せないということを数字で確認しているだけw 基本再生算数はウイルスの特徴で決まるもので 実行再生算数は更に環境や人の行動の影響で変化するものだと理解している 実行再生算数を減らすようにするには どう行動したらいいか かつ経済活動を出来るだけ下げずに 西浦モデル、一度感染日ごとの人数を最尤法で確定して、そこからベイズを回してRtを求めてるんだよなあ。 本来は、日々の感染者数も確率的に揺れがあるはず。だからRtの誤差幅は発表よりも多く見積もるべきかもしれん。 【新型コロナ】 東京0.6%、東北6県0.4%陽性・・・抗体検査1000人実施 [影のたけし軍団★] https://asahi.5ch.net/test/read.cgi/newsplus/1589502801/ 複数の検査キットの性能評価と感染状況の確認が目的でしたが、東京都で献血した500人のうち3人、 東北6県で献血した500人のうち2人がいずれかの検査キットで陽性と判定されました。 満員電車など人との接触の多い東京とそうでない東北で陽性率に有意差はあるか? >>513 検査キットの特異度が99.5%以上あることが検証されました。 ってだけの話じゃね? まぁ、東北6県と東京が同じってことは、特異度はそれより 高くはないってことだろうね。 言い換えると、真の陽性率はわかりません、ってこと。 >>513 陽性率の確率分布を一様分布にすると事後分布は https://i.imgur.com/YF6m869.png なるけど、重なりの部分の面積が差がないことの度合いを示していると考えていいかな? これは、味気がないな。 > 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 Ratio: 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 ------------------------------------------------------ 4/2の時点で感染者数を6000くらいと見積もってるね。>アヒル本の人 実数の三倍程度。いい線かもしれない。 int<lower = K> upper_bound; ↓ int upper_bound; にしてもエラーがでる。 stan_dataで upper_bound = 147 にすると動くけど、何をやってんのか自分でもよくわからん。 >>525 名著として名高いStanの入門書 StanとRでベイズ統計モデリング (Wonderful R) (日本語) 単行本 ? 2016/10/25 表紙の色がアヒルのくちばしの色に似ているかららしい。 >>524 自己解決 select ↓ dplyr::select でちゃんと動作した R に EpiEstimというパッケージがあって、再生産数を算出する関数が搭載されている。 結局、infecterとinfecteeが発症するまでの期間serial intervalの分布をどう設定するかで結果が変わるみたいだなぁ。 Rのヘルプファイルを解読中。 Rのヘルプファイルは不親切設計で有名(理解できている人の備忘録みたいな性格だから)。 ここに居る人達って何者なの 学部の知識超えてるよね 統計でご飯食べてる人たち? >>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.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 正規分布で近似してもよさそうな感じだな。 >>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 別の論文だと対数正規分布がフィットすると西浦氏は記載している。 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項分布で実線がワイブル分布 “頭脳王”東大生・河野玄斗 基本的な数学でコロナウイルス検査を全員にしても意味がないことを証明してみた https://www.youtube.com/watch?v=jMIScCb04qs 実効再生産数を感染者数の推移から推定する数理的原理をキチンと解説した本、または論文誰か知りませんか? 勉強してみたい。 >>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/ Rで引数なしで関数を実行させようとすると、ソースが表示される。 その関数が呼び出した関数のソースを次々に辿っていけば内部計算がわかるので、そこから原理がつかめるかも(俺には無理なのでしないけど) 内部関数のときは:が3つとか、パッケージ名:::関数、パッケージ名:::関数.default(例、t検定のソース表示は stats:::t.test.default )で表示される。 EpiEstim::estimate_R EpiEstim:::process_si_data EpiEstim:::process_config_si_from_data coarseDataTools::dic.fit.mcmc >>529 感染者に0が続くと再生産数の信頼区間幅がどんどん広くなってくる。 まあ、疫病用のソフトウェアと理解しておこう。 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))) >>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 基本コンセプトはこれだろうな。 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. 公表された感染者数で計算するために発症から診断/公表までの時間も考慮した点が西浦モデルの優れた点ではないかと感じている。 >>537 ありがとう。 読んでみます。 でも"a new frame work"っていうならもっと大本のこの道の研究者なら「Rtの推定はコレ」みたいな標準的なものがあるんですかね? どなたかご存知ないですか? >>540 おぉ、こっちがRで採用されてる推定法なんですね。 読んでみます。 結局、これが核心部分 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 ガンマ分布で近似するのが定石らしい。 >>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 知ってる人いたら教えてもらおうと。 まだ論文読んでないので書いてあるかもですけど。 >>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 (cはcohortの頭文字)として言及されているから、まあ、理論疫学者の間ではcommonな定義なんだろうと思う。 >>546 >4/1〜4/30までのデータx1,x2,‥,x30 x1,...,x30が感染者数なら非負整数で与えられるから、CIを考える必要があるかな? 集計ミスで19人が117人とかあったらしいから、信頼区間を考える必要があるのかもしれないとは思うけど。 確かに、発症日に関してはファジーな部分があるとは思う。 いつから熱がありましたか?味がわからなくなりましたか?と問われても1〜3日位の幅はでるだろう。 RのEpiEstim::estimate_Rをヘルプファイルの実例を使って走らせてみた。 serial intervalの分布をデータから推定させるのに infecterの発症日の下限・上限 infecteeの発症日の下限・上限 を設定する項目が(EL,ER,SL,SR)があった。 このデータに合致する分布関数(ガンマ、ワイブル、対数正規分布が指定可能)のパラメータを算出して計算させているみたい。 西浦モデルではワイブル分布で固定。 潜伏期間にも変動があるから、誰がinfecterで誰がinfecteeかを決定するのも難しいだろうとは思う。 後から発症した人間が感染源というのもありうるし。 >>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の中の立方体の体積がλになるように真ん中にとるんだろうなぁと思うんですけど。 職場でも最初に発症した人が感染源のように扱われるけど 潜伏期間の分布を考えたら断定はできない。 COVID19の潜伏期間の論文 https://www.nejm.org/doi/full/10.1056/NEJMoa2001316 から、 潜伏期間は対数正規分布で近似できてそのパラメータは #--- incubation period --- # from Li et al NEJM 2020 # lognormal mean = 5.2 ln.par1 = 1.434065 ln.par2 = 0.6612 という。 ある人物Xが新型コロナ肺炎に罹患したとする。 行動調査によって発症前にキャバクラに行っており接客したキャバ嬢がX発症の2日後に発症していたことがわかった。 Xがキャバ嬢から移された確率を求めよ。 >>549 Highest Density Probability Intervalを求めればいいんじゃないの? >>550 正確にはキャバ嬢がXより先に感染していた確率だな。 >>554 ただ、それだとそもそもモードに近いとこをやってます。 信頼区間は密度関数を横に切るのではなく両裾を縦に切ってハズレ部分が1-λになるようにするので少しイメージが違うしきがします。 モードなのかメジアンなのかの違いです。 いずれにせよ、こうやればいいという拡張のための俺様ルールを設定するのはいくらでもできますが、統計の話なのでそんな俺様ルールについて語っても意味ありません。 >>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の方が幅が小さい。 ■ このスレッドは過去ログ倉庫に格納されています
read.cgi ver 07.5.4 2024/05/19 Walang Kapalit ★ | Donguri System Team 5ちゃんねる