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

廿TT

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

ワイブル分布のパラメータ推定(Stan vs survreg)

R 生存時間分析 Stan

※初公開時は Stan のコードが間違ってたせいでパラメータがうまく求まっていなかった。(修正:7/28)

ワイブル分布のパラメータ推定は意外とむずかしくって、打ち切りが多いときとかパラメータの数が多いときはシェイプパラメータが変なところに飛んでっちゃうことがよくある。

サンプルサイズが小さくて最尤法ではパラメータの値が不安定になるとき、MCMC 使ってベイズ推定するといいといううわさ( MCMC(Markov Chain Monte Carlo、マルコフ連鎖モンテカルロ法)について - iAnalysis 〜おとうさんの解析日記〜 )なので、はじめて Stan をさわってみた。

Stan のコードは、
http://blog.gepuro.net/archives/105
から丸写しさせていただく。すみません。

data {
  int<lower=1> J; // number of data
  int status[J];
  real<lower=0> tt[J];
}
parameters {
  real<lower=0> m; 
  real<lower=0> eta;
}
model {
  for(j in 1:J){
    if(status[j] == 1){
      increment_log_prob(weibull_log(tt[j], m, eta));
    }
    if(status[j] == 0){
      increment_log_prob(weibull_ccdf_log(tt[j], m, eta));
    }
  }
}

R によるシミュレーションのコードはこう。中央値をパラメータの推定値とした。

library(survival)
library(rstan)

n <- 100
it <- 100

m1 <- m2 <- eta1 <- eta2 <- numeric(it)

stanmodel <- stan_model("weibull.stan")

for(i in 1:it){
  ###データ生成###
  s <-rweibull(n,shape=.5, scale=5) #生存時間
  flag <- s > 10
  s[flag] <- 10 #定時打ち切り
  d <- as.numeric(!flag)
  data4stan <- list(J = n, tt = s, status = d)
  ################
  fit <- sampling(stanmodel, data=data4stan)
  ex <- extract(fit)
  
  m1[i] <- median(ex$m)
  eta1[i] <- median(ex$eta)
  
  res1 <-survreg(Surv(s,d)~1,dist="weibull")
  m2[i] <-1/res1$scale
  eta2[i] <-exp(res1$coefficients)
}
#
boxplot(data.frame("stan"=m1,"survreg"=m2),main="shape parameter")
abline(h=0.5,lty=2)
boxplot(data.frame("stan"=eta1,"survreg"=eta2),main="scale parameter")
abline(h=5,lty=2)

f:id:abrahamcow:20150728103341p:plain

f:id:abrahamcow:20150728103353p:plain

点線がシミュレーションで仮定した真値。

> mean(m1)
[1] 0.4991982
> mean(m2)
[1] 0.5019523
> mean(eta1)
[1] 5.622107
> mean(eta2)
[1] 5.240672

どうも最尤法のほうが若干優秀……? もっとパラメータを増やしたら Stan の真価が発揮されるかもしれない。