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

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

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

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

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

まぁそれなら
感度+特異度>1
のときは検査数を上げていいなら望むだけ市中感染率を正しく推定できる
ですな。
なので上先生の勝ち。
東大生の負け。

r=x/nとして
Θ = (q + r - 1)/(p + q - 1)

>>85
プログラムを組んで計算させてみた（有病率の事前分布は一様分布を仮定)

感度0.7 特異度0.9として
100人に１人陽性だった場合の有病率の期待値、中央値、最頻値、９５％信頼区間
> sim(100,1)
mean median mode lower upper
0.0196 0.0166 0.0101 0.0004 0.0464

1000人に１0人陽性だった場合の有病率の期待値、中央値、最頻値、９５％信頼区間
> sim(1000,10)
mean median mode lower upper
0.0110 0.0107 0.0099 0.0050 0.0175

10000人に１00人陽性だった場合の有病率の期待値、中央値、最頻値、９５％信頼区間
> sim(10000,100)
mean median mode lower upper
0.0101 0.0101 0.0100 0.0082 0.0121

100000人に１000人陽性だった場合の有病率の期待値、中央値、最頻値、９５％信頼区間
> sim(100000,1000)
mean median mode lower upper
0.0100 0.0100 0.0100 0.0094 0.0106

>>85
感度+特異度>1　の条件は必要？
感度0.4　特異度0.5でシミュレーションしてみたけど
ｎを増やせば信頼区間が狭くなっていく


> sim(100,1,sen=0.4,spc=0.5)
mean median mode lower upper
0.0196 0.0166 0.0102 0.0004 0.0463
> sim(1000,10,sen=0.4,spc=0.5)
mean median mode lower upper
0.0110 0.0107 0.0099 0.0050 0.0176
> sim(10000,100,sen=0.4,spc=0.5)
mean median mode lower upper
0.0101 0.0101 0.0101 0.0082 0.0121
> sim(100000,1000,sen=0.4,spc=0.5)
mean median mode lower upper
0.0100 0.0100 0.0100 0.0094 0.0106

>>89
ああ、p+q<1だと逆にあてにならなすぎて逆張りしてθが推定できるんだな。

プログラムにバグがあったので修正というより、stanでのMCMCに変更。

> PCR(100,1)
mean lower upper
0.01497496587 0.00004299248 0.04480331292
> PCR(1000,10)
mean lower upper
0.001642898836 0.000001299175 0.004982795701
> PCR(10000,100)
mean lower upper
0.000174692810 0.000000052897 0.000560041026
> PCR(10000,1000)
mean lower upper
0.003876080828 0.000002574352 0.010176384444

まあ、ｎを増やすほど信頼区間が狭くなるという結論には変わりない。

【富山県最強伝説】新型コロナウイルスPCR検査件数　54人　陽性0人
https://asahi.5ch.net/test/read.cgi/newsplus/1584811696/

ある集団から54人を無作為に選んでPCR検査したら陽性0であった。PCR検査の感度0.7　特異度0.9としてこの集団の有病率の期待値と95%信頼区間を求めよ。

>>92
> PCR(54,0)
mean lower upper
0.0282436750 0.0000053869 0.0839653803

日本人の平均身長を推測するのにその値は１〜２ｍの間であるという弱情報事前分布は合理的。
現時点での新型コロナの有病率は0.1未満の一様分布という弱情報事前分布として
【富山県最強伝説】新型コロナウイルスPCR検査件数　54人　陽性0人
ある集団から54人を無作為に選んでPCR検査したら陽性0であった。感度0.7　特異度0.9としてこの集団の有病率の期待値と9信頼区間を推測する。
事前分布のパラメータを変えるとstanだとコンパイルが必要になるのでjagsでプログラムを組んでみた。



# 感度SEN, 特異度SPCの検査でN人中X人が陽性であったときの推定有病率prevalence
# 弱情報事前分布はprevalence ~ dunif(0,UL) , UL:上限
library(rjags)
PCRj <- function(N,X,UL=1,SEN=0.7,SPC=0.9,verbose=TRUE){ # UL:upper limit of dunif(0,UL)
modelstring=paste0('
model
{
x ~ dbin(p,n)
p <- prev*sen + (1-prev)*(1-spc)
prev ~ dunif(0,',UL,')
}
')
if(verbose & UL!=1) cat(modelstring)
writeLines(modelstring,'TEMPmodel.txt')
dataList=list(n=N,x=X,sen=SEN,spc=SPC)
jagsModel = jags.model( file="TEMPmodel.txt" ,data=dataList,quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples( jagsModel ,
variable=c("prev","p"), n.iter=1e6, thin=10)
js=as.matrix(codaSamples)
BEST::plotPost(js[,'prev'],xlab='prevalence',showMode = TRUE) ; lines(density(js[,'prev']),col='skyblue')
round(c(mean=mean(js[,'prev']),HDInterval::hdi(js[,'prev'])[1:2]),10)
}

実行結果
> PCRj(54,0,UL=0.1)

model
{
x ~ dbin(p,n)
p <- prev*sen + (1-prev)*(1-spc)
prev ~ dunif(0,0.1)
}
|**************************************************| 100%
mean lower upper
0.0245104429 0.0000003782 0.0703606657

統計のことははぜんぜんわからんが、確率論的には
検査陽性率の期待値＝有病率×感度＋(1-有病率）×（1ー特異度）
なんだから、
有病率<<1なら、検査陽性率の期待値≒有病率×感度＋(1- 特異度)

っちゅうことで、感度７０％で特異度が100％なら検査陽性率の３割増し
が有病率だとみなせばよろしいんでしょ。
一方、特異度が90%しかなかったりすると陽性率の期待値が10%も
水増しされちゃうから、有病率の推定が大幅に困難になる。

つまり、特異度がよほど高くなければ(99%とかね）、有病率が数％以下
の状況でランダム検査しても偽陽性が真陽性を上回って混乱をきたす。
（見かけ上致死率は下がる、ってのが南朝鮮の状況か）

>>95
同感。
有病率0.1%で特異度0.9なら偽陽性だらけになるんだよね。

>>95-96
疾病率＝陽性率、すなわち無作為抽出ならその算段も成立するかもね。
現行の制度では推定される市中感染率が1/10000程度で陽性率が5%ほどらしいから現行制度下での検査はうまく行ってますね。
ただ感度＋特異度が1.7位あるので検査数増やした方が有益である事は間違いがない。

>>95
誤解や、誤解を引き起こしかねない内容があったので、勝手に補足させていただきます。

>>っちゅうことで、感度７０％で特異度が100％なら検査陽性率の３割増し
>>が有病率だとみなせばよろしいんでしょ。
この場合、検査陽性率の期待値＝有病率×感度　なのだから、1/(0.7)=1.4285...
3割増しではなく、4割強増しと言うべき。

>>一方、特異度が90%しかなかったりすると陽性率の期待値が10%も
>>水増しされちゃうから、有病率の推定が大幅に困難になる。

感度70、特異度100で、有病率　1%、0.1%、0.01%、0.001%の時、10万人に検査したときの
検査陽性者数は700,70,7,0.7人です。　有病率に比例して増減し、これらは、はっきり見極めできます。
一方、感度70、特異度90で、同じ事をすると、それぞれ、10600,10060,10006,10000.6人です。
有病率に応じて差はあるのですが、常に偽陽性が約10000人いて、誤差を考えると、見極めは困難です。

「陽性率の期待値が10%も水増しされちゃうから」と書かれていますが、これは、
検査対象者10万人に対する約10%=1万人が常に水増しされているのであって、
陽性と判定される人の数が10%水増しされていると誤解しないよう、補足しておきます。

>>98
心配性だねｗ
3割が4割でもたいして変わらんがな。
陽性率の期待値が10%水増しってのも、期待値の10%じゃなくて、
期待値そのものが10%上昇するんだってことくらい、元の式見りゃ
自明だしね。

そもそも感度や特異度の具体的な数値はよくわかんないんだから、
具体的な数字にこだわってもしょうがない。
有病率が低いときの、有病率と陽性率と特異度の関係がわかれば
よろし。

こういう計算が必要になる。

事前分布を選択する（例. 有病率は高々10%として(0.0.1]の一様分布とする)、
陽性確率は真陽性確率と偽陽性確率の和、
陽性数はこの確率で二項分布、
以上を実際に得られた検査数と陽性数から最尤値となるパラメータとして有病率の分布を出して期待値や信頼区間を出す。
手計算では無理。

感度0.7　特異度0.9でコンピュータに計算させると

陽性率が30％なら有病率の推測値は

> PCRj(100,30)
mean lower upper
0.3396123 0.1994400 0.4964517

> PCRj(1000,300)
mean lower upper
0.3343576 0.2876152 0.3832246
> PCRj(10000,3000)

mean lower upper
0.3333387 0.3186388 0.3481056
> PCRj(100000,30000)
mean lower upper
0.3332425 0.3286774 0.3380952

と　期待値も信頼区間もそれらしい値になるが、


陽性率が1%なら有病率の推測値は偽陽性が増えまくるので陽性率と有病率が乖離する、
> PCRs(100,1)
mean lower upper
0.01497496587 0.00004299248 0.04480331292
> PCRs(1000,10)
mean lower upper
0.001642898836 0.000001299175 0.004982795701
> PCRs(10000,100)
mean lower upper
0.000174692810 0.000000052897 0.000560041026
> PCRs(100000,1000)
mean lower upper
0.000016780530193 0.000000002738145 0.000051489256688


陽性率が60%なら、今度は偽陰性が増えまくるので偽陰性が増えまくるので陽性率と有病率が乖離する。
> PCRs(100,1)
> PCRs(100,60)
mean lower upper
0.8280581 0.6766480 0.9764282
> PCRs(1000,600)
mean lower upper
0.8335023 0.7826355 0.8825766
> PCRs(10000,6000)
mean lower upper
0.8334634 0.8187334 0.8492064
> PCRs(100000,60000)
mean lower upper
0.8332873 0.8289242 0.8379535

Rとstanでベイズ統計ができるなら、以下のコードで実行可能。
パッケージ rstan　と　BEST(信頼区間グラフ描出用)が必要


library("BEST")
library("rstan")
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
options(scipen = 5)

model.string='
data{
int n; // sample size
int x; // number of positive test
real<lower=0,upper=1> sen; // sensitivity 0.7
real<lower=0,upper=1> spc; // specificity 0.9
real<lower=0,upper=1> ul; // uniform(0,ul)
}

parameters{
real<lower=0,upper=1> prev; // prevalence
}

transformed parameters{
real<lower=0,upper=1> p;
p = prev*sen + (1-prev)*(1-spc) ; // probability of positive test result
}

model{
x ~ binomial(n,p);
prev ~ uniform(0,ul);
}

'
writeLines( model.string , con='model.stan' )
corona1.model=stan_model('model.stan')
# saveRDS(corona1.model,file='corona1.rds')
# corona1.model=readRDS('corona1.rds')

PCRs <- function(N=1000,X=10,UL=1,SEN=0.7,SPC=0.9,verbose=FALSE,...){
data = list(n=N,x=X,sen=SEN,spc=SPC,ul=UL)
fit.corona = sampling(corona1.model, data=data,
seed=1234,control=list(adapt_delta=0.99),...)
if(verbose) print(fit.corona, prob=c(0.025,0.5,0.975),pars=c('prev'),digits=8)
ms=rstan::extract(fit.corona)
BEST::plotPost(ms$prev,showMode = T,xlab='prevalence') ; lines(density(ms$prev),col='skyblue')
c(mean=mean(ms$prev),HDInterval::hdi(ms$prev)[1:2])
}

>>102(補足)

非対称分布の信頼区間計算にパッケージHDIntervalも必要。

>>76
偽陽性率が現時点での有病率を大きく上回るから東大生の言い分が正しい。
有病率が３割程度になれば上先生の言い分が正しい。

>>104
1000人調べたときの検査陽性率と推定陽性率をグラフにしてみた。
灰色直線は検査陽性率＝推定陽性率の直線
検査陽性率が低いときは過小評価、高いときは過大評価する。
https://i.imgur.com/wmOAj5i.png

何言ってんの？
市中感染率をいくらでも正しく推定できるかなんでしよ？
結論のために問題変えるなよ。
政治の話をここに持ち込むなよ。

>>104
上の言い分も一理あるし、東大生の言い分も一理あるｗ

特異度が９０％を越える高い値だという前提があればこうなる。

１）検査陽性率が低く（数％以下）、特異度がほぼ100%でないの
　なら、市中感染率を推定するのは難しい。とはいえ、上限は
　（検査陽性率/感度）程度で抑えられる。つまり10%以下くらいの
　ことは言えるが、0%かもしれない。

２）一方、検査陽性率が高ければ（数十％以上）、下限は検査
　陽性率程度と見込めるが、市中感染率は感度に依存して大きく
　変化する。つまり、数十％以上とは言えるが100%近いかどうか
　までは不明。

ってことで、検査陽性率からある程度市中感染率の目安は立つが、
それがどこまで意味があるとみなせるかは疑問。TPO次第か。

あと、補足すると、感度と特異度が正確にわかっているのなら、
統計学的に市中感染率を推定することはある程度可能だけど、
実際はそうではないから、上様の言い分には意味がない。

>>105
感度と特異度を変化させて、検査陽性率と推定有病率の関係をグラフにしてみた。
https://i.imgur.com/YEtcSfn.png

だから自分の言いたい結論が先に決まっててそれに合わせて好き勝手に問題読み替えてるだけじゃん？
そんな考え方しかできないなら理系板に来んなよ。

>>108
感度も特異度も定数でなく何らかの分布に従うパラメータとしてモデルを組めばいいだけの話だろ。

感度が最頻値0.7 標準偏差0.05のベータ分布β(58.229, 25.527 )
特異度が最頻値0.9 標準偏差0.05のベータ分布β(36.172, 4.908)
に従うと過程して

model
{
x ~ dbin(p,n)
p <- prev*sen + (1-prev)*(1-spc)
sen ~ dbeta(sn[1],sn[2])
spc ~ dbeta(sp[1],sp[2])
prev ~ dunif(0,ul)
}
こういうモデルでMCMCすれば可能。

実行結果

> PCRj2(1000,10) # 陽性率1％で有病率を推定
mean lower upper
0.001667604624 0.000000053909 0.004956969423

> PCRj2(1000,300) # 陽性率30％で有病率を推定
0.33414 0.28797 0.38253

> PCRj2(1000,600) # 陽性率60％で有病率を推定
mean lower upper
0.83296 0.78428 0.88496

β分布のパラメータを出すRスクリプト

# a,b to Mode,mean,variance
ab2Mmv<-function(a,b){
M<-(a-1)/(a+b-2)
m<-a/(a+b)
v<-a*b/((a+b)^2*(a+b+1))
cat('Mode =',M,'mean =',m,'variance =',v,'\n')
invisible(c(Mode=M,mean=m,variance=v))
}

# Mode,kappa to mean,variance
Mk2mvab= function( mode , kappa ) {
# if ( mode < 0 | mode > 1) stop("must have 0 <= mode <= 1")
# if ( kappa <2 ) stop("kappa must be >= 2 for mode parameterization")
a = mode * ( kappa - 2 ) + 1
b = ( 1.0 - mode ) * ( kappa - 2 ) + 1
m=a/(a+b)
v=m*(1-m)/(a+b+1)
return( c( mean=m , variance=v,a=a,b=b ) )
}

# Mode,variance to a,b
Mv2ab = function(mode,vari){
f=function(kappa) Mk2mvab(mode,kappa)[2] - vari
kappa=uniroot(f,c(2,10000))$root
ab=Mk2mvab(mode,kappa)[c('a','b')]
ab2Mmv(ab[1],ab[2])
return(ab)
}
(sn=Mv2ab(0.7,0.05^2))
curve(dbeta(x,sn[1],sn[2]),bty='l')
(sp=Mv2ab(0.9,0.05^2))
curve(dbeta(x,sp[1],sp[2]),bty='l')

上記の準備をして以下で実行

PCRj2 <- function(
N,X,
UL=1,
SEN=0.7,
SPC=0.9,
SD=0.05,
print=TRUE){
# UL:upper limit of dunif(0,UL)
library(rjags)
library(BEST)
sn=Mv2ab(SEN,SD^2)
sp=Mv2ab(SPC,SD^2)

modelstring=paste0('
model
{
x ~ dbin(p,n)
p <- prev*sen + (1-prev)*(1-spc)
sen ~ dbeta(sn[1],sn[2])
spc ~ dbeta(sp[1],sp[2])
prev ~ dunif(0,ul)
}
')
writeLines(modelstring,'TEMPmodelj.txt')
dataList=list(n=N,x=X,ul=UL,sen=SEN,spc=SPC,sn=sn,sp=sp)
jagsModel = jags.model( file="TEMPmodelj.txt" ,data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples( jagsModel ,
variable=c("prev","p","sen","spc"), n.iter=1e5, thin=5)
js=as.matrix(codaSamples)
if(print){
BEST::plotPost(js[,'prev'],xlab='prevalence',showMode = TRUE)
lines(density(js[,'prev']),col='skyblue')}
re=c(mean=mean(js[,'prev']),HDInterval::hdi(js[,'prev'])[1:2])
return(re)
}

options(digits = 5)
options(scipen = 5)

PCRj2(1000,10) # 陽性率1％で有病率を推定
PCRj2(1000,300) # 陽性率30％で有病率を推定
PCRj2(1000,600) # 陽性率60％で有病率を推定

>>99
期待値というのは、無次元量ではない。観測値とか物理量と同じように単位をつけて議論できる量。
従って「期待値が10%増える」等という言葉があれば、期待値が1.1倍になるのだろうと感じるのが普通。
そのような性質を持つ期待値に対し、「10%増える」と表現し、
「期待値の値そのものが、0.1増えることを意味している」
と説明しなければならないならば、やはり誤解を招きやすい表現だと思う。

今回の期待値は比率であり、無次元量であったから、「10%」と言うのが、
どちらの意味としても、通用したため発生したとは言えるが、読み手の立場に立った表現を望む。

似た議論に、選挙時の投票率がある。前回の投票率が40%。今回の投票率が50%だとする。
「前回に比べ、今回は10%増えました」
「前回に比べ、今回は25%増えました」
どちらも、言い得る表現。聞き手の混乱を避けるため、前者の意味で使う場合、
「10%ポイント増えました」とコメントするのを最近聞くようになった。
私にはよい傾向と感じるが、中には、違いは何かとか、混乱の源の存在さえ意識していない人もいるようだ。

「３割増も４割強増も大した差ではない」には、「式が違っても結果が誤差範囲なら問題ない」
という考えが背景に見える。そのような方が、混乱を引き起こしかねない表現を用いた。
だから、補足した。果たして本当に杞憂だったのだろうか？

富山では62人PCR検査して陽性0人（3月22日までの集計)有病率を推定とその信頼区間を推定したい。

http://www.pref.toyama.jp/cms_pfile/00021629/01366377.pdf

PCR検査の感度は最頻値0.6標準偏差0.１、特異度は最頻値0.9標準偏差0.05のベータ分布（正規分布は負になったり1を超えるので不適）、
有病率は一様分布として、推定される有病率の期待値と95％を計算せよ。

図示するとこんな感じ。
https://i.imgur.com/Ip6gSCa.png

stanのモデルのスクリプトはこれ
sn,spはβ分布のパラメータ、その計算法は既述

data{
int n; // sample size
int x; // positive test result
real<lower=0,upper=1> ul; // uniform(0,ul)
real<lower=0> sn[2]; // sen ~ beta(sn[1],sn[2])
real<lower=0> sp[2]; // spc ~ beta(sp[1],sp[2])
}

parameters{
real<lower=0,upper=1> prev; // prevalence
real<lower=0,upper=1> sen; // sensitivity
real<lower=0,upper=1> spc; // specificity
}

transformed parameters{
real<lower=0,upper=1> p;
p = prev*sen + (1-prev)*(1-spc) ; // probability of positive test result
}

model{
x ~ binomial(n,p);
prev ~ uniform(0,ul);
sen ~ beta(sn[1],sn[2]);
spc ~ beta(sp[1],sp[2]);
}

>>116
ここで問題
感度特異度の分布はそのベータ分布として
何人陰性が続けば95%信頼区間の上限が0.05を下回るか？

>>117
感度、特異度の分布？？？

>>118
何でも確率変数にするのがベイズ推計。
ｐ値の分布すら考えるぞ。

PCR検査の感度は最頻値0.6標準偏差0.１、特異度は最頻値0.9標準偏差0.05のベータ分布を事前分布にしたけど、
事後分布はstanによるMCMCで
感度は期待値0.57　95％信頼区間は[0.37,0.77]
特異度は期待値0.96 95%信頼区間は[0.91,0.99]
とコンピュータが計算してくれる。

>116のように弱情報事前分布を設定することで事後分布は次のように描ける。

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

>>54
いや、特異度の事前分布を設定することで事後分布をMCMCで求めることができる。
>116の設定での結果が>121

>>122
結局事前分布の設定次第ってことはないの？

>>123
日本人の平均身長を推測するのにその値は１〜２ｍの間であるという弱情報事前分布は合理的。
感度特異度の分布に正規分布を使うのはアホ。
負になったり、１を超えたりするから。

>>123
感度を0.4-0.8の一様分布、特異度を0.8-1.0の一様分布にしても有病率の推定値は
> round(re$mci,5)
mean lower upper
0.02827 0.00000 0.08592
であまり変わらないね。

sensitivity ~ N(m=0.6,sd=0.1) specificity ~ N(m=0.9, sd=0.05)
にしても推測有病率は平均3％弱で　95%CIは0-8％とあまり分布の形にはよらないね。
mean lower upper
0.026841384 0.000000153 0.081071379

確率だと定義域が0-1で計算しやすいのでβ分布を使うことが多い。

>>119
そう？
統計の推定の理論で推計する母数は確率変数ではないと習ったけど？

>>127
時代は頻度主義統計からベイズ統計だよ。

頻度主義統計でも最尤推定では
データを固定してパラメータを動かすだろ。

>>127
階層ベイズモデルを扱ったことないの？
>112は簡単な実例。

>>128 >>130
何を持ってベイズ統計っていってんのか知らん。
pcr検査の感度とは被験者が感染者である場合の検査結果が陽性となる条件付き確率でしょ？
条件付き確率の分布ってどういうことよ？
確率がまた確率になるってなんの話してんの？
変数Xの平均とか分散とかは統計学においては推定すべき定数であって確定値。
それの分散なんて数学的に意味不明。
一体どこの統計学の教科書にそんなデタラメ書いてあんの？

まちがえた。
確率の平均がまた確率変数になるってどういうことよ、ね。
式でかけば確率変数Xの平均E(X)の分散ってなんの話ってことになる。
確率変数Xはある標本空間上の関数だけどE(X)は実数だよ？

>>131
ベイズ階層モデルも組めない奴とは議論にならんね。
分散の事前分布に逆ガンマ分布でなく半コーシーを使う方がいいとかいう議論も理解できんだろ。

>>131
ベータ分布は定義域が[0,1]で二項分布の確率の確率密度関数としてベイズ階層モデルでは頻用されるよ。

ベイズ階層モデルを使わずにこの計算できるならやってみてくれ。
020/3/24 11:00時点で検査人数での陽性率は171/2013であるという。
新型コロナ肺炎のPCR検査の感度は５〜７割、特異度は９割前後らしい。幅をもたせた値を使って検査をうけたグループの有病率を計算せよ。

>>133
まぁ言わんとする事はもちろんわかるし伝わるけど、疫学だから数学やってる人間がなんとなく伝わるではダメだろ？
数学だけの話ではなく、疫学は実社会とキチンと繋がってるんだから？
統計学ではあくまで検定する母数は定数。
それは確率モデルでは実数値であり、定数。
そして統計データを確率変数に割り当てる。
当然それらの確率変数は一つの測度空間の一つしかない確率変数であり、平均も分散もひとつしかない定数値。
それらをいっぱい考えてどうこう言ってるんだろうとは思うけどそんなの統計学や疫学の一般的な考えにはない。
何故なら現実世界はひとつしかなく、確率変数に対応している統計量も一個しかない。
もちろん母数がめちゃめちゃ大きい統計量で例えば10000個のデータを100こずつ切って100個の統計量を100の世界からとってきたなんて考えが無理クリできなくはないが、そんな考え方は普通しない。
それはあくまで100個ずつに区切られた10000個の一つの世界の確率変数としか扱わない。
そういうオリジナルな考えで捉えたいならそれは勝手だけど、それならそれで話の中で明示しないとダメ。
数学の世界なら言わずもがなの話は言わなくてもエスパーしてもらえても、疫学、統計学の世界では実社会とつながる話だからダメ。

>>135
能書きいいから、

ベイズ階層モデルを使わずにこの計算できるならやってみてくれ。

020/3/24 11:00時点で検査人数での陽性率は171/2013であるという。
新型コロナ肺炎のPCR検査の感度は５〜７割、特異度は９割前後らしい。幅をもたせた値を使って検査をうけたグループの有病率を計算せよ。

>>136
こういう判断が現実には必要。
検査特性を無視して単純な割り算だと検査を受けた人の有病率は8.5%弱になるけどこれは過大評価か過小評価か？

検査感度が5-7割、特異度が9割前後なら
検査陽性率=有病率とすると常に過大評価かどうか気になったので陽性数を変化させて計算してみた。
検査感度はmode=0.6,sd=0.1　特異度はmode=0.9,sd=0.05のベータ分布に設定してJAGSでベイズ階層モデルをたてて計算。


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

陽性率が20%未満のときは過大評価、それ以上のときは過小評価である、という結論になった。

ベイズ統計を理解できている人の検証希望。

>>138
プログラムの練習がてらに、
MCMCのアルゴリズムの異なるstanでベイズ階層モデルを組んで検証。
当然ながら、同様の結果。　検査陽性率が２０％を境に過大評価と過小評価が入れ替わる。

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

>>136(自己レス)

今日の都の発表で(171+41)/(2013+89) に検査陽性率が増えたので再計算。

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

>>138
サンプリング回数を増やしてグラフを完成。

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

>>140
89は検査数で検査人数は74という。
計算し直すと

> PCRj2(N,r,SEN=0.6,SD1=0.1,SPC=0.9,SD2=0.05,N.ITER=5e5)
|**************************************************| 100%
mean lower upper
0.05720165 0.00000015 0.1332385

41/74の推測有病率は

mean lower upper
0.8121975 0.5957315 0.9999992

オリンピック延期発表後の検査陽性率は88/169で52%だが、
PCR検査の感度と特異度がはっきりしないので、検査陽性率をこの集団の有病率とするのは正しくない。
88/169のときの感度・特異度と推定有病率の関係をグラフにしてみた。
https://i.imgur.com/iQC88tZ.png
感度0.6、特異度0.9のときの推定有病率は85%で陽性率からの憶測は過小評価といえる。

検査陽性率が小さいときには、実際の有病率より過大評価してるし、
検査陽性率が高いときは、過小評価してるだろうってことでしょ。
そのくらいは定性的に理解できる。

>>145
どこが境目かは直感じゃわからんね。

そりゃ感度や特異度次第だからな。

まあ、数％と数十％では違うんだということがわかればいいんじゃね？
境目なんかどうでもいいでしょ。

陽性率が15%でこれを有病率の推測値に使うのは過大評価なのか過小評価がわからんのはまずいね。

オリンピック延期決定以後の検査数と陽性数
subjects=c(74,95,87)
positives=c(41,47,40)
PCRs3(subjects,positives,iter=10000,warmup=1000)
として、
感度・特異度を考慮した推定有病率は
mean lower upper
0.77417 0.56756 0.99944
>
日々の陽性数が二項分布に従うとして計算。

志村けんみたいな有名人がコロナに感染してることから日本全体のコロナ感染者数を推定してみる。

まず、日本の有名人が1000人いるとしよう。
つぎに、日本でコロナに感染していない確率をxとしよう。
すると、有名人1000人が一人も感染していない確率は、xの1000乗となる。これをyとおこう。すると、有名人が一人でも感染している確率は(1-y)となる、これをzとおこう。
まとめると以下の関係がなりたつ。

・コロナに感染しない確率:x
・有名人が一人もコロナに感染しない確率:y=x^1000
・有名人が一人でもコロナに感染している確率:z = 1-y

志村は感染したわけなので、以下、2つのケースにわける

ケース1: zが10%のとき
z=0.1, 故にy = 1-0.1=0.9
故にx = y^0.001よりx=0.9^0.001=0.999894
これがコロナに感染していない確率なので、
コロナ感染確率は、1-0.999894=0.000106
よって日本のコロナ感染者数は推定
120,000,000*0.000106=12,720人

ケース2:zが50%のとき
ケース1と同様の計算で、
日本のコロナ感染者数は推定
120,000,000*0.000693=83,160人

検査が少ないから感染者増が緩やか？数学的に検証してみた
http://agora-web.jp/archives/2045047.html

主な関係国について、新型コロナ感染者数の片対数グラフがある。
http://agora-web.jp/cms/wp-content/uploads/2020/03/WS000876.jpg
FT.comより

感染者数の伸びが日本は緩やかと解釈するのが普通だが、検査が少ないからとする解釈もある。本当はどうなのか計算してみる。

結論を先に書くと、検査が多いか少ないかは関係ない。

対数をとると係数（感染者の発見率）は定数項になり、今回の片対数グラフの整理法の前提としてキャンセルされる
日本が展開しているのは患者認定の精度上昇であり、医療リソースの効果を最大化して死者数を低く抑えている要因の一つといえる

キャバクラ客100人から無作為に5人から検体を採取してこの検体を混合攪拌してコロナ検査したところ陽性であった。

（１）100人のキャバクラ客の陽性数の期待値と95％信頼区間を求めよ。
（２）PCR検査の感度0.6、特異度0.9として100人のキャバクラ客の感染数の期待値と95%信頼区間を求めよ。

>>151
> m=1000 # 有名人の人数
> n=1.268e5 # 日本の人口
> x=0:n # 感染者数:x, 非感染数:n-x
> pmf=1- chooseZ(n-x,m)/chooseZ(n,m) # 1 - (ｍ人全員非感染の確率)
> pdf=pmf/sum(pmf) # 確率密度関数化して
> (E=sum(x*pdf)) # 期待値を計算
Big Rational ('bigq') :
[1] 63590201/1002
> as.numeric(E)
[1] 63463.27

6万3000人と計算された。

>>155
良く分からんが、ありがとう。
こちとら高校レベルの確率の知識しかないもんで。

>>156
n(=10)人の中にi人の感染者がいるとき無作為にm(=2)人を選ぶ。
選ばれた2人の中に少なくとも一人の感染者がいる確率をP[x]として、
n個からr個選ぶ組み合わせの数をChoose(n,r)で表すと

P[xi]=1- choose(10-x,2)/choose(10,2)

xを0から10まで変化させて、

Σx*P[x]/(ΣP[x])で

期待値が求まる。

タイプミス修正

P[x]=1- choose(10-x,2)/choose(10,2)

何このアホスレ？

有名人の数を増やしてみても同様の結果になった。

> # 有名人が感染
> library(gmp)
> m=18200 # 有名人の数(桜を見る会参加人数)
> n=1.268e5　# 日本の人口
> x=0:n # 感染者数:x, 非感染数:n-x
> pmf=1- chooseZ(n-x,m)/chooseZ(n,m) # 1 - (ｍ人全員非感染の確率)
> pdf=pmf/sum(pmf) # 確率密度関数化して
> (E=sum(x*pdf)) # 期待値を計算
Big Rational ('bigq') :
[1] 1154070201/18202
> as.numeric(E) # E=63463.27 (m=1000) , E=1154070201/18202=63403.48(m=1.268e5)
[1] 63403.48

>>160
なんだってー。直感に反するな

>>161
総人口100人として有名人の数を1〜100人まで変化させて、有名人に感染者がいたときの100人中の感染者の数をグラフにすると

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

有名人の数を変化さえても期待値にさほどの変化はない。

>>162
感染者が1名以上という条件だと、
有名人の割合が一定以上になると飽和するんだな。

Ax: x人の感染者がいる(x=0〜n)という事象
B:最低一人の感染陽性判定という事象
Pr[Ax|B]=Pr[B|Ax]Pr[Ax]/Pr[B]
Pr[Ax]:事前確率
Pr[B|Ax]:尤度
Pr[B]:周辺尤度(規格化定数)

求めたい期待値Eは
Σ(x*Pr[Ax|B])/ΣPr[Ax|B] = Σ(x*Pr[B|Ax]Pr[Ax])/Σ(Pr[B|Ax]Pr[Ax])
Pr[Ax]がxにかかわらず定数であれば
E=Σ(x*Pr[B|Ax])/Σ(Pr[B|Ax])

事前確率分布を一様分布と仮定しての計算
つまり、感染者が1人の確率も50人の確率も100人の確率,....も一定という前提での計算。

>>163
そうみたいですね。

> data.frame(有名人=1:10,期待値=sapply(1:10,function(x) fn(100,x)$mean))
有名人 期待値
1 1 67.00000
2 2 62.75000
3 3 60.20000
4 4 58.50000
5 5 57.28571
6 6 56.37500
7 7 55.66667
8 8 55.10000
9 9 54.63636
10 10 54.25000
> data.frame(有名人=1:10*10,期待値=sapply(1:10*10,function(x) fn(100,x)$mean))
有名人 期待値
1 10 54.25000
2 20 52.31818
3 30 51.59375
4 40 51.21429
5 50 50.98077
6 60 50.82258
7 70 50.70833
8 80 50.62195
9 90 50.55435
10 100 50.50000

「有名人」を「wikに載ってる人」と定義し
その数を１００００人としてそのうち４人（志村、藤浪、長坂、伊藤隼人）
感染したとしても結果は変わらない

昨日の東京のコロナ陽性者は87人検査して63人陽性であったという。
検査の感度0.6　特異度0.9と仮定して、87人中に感染者は何人と推定されるか？

真陽性率＝感度=0.6
偽陽性率＝１−特異度=0.1

87人中の感染者数をxとすると

陽性者数= 感染者数*真陽性率 + 非感染者数*偽陽性率

63=x*0.6+(87-x)*0.1

これを解くとあり得ない答になる。

>>166
総人口ｎ人、有名人ｍ人、そのうち感染者ｋ人とすると
ｎ人中の感染者の期待値は
x =　0 〜 nとして 、xCkはｘ人からｋ人選ぶ組み合わせの数を表す

Σ(x*(xCk/nCm))/Σ(xCk/nCm)　= =Σ(x*(xCk))/Σ(xCk)
となるのでｍの値には依存しない。

ｎ

　

>>168
するとこの計算で出てくる推定感染者数6万人って値は意味ない感じですか？

>>167

陽性者数が87人中63人になるような感度と特異度を最小二乗法で求めると。

> (opt=optim(c(0.6,0.9,63),nazo,method='CG'))
$par
[1] 0.916014625 0.779617519 63.002729987

学術の巨大掲示板群 - アルファ・ラボ
ttp://x0000.net

数学 物理学 化学 生物学 天文学 地理地学
IT 電子 工学 言語学 国語 方言 など

>>170
初期値に依存するから意味のないスクリプトであると判明したので撤回します。

>>169
単なる数字の遊びだろうね。

>>169
前提となっているのが、
日本人1億2680万人いるとして
日本人の感染者数が１人である確率も1億人である確率も同じと、一様分布を仮定しているのが現実離れしている。
よって現実的には意味がない。

https://youtu.be/WUMN_71p3Js?t=56

専門家会議がモデルを出したから議論してくれ

>>175
Reed Frost モデルかな？
何を使ったかには言及がなかった。

Reed -Frostはパラメータが１個ですむから推定しやすいんだろう。

>>76
54119人という値になった。
計算プログラムは以下の通り。

# width of 99% confidence interval when 1000 subjects are examined
p2w <- function(
prevalence,
subjects=1000,
sensitivity=0.6,
specificity=0.9,
conf.level=0.99){
# prevalence -> width of 99% confidence interval
n=subjects
p=prevalence*sensitivity+(1-prevalence)*(1-specificity) # positive rate=prev*TP+(1-prev)*FP
q=1-p
2*qnorm(1-(1-conf.level))*sqrt(p*q/n) # width of 99%CI

}

p2w=Vectorize(p2w)
prevalence=seq(0,1,by=0.01)
plot(prevalence,p2w(prevalence),bty='l',type='l',lwd=2,ylab='99%CI width',
main='subjects:1000\nsensitivity:0.6\nspecificity:0.9')
optimize(p2w,c(0,1),maximum=TRUE)
#
sj2w <- function(subjects){ # subjects -> maximum 99%CI width & its prevalence
optimize(function(prev) p2w(prev,subjects),c(0,1),maximum = TRUE)
}
# at how many subjects 99%ci width equals 0.01
uniroot(function(x,u0=0.01) sj2w(x)$objective-u0,c(1000,100000))

上って、灘校東大理IIIの超秀才のはずなのに、なんで
あんなに頭の悪い発言ばかりしてんの？

変な宗教にでも取り憑かれて理性が狂わされてるのかな？

日本人1億2680万人からX人を無作為に抽出してPCR検査して、感染者数(≠検査陽性者数)を信頼区間99%誤差±1%で検定したい。
PCR検査は感度0.6，特異度0.9とする。

何人を抽出すれば十分といえるか？

54000人程度になったけど、あってる？

>>179
超秀才は理Iに行くんじゃないの？

>>180
感染率の程度、感度・特異度の値の精度の言及無しに出された結論に、ほとんど説得力は無い。

>>182
感度 beta(13.6991,9.4661)でmode 0.6 sd=0.1
特異 beta(36.172,4.908) でmode 0.9 sd=0.05
でベイズの階層モデルを組んでみるかな。

>>183
そのβ分布を弱情報事前分布に設定して、乱数発生させて計算すると

54000人で99%信頼区間の幅の分布は

> summary(s2w(54000))
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.008144 0.009912 0.009981 0.009927 0.010005 0.010011

となるから、まあ、概ねあっていると思うな。