廿TT

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

Becker et al. (1991) method of non-parametric back-projection を試す

A Method of Non-Parametric Back-Projection and Its Application to AIDS Data - PubMed のテストです。

背景

感染症は潜伏期間があるために発病日がわかっても感染日はわからない。

そこで時刻 t における発症者の人数  y_t から時刻 t に感染した人の数  N_t を見積もることがしたい。

潜伏期間の分布はあらかじめわかっているものとして、潜伏期間が t+1 である確率を f_t と置く。

 N_t の期待値を  \lambda_t として、感染がポアソン過程に従っていることにすると、

 \displaystyle y_t \sim \operatorname{Poisson}\left(\sum_{i=1}^{t} \lambda_i f_{t-i+1}\right)

となる。

たとえば:

時刻tに感染して、時刻tに発症する人の数の期待値は \lambda_t f_1

時刻t-1に感染して、時刻tに発症する人の数の期待値は \lambda_{t-1} f_2

時刻t-2に感染して、時刻tに発症する人の数の期待値は \lambda_{t-2} f_3

なので、これを足し合わせたものが y_t の期待値となる。

ただし、このままだと対数尤度を計算するときに log の中に和が入っている形となってしまい、計算しにくい。

そこで生成モデルを次のように改めて考える。

 \displaystyle s_{td} \sim \operatorname{Poisson}( \lambda_d f_{t-d+1}), \quad d=1, \ldots, t

 \displaystyle y_t = \sum_{d=1}^{t} s_{td}

これは最初の式とまったく等価。

未観測の s_{td} が得られているとして、対数尤度をかくと、

\displaystyle l(\lambda_t) = \sum_{d=1}^t \{ s_{td} \log(\lambda_t f_{t-d+1}) -\lambda_t f_{t-d+1} \}

となる。

 \lambda_t に関して微分して 0 と置いて解くと、

\displaystyle \hat \lambda_t =\left( \sum_{d=1}^t s_{td} \right) / \left( \sum_{d=1}^t f_{t-d+1}\right)

これが EM アルゴリズムの M ステップ。

E ステップは s_{td} の条件付き期待値を求める。

ポアソン分布の和は和の値で条件つけてやると多項分布になるので、

 \displaystyle E[s_{td}] = y_t \lambda_t f_{t-d+1} / \left(\sum_{d=1}^{t} \lambda_t f_{t-d+1}\right)

となる。

Rによるシミュレーション

シミュレーションでは t=10,11,..20 \lambda_t = 100, 他の場合  \lambda_t = 0 とした。

f:id:abrahamcow:20200518134058p:plain

縦棒が発症者の推移の観測値、点線が推定された分布の平均。

f:id:abrahamcow:20200518134143p:plain

縦棒が未観測の感染者の推移、丸が推定された  \hat \lambda_t で、ちょっとガタガタしてるけどまあまあうまく求まっていることがわかる。

ただ、EMアルゴリズムをどこで止めるかで結果がだいぶかわる。

上記のは対数尤度の増加が0.1未満になったら収束と判定した。

収束判定を厳しくして 0.001 とかにしてみるとこんな感じ。

f:id:abrahamcow:20200518134708p:plain

f:id:abrahamcow:20200518134648p:plain

尤度は大きくなっていくんだけど  \hat \lambda_t はどんどんがたがたしてくる。

実用上はなんらかの方法で正則化とかしたほうがいい気がする。

genefun_c <- function(maxrd,maxed,mined,lambda,m,eta){
  np <- rpois(1,(maxed-mined)*lambda)
  ti <- sort(runif(np,mined,maxed))
  si <- rweibull(np,m,eta)
  rd <- ceiling(ti+si)
  id <- ceiling(ti)
  infection <- sapply(1:maxrd,function(t)sum(id==t))
  obs <- sapply(1:maxrd,function(t)sum(rd==t))
  data.frame(obs,infection)
}

lambdaest <- function(count_obs, ft, maxit=1000,tol=0.1){
  len <- length(count_obs)
  lamhat <- rep(mean(count_obs),len)
  den <-rev(cumsum(ft))
  mu <- sapply(1:len, function(t)sum(lamhat[1:t]*f[t:1]))
  ll <- numeric(maxit)
  ll[1] <- sum(dpois(count_obs,mu,log = TRUE))
  for (i in 2:maxit) {
    #E stpp
    yp <- lapply(1:len, function(t)count_obs[t]*lamhat[1:t]*ft[t:1]/sum(lamhat[1:t]*ft[t:1]))
    #M step
    num <- sapply(1:len,function(t)sum(sapply(yp, function(x)x[t]),na.rm = TRUE))
    lamhat <- (num)/(den)
    mu <- sapply(1:len, function(t)sum(lamhat[1:t]*ft[t:1]))
    ll[i] <- sum(dpois(count_obs,mu,log = TRUE))
    if(ll[i]-ll[i-1]<tol){
      break
    }
  }
  return(list(lambda=lamhat,mu=mu,it=i,ll=ll[1:i]))
}

set.seed(123)
dat <- genefun_c(40,20,10,100,2,4)
f <- diff(pweibull(0:40,2,4))
lambda_out <- lambdaest(dat$obs, f, maxit  = 100)

print(lambda_out$it)

plot(lambda_out$ll,type="l", ylab="log-likelihood")

plot(dat$obs,type="h",ylim=c(0,max(lambda_out$mu,dat$obs)),ylab = "symptom onset")
lines(lambda_out$mu,lty=2)
plot(lambda_out$lambda,ylim=c(0,max(lambda_out$lambda,dat$infection)),ylab = "infection")
points(dat$infection,type = "s",lty=2)