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

[過去ﾛｸﾞ] 数学統計に詳しい人が語るコロナウイルス (1002ﾚｽ)
上下前次 1-新
通常表示 512ﾊﾞｲﾄ分割ﾚｽ栞

このｽﾚｯﾄﾞは過去ﾛｸﾞ倉庫に格納されています｡
次ｽﾚ検索歴削→次ｽﾚ栞削→次ｽﾚ過去ﾛｸﾞﾒﾆｭｰ

68: １３２人目の素数さん [sage] 2020/03/16(月) 06:09:25.25 ID:CVVw1pKV

国内で患者数が大幅に増えたときに備えた医療提供体制の確保について
今後、国内で新型コロナウイルス感染症患者数が大幅に増えたときに備え、各都道府県、保健所設置市及び特別区（以下「都道府県等」という。）における外来を受診する患者数等について、以下の数式を用いて計算いただき、ピーク時の医療需要の目安としてご活用の上、必要な医療提供体制を確保していただくようお願いいたします。
（１）（ピーク時において1 日あたり新たに新型コロナウイルス感染症を疑って外来を受診する患者数）＝（0-14 歳人口）×0.18／100＋（15-64 歳人口） ×0.29／100＋（65 歳以上人口） ×0.51／100
（２）（ピーク時において1 日あたり新型コロナウイルス感染症で入院治療が必要な患者数）＝（0-14 歳人口）×0.05／100＋（15-64 歳人口）×0.02／ 100＋（65 歳以上人口） ×0.56／100
（３）（ピーク時において1 日あたり新型コロナウイルス感染症で重症者として治療が必要な患者数）＝（0-14 歳人口）×0.002／100＋（15-64 歳人口） ×0.001／100＋（65 歳以上人口） ×0.018／100
注１）ピーク時は、各都道府県等において疫学的関連性が把握できない程度に感染が拡大した時点から概ね３か月後に到来すると推計されている。ただし、公衆衛生上の対策を行うことにより、ピークが下がるとともに後ろ倒しされる。
注２）重症者とは、集中治療や人工呼吸器を要する管理が必要な患者を指す。
注３）当該計算式は、都道府県等の単位以下における医療提供体制を確保するためのものであるとともに、各都道府県等によってピークを迎える時期が異なるため、全国の人口を用いて計算することや単純に各自治体が算出するピークの数値を足し合わせることは、不適切な取扱いとなることに留意いただきたい。なお、当該計算式については、今後新たな知見等により変更される可能性がある。
注４）実際には、ピーク時に至るまでの日々の患者数の増加はばらつきがあり、増加曲線は推計通りの形にならない可能性が高いため、現実の患者の発生動向も踏まえて適切に体制を確保することが必要。
注５）当該計算式については、今後新たな知見等により変更される可能性がある。

http://rio2016.5ch.net/test/read.cgi/math/1582910321/68

94: １３２人目の素数さん [sage] 2020/03/23(月) 03:28:01.22 ID:uvHIelYA

日本人の平均身長を推測するのにその値は１〜２ｍの間であるという弱情報事前分布は合理的。
現時点での新型コロナの有病率は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

http://rio2016.5ch.net/test/read.cgi/math/1582910321/94

101: １３２人目の素数さん [sage] 2020/03/24(火) 08:34:40.14 ID:TnHQvRcs

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

http://rio2016.5ch.net/test/read.cgi/math/1582910321/101

102: １３２人目の素数さん [sage] 2020/03/24(火) 08:42:02.80 ID:TnHQvRcs

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])
}

http://rio2016.5ch.net/test/read.cgi/math/1582910321/102

上下前次 1-新書関写板覧索設栞歴

あと 890 ﾚｽあります
ｽﾚ情報赤ﾚｽ抽出画像ﾚｽ抽出歴の未読ｽﾚ AAｻﾑﾈｲﾙ

ぬこの手ぬこTOP 0.020s