読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

SIR モデルと非定常ポアソン過程

これを読んで、疫学の人は(最尤法の枠組みの中で)非定常ポアソン過程を使ったりはしないんだろうか、と思った。

この記事では最尤法の枠組みの中で非定常ポアソン過程を推してみたい。

非定常ポアソン過程のなんたるかについては、
http://ebsa.ism.ac.jp/ebooks/sites/default/files/ebook/1868/pdf/vol2_ch3.pdf
などを参照。

SIR モデルのなんたるかについては前に一回書いたことがある(感染症のモデル(SIRモデル)に入門した - 廿TT)が、メモ的な羅列をしてるだけなので、ウィキペディアの記事かなにかを読んだほうがいいと思う。

MCMC を使った推定は、

に具体例がある。

微分方程式モデルをベイズモデル化する利点は、ぼくの理解だと、

  1. 解ががちゃがちゃした感じになる微分方程式でもパラメータが推定できる
  2. 状態変数に確率的な変動を乗せられるからよりデータにフィットする

の二点くらいだろうか。

一点目:SIR モデルから派生したモデルはいっぱいあって、その中には解ががちゃがちゃした感じになるのもあるだろうけど、シンプル SIR モデルなら最尤法でも推定は難しくない(はず)。

二点目:状態変数に確率的な変動を乗せるとその分モデルが柔らかくなるけど複雑になって定性的な理解が難しくなる。非定常ポアソン過程は動きが固いけど微分方程式をシンプルにそのまま扱える。

分析対象

データは2015年韓国におけるMERSの流行 - Wikipedia
からとって来ました。

MERS2015.csv · GitHub に置いておきます。

一応、スクレイピングしたコードも書いておきます(R です)。

library(rvest)
library(dplyr)
library(tidyr)
url1 <-"https://ja.wikipedia.org/wiki/2015%E5%B9%B4%E9%9F%93%E5%9B%BD%E3%81%AB%E3%81%8A%E3%81%91%E3%82%8BMERS%E3%81%AE%E6%B5%81%E8%A1%8C"
html1 <-read_html(url1)
tab1 <-html_table(html1)
tmp1 <-sub("人.+","",tab1[[2]]$`感染者数[32]`) 
tmp2 <-sub("人","",tmp1)
len <-length(tmp2)
Infects <-as.integer(tmp2[-c(len-1,len)])
tmp1 <-sub("^.+年","",tab1[[2]]$日付[-c(len-1,len)])
tmp2 <- sub("月","/",tmp1)
tmp3 <- sub("日","",tmp2)
Date <-as.Date(paste("2015",tmp3,sep="/"))
df1 <- data_frame(Date,Infects)
df2 <- data_frame(Date=seq.Date(as.Date("2015/05/20"),as.Date("2015/07/28"),by="1 day"))
df3 <- left_join(df2,df1,by="Date") 
out <-df3$Infects
for(i in 1:length(out)){
  if(is.na(out[i])){
    out[i] <- tmp
  }else{
    tmp <- out[i]
  }
}
df3$Recovered <-out
dat <- select(df3,Date,Recovered)
write.csv(dat,"MERS2015.csv",row.names = FALSE)

ウィキペディアではデータは日別で「感染者数」として記録されているが、用語に注意。

「感染者数」は(おそらく)病院に来た人の数なので、SIR モデルの R に対応する。

f:id:abrahamcow:20170422104814p:plain

グレーが未観測、白が観測される変数。

推定

非定常ポアソン過程の尤度関数については、abrahamcow.hatenablog.com に書いたことがある。

非定常ポアソン過程の累積強度関数に SIR モデルの R を入れたのが今回扱うモデルとその尤度関数です。

{\begin{aligned}{\frac  {d}{dt}}S(t)&=-\beta S(t)I(t)\\{\frac  {d}{dt}}I(t)&=\beta S(t)I(t)-\gamma I(t)\\{\frac  {d}{dt}}R(t)&=\gamma I(t)\end{aligned}}

未知パラメータは観戦可能者の初期状態 S(0)\beta\gamma です。

I(0)=1R(0)=0 としました。

R については解析的には解けないので、deSolve パッケージの ode 関数を使って解きます。

対数尤度関数の最大化は optim にまかせます。

dat <-read.csv(file = "MERS2015.csv",stringsAsFactors = FALSE)

library(deSolve)
SIR <- function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
    dS <- - beta*S*I
    dI <- beta*S*I - gamma*I
    dR <- gamma*I
    list(c(dS, dI, dR))
  })
}

dR <-diff(dat$Recovered)
len <- length(dR)

#対数尤度関数
llfunc <- function(parms){
  parms <- exp(parms)
  out <-ode(y=c(S=parms[1],I=1,R=0),times=0:len, func=SIR,
            parms=parms[2:3])
  pred =diff(out[,"R"])
  sum(dR*log(pred)-pred)
}

opt <- optim(par=c(5,beta=-1,gamma=-1),llfunc,control = list(fnscale=-1), hessian = TRUE)

当てはめプロット

library(cowplot)
p1 <-ggplot(dat,aes(x=as.Date(Date),y=Recovered))+
  geom_point()+geom_line()+xlab("")+
  scale_x_date(date_labels ="%m/%d")

#print(p1)

out <-ode(y=c(S=exp(opt$par[1]),I=1,R=0),times=0:(len+1), func=SIR,
          parms=exp(opt$par))

p1 + geom_line(aes(y=out[-1,"R"],colour="2"),size=1)+
  theme(legend.position = "none")

観測値が黒、推定された R が赤です。

f:id:abrahamcow:20170408061003p:plain

当てはまりはいまいちだなーと思った。

けれども、SIR モデルがそのままうまく当てはまってしまうとしたら、それは「感染症対策をなりゆきに任せた」ということになるのではないでしょうか。

モデルが当てはまらないのは、良いことと思うべきではないでしょうか。

感染症の研究してる人はそのへんどう考えてるんだろう。提案モデルが素晴らしくって感染症の伝播を完璧に予測して、感染症の拡大を防いでいたとしたら、振り返ってデータを見ると予測が外れたことになってしまう?)


S、I、R のそれぞれの動きはこんな風でした。

f:id:abrahamcow:20170408061301p:plain


考察を書きたいけど、ちょっと疲れてきたので、今回はここで終わります。

とりあえず、こんなふうにすれば推定ができますという話でした。

続くかも。

微分方程式で数学モデルを作ろう

微分方程式で数学モデルを作ろう

付録

tikz 練習

\documentclass[dvipdfmx]{standalone}
\usepackage{tikz}
\usetikzlibrary{shapes, positioning,arrows,automata}
\begin{document}
\begin{tikzpicture}[->, auto,thick]
    \node[state, fill=gray, text=white] (S) {$S(t)$};
    \node[state, fill=gray, text=white][right = of S] (I)  {$I(t)$};
    \node[state][right = of I] (R)  {$R(t)$};
    \path[->](S) edge  node {$\beta$}(I);
    \path[->](I) edge  node {$\gamma$}(R);
\end{tikzpicture}
\end{document}