諦めていた。
mikuhatsune.hatenadiary.com
しかし、なんか出来そうになった。
mikuhatsune.hatenadiary.com
（型の）コロナウイルスの流行具合から、相互に免疫があるものとして、新型肺炎の今後数年の流行がどうなっていくかをSEIRSモデルで、かつ2種類のウイルス（two strain）がいるモデルでシミュレーションした話。
最初の記事の16コンパートメントの推移は式が間違っていたようなのでがんばって修正した。
普通にRでを足し合わせた場合の結果はこんな感じになる。
のコンパートメントは
、
のコンパートメントは
とすると、図になる。
さて、rstanを使ってパラメータの変動を考慮しつつSEIRS two strain モデルをやってみる。
一生懸命コードを書くだけであるが、16コンパートメントのsimplex を時系列T 分だけ定義する。時間依存パラメータ（ここでは）も、時間依存でないパラメータたちも全部、ひとまず時間依存パラメータとして（つまりすべての
で一定）として、パラメータの数と同じ長さをもつ
vectorをもつvector で定義する。
こんな感じでできた。
# SEIRS をシミュレーション # mu beta nu x12 x21 gamma delta1 delta2 library(stringr) labs <- c(outer(paste0(lab, 1), paste0(lab, 2), paste0)) times <- 0:200 y <- rep(0, 16) y[c(1, 2, 5)] <- c(0.997, 0.002, 0.001) Y <- y p <- c(1/80, 2.2, 1/3, 0.78, 0.51, 1/5, 1/45, 1/45) for(i in 2:length(times)){ dydt <- rep(0, length(y)) I1 <- (y[3]+y[7]+y[11]+y[15]) I2 <- (y[9]+y[10]+y[11]+y[12]) dydt[1 ] = -p[2]*(I1+I2)*y[1] + p[7]*y[4] + p[8]*y[13] + p[1] - y[1]*p[1]; dydt[2 ] = -(p[1]+p[3]+(1-p[4])*p[2]*I2)*y[2] + p[2]*I1*y[1] + p[8]*y[14]; dydt[3 ] = -(p[1]+p[6]+(1-p[4])*p[2]*I2)*y[3] + p[3]*y[2] + p[8]*y[15]; dydt[4 ] = -(p[1]+p[7]+(1-p[4])*p[2]*I2)*y[4] + p[6]*y[3] + p[8]*y[16]; dydt[5 ] = -(p[1]+p[3]+(1-p[5])*p[2]*I1)*y[5] + p[2]*I2*y[1] + p[7]*y[8]; dydt[9 ] = -(p[1]+p[6]+(1-p[5])*p[2]*I1)*y[9] + p[3]*y[5] + p[7]*y[12]; dydt[13] = -(p[1]+p[8]+(1-p[5])*p[2]*I1)*y[13] + p[6]*y[9] + p[7]*y[16]; dydt[6 ] = -(2*p[3]+p[1])*y[6] + p[2]*((1-p[4])*I2*y[2] + (1-p[5])*I1*y[5]); dydt[7 ] = -(p[1]+p[3]+p[6])*y[7] + (1-p[4])*p[2]*I2*y[3] + p[3]*y[6]; dydt[8 ] = -(p[1]+p[3]+p[7])*y[8] + (1-p[4])*p[2]*I2*y[4] + p[6]*y[7]; dydt[10] = -(p[1]+p[3]+p[6])*y[10] + (1-p[5])*p[2]*I1*y[9] + p[3]*y[6]; dydt[14] = -(p[1]+p[3]+p[8])*y[14] + (1-p[5])*p[2]*I1*y[13] + p[6]*y[10]; dydt[11] = -(2*p[6]+p[1])*y[11] + p[3]*(y[7]+y[10]); dydt[12] = -(p[1]+p[6]+p[7])*y[12] + p[3]*y[8] + p[6]*y[11]; dydt[15] = -(p[1]+p[6]+p[8])*y[15] + p[3]*y[14] + p[6]*y[11]; dydt[16] = -(p[1]+p[7]+p[8])*y[16] + p[6]*(y[12]+y[15]); Y <- rbind(Y, y+dydt) y <- y + dydt } i1 <- rep(1:4, 4) i2 <- rep(1:4, each=4) i12 <- list(i1, i2) Y0 <- c(list(Y), mapply(function(z) t(apply(Y, 1, tapply, z, sum)), i12, SIMPLIFY=FALSE)) cols <- sample(na.omit(str_extract(colors(), "^(?!.*gr[ea]y).*$")), length(labs)) lwd <- 3 par(mfrow=c(1, 3), las=1, mar=c(4, 5, 2, 2), cex.lab=1.5) matplot(Y0[[1]], type="l", lty=1, lwd=lwd, xlab="Time", ylab="Proportion", col=cols) legend("topright", legend=labs, ncol=4, col=cols, pch=15, xpd=TRUE, bty="n", cex=1.5) ci <- tapply(sample(seq(labs), 8), rep(1:2, each=4), c) for(i in 2:3){ matplot(Y0[[i]], type="l", lty=1, lwd=lwd, col=cols[ci[[i-1]]], xlab="Time", ylab="Proportion") legend("topright", legend=sprintf("%s%d", lab, i-1), ncol=4, col=cols[ci[[i-1]]], pch=15, xpd=TRUE, bty="n", cex=1.5) } # beta をb(t) にしたうえでSEIRS library(rstan) rstan_options(auto_write=TRUE) options(mc.cores=parallel::detectCores()) code <- " data { int<lower=1> T; // 時間の数 int<lower=1> nC; // コンパートメントの数 real TS[T]; // 各時刻 real Y0[nC]; // Y の初期値 SEIR real m01[2]; // mu の範囲 real n01[2]; // nu の範囲 real x01[2]; // x12 の範囲 real x02[2]; // x21 の範囲 real g01[2]; // gamma の範囲 real d01[2]; // delta1 の範囲 real d02[2]; // delta2 の範囲 real p01[2]; // phi の範囲 real f01[2]; // f の範囲 real R01[2]; // R0 の範囲 } parameters { real<lower=m01[1], upper=m01[2]> mu; real<lower=n01[1], upper=n01[2]> nu; real<lower=x01[1], upper=x01[2]> x12; real<lower=x02[1], upper=x02[2]> x21; real<lower=g01[1], upper=g01[2]> gamma; real<lower=d01[1], upper=d01[2]> delta1; real<lower=d02[1], upper=d02[2]> delta2; real<lower=p01[1], upper=p01[2]> phi; real<lower=f01[1], upper=f01[2]> f; real<lower=R01[1], upper=R01[2]> R0; } transformed parameters { simplex[16] y[T]; vector<lower=0>[8] p[T]; real<lower=0> I1[T]; real<lower=0> I2[T]; y[1] = to_vector(Y0); // Initialize for(i in 1:T){ p[i][1] = mu; p[i][2] = gamma * R0 * (f/2*cos(2*pi()/52*(TS[i]-phi)) + (1-f/2)); // beta p[i][3] = nu; p[i][4] = x12; p[i][5] = x21; p[i][6] = gamma; p[i][7] = delta1; p[i][8] = delta2; } for(i in 1:(T-1)){ // S1S2 E1S2 I1S2 R1S2 1-4 // S1E2 E1E2 I1E2 R1E2 5-8 // S1I2 E1I2 I1I2 R1I2 9-12 // S1R2 E1R2 I1R2 R1R2 13-16 I1[i] = (y[i][3]+y[i][7] +y[i][11]+y[i][15]); I2[i] = (y[i][9]+y[i][10]+y[i][11]+y[i][12]); y[i+1][1 ] = y[i][1 ] + -p[i][2]*(I1[i]+I2[i])*y[i][1] + p[i][7]*y[i][4] + p[i][8]*y[i][13] + p[i][1] - y[i][1]*p[i][1]; y[i+1][2 ] = y[i][2 ] + -(p[i][1]+p[i][3]+(1-p[i][4])*p[i][2]*I2[i])*y[i][2] + p[i][2]*I1[i]*y[i][1] + p[i][8]*y[i][14]; y[i+1][3 ] = y[i][3 ] + -(p[i][1]+p[i][6]+(1-p[i][4])*p[i][2]*I2[i])*y[i][3] + p[i][3]*y[i][2] + p[i][8]*y[i][15]; y[i+1][4 ] = y[i][4 ] + -(p[i][1]+p[i][7]+(1-p[i][4])*p[i][2]*I2[i])*y[i][4] + p[i][6]*y[i][3] + p[i][8]*y[i][16]; y[i+1][5 ] = y[i][5 ] + -(p[i][1]+p[i][3]+(1-p[i][5])*p[i][2]*I1[i])*y[i][5] + p[i][2]*I2[i]*y[i][1] + p[i][7]*y[i][8]; y[i+1][9 ] = y[i][9 ] + -(p[i][1]+p[i][6]+(1-p[i][5])*p[i][2]*I1[i])*y[i][9] + p[i][3]*y[i][5] + p[i][7]*y[i][12]; y[i+1][13] = y[i][13] + -(p[i][1]+p[i][8]+(1-p[i][5])*p[i][2]*I1[i])*y[i][13] + p[i][6]*y[i][9] + p[i][7]*y[i][16]; y[i+1][6 ] = y[i][6 ] + -(2*p[i][3]+p[i][1])*y[i][6] + p[i][2]*((1-p[i][4])*I2[i]*y[i][2] + (1-p[i][5])*I1[i]*y[i][5]); y[i+1][7 ] = y[i][7 ] + -(p[i][1]+p[i][3]+p[i][6])*y[i][7] + (1-p[i][4])*p[i][2]*I2[i]*y[i][3] + p[i][3]*y[i][6]; y[i+1][8 ] = y[i][8 ] + -(p[i][1]+p[i][3]+p[i][7])*y[i][8] + (1-p[i][4])*p[i][2]*I2[i]*y[i][4] + p[i][6]*y[i][7]; y[i+1][10] = y[i][10] + -(p[i][1]+p[i][3]+p[i][6])*y[i][10] + (1-p[i][5])*p[i][2]*I1[i]*y[i][9] + p[i][3]*y[i][6]; y[i+1][14] = y[i][14] + -(p[i][1]+p[i][3]+p[i][8])*y[i][14] + (1-p[i][5])*p[i][2]*I1[i]*y[i][13] + p[i][6]*y[i][10]; y[i+1][11] = y[i][11] + -(2*p[i][6]+p[i][1])*y[i][11] + p[i][3]*(y[i][7]+y[i][10]); y[i+1][12] = y[i][12] + -(p[i][1]+p[i][6]+p[i][7])*y[i][12] + p[i][3]*y[i][8] + p[i][6]*y[i][11]; y[i+1][15] = y[i][15] + -(p[i][1]+p[i][6]+p[i][8])*y[i][15] + p[i][3]*y[i][14] + p[i][6]*y[i][11]; y[i+1][16] = y[i][16] + -(p[i][1]+p[i][7]+p[i][8])*y[i][16] + p[i][6]*(y[i][12]+y[i][15]); } I1[T] = (y[T][3]+y[T][7] +y[T][11]+y[T][15]); I2[T] = (y[T][9]+y[T][10]+y[T][11]+y[T][12]); } model { mu ~ uniform(m01[1], m01[2]); nu ~ uniform(n01[1], n01[2]); x12 ~ uniform(x01[1], x01[2]); x21 ~ uniform(x02[1], x02[2]); gamma ~ uniform(g01[1], g01[2]); delta1 ~ uniform(d01[1], d01[2]); delta2 ~ uniform(d02[1], d02[2]); phi ~ uniform(p01[1], p01[2]); f ~ uniform(f01[1], f01[2]); R0 ~ uniform(R01[1], R01[2]); } " m0 <- stan_model(model_code=code) times <- 0:200 y0 <- rep(0, 16) y0[c(2, 5)] <- c(0.002, 0.001) y0[1] <- 1-sum(y0) standata <- list(T=length(times), TS=times, Y0=y0, nC=length(y0), m01=1/c(81, 79), n01=1/c(14, 3), x01=c(0, 1), x02=c(0, 1), g01=1/c(14, 3), d01=1/c(100, 25), d02=1/c(100, 25), f01=c(0.18, 0.25), p01=c(1.5, 2.4), R01=c(2, 3)) standata <- list(T=length(times), TS=times, Y0=y0, nC=length(y0), m01=1/c(81, 79), n01=1/c(14, 3), x01=c(0.5, 1), x02=c(0, 0.5), g01=1/c(14, 3), d01=1/c(50, 20), d02=1/c(30, 25), f01=c(0.18, 0.25), p01=c(1.5, 2.4), R01=c(3, 5)) fit <- sampling(m0, standata, iter=1000, warmup=0, chain=4) 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(function(z) cbind(c(0, y0[z], 0), apply(ex$y[,,z], 2, quantile, cia)), seq(y0), SIMPLIFY=FALSE), along=3) lab <- c("S", "E", "I", "R") labs <- c(outer(paste0(lab, 1), paste0(lab, 2), paste0)) names(y0) <- labs # 透明にしながら実線も描けるようにする cols256 <- list(S1S2='lightgreen', E1S2='yellow', I1S2='orange', R1S2='blueviolet', S1E2='yellow2', E1E2='yellow3', I1E2='orange2', R1E2='cadetblue', S1I2='orangered', E1I2='orangered2', I1I2='orangered3', R1I2='deepskyblue', S1R2='powderblue', E1R2='royalblue', I1R2='skyblue', R1R2='pink') cols256 <- lapply(lapply(cols256, col2rgb), c) 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, 2, 2), cex.lab=1.5, las=1) matplot(m[2,,], type="n", col=colsa, lty=1, lwd=3, xlab="Time (days)", ylab="Proportion of S/E/I/R patients") #legend(0, 1, names(y0), ncol=length(y0), col=cols, pch=15, yjust=-0.1, xpd=TRUE, bty="n", cex=1.6) legend("topright", legend=names(y0), ncol=4, col=cols, pch=15, yjust=-0.1, xpd=TRUE, bty="n", cex=1.4) #legend(par()$usr[2], 1, names(y0), col=cols, pch=15, xpd=TRUE, bty="n", cex=1.6) 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) } ### library(googleVis) library(stringr) gmat <- cbind.data.frame(Time=times, head(m[2,,], -1)) colnames(gmat) <- c("Times", labs) url <- "https://science.sciencemag.org/content/early/2020/04/14/science.abb5793" webpage <- "Science. 2020 Apr 14. pii: eabb5793. doi: 10.1126/science.abb5793." headertxt <- sprintf("<span><a href='%s' target='_blank'>%s</a></span> on %s", url, webpage, Sys.Date()) # googlevis 用のプロットオプション ops <- list( fontSize=12, height=600, width=600, lineWidth=1, vAxis="{title:'Proportion', fontSize: 20, titleTextStyle: {fontSize: 20, bold: 'true', italic: 'false'}, viewWindow: {min: 0, max: 1}, gridlines: {color: '#00F', minSpacing: 100}, logScale: 'false'}", sizeAxis="{minValue: -80, maxValue: 75, maxSize: 15}", hAxis="{title:'Time', fontSize: 20, titleTextStyle: {fontSize: 20, bold: 'true', italic: 'false'}, viewWindow: {min: 0, max: 205}}", explorer="{actions: ['dragToZoom', 'rightClickToReset']}", focusTarget="category", selectionMode="multiple", legend="{position: 'right', textStyle: {color: 'black', fontSize: 10}, pageIndex: 1}" ) gvis <- gvisLineChart(gmat, options=ops) gvis$html$caption <- NULL gvis$html$footer <- str_replace(gvis$html$footer, "<span> \n .+\n .+\n</span>", headertxt) plot(gvis)
Science. 2020 Apr 14. pii: eabb5793. doi: 10.1126/science.abb5793. on 2020-04-23
