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

廿TT

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

離散時間データからのワイブル分布のパラメータの最尤推定

尤度

イベント発生が区間 (L_i,R_i] に起こったことがわかっており、イベントが発生した時間そのものはわからない状況を考える。

このようなデータを区間打ち切り(interval censored)データと呼ぶ。

尤度は、

 \displaystyle \prod_{i=1}^n \left[ \{1-F(L_i)\} - \{ 1-F(R_i) \} \right]

である。

シミュレーション

R の survival パッケージでは Turnbull のアルゴリズムを用いて区間打ち切りデータの生存関数のノンパラメトリックな推定量を求めることができる。

パラメトリックな推定は自分で尤度を書いてやる必要がある。

f:id:abrahamcow:20150927153616p:plain

図は Turnbull のアルゴリズムを用いたノンパラメトリックな経験生存関数にワイブル分布の生存関数を重ねてプロットしたもの。

###データ生成###
  set.seed(1234)
  z <- rbinom(100,1,.5)
  Z <- cbind(1,z)
  Beta <- c(1,10)
  y <-rweibull(100,5,Z%*%Beta)
  y <- floor(y)
###############

ll <- function(par){
  a <- par[1]
  Beta <- par[2:3]
  b <- Z%*%Beta
  sum(log(
    pweibull(y,a,b,lower.tail = FALSE)-
      pweibull(y+1,a,b,lower.tail = FALSE)
  ))
}
res <-optim(c(1,1,1),ll,control = list(fnscale=-1))
library(survival)
sf1 <-survfit(Surv(y,y+1,type="interval2")~z)
plot(sf1,col=c("tomato","royalblue"),lwd=2)
curve(pweibull(x,res$par[1],res$par[2],lower.tail = FALSE)
      ,add=TRUE,n=1000,col="tomato",lty=2, lwd=2)
curve(pweibull(x,res$par[1],sum(res$par[2:3]),lower.tail = FALSE),
      add=TRUE,n=1000,col="royalblue",lty=2,lwd=2)
legend("topright",c("z=0","z=1"),lty=1,col=c("tomato","royalblue"),lwd=2)
> res
$par
[1] 4.7095238 0.9871803 9.8832596

$value
[1] -140.8783

$counts
function gradient 
     242       NA 

$convergence
[1] 0

$message
NULL