2014年度マイコースプログラム

4回生K.K. @統計遺伝学分野

9月12日(金)午前

リード・フロストモデルによるシミュレーション

一応ネットでも調べておくか…と思って「リード・フロストモデル」と検索したら、自分の昨日の記事が4番目くらいに出てきて唖然とした。
ネット検索において、価値のある情報が上位にくるわけではないというのはたまに言われることだが、それを実感した。

リード・フロストモデルとは、各期間内において、

(1)集団内の人々は感染・未感染を問わずランダムにまじりあう。
(2)感染者と未感染者が接触して伝播する一定の確率が存在する。
(3)一度感染した人はもう感染しない。
(4)その集団は他の集団から隔離されている。

という仮定のもと、

{C_{t+1}=S_t(1-(1-p)^{C_t})}

{C_{t}}:時刻tにおける感染者数、{C_{t+1}}:時刻t+1における感染者数、{S_t}:時刻tにおける未感染者、{p}:ある1人の感染者が1世代時間にある1人の未感染者に感染を伝播する確率)

という式を用いるモデルである。

時間の単位は「世代時間」で、感染してから人に伝播するまでの時間である。

Ct <- function(t,N,n,p){
if(t==0){
return(0)
}else if(t==1){
return(1)
}else{
x <- 0
y <- 1 
for(i in 1:(t-1)){
x <- x+Ct(i,N,n,p)
y <- (N-x)*(1-(1-p)^y)
}
return(y)
}}

Rt <- function(t,N,n,p){
if(t==0){
return(0)
}else{
x <- c()
for(i in 0:(t-1)) x <- c(x,Ct(i,N,n,p))
sum(x)
}}

Nt <- function(t,N,n,p){
N-Rt(t,N,n,p)-Ct(t,N,n,p)
}

N <- 100
n <- 1
p <- 0.04

x <- 0:15

y <- c()
for(i in x)
y[i+1] <- Ct(i,N,n,p)

z <- c()
for(i in x)
z[i+1] <- Rt(i,N,n,p)

w <- c()
for(i in x)
w[i+1] <- Nt(i,N,n,p)

plot(x,y,xlim=c(0,15),ylim=c(0,100),xlab="時間",ylab="人数",main="p=0.04のとき",type="l",lty=1)

par(new=T)
plot(x,z,xlim=c(0,15),ylim=c(0,100),ann=F,type="l", lty=2) 

par(new=T)
plot(x,w,xlim=c(0,15),ylim=c(0,100),ann=F,type="l", lty=3) 

text(8,90,"免疫獲得者")
text(1.2,80,"感受性ある者")
text(7,20,"感染者")

f:id:kouri_don:20140912113932p:plain

f:id:kouri_don:20140912120429p:plain

一応、「ロスマンの疫学」に書いてあった通りのシミュレーションはできた。

せっかくなので午後は自分なりのアレンジを加えてちょっとやってみる。

以上