20/05/18 20:05:42 1P7V5xJn.net
>>550(自答)
#
# 人物ｄが発症してdelay日後に濃厚接触したキャバ嬢ｃが発症
# ｃの感染がｄより先行していた確率は？
rm(list=ls())
stancode=
"
data{
real onset_delay;
real ln_par1;
real ln_par2;
}
parameters{
real <lower=0> d_incubation;
real <lower=0> c_incubation;
}
transformed parameters{
real infection_delay = onset_delay + d_incubation - c_incubation;
}
model{
d_incubation ~ lognormal(ln_par1,ln_par2);
c_incubation ~ lognormal(ln_par1,ln_par2);
}
"
model=stan_model(model_code = stancode)
#--- incubation period ---
# from Li et al NEJM 2020
# lognormal mean = 5.2
ln_par1 = 1.434065
ln_par2 = 0.6612
fn.stan <- function(delay, print=FALSE, ...){
dataList=list(onset_delay=delay,ln_par1=ln_par1,ln_par2=ln_par2)
fit=sampling(model,data=dataList, ...)
ms=rstan::extract(fit)
if(print) BEST::plotPost(ms$infection_delay,compVal=0,xlab='infection delay')
mean(ms$infection_delay < 0)
}
fn.stan(2,print=T,iter=5000,warmup=1000)
onset_delays=0:20
y=sapply(onset_delays,fn.stan)
plot(onset_delays,y, ylab='Pr[ Infected Later ])')
２日後の発症だと
> fn.stan(2,print=T,iter=5000,warmup=1000)
[1] 0.2945
３割くらいはあとから発症した方が先に感染していた可能性がある。