シグモイド型成長曲線のパラメータを閉じた形で最尤推定する

モデル

次のような確率モデルを考える.

$\displaystyle \Delta I_t \sim {\mathrm Binomial}( S_{t-1}, bI_{t_1})$

ここで $\Delta$ は差分を意味する記号であり, $\Delta I_t =I_t - I_{t-1}$ を表す.
また ${\mathrm Binomial}( S_{t-1}, p)$ は試行回数パラメータ $S_t$ , 成功確率パラメータ $p$ の二項分布を表し, 記号 $\sim$ （実は読み方を知らない. 「チルダ」でいいのか？）は左辺の $\Delta I_t$ が右辺の二項分布に従う確率変数であることを表す.

イメージしやすくするため, $t$ は時間を表す添字であり, $S_t$ は $t$ 時点で, ある感染症に感染するリスクがある集団の人数, $I_t$ は $t$ 時点で病気にかかっていて人に病気を移す可能性がある集団の人数だと思ってほしい.

（ただし今日のこの日記の目的は感染症の数理モデルを解説することではないので, 過度な期待はしないほうがよい.）

このモデルは「感染せずに残っている $S_{t-1}$ 人が, 確率 $bI_{t_1}$ で感染者になる」というもので, 感染者が1人から2人に増えたら, 感染する確率も2倍になるだろうという考えで作られた.

一方, 次の微分方程式を考える.

$\displaystyle \frac{d}{dt} S(t) = -b S(t)I(t)$
$\displaystyle \frac{d}{dt} I(t) = b S(t)I(t)$

総人口 $N$ を一定としたので, $S(t)=N-I(t)$ であり, 組の方程式を一つにまとめることができる.

$\displaystyle \frac{dI}{dt}I(t) = b (N-I(t))I(t)$

これを解くと, 次の関数が得られる.

$\displaystyle I(t)=\frac{N}{1+\frac{N-I_0}{I_0}\exp(-Nbt)}$

微分方程式によるモデルと二項分布のモデルは違うものであるが, パラメータ $b=0.15/N$ , 0 時点めの感染者数 $I_0=1$ , 総人口 $N=1000=S_t+I_t$ と数字を適当に決め, 両方プロットしてみよう.

黒い実線が乱数による $I_t$ , 青い点線が微分方程式の解である.

累積でなく日毎の新規感染者数 $Delta I_t$ に関心があることも多いだろう. 微分方程式の解を微分して,

$\displaystyle \frac{dI(t)}{dt} = \frac{N\frac{N-I_0}{I_0}\exp(-Nbt)}{(1+\frac{N-I_0}{I_0}\exp(-Nbt))^2}$

もプロットしてみる.

黒い点が乱数による $\Delta I_t$ , 青い実践が微分方程式の解によるものである.

この場合は両者でそこそこ近い値が得られている（運がわるいと大きくずれることもある）.

最尤推定

時刻 0 の感染者数と総人口( $=S_t+I_t$ )はわかっているものとする. パラメータ $b$ をデータから推定したい.

そこで最尤法を用いることにする. 対数尤度関数はモデルから直ちに,

$\displaystyle l(b) = \sum_{t=1}^T [ x_t \log (b I_{t-1}) - S_{t} \log (1-b I_{t-1}) ]+ {\mathrm const.}$

と書ける. 記号が煩雑になるのをさけるため, $x_t = \delta I_t$ とした. この $x_t$ ( $t=1,2,\ldots,T$ ) が観測されるデータである.
またここでの ${\mathrm const.}$ はパラメータ $b$ が登場しない項という意味である.

対数尤度関数を $b$ で微分した式を0とおいて解けば最尤推定量が求まるが, 計算を楽にしたいので, 一旦 $b_t = b I_{t-1}$ とおいて $b_t$ で微分することにする.

$t$ をどこか一時点に固定して考えると, $b_t$ で微分することにより, $t$ 時点め以外の項は消える.

$\displaystyle d l(b) / db_t= x_t /b_t - (S_{t-1} - x_t )/ (1-b_t)$

$S_t =S_{t-1} - x_t$ であるから,

$\displaystyle d l(b) / db_t= x_t /b_t - S_{t} / (1-b_t)$

である.

$d l(b) / db = (d l(b) / db_t) (db_t/db)$ であり, $db_t/db = I_{t-1}$ なので,

$\displaystyle d l(b) / db = [ x_t /b_t - S_{t} / (1-b_t)] I_{t-1}$

である. これを0とおいて解きたいが $t=1,2,\ldots,T$ と動くので,

$\displaystyle x_t /b_t - S_{t} / (1-b_t)= 0\\ \displaystyle x_t (1-b_t) - S_{t} b_t]= 0 \\ \displaystyle b_t( x_t + S_{t})= x_t \\ \displaystyle bI_{t-1}( S_{t-1})= x_t$

という方程式が複数できる. 全部足し合わせて, 最尤推定量は

$\displaystyle \hat b= \sum_{t \in \tau}\frac{x_t}{I_{t-1}( S_{t-1})}$

である. ここで, $\tau$ は $\{1,2,\ldots,T\}$ から $I_{t-1}=0$ または $S_{t-1}=0$ になる場合を除いた集合を指す.
（0=0という方程式には情報がないので）

アニメーション

パラメータ $b=0.15/N$ , 0 時点めの感染者数 $I_0=1$ , 総人口 $N=1000(=S_t+I_t)$ と数字を適当に決め, 時点 $t$ までの情報で $b$ を最尤推定し, 微分方程式の解に代入してプロットしてみよう.

成長曲線の推定あるあるとして, このように単純なモデルでも, 成長曲線がサチってくる点（変曲点）を超えるまで未来の予測はなかなかうまくいかないというのがある.
（本当は流行初期に未来のことを知りたいのに…）

アニメが苦手な人のために, 時間を色の濃淡に変えた2次元の図も載せておく.

パラメータの推定値の時点ごとの履歴は次のようであった.

おしまい.

コード

Rです.

library(ggplot2)
library(gganimate)

simSigmoid <- function(Tim,N,I0,b,seed=1){
  set.seed(seed)
  xt <- integer(Tim) #Delta I_t
  St <- integer(Tim)
  It <- integer(Tim)
  St[1] <- N-I0
  It[1] <- I0
  for(t in 1:(Tim-1)){
    xt[t] <- rbinom(1,St[t],b*It[t])
    It[t+1] <- It[t] + xt[t]
    St[t+1] <- St[t] - xt[t]
  }
  data.frame(time=0:(Tim-1),St=St,It=It,xt=xt)
}

estparam <- function(St,It,xt){
  flag <- St!=0 & It!=0
  b_est = sum(xt[flag])/sum(St[flag]*It[flag])
  return(b_est)
}

Sigmoid <- function(t,N,I_0,b){
  N/(1+((N-I_0)/(I_0))*exp(-N*b*t))
}

dSigmoid <- function(t,N,I_0,b){
  A <- N*b
  B <- ((N-I_0)/(I_0))
  (N*B*A*exp(-A*t))/(1+B*exp(-A*t))^2
}

dat <- simSigmoid(100, 1000, 1, 0.15/1000, seed=123)

ggplot(dat, aes(x=time, y=It))+
  geom_line()+
  stat_function(fun = Sigmoid, args=list(N=1000,I_0=1,b=0.15/1000), 
                colour="royalblue", linetype=2)+
  labs(y="Infection(cumulative)")+
  theme_bw(16)
ggsave("fixcum.png")

ggplot(dat, aes(x=time, y=xt))+
  geom_point()+
  stat_function(fun = dSigmoid, args=list(N=1000,I_0=1,b=0.15/1000), 
                colour="royalblue")+
  labs(y="Infection(by date)")+
  theme_bw(16)
ggsave("fix.png")

bhat <- sapply(1:100, function(t)estparam(dat$St[1:t],dat$It[1:t],dat$xt[1:t]))
predm <- sapply(1:100, function(t)dSigmoid(1:100, 1000, 1, bhat[t]))

df <- data.frame(time=1:100, xt = dat$xt,
           pred=as.vector(predm), time_dummy=rep(1:100,each=100))

ggplot(df, aes(x=time, group=time_dummy))+
  geom_line(aes(y=pred, colour=time_dummy), alpha=0.3)+
  geom_point(aes(y=xt))+
  labs(y="Infection(by date)", colour="time")+
  theme_bw(16)
ggsave("inf.png")

ggplot(df, aes(x=time, group=time_dummy))+
  geom_vline(aes(xintercept = time_dummy), linetype=2)+
  geom_point(aes(y=xt))+
  geom_line(aes(y=pred), colour="royalblue")+
  labs(y="Infection(by date)")+
  theme_bw(16)+
  transition_manual(time_dummy)

anim_save("inf.gif")

dfb <- data.frame(time=1:100, estimates=bhat)
ggplot(dfb,aes(x=time, y=estimates))+
  geom_hline(yintercept = 0.15/1000, linetype=2)+
  geom_point()+
  theme_bw(16)
ggsave("est.png")