読んだ。
A Simulation on Potential Secondary Spread of Novel Coronavirus in an Exported Country Using a Stochastic Epidemic SEIR Model. - PubMed - NCBI
COI：筆者はこの著者とは直接の関係はないので、純粋に統計解析のツッコミです。
こんなツイーヨを観測した。
A Simulation on Potential Secondary Spread of Novel Coronavirus in an Exported Country Using a Stochastic Epidemic SEIR Model. - PubMed - NCBI https://t.co/P2FQHeJkcX
— 岩田健太郎 Kentaro Iwata (@georgebest1969) 2020年4月7日
やっていることとしては、SEIRモデルでもし日本が
SEIRモデルは、パラメータ
として、これを
しかし、 の部分は定義されていない（というか不要）なパラメータ
がついていて、正しくは
である。
は[0.01, 0.02], [0.04, 0.06], [0.1, 0.2], [0.4, 0.6], [0.8, 1]
は[1/4, 1/2, [1/10, 1/4], [1/14, 1/10]
は[1/2, 1], [1/7, 1/2], [1/10, 1/7]
の5*3*3=45 パターンについて検討し、再生産数についても推定する。
これ、 のステップは0.2 刻みで、100日までシミュレーションした、と書いてあるが、離散的にやっているのに
とはどういうことなん…?よくわからん。
rstanでSEIRモデルをシミュレーションするには、ここらへんを参考にスクリプトをパクる。
観測したデータはないので、model ブロックでのサンプリングはない。代わりに、parametersブロックでbeta, delta, nu のそれぞれが（一様分布から）シミュレーション的にサンプリングされ、そのパラメータと、Y0に入力する初期値（,
,
,
のこと）に従って、rstan内では
integrate_ode関数がなんやかんや微分方程式を解きながら数値解をy_hatに格納してくれる、ということになっている。
statmodeling.hatenablog.com
Contemporary statistical inference for infectious disease models using Stan. - PubMed - NCBI
functions { real[] sir( real t, real[] y, real[] theta, real[] x_r, int[] x_i ) { real dydt[4]; dydt[1] = - theta[1]*y[1]*y[3]; dydt[2] = theta[1]*y[1]*y[3] - theta[2]*y[2]; dydt[3] = theta[2]*y[2] - theta[3]*y[3]; dydt[4] = theta[3]*y[3]; return dydt; } } data { int<lower=1> T; // 時間の数 real T0; // time 0 real TS[T]; // 各時刻 real Y0[4]; // Y の初期値 SEIR real b01[2]; // beta の範囲 real d01[2]; // delta の範囲 real n01[2]; // nu の範囲 } transformed data { real x_r[0]; int x_i[0]; } parameters { real<lower=b01[1], upper=b01[2]> beta; real<lower=d01[1], upper=d01[2]> delta; real<lower=n01[1], upper=n01[2]> nu; } transformed parameters { real y_hat[T,4]; real theta[3]; theta[1] = beta; theta[2] = delta; theta[3] = nu; y_hat = integrate_ode(sir, Y0, T0, TS, theta, x_r, x_i); } model { // 何もしない }
さて、論文でいうところの最悪のシナリオである、 [0.8, 1],
[1/14, 1/10],
[1/10, 1/7] をやってみる。
論文ではFigure 1 のはずだが、なんか全然違う。
の幅と日付かと思ったが、それでも
と
の数が全然違うので、なんかスクリプトがおかしいのだろう。ちなみに
だった。
は
の多さを表すので、これが大きいと最悪のシナリオになるのは当然だと思うが、論文っぽさを出すなら
[0.001, 0.002] くらいじゃないかと思った。ただしこれだと
で何もしなくても勝手に感染は収束する。
ベストなシナリオだと言っている [0.01, 0.02],
[1/4, 1/2],
[1/2, 1] だが、これだと
はほとんど減らず、
も逆にまったく増えないので、感染は爆発しない。はずだが、自分がやるとまたおかしい。Figure 2 に相当するはずだった。どちらかというとfigure 4 に近い。
本当はスクリプトをもっと精査したかったが、ぶっちゃけて言えば生データを使っているわけでなくパラメータいじっただけなんだからスクリプトくらいsupple で出したら?って思った。
rstanでパラメータを与えた時の、モデルの推定範囲が出せるようになったのでよかった。
library(rstan) rstan_options(auto_write=TRUE) options(mc.cores=parallel::detectCores()) code <- " functions { real[] sir( real t, real[] y, real[] theta, real[] x_r, int[] x_i ) { real dydt[4]; dydt[1] = - theta[1]*y[1]*y[3]; dydt[2] = theta[1]*y[1]*y[3] - theta[2]*y[2]; dydt[3] = theta[2]*y[2] - theta[3]*y[3]; dydt[4] = theta[3]*y[3]; return dydt; } } data { int<lower=1> T; // 時間の数 real T0; // time 0 real TS[T]; // 各時刻 real Y0[4]; // Y の初期値 SEIR real b01[2]; // beta の範囲 real d01[2]; // delta の範囲 real n01[2]; // nu の範囲 } transformed data { real x_r[0]; int x_i[0]; } parameters { real<lower=b01[1], upper=b01[2]> beta; real<lower=d01[1], upper=d01[2]> delta; real<lower=n01[1], upper=n01[2]> nu; } transformed parameters { real y_hat[T,4]; real theta[3]; theta[1] = beta; theta[2] = delta; theta[3] = nu; y_hat = integrate_ode(sir, Y0, T0, TS, theta, x_r, x_i); } model { // 何もしない } " # コンパイルしておく m0 <- stan_model(model_code=code) y0 <- c(S=999, E=1, I=0, R=0) # 初期値 N <- sum(y0) times <- seq(0, 100, by=1)[-1] # 時刻 0.2 刻みはどういうこと? # 各パラメータを振る範囲 lim <- list(b01=list(c(0.001, 0.002), c(0.04, 0.06), c(0.1, 0.2), c(0.4, 0.6), c(0.8, 1)), d01=list(c(1/4, 1/2), c(1/10, 1/4), c(1/14, 1/10)), n01=list(c(1/2, 1), c(1/7, 1/2), c(1/10, 1/7)) ) # for で回したい時これをindex にする b <- 1 d <- 3 n <- 3 standata <- list(T=length(times), T0=0, TS=times, Y0=y0, b01=lim$b01[[b]], d01=lim$d01[[d]], n01=lim$n01[[n]]) fit <- sampling(m0, standata, iter=200, warmup=0, chain=2, seed=1234) ex <- extract(fit, pars=head(fit@model_pars, -1)) alpha <- 0.05 cia <- c(alpha/2, 0.5, 1-alpha/2) m <- abind::abind(mapply(cbind, 0, mapply(function(z) apply(ex$y_hat[,,z], 2, quantile, cia), seq(y0), SIMPLIFY=FALSE), SIMPLIFY=FALSE), along=3) m <- abind::abind(mapply(function(z) cbind(c(0, y0[z], 0), apply(ex$y_hat[,,z], 2, quantile, cia)), seq(y0), SIMPLIFY=FALSE), along=3) parms <- apply(ex$theta, 2, quantile, cia) R0s <- ex$beta / ex$nu # 透明にしながら実線も描けるようにする cols256 <- list(S=c(244, 0, 255), E=c(255, 102, 5), I=c(0, 255, 0), R=c(0, 255, 255)) cols <- mapply(function(z) rgb(z[1], z[2], z[3], maxColorValue=256), cols256) colsa <- mapply(function(z) rgb(z[1], z[2], z[3], 20, maxColorValue=256), cols256) par(mar=c(5, 5, 3, 5), cex.lab=1.5) matplot(m[2,,], type="n", col=colsa, lty=1, lwd=3, xlab="Time (days)", ylab="Number of S/E/I/R patients") #legend(0, N, names(y0), ncol=length(y0), col=cols, pch=15, yjust=-0.1, xpd=TRUE, bty="n", cex=2) legend(par()$usr[2], N, names(y0), col=cols, pch=15, xpd=TRUE, bty="n", cex=2) for(i in seq(dim(m)[3])){ lines(c(0, times), m[2,,i], col=cols[i], lwd=4) polygon(c(c(0, times), rev(c(0, times))), c(m[1,,i], rev(m[3,,i])), col=colsa[i], border=NA) } # SEIR のパラメータ txt <- lapply(list(b1=lim$b01[[b]][1], b2=lim$b01[[b]][2], d1=lim$d01[[d]][1], d2=lim$d01[[d]][2], n1=lim$n01[[n]][1], n2=lim$n01[[n]][2]), round, 3) sub <- as.expression(substitute(beta~"["*b1*","~b2*"],"~delta~"["*d1*","~d2*"],"~nu~"["*n1*","~n2*"]", txt)) text(mean(times), par()$usr[4], sub, xpd=TRUE, pos=3, cex=1.4) # R0 R0ci <- round(quantile(R0s, cia), 3) subR0 <- as.expression(substitute(R[0]==r1~"["*b1*","~b2*"]", list(r1=R0ci[2], b1=R0ci[1], b2=R0ci[3]))) text(par()$usr[2], par()$usr[3], subR0, xpd=TRUE, adj=c(0.7, 4), xpd=TRUE, cex=1.2)
