廿TT

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

区間打ち切りデータからの新型コロナウイルスの潜伏期間の推定(rstan版)

今日の川柳

追記:計算ミスがあったので修正しました(4月14日)

Stanのコードは GitHub - abikoushi/Abe2020: This is supplement scripts for Abe (2020) に上がってます。

背景の整理

新型肺炎COVID-19 の潜伏期間をrstanで推定する - 驚異のアニヲタ社会復帰の予備 で紹介されている論文

Incubation period of 2019 novel coronavirus (2019-nCoV) infections among travellers from Wuhan, China, 20-28 January 2020. - PubMed - NCBI

について書きます。以下,この論文を Backer (2020) と略します。

分析対象となるデータには,感染者が新型コロナウイルスのある地域への暴露を開始した日付と終了した日付,発症した日付が記録されている。

ここから新型コロナウイルスの潜伏期間を推定するのが目的。

感染した日はわからないが,暴露を開始した日付と終了した日付の間のどこかに落ちる。

患者  i が発症した日を原点 0 として, 暴露が終了した日を  t_i 日前とする。 暴露が終了した日からさかのぼって  s_i 日前が感染した日だとする。

f:id:abrahamcow:20200321143627p:plain

この記法のもとで,もし感染した日がわかっている完全な観測が得られた場合, 発症までの待ち時間の確率密度は  f(s_i+t_i) である。

Turnbull のアルゴリズム(R の survival パッケージ)を用いた新型コロナウイルスの潜伏期間の推定 - 廿TT では  s_i は未観測なので積分消去した。

感染者  i新型コロナウイルスのある地域へ暴露していた期間を  u_i とすると,ある患者 i の観測についての尤度は,

\int_0 ^{u_i} f(s_i+t_i) \,ds_i

である。 s_i+t_i をあらためて x_i とおくとこの積分

\int_{t_i} ^{u_i+t_i} f(x_i) \,dx_i

と書き換えられる。

これは,生存時間分析の分野で区間打ち切りとよばれる観測と同じ。

Backer (2020) では積分を明示的に計算しないで, 未観測の s_i をサンプリングする方法をとっているが本質的にはやってることは同じになるはず。

この投稿では積分を明示的に計算する方法で求めた結果を公開する。

実装上の注意

区間打ち切りデータには3種類の観測が存在する。

1つ目は区間が一点になっている場合,つまり u_i が0のとき。

この場合尤度を構成する因子は確率密度そのものになる。

2つ目は区間が普通の区間のとき, つまり u_i が有限の正の実数のとき。

尤度を構成する因子の積分を計算すると

\int_{t_i} ^{u_i+t_i} f(x_i) \,dx_i = F(u_i+t_i) - F(t_i)

になる。ただし F(x) は分布関数。実装上はこの引き算をlog_diff_exp関数を使って計算して,アンダーフローやオーバーフローを避けている。

3つ目は区間の右端が不明のとき。

この場合は無限大まで積分することになるので,尤度を構成する因子は

 \lim _ {u_i \to \infty } \int_{t_i} ^ {u_i+t_i} f(x_i) \,dx_i = 1 - F(t_i)

になる。

モデル比較

最尤法の場合

パラメータの数が同じなので,対数尤度の大きいモデルを選ぶことが AICBIC でモデルを選ぶのと同じになる。

rstanではoptimizing関数を使うことで最尤推定やMAP推定が行える。

モデル 対数尤度
ワイブル分布 -34.99795
ガンマ分布 -34.77475
対数正規分布 -34.7193

対数正規分布がもっとも良い。supplementにあるデータを使っているが、Backer (2020)とちょうど逆の結果になった。

ただ、対数正規分布とガンマ分布は五十歩百歩というところ。

MCMCによるベイズ推定の場合

パラメータの数が同じなので単に対数尤度の平均を計算するだけでいい気がするが, ここでは Backer (2020) と同じく loo パッケージを使って looic を計算する。

モデル looic
ワイブル分布 73.86812
ガンマ分布 73.12318
対数正規分布 73.47403

ついでにWAICも載せておく。loo パッケージを使った。

モデル WAIC
ワイブル分布 73.84497
ガンマ分布 73.12318
対数正規分布 73.4934

潜伏期間の予測分布

最尤法の場合

正規分布により近似した予測分布のヒストグラムを次に示す。

f:id:abrahamcow:20200414215837p:plain

予測区間を次に示す。

2.5% 50% 97.5%
Weibull 2.59 6.73 11.58
gamma 2.81 6.59 13.72
log-normal 2.89 6.52 15.78

ベイズ法の場合

MCMCにより近似した予測分布のヒストグラムを次に示す。

f:id:abrahamcow:20200414215819p:plain

予測区間を次に示す。

2.5% 50% 97.5%
Weibull 2.54 6.97 12.48
gamma 2.91 6.51 12.59
log-normal 2.65 6.74 18.39

R のコード

loo パッケージ使うのはじめてだしよくわかっていません。

library(tidyverse)
library(readxl)
library(cowplot)
library(rstan)
library(loo)
rstan_options(auto_write = TRUE)

mod_wei <- stan_model("weibull.stan")
mod_gam <- stan_model("gamma.stan")
mod_lnorm <- stan_model("lognormal.stan")

#データは論文のSuppelemental materialsにある
data_raw <- read_tsv(file = "~/Downloads/20-00062_BAKER_SupplementMaterial/BAKER_Supplementary material S1_data.tsv")

data2 <- data_raw %>%
  mutate(exposure = as.integer(as.Date(exposure_end,format="%m/%d/%Y")-as.Date(exposure_start,format="%m/%d/%Y")),
         incubation = as.integer(as.Date(symptom_onset,format="%m/%d/%Y") - as.Date(exposure_end,format="%m/%d/%Y"))) %>% 
  dplyr::filter(incubation > 0)

dat4stan <- list(
  N=nrow(data2),
  d=ifelse(is.na(data2$exposure),0,ifelse(data2$exposure==0,1,2)),
  tE=ifelse(is.na(data2$exposure),0,data2$exposure),
  tI=data2$incubation
)

fit_wei <- optimizing(mod_wei,dat4stan,seed=1,draw=10000)


fit_gam <- optimizing(mod_gam,dat4stan,seed=1,draw=10000)


fit_lnorm <- optimizing(mod_lnorm,dat4stan,seed=1,draw=10000)

fit_wei$value
fit_gam$value
fit_lnorm$value

pred_wei <- fit_wei$theta_tilde[,"pred"]
pred_gam <- fit_gam$theta_tilde[,"pred"]
pred_lnorm <- fit_lnorm$theta_tilde[,"pred"]

dists <- rev(c("Weibull","gamma","log-normal"))
dfpred <- data.frame(dist=rep(factor(dists,levels = dists),each=10000),
                     pred=c(pred_wei,pred_gam,pred_lnorm))


ggplot(dfpred,aes(x=pred),size=1)+
  geom_histogram(bins = 100,position = "identity")+
  facet_grid(dist~.)+
  theme_cowplot(24)+
  theme(strip.text.y = element_text(angle=0),strip.background = element_blank())

round(t(simplify2array(tapply(dfpred$pred, dfpred$dist, quantile, prob=c(0.025,0.5,0.975)))),2)

####
#ここからMCMC
smp_wei <- sampling(mod_wei,dat4stan,seed=1)
traceplot(smp_wei)
summary(smp_wei)$summary


smp_gam <- sampling(mod_gam,dat4stan,seed=1)
traceplot(smp_gam)
summary(smp_gam)$summary


smp_lnorm <- sampling(mod_lnorm,dat4stan,seed=1)
traceplot(smp_lnorm)
summary(smp_lnorm)$summary

loo_wei <- loo(extract_log_lik(smp_wei))
loo_gam <- loo(extract_log_lik(smp_gam))
loo_lnorm <- loo(extract_log_lik(smp_lnorm))

loo_wei$estimates[3,1]
loo_gam$estimates[3,1]
loo_lnorm$estimates[3,1]

waic_wei <- waic(extract_log_lik(smp_wei))
waic_gam <- waic(extract_log_lik(smp_gam))
waic_lnorm <- waic(extract_log_lik(smp_lnorm))

waic_wei$estimates[3,1]
waic_gam$estimates[3,1]
waic_lnorm$estimates[3,1]

pred_wei <- extract(smp_wei)$pred
pred_gam <- extract(smp_gam)$pred
pred_lnorm <- extract(smp_lnorm)$pred

dists <- rev(c("Weibull","gamma","log-normal"))
dfpred <- data.frame(dist=rep(factor(dists,levels = dists),each=4000),
           pred=c(pred_wei,pred_gam,pred_lnorm))

ggplot(dfpred,aes(x=pred))+
  geom_histogram(bins = 100)+
  facet_grid(dist~.)+
  theme_cowplot(24)+
  theme(strip.text.y = element_text(angle=0),strip.background = element_blank())

round(t(simplify2array(tapply(dfpred$pred, dfpred$dist, quantile, prob=c(0.025,0.5,0.975)))),2)