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

東大数学科卒の元官僚はこう分析してるが、お前らはどうなると思う？

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

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

こうした推計をすると、今の感染者は氷山の一角だと思われるが、今後の増加ペースはどうなるだろうか。
新型コロナウイルスの検査は簡単に行えるので、今後、日本での感染者数は増えていくだろう。
ある時点ではそれがネズミ算的に増えるかのように思える局面もあるだろうが、
筆者の推計が正しければ、現時点ではせいぜい数百人が一つのメドだ。

抗体検査キットの評価をしたら特異度はみな5/5だったらしいけど、サンプルが
５つだけって意味だとすると、せいぜい90％以上くらいのことしか言えないな。
特異度95%だとしたら、５％が陽性っていわれてもねぇｗ

>>436
事前分布にかなり影響をうけるが、

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

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

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

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

有病率の期待値は2.3%、最頻値は0.31%　少数データなので信頼区間が広い。
有病率が平均値50%の一様分布というのは現実離れした分布だとは思う。

>>439
すべてが陰性に出るイカサマキットなら特異度が100％
>18参照

>55人中5人と147人中7人じゃ、有意差ないだろ。

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

> fn(4)
[1] 0.045678

>>440（追記)
>有病率が平均値50%の一様分布というのは現実離れした分布だとは思う。
オーストリアのデータ0.3%を参考に有病率の事前分布が0.2-0.4％
β(10.865, 3279.34)に相当として、MCMCしてみた。
感度と特異度の事前分布は前回と同じ。

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

さほど、有病率が高いという結論は引き出せないな。

自己流な検定

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

モンテカルロ法で検証

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

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

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

EXCEL VBAソースコード概略

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

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

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

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

Next

β分布と乱数発生で比較。

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

a=0.5 ; b=0.5
r1=5 ; r2=7 ; n1=55 ; n2=147
layout(1)
layout(matrix(1:2,2))
curve(dbeta(x,a+r2,b+n2-r2),0,0.3,bty='l',ann=F,lwd=2)
curve(dbeta(x,a+r1,b+n1-r1),col=2,add=T,lwd=2)
legend('center',bty='n',legend=c('general','medical'),lwd=2,col=1:2)
k=1e7
general=rbeta(k,a+r2,b+n2-r2)
medical=rbeta(k,a+r1,b+n1-r1)
BEST::plotPost(medical/general,compVal = 1)

>>441
すべてが陰性に出るなら陽性は０人だろｗ

>>446
偽陽性率０％だから特異度は１００％

>>447
だから、陽性が何人か出てんだから感度0%のインチキじゃなかろうってこと。

最近の報道

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

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

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

問題

1試薬は1回しか検査できないとして
これがイカサマキットであることを証明する手段はあるか？

3月の宿題で(1)のみ正解の数弱@shukudai_sujaku

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

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

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

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

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

https://twitter.com/shukudai_sujaku
https://twitter.com/5chan_nel (5ch newer account)

高橋洋一(統計数理研究所→大蔵省)
スウェーデンの感染症対策を解説
https://youtu.be/2bsYRuorwMI?t=192

https://toyokeizai.net/sp/visual/tko/covid19/
のデータを使って
P=14895/153581 　#　2020/05/04
PCR検査の感度30-70％のモデルM1と感度70-90％のモデルM2のどちらが信憑性があるか、ベイズファクターで計算してみる。
M1は感度が最頻値60%標準偏差10%、M2は最頻値80%標準偏差10%のベータ分布に設定
特異度はいずれも最頻値95%標準偏差2.5%に設定し、有病率は一様分布を仮定
陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする。
事後確率分布は
https://i.imgur.com/81bK4KE.png
陽性率P=14895/153581=0.09698465での事後確率分布の密度比(Savage-Dickey density ratio)でベイズファクターを出すと
> d1/d2 # Savage-Dickey densiti ratio = BF12
[1] 1.007722

まあ、ちょっぴり、感度30-70％のモデルの方がいいかも、という結果。
陽性数/検査数の時系列データでもあればもう少し差がでるかもしれん。
東京都のデータで計算させようかと思ったが、東京都は検査人数を隠蔽しているので使いものにならない。

訂正

陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする
↓
陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする

3月の宿題で(1)のみ正解の数弱@shukudai_sujaku

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

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

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

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

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

https://twitter.com/shukudai_sujaku
https://twitter.com/5chan_nel (5ch newer account)

>>449
抗体のはいったサンプルを２，３０個試してみりゃわかんじゃないの？
いくらなんでも感度が３％じゃつかえねぇ。

イカサマキットの感度特異度の事前分布を一様分布に設定して
抗体のはいったサンプルを２０個で全部陰性であったので３０個試したら全部陰性であったとすると
感度・特異度の事後分布は

https://i.imgur.com/6ZYn34k.png

ところが、有病率33/1000であったときに無作為に２０人および３０人を選んで全部陰性であっても
感度・特異度の事後分布は
https://i.imgur.com/jWg4RgH.png

>>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%以下じゃつかいもんにならんわなぁ。インチキで確定。

https://www.buzzfeed.com/jp/naokoiwanaga/covid-19-antibody-test
大阪市立大学が一昨年の血液を使って検査したところ、偽陽性は1人もいなかったそうだ。

>>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
５０人で０でしょ？
だったら、特異度として保証されるのはせいぜい９８％程度じゃん。
１％の陽性率に対して、それでは不十分。

>>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%
ってコイントス変わらんな。

結局のところ、断定的な結論は出せないねということを数字で確認しているだけｗ

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以外は非対称。

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

Wilson法は値が０や１に近くても信頼できるという人がいるけど、どうやって検証するのかはよくわからん。

>>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)
"

どの週においても両群で症状改善率は不変
症状改善率はどちらも一様分布を事前分布に仮定して

症状改善率の差δ=:０というモデルとδ≠０というモデルでのベイズファクター（周辺尤度比＝エビデンス比）を求めると
> (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)

確率分布を図示すると

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

100例もやって差がでないのなら、はっきり言って効き目なしとかわらん。
むしろ副作用がこわい。

アビガンはどうなんだろうね。なぜか国内のまとまったデータがまだ出てこない。

まあ、軽症者が８割とかだと、早期投与の効き目を判定するには相当の
人数に使わないとわかんないだろうけど。

死亡率に関してはベイズファクターで推論できた。

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はレムデシじゃなくてアビガンについての話だかんね。
レムデシは重症者向けらしいけど、たぶん駄目だろ。

【政府が隠す３．１１】　地震波形が、核実験と同じ
http://rio2016.5ch.net/test/read.cgi/earth/1575170033/l50

5/10夜『NHKスペシャル』「新型コロナウイルス　出口戦略は」で
厚労省新型コロナ対策班西浦博北大教授と氏が分析中のノートPC画面が
映されていた。氏の使用解析ソフトウェア名は何というのでしょうか？

西浦先生はRじゃないかな。
今週ぐらいに、実効再生産数を算出するRのコードを公開できるかもと言ってた。

火曜にはニコ生があるから、そこでたっぷりと聞けるはず。
https://sp.live.nicovideo.jp/watch/lv325833316

>>488
>11に数式があるけど

>>489
いや、与えられた何日分かの新規患者数から再生産数とか回復率とかを推定するための数式じゃないの？

>>490
SEIRより単純なReed-Frostモデルだよな
https://en.m.wikipedia.org/wiki/Reed%E2%80%93Frost_model

>>489
自分で計算してみた人ならわかるけど、日々の再生算数を出すにはA「感染してから他人にうつすまでの日数の分布」が必要。西浦先生はこれをまず求めていて、その結果は割合と世界中で使われている。
加えてB「感染してから発表の数字に載るまでの日数の分布」も必要で、発表の数字をBで逆畳み込み積分してC「感染推定日毎の数字」に変換し、CをAで畳み込み積分したものでCを割れば日々のRtが求められる。
とは言うものの、そのままやると誤差で数字が発散するし、逆畳み込みは一意ではないのでどう推定するかは色々と考えなくてはならない。多分日本独自の事情もコードに入れなきゃいけない。大変だろうなあとは思う。

>>492
レスありがとうございます。
確かにSEIRモデルやSIRIモデルだと
「感染してから発表の数字に載るまでの日数の分布」って考慮されてないね。

>>488
公開されたぞー。StanとRだ。

https://nbviewer.jupyter.org/github/contactmodel/COVID19-Japan-Reff/blob/master/scripts/C.%20Calculating%20the%20Rt%20in%20Stan.ipynb

ウイルス感染のモデルも必要だと思うけど
社会経済活動の影響や
自粛解除した後の死者増加や
自粛継続した場合の法人の倒産や
治療薬やワクチンができた後の社会の経済や文化面の回復
そこまで考慮してどういう政策を選択したら
損失を最小化できるかという問題は解ける？

制約付きの最小化問題だと思う
数値化するのが難しい考慮要素もあるからその辺は仮の値とかで

>>495
モデルを精密複雑化すればするほど、必要なパラメータが増えて
それを弱情報事前分布にすると事後分布の信頼区間が広がるんだなぁと思う。

>>494
これにある、
JapaneseDataCOVID19
ってどこかからダウンロードできるんだろうか？

>>497

https://github.com/contactmodel/COVID19-Japan-Reff/tree/master/data

>>498
ありがとうございます。

親ディレクトリ（フォルダ）探せばよかったんだね。

https://nbviewer.jupyter.org/github/contactmodel/COVID19-Japan-Reff/tree/master/

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

再生算数を０〜１０人の一様分布にすると、収束しない。

> 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

どうも、こういう境地だなぁ。

断定的な結論は出せないということを数字で確認しているだけｗ

基本再生算数はウイルスの特徴で決まるもので
実行再生算数は更に環境や人の行動の影響で変化するものだと理解している
実行再生算数を減らすようにするには
どう行動したらいいか
かつ経済活動を出来るだけ下げずに

西浦モデル、一度感染日ごとの人数を最尤法で確定して、そこからベイズを回してRtを求めてるんだよなあ。
本来は、日々の感染者数も確率的に揺れがあるはず。だからRtの誤差幅は発表よりも多く見積もるべきかもしれん。

コードを解説してる人
http://mikuhatsune.hatenadiary.com/entry/2020/05/13/205824

再生産数の事前分布を色々かえて事後分布を出してみた。

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

>>504
切断分布だと収束しないみたいなので、
一様分布[0,10]に近そうな正規分布[5,3]
https://i.imgur.com/h8vMZUM.png
を事前分布にして走らせてみた。

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

【新型コロナ】 東京0.6％、東北6県0.4％陽性・・・抗体検査1000人実施 [影のたけし軍団★]
https://asahi.5ch.net/test/read.cgi/newsplus/1589502801/

複数の検査キットの性能評価と感染状況の確認が目的でしたが、東京都で献血した500人のうち3人、
東北6県で献血した500人のうち2人がいずれかの検査キットで陽性と判定されました。

満員電車など人との接触の多い東京とそうでない東北で陽性率に有意差はあるか？

>>513
検査キットの特異度が99.5%以上あることが検証されました。
ってだけの話じゃね？

まぁ、東北６県と東京が同じってことは、特異度はそれより
高くはないってことだろうね。

言い換えると、真の陽性率はわかりません、ってこと。

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

再生産数を計算するRのプログラムあったんだな

https://www.rdocumentation.org/packages/EpiEstim/versions/2.2-1/topics/estimate_R

>>509
後者の分析、RStan本の訳者が今やっているようだ。

https://twitter.com/hankagosa/status/1261430169283125248
https://twitter.com/5chan_nel (5ch newer account)

>>519

https://twitter.com/hankagosa/status/1261426113374416897
https://twitter.com/5chan_nel (5ch newer account)

>>520
アヒル本の著者だね。俺も読んだ。

4/2の時点で感染者数を６０００くらいと見積もってるね。＞アヒル本の人
実数の三倍程度。いい線かもしれない。

https://github.com/contactmodel/COVID19-Japan-Reff/blob/master/nishiura_Rt%E4%BC%9A%E8%AD%B0_12May2020.pdf

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
でちゃんと動作した

>>526
thx

R に　EpiEstimというパッケージがあって、再生産数を算出する関数が搭載されている。
結局、infecterとinfecteeが発症するまでの期間serial intervalの分布をどう設定するかで結果が変わるみたいだなぁ。

Rのヘルプファイルを解読中。
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（例、ｔ検定のソース表示は　stats:::t.test.default　）で表示される。

EpiEstim::estimate_R
EpiEstim:::process_si_data
EpiEstim:::process_config_si_from_data
coarseDataTools::dic.fit.mcmc