2重区間打ち切りされた観測についての尤度
あらまし
Reich et al. (2009) は感染症の潜伏期間を推定する際によく参照される文献で2重区間打ち切り(doubly interval censored)のある観測について論じている.
Reich さんたちは Reich et al. (2009) の考察に基づく R のパッケージ coaseDataTools
も提供しており,このパッケージでは2重区間打ち切りされた観測に対してワイブル分布,ガンマ分布,対数正規分布などを当てはめることができる.
coaseDataTools
の中身を見てみると尤度の中の積分は数値積分で評価している様子だったが,ウルフラム・アルファ を使って計算してみたらワイブル分布,ガンマ分布,対数正規分布についてはこの積分が閉じた形(といっても特殊関数は使う)で書けたのでここで報告する.
ガンマ関数などの特殊関数はそれ用の数値計算のプログラムがうまく作られていたりするので特殊関数を使って書けると計算が速くなることある(もともとそこまで時間のかかる計算ではなかったので計算を速くしたいという意欲はそんなにはないかもしれないが……).
2重区間打ち切り
準備としてまず観測の枠組みを説明する.
感染の起こった時点に対応する確率変数を
潜伏期間は
このとき尤度は
と書ける.
ここで
このようにイベントが起こったことはわかっているがイベントの時点が幅を持って観測される状況を区間打ち切りと呼ぶ.
この場合,
と書ける.
さらに,感染した時点
Reich et al. (2009) はこのような観測を2重区間打ち切り(doubly interval censored)と呼んだ.
と書ける. この記事ではこの尤度を具体的に計算してみる.
(ところで,感染の打ち切りは発生するだろうなと思うけど発症した日は大体わかりそうな気がする.現実のどういう場面で発症日が打ち切りになるのかいまいちわかってないのでだれか教えてください.)
1-F(y) の積分
準備:先の積分について具体的に計算する前に,少しだけ記号を準備しておく.
分布
あとちょっと面倒を避けるために
打ち切りデータを扱う生存時間分析の文脈だと生存関数
であり,
が成り立つ.なぜかというと,部分積分により,
が成り立つからである.最後の変形で第1項が 0 になるのは,
であり,
であることからわかる.
ここから生存関数を
とすると,
この確率密度を持つ分布を Ross (1995) "Stochastic Processes" にならい,均衡分布(equibillium distribution)と呼ぶことにする.
均衡分布
ガンマ分布,ワイブル分布,対数正規分布について,均衡分布の分布関数を求めてみたところ以下に示すような形になった.
興味のある方はウルフラム・アルファ を開いて integral[ Gamma(a, x/b) , x]
などと打ち込んでみてほしい.
分布関数を使って書けたので(この記事では扱わないが)Stan などの確率的プログラミング言語でも実装しやすいのではないかと思う.
ガンマ分布の場合
とする. 形状パラメータ
である.
正規化されたガンマ関数はガンマ分布の分布関数で尺度パラメータ
eqgamma <- function(x, shape, scale){
pgamma(x/scale, shape+1)+(x/scale)*pgamma(x/scale, shape, lower.tail=FALSE)/shape
}
ワイブル分布の場合
形状パラメータ
R での実装例は:
eqweibull <- function(x, shape, scale){
pgamma((x/scale)^(shape), 1/shape)
}
対数正規分布
とすると,対数平均パラメータ
R での実装例は:
erf <- function(x){
2*pnorm(sqrt(2)*x) - 1
}
erfc <- function(x){
2*pnorm(-sqrt(2)*x)
}
eqlnorm <- function(x, meanlog, sdlog){
ifelse(
x < 0,
0,
0.5*(1+x*erfc((log(x) - meanlog)/(sdlog*sqrt(2)))*exp(-(meanlog + (sdlog^2)/2)) - erf((meanlog + sdlog^2 - log(x))/(sdlog * sqrt(2))))
)
}
2重区間打ち切りされた観測についての尤度
再び2重区間打ち切りされた観測についての尤度に戻る.
これまでの結果を使って,2重区間打ち切りされた観測についての尤度を均衡分布で表してみる.
積分区間に注意して3つの場合にわける.
感染が起こった時点は少なくとも発症よりは前のはずなので,
独立同分布の仮定のもとで
coaseDataTools
との比較
coaseDataTools
パッケージのデモ用に入っているインフルエンザの潜伏期間のデータで,最尤推定の計算時間を比較してみる.
R のコードは少し長くなるので最後にまとめおく.
期待通り速くなり,結果はおおよそ一致した.
#上が coaseDataTools で実装されている尤度, 下は僕の実装
#ガンマ分布
> print(bm_g)
test
1 flu_opt_g0 <- optim(c(0, 1), loglikhd, dat = fluA.inc.per, dist = "G")
2 flu_opt_g <- optim(c(0, 1), lpG, dat = fluA.inc.per)
replications elapsed relative user.self sys.self user.child sys.child
1 100 36.063 5.408 35.949 0.117 0 0
2 100 6.668 1.000 6.510 0.027 0 0
> print(all.equal(flu_opt_g0$par, flu_opt_g$par))
[1] TRUE
#ワイブル分布
> print(bm_w)
test
1 flu_opt_w0 <- optim(c(0, 1), loglikhd, dat = fluA.inc.per, dist = "W")
2 flu_opt_w <- optim(c(0, 1), lpW, dat = fluA.inc.per)
replications elapsed relative user.self sys.self user.child sys.child
1 100 31.457 7.398 31.174 0.160 0 0
2 100 4.252 1.000 4.238 0.011 0 0
> print(all.equal(flu_opt_w0$par, flu_opt_w$par))
[1] TRUE
#対数正規分布
> print(bm_l)
test
1 flu_opt_l0 <- optim(c(0, 1), loglikhd, dat = fluA.inc.per, dist = "L")
2 flu_opt_l <- optim(c(0, 1), lpL, dat = fluA.inc.per)
replications elapsed relative user.self sys.self user.child sys.child
1 100 52.890 5.015 52.598 0.206 0 0
2 100 10.547 1.000 10.526 0.024 0 0
> print(all.equal(flu_opt_l0$par, flu_opt_l$par))
[1] TRUE
せっかくなので生存関数もプロットしておく.
c3 <- c("royalblue", "orangered", "forestgreen")
curve(pweibull(x, exp(flu_opt_w$par[1]), exp(flu_opt_w$par[2]), lower.tail = FALSE),
xlim=c(0,5), col=c3[1], lty=2, ylab="CCDF", xlab="incubation period")
curve(pgamma(x, exp(flu_opt_g$par[1]), scale = exp(flu_opt_g$par[2]), lower.tail = FALSE),
add=TRUE, col=c3[2], lty=3)
curve(plnorm(x, flu_opt_l$par[1], sdlog = exp(flu_opt_l$par[2]), lower.tail = FALSE),
add=TRUE, col=c3[3], lty=4)
legend("topright", legend = c("Weibull", "gamma", "lognormal"), lty=2:4, col=c3, lwd=1)
選んだモデルの範囲だとどれもそれほど大きな違いはない.
横軸の単位は日だと思うが,ウイルスは変異とかもあるしこの潜伏期間が文字通りに解釈できるかはわからない(念のための注意).
以下,R のコード:
#install.packages("coarseDataTools")
library(coarseDataTools)
library(dplyr)
library(ggplot2)
library(rbenchmark)
#gamma dist
eqgamma <- function(x, shape, scale){
pgamma(x/scale, shape+1)+(x/scale)*pgamma(x/scale, shape, lower.tail=FALSE)/shape
}
meangamma <- function(shape, scale){
scale*shape
}
#weibull dist
eqweibull <- function(x, shape, scale){
pgamma((x/scale)^(shape), 1/shape)
}
meanweibull <- function(shape, scale){
scale*gamma(1+1/shape)
}
#lognormal dist
erf <- function(x){
2*pnorm(sqrt(2)*x) - 1
}
erfc <- function(x){
2*pnorm(-sqrt(2)*x)
}
eqlnorm <- function(x, meanlog, sdlog){
ifelse(
x < 0,
0,
0.5*(1+x*erfc((log(x) - meanlog)/(sdlog*sqrt(2)))*exp(-(meanlog + (sdlog^2)/2)) - erf((meanlog + sdlog^2 - log(x))/(sdlog * sqrt(2))))
)
}
meanlnorm <- function(mu, sigma){
exp(mu+sigma^2/2)
}
#doubly interval censored
evaluatelp0 <- function(EL, ER, SL, SR, par, eqcdf, meanf){
ll = 0
mu = meanf(par)
if(ER <= SL){
ll = log(eqcdf(SR-ER, par) - eqcdf(SL-ER, par) - (eqcdf(SR-EL, par) - eqcdf(SL-EL, par))) + log(mu)
}else if (SL < ER & ER < SR){
ll = log((ER - SL) + mu*(eqcdf(SR-ER, par) - (eqcdf(SR-EL, par) - eqcdf(SL-EL, par))))
}else if(SR <= ER){
ll = log((SR - SL) - mu*(eqcdf(SR-EL, par) - eqcdf(SL-EL, par)) )
}
return(ll)
}
# single interval
evaluatelp1 <- function(EL, ER, SL, SR, par, cdf){
ll = log(cdf(SR-EL, par) - cdf(SL-EL, par))
return(ll)
}
# no censored
evaluatelp2 <- function(EL, ER, SL, SR, par, logpdf){
logpdf(SL-EL, par)
}
eval_lp <- function(par, EL, ER, SL, SR, type, dist){
if(dist=="weibull"){
eqcdf <- function(x, par){
eqweibull(x, shape = exp(par[1]), scale = exp(par[2]))
}
cdf <- function(x, par){
pweibull(x, shape = exp(par[1]), scale = exp(par[2]))
}
logpdf<- function(x, par){
dweibull(x, shape = exp(par[1]), scale = exp(par[2]), log=TRUE)
}
meanf <- function(par){
meanweibull(exp(par[1]), exp(par[2]))
}
}else if(dist == "gamma"){
eqcdf <- function(x, par){
eqgamma(x, shape = exp(par[1]), scale = exp(par[2]))
}
cdf <- function(x, par){
pgamma(x, shape = exp(par[1]), scale = exp(par[2]))
}
logpdf<- function(x, par){
dgamma(x, shape = exp(par[1]), scale = exp(par[2]), log=TRUE)
}
meanf <- function(par){
meangamma(exp(par[1]), exp(par[2]))
}
}else if(dist == "lognormal"){
eqcdf <- function(x, par){
eqlnorm(x, meanlog = par[1], sdlog = exp(par[2]))
}
cdf <- function(x, par){
plnorm(x, meanlog = par[1], sdlog = exp(par[2]))
}
logpdf<- function(x, par){
dlnorm(x, meanlog = par[1], sdlog = exp(par[2]), log=TRUE)
}
meanf <- function(par){
meanlnorm(par[1], exp(par[2]))
}
}else{
warning("this distribution is not implemented")
}
N <- length(EL)
ll = 0
for(i in 1:N){
if(type[i]==0){
ll = ll + evaluatelp0(EL[i], ER[i], SL[i], SR[i], par, eqcdf = eqcdf, meanf = meanf)
}else if(type[i]==1){
ll = ll + evaluatelp1(EL[i], ER[i], SL[i], SR[i], par, cdf = cdf)
}else if(type[i]==2){
ll = ll + evaluatelp2(EL[i], ER[i], SL[i], SR[i], par, logpdf = logpdf)
}else{
warning("this censor type is not implemented")
}
}
return(ll)
}
data(fluA.inc.per)
#head(fluA.inc.per)
lpW <- function(par, dat){
with(dat, -eval_lp(par, EL, ER, SL, SR, type = type, dist = "weibull"))
}
lpG <- function(par, dat){
with(dat, -eval_lp(par, EL, ER, SL, SR, type = type, dist = "gamma"))
}
lpL <- function(par, dat){
with(dat, -eval_lp(par, EL, ER, SL, SR, type = type, dist = "lognormal"))
}
bm_w <- benchmark(flu_opt_w0 <- optim(c(0,1), loglikhd, dat=fluA.inc.per, dist = "W"),
flu_opt_w <- optim(c(0,1), lpW, dat=fluA.inc.per), order = NULL)
print(bm_w)
print(all.equal(exp(flu_opt_w0$par),exp(flu_opt_w$par)))
bm_g <- benchmark(flu_opt_g0 <- optim(c(0,1),loglikhd,dat=fluA.inc.per,dist = "G"),
flu_opt_g <- optim(c(0,1), lpG, dat=fluA.inc.per), order = NULL)
print(bm_g)
print(all.equal(exp(flu_opt_g0$par),exp(flu_opt_g$par)))
bm_l <- benchmark(flu_opt_l0 <- optim(c(0,1),loglikhd, dat=fluA.inc.per,dist = "L"),
flu_opt_l <- optim(c(0,1), lpL, dat=fluA.inc.per), order = NULL)
print(bm_l)
print(all.equal(exp(flu_opt_g0$par),
exp(flu_opt_g$par)))
Discussion