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

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

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

10: １３２人目の素数さん [sage] 2020/03/10(火) 11:16:04 ID:H1fx2jVB

感度70％特異度90％で
有病率と陽性的中率・陰性的中率の関係をグラフにしてみた。
https://i.imgur.com/HGqPv2y.jpg

陽性的中率が0.8になるのは有病率が0.63のとき

そのRのコードはこれ。
rm(list=ls())

pr2pv <- function( # prevalence to predicative value
pr    ,# prevalence
sn=0.7, # sensitibity=TP/(TP+FN)
sp=0.9) # specificity=TN/(TN+FP)
{
N=1 # polutaion million, billion,or any proper unit
si=pr*N # sick population
he=(1-pr)*N # healthy population
TP=si*sn
FN=si*(1-sn)
TN=he*sn
FP=he*(1-sn)
PPV=TP/(TP+FP)
NPV=TN/(TN+FN)
PV=c(PPV=PPV,NPV=NPV)
return(PV)
}

prev=seq(1e-7,1,length.out = 1000)
plot(prev,sapply(prev, function(x) pr2pv(x)['PPV']),bty='l',type='l',
ylab='predicative vale',xlab='prevalence(log)',main='sensitity=0.7,specificity=0.9',log='x',lwd=2)
lines(prev,sapply(prev, function(x) pr2pv(x)['NPV']),lty=3,lwd=2)
legend('center',bty='n',legend=c('Posivive Predicative Value','Negative Predicative Value'),lty=c(1,3),lwd=2)

abline(h=0.8,col='gray')
uniroot(function(x) pr2pv(x)['PPV']-0.8, c(0,1))

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

11: １３２人目の素数さん [sage] 2020/03/10(火) 12:10:41 ID:H1fx2jVB

"
SEIR MODEL
dS(t)/dt = mu*(N-S) - b*S(t)*I(t)/N - nu*S(t)
dE(t)/dt = b*S(t)I(t)/N - (mu+sig)*E(t)
dI(t)/dt = sig*E(t) - (mu+g)*I(t)
dR(t)/dt = g*I(t) - mu*R + nu*S(t)
mu:自然死亡率 b:感染率(S->I)
nu:ワクチン有効率(S->R) sig:発症率(E->I),g:回復率(I->R)
"
SEIRモデルのパラメータ

SEIR2 <- function(
# Parameters
contact_rate = 10,                     # number of contacts per day
transmission_probability = 0.01,       # transmission probability
beta = contact_rate * transmission_probability,  # tranmission rate
infectious_period = 20,                # infectious period
gamma = 1 / infectious_period,        # Prob[infected -> recovered]
latent_period = 5,                     # latent perior
sigma = 1/latent_period,               # The rate at which an exposed person becomes infective
mu = 0,                                # The natural mortality rate
nu = 0  ,  # vaccination moves people from susceptible to resistant directly, without becoming exposed or infected.

Ro = beta/gamma, # Ro - Reproductive number.

# Initial values for sub-populations.
s = 99,          # susceptible hosts
e = 0,           # exposed hosts
i = 1,           # infectious hosts
r = 0,           # recovered hosts
# Compute total population.
N = s + i + r + e,
# Output timepoints.
timepoints = seq (0, 365, by=0.5),
...
)

有病率を１％とすると、3000人にクルーズ船でも100人の屋形船でも感染者のピークは変わらないな。
同一時間あたりのcontact_rateとtransmission_probabilityが宴会での方が高いからだろうな。

パラメータを変えてグラフを書いてみた。

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

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

17: １３２人目の素数さん [sage] 2020/03/11(水) 06:41:46 ID:hVKkfTiV

PCR検査の感度を0.7、特異度を0.9とする。
広島県で第一号の感染発見例は
県の検査で１回陰性、病院の検査で２回陽性、症状軽快した現時点で陰性（何回やったか報道がないので１回陰性とする）であるという。
ここで問題：
検査前の感染確率の分布が一様分布であると仮定して、
現在患者が感染している確率とその95%CIを計算してみた。
"

https://i.imgur.com/4CzTEWA.jpg

sn=0.7 # sensitivity
sp=0.9 # specificity
plus=2 # how many positive result?
minus=2 # how many negative result?
n=1e7 # how large the simulation
p0=runif(n,0,1)
oz0=p0/(1-p0) # prob -> odds
pLR=sn/(1-sp) # TP/FP
nLR=(1-sn)/sp # FN/TN
oz1=oz0*pLR^plus*nLR^minus # Bayesian formula
p1=oz1/(1+oz1) # odds -> prob
BEST::plotPost(p1,showMode =T) # show mode instead of mean
BEST::plotPost(p1,showMode =F)
HDInterval::hdi(p1) # Highest Density Interval
quantile(p1,c(.025,0.5,.975)) # 95%CI by quantile
summary(p1) # mean, median
MAP <- function(x) {
dens <- density(x)
mode_i <- which.max(dens$y)
mode_x <- dens$x[mode_i]
mode_y <- dens$y[mode_i]
c(x=mode_x, y=mode_y)
}
MAP(p1)['x'] # show mode

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

26: １３２人目の素数さん [sage] 2020/03/11(水) 14:21:53 ID:hVKkfTiV

>>17
数値を変化させてグラフ化できるように関数化
事前分布は一様分布でなくJeffereysを採用
必要に応じて指定

"
PCR検査の感度を0.7、特異度を0.9とする。
広島県で第一号の感染発見例は
県の検査で１回陰性、病院の検査で２回陽性、症状軽快した現時点で陰性（何回やったか報道がないので１回陰性とする）であるという。

検査前の感染確率の分布が一様分布であると仮定して、
現在患者が感染している確率とその95%CIを計算してみた。
"
PCR2prob <- function(
sn=0.7, # sensitivity
sp=0.9, # specificity
plus=2, # how many positive result?
minus=2, # how many negative result?
n=1e5,
p0=rbeta(n,0.5,0.5), # prior : Jeffreys' distribution
print=TRUE) # how large the simulation
{
oz0=p0/(1-p0) # prob -> odds
pLR=sn/(1-sp) # TP/FP
nLR=(1-sn)/sp # FN/TN
oz1=oz0*pLR^plus*nLR^minus # Bayesian formula
p1=oz1/(1+oz1) # odds -> prob
if(print & length(p0)>1){ # p0 ~ some distribution
BEST::plotPost(p1,showMode =T) # show mode instead of mean
print(HDInterval::hdi(p1)) # Highest Density Interval
print(quantile(p1,c(.025,0.5,.975))) # 95%CI by quantile
print(summary(p1)) # mean, median
dens = density(p1) # print mode
mode_i = which.max(dens$y)
print(c(Mode=dens$x[mode_i]))
}
if(length(p0)==1) print(p1) # when p0 is point-designated
invisible(p1)
}

PCR2prob()
PCR2prob(p0=rbeta(1e5,1,1)) # p0 ~ uniform distiribution
PCR2prob(p0=0.5) # point probability
PCR2prob(minus=3) # one more negative result
PCR2prob(minus=4) # two more negative result

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

68: １３２人目の素数さん [sage] 2020/03/16(月) 06:09: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

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

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

ぬこの手ぬこTOP 0.023s