Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2015-09-06

MikuHatsune2015-09-06

世界一わかりやすいワクチン講義(誇大広告)

炎上案件を見かけた。本当は祭りになっているときに早めに投下して鮮度を大事にしたかったけど、個人的用事で忙しいこととシミュレーションプログラムを作るのに時間がかかったこともあり、この記事の鮮度はもう0よッ…

おそらく一般人が全然知らないワクチンの効能として、ワクチンを打った自分自身が病気にかかりにくくなる、というのは当然のこととして、他人が病気にかかりにくくなり、集団全体として大流行を防ぐことができる、というのがある。これをSIR モデルという数理モデル情報科学的にR によるシミュレーションで示そう。

 

教室に規則正しく並べた机のごとく、nr^2の大きさの適当な正方形(長方形でもいいけど)を考える。1セルが1人の人に相当するし、正方形全体がとある集団、村とか都市とか、そういうのに相当する。

ここに、外から流入してきたとか、自然発生したとか、時刻0 の時点で何人かが感染症発症したとする(start_infect)。この感染症は、とある感染者1人の周囲8セルに確率論的にp¥hspace{3}(0<p¥leq1)の確率で感染する(infect_prob)。感染力が強ければpは1に近いし、弱ければ0に近い。

ここで、ワクチンを打つことで感染確率が下がるというのが個人に対するワクチンの効能である。これをp_v¥hspace{3}(0<p_v¥leq1)とし、1に近いと感染確率が例えば95% 下がる、という感じにする。このようなワクチンを、nr^2 の人口のうち vaccine_prop の割合で接種されているとする。とすると、vaccine_n 人だけ集団内にワクチン接種された人がいて、それが誰かは v の行列である。

 

ここで、SIR モデルが何かというと、適当にググればいくらか説明がある。これは微分方程式として解が求められるが、世界一わかりやすい()を標榜する記事としてはシミュレーション可視化する。

感染した人は i_day で決められた期間、感染力を持つ。現実にはこの期間は症状がなかったり、症状があったりと感染症によって異なる。感染した人は感染したら二度と感染しないとすれば、R に移行する。もちろんこのあたりは何回も反復感染するとか、潜伏期間と発病期間が違うとかいろいろやればいいが、面倒なので割愛する。

 

シミュレーションに必要な条件定義は終わったので、適当に数値をいれてシミュレーションしよう。イメージとしては、インフルエンザシーズンを考えている。100^2=10000人の集団に、突発的に20人のインフルエンザ発症した、もしくは持ち込まれたとする。インフルエンザ自体の感染力は、隣の人に対してp=0.7の確率で起こるとする。

ワクチンを打ったとすれば、インフルエンザ感染確率を90% 低下させる、とする。つまり、ワクチンを打った人は、隣の人からインフルエンザにかかる確率が 0.7*(1-0.9)=0.07 となる。

インフルエンザにかかってしまったとしたら、3日間は周囲の人に感染させる危険があるとする。ただし、感染してから4日目以降は治療もしくは勝手に治癒したとして、しかもこのシーズン中には二度と感染しないし感染させないとする(AとB型で2回かかるじゃん? という質問。残念、これはシミュレーションでしたァ!!)

 

ここで、この集団内でワクチン接種率が非常に高い状況を考えよう。この集団では95% もの人たちがワクチン接種をしているとすると、集団の10% くらいは感染してしまうが、全体としてはほとんどの人が感染せずに感染収束する。また、ワクチンを接種していない人も20% 程度の人の感染で収まる。

f:id:MikuHatsune:20150907153024g:image

というわけで

ワクチン接種しなくても感染したことがない

という人は、免疫力があるから、というのもそうだろうけど、実はまわりの人が感染させないようにしているからかかっていない、ということを知ってもいいんじゃないかな。

 

次に、ワクチン接種率が少し低い、70% の状況を考える。このとき、感染している人の赤と黄色の波がどんどん広がって収集がつかなくなっている。最終的には80% もの人が感染したことになる。

30日前後の時点で、感染力をもつ感染者数は全体の6% 程度存在する。一般人口において実際にこれくらいの割合がインフルエンザになると、病院の救急外来はパンクすると思ってもらってかまわない。

f:id:MikuHatsune:20150907153025g:image

 

ワクチン接種していない人を考える。ワクチン接種しなかった30% の人はほぼ99% 感染してしまっている。

ワクチン接種した人を考える。ワクチン接種した人の70% 近くが、ワクチン接種をしたにも関わらず感染してしまっている。ワクチン接種により感染を防ぐことができる確率が90% だったが、70% は感染してしまったので実質 30% にまで低下してしまっている。

というわけで

ワクチン接種しても感染症にかかる

という人は、ワクチンがそもそも効いていない(効きにくい)、というのもそうだろうけど、まわりの人から感染するチャンスが多すぎて守られていないから、ということを知ってもいいんじゃないかな。

 

というわけで、できることならば、接種割合は100% にしたい、というのが数理モデルからの結論である。

しかしながら、ワクチン接種で副作用あるじゃん、だから接種したくない、という気持ちも分からないわけではない。これについては、ワクチン接種による副作用が出る確率と、それによる健康被害(DALY や金額など) をシミュレーションに組み込めば、まあ計算できないこともないけれども今回は面倒なのでやらない。また、unwise decision という単語があったけれども、このような、ワクチンを接種しない、という集団が、大流行を阻止するにあたってどれだけ許容されるか、もシミュレーションでは分かる。今回は95% の接種割合のときには大流行を阻止できていたので、逆に言えば5% は接種していない集団がいてもどうにかなりうる。

ワクチンを提供する側から言えば、今回のパラメータワクチンの効能と感染確率がある。いいワクチンを作れば90% の効果から99% にできるだろうし、感染確率はたとえば早期治療介入とか、ベタだが手洗い・うがいの徹底で周辺の8マスへの感染確率を下げることが大事である。

nr <- 150 # 自乗が人口
i_day <- 3 # 感染力のある日数
obs_time <- 100 # 観察期間
start_infect <- 20 # 最初の感染者数
infect_prob <- 0.7 # 感染力
vaccine_power <- 0.9 # ワクチンによって感染力がどれくらい下がるかの効果
vaccine_prop <- 0.950 # ワクチン接種している人の割合

mat <- mat0 <- v <- diag(0, nr+2)
mat[c(1, nr+2),] <- 1
mat[,c(1, nr+2)] <- 1
idx0 <- which(c(mat) == 1)   # 周辺のいらない項
idx1 <- seq((nr+2)^2)[-idx0] # 実際の人の項
mat[idx0] <- 0

vaccine_n <- floor(vaccine_prop*nr^2) # ワクチン接種している人数
v[sample(idx1, vaccine_n)] <- 1 # ワクチン接種した人


###

sir <- list(s=v, i=mat, r=mat, v=mat)
idxr <- which(sir$r > 0) # 既感染の id
sir$i[sample(setdiff(idx1, which(v == 1)), start_infect)] <- 1 # 最初の感染者
res <- matrix(0, nr=obs_time, nc=length(sir))
colnames(res) <- names(sir)
b <- matrix(0, nr=obs_time, nc=5)
xl <- c(1, obs_time)
yl <- c(0, 1)
library(animation)
#saveGIF(
for(j in seq(obs_time)){
  infected <- setdiff(which(sir$i > 0), idxr) # 現感染者の id
  idx8 <- unlist(cbind(mapply(function(x) infected+x, -1:1), cbind.data.frame(mapply(function(k) mapply(function(j) infected+(nr+j)*k, 1:3), c(1, -1), SIMPLIFY=FALSE))), use.name=FALSE) # 現感染者の周囲8人
  idx8 <- idx8[!(idx8 %in% union(idx0, c(infected, idxr)))] # 周辺項と既感染者を除く
  a <- ifelse(
    tapply(mapply(function(x)
    sample(0:1, size=1, prob=c(1, 0)+c(-1, 1)*infect_prob*ifelse(x==1, 1-vaccine_power, 1)),
    v[idx8]), idx8, sum)>0, 1, 0) # 現感染者から周囲に感染したか

  sir$i <- ifelse(sir$i>0, sir$i+1, 0) # 感染している人に感染時間を足す
  sir$i[as.numeric(names(a))[a>0]] <- 1 # 新たに感染した人

  sir$r[sir$i > i_day] <- 1 # 既感染
  idxr <- which(sir$r > 0) # 既感染の id 

  sir$v[which((sir$i>0) * (v>0) > 0)] <- 1 # ワクチン打ったのに感染した人

  k <- mapply(function(x) ifelse(x[idx1]>1, 1, x[idx1]), sir) # 人数の計算

  k[,"i"] <- replace(k[,"i"], k[,"r"]==1, 0) # 現感染者と既感染者をかぶらないように
  res[j,] <- colMeans(k)
  b[j,] <- c(sum(res[j,c("i","r")]), res[j,c("i","r")], mean(k[k[,"s"]==1,"v"]), mean(k[k[,"s"]==0, "r"]))
  flag <- rowSums(sweep(k[,-4] ,2, 2^(0:2), "*"))
  cols <- c("white", "lightgreen", "red", "yellow", "black", grey(0.8))
  # ワクチンを打たずに感染してない人
  # ワクチンを打ってて感染してない人
  # ワクチンを打たずに感染した人
  # ワクチンを打ってて感染した人
  # ワクチンを打たずに既感染の人
  # ワクチンを打ってて既感染の人
  mat0[idx1] <- flag
  mat0[1, 1:length(cols)] <- seq(cols)-1
  xy <- seq(nr+2)
  par(mfrow=c(1, 2))
  image(xy, xy, mat0, col=cols, frame=F, axes=F, xlab="", ylab="")
  title(paste("Vaccination", round(vaccine_prop, 3)))
  polygon(c(0.5,1.5,1.5,0.5), rep(c(0.5, length(cols)+0.5), each=2), col="white", border=NA)
  matplot(b[seq(j),], type="l", las=1, lwd=3, xlim=xl, ylim=yl, lty=1, xlab="Time", ylab="Proportion")
  legend("topleft", legend=c("総感染者", "現感染者", "既感染者", "ワクチン感染者", "非ワクチン感染者"), lwd=3, col=seq(ncol(b)))
}

#, movie.name = "brownian_motion.gif", interval = 0.1, nmax = 30, ani.width = 800, ani.height = 400)

 

R初心者R初心者 2017/04/25 14:58 はじめまして. 非常に勉強になりました. 麻疹や水痘など病院内感染の講義でこのスクリプトを応用したgifアニメを使用できればと考えているのですが, Rの初心者でgifアニメ出力のやり方がわかりません. もしよろしければ教えていただけないでしょうか.

MikuHatsuneMikuHatsune 2017/04/26 10:23 院内感染講義なら麻疹や水痘といった個別の感染症より手洗い標準防護策で感染伝播を防ごう!! っていうほうがまあ、はい。

どれくらいR 初心者なのかわかりませんがすくなくともR はインストールされているものと仮定します。
R コンソールに上記スクリプトをコピペすると、library(animation) 以外はエラーなく動くはずです。
gif 化したければ、
install.package(animation)
でanimation という追加アドインをネットから入手して、(自動で)インストールします。
その後、saveGIF, movie.name=... の2行頭にある # を消します。
再度コピペして実行すると、GIF ができます。実行しているOS 環境はわかりませんが、windows, max, linux いずれもtemp という中間ファイルフォルダができる(コンソールを読むとどこかに書いてあると思います)ので、適当に名前をつけかえてファイルを移動します。

GIF でなく静止画像が欲しければ、GIF を適当なソフトでばらしてもいいですが、
png(sprintf("保存先フォルダ%03d.png", j))
というコマンドをpar(mfrow=c(1,2)) の前におき、
def.off()
というコマンドをlegend("topleft", ...) の後に置くと、保存先フォルダ001.png, ... という連番でpng が保存されます。

スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト


画像認証

トラックバック - http://d.hatena.ne.jp/MikuHatsune/20150906/1441521134