諦めていた。
mikuhatsune.hatenadiary.com
しかし、なんか出来そうになった。
mikuhatsune.hatenadiary.com
（
型の）コロナウイルスの流行具合から、相互に免疫があるものとして、新型肺炎の今後数年の流行がどうなっていくかをSEIRSモデルで、かつ2種類のウイルス（two strain）がいるモデルでシミュレーションした話。
最初の記事の16コンパートメントの推移は式が間違っていたようなのでがんばって修正した。
普通にRで
を足し合わせた場合の結果はこんな感じになる。
のコンパートメントは
、
のコンパートメントは
とすると、図になる。
さて、rstanを使ってパラメータの変動を考慮しつつSEIRS two strain モデルをやってみる。
一生懸命コードを書くだけであるが、16コンパートメントのsimplex を時系列T 分だけ定義する。時間依存パラメータ（ここでは
）も、時間依存でないパラメータたちも全部、ひとまず時間依存パラメータとして（つまりすべての
で一定）として、パラメータの数と同じ長さをもつvectorをもつvector で定義する。
こんな感じでできた。
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)
}
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("topright", legend=names(y0), ncol=4, col=cols, pch=15, yjust=-0.1, xpd=TRUE, bty="n", cex=1.4)
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())
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)
LineChartID105e2da79c06
Science. 2020 Apr 14. pii: eabb5793. doi: 10.1126/science.abb5793. on 2020-04-23
