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

廿TT

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

Laslett (1982) の尤度関数とワイブル分布を仮定した最尤推定のシミュレーション

Laslett (1982) の尤度関数

Laslett, G. M., (1982) The Survival Curve Under Monotone Density Constraints With Application to two-Dimensional Line Segment Processes. Biometrika, 69: pp. 153-160
JSTOR: An Error Occurred Setting Your User Cookie
ではこんなことが書いてある(と思う)。

 生存時間分析の分野では,観測されたデータに打ち切りがある場合の最尤推定について研究がなされている.
まず,観察期間 (o,o+t) 内のランダムな時間に到着した m 人の患者の生存時間分布の推定を考える.
観測されるのは変数の組
(Z_i=z_i, D_i=d_i)
↑(i=1,...,m)
である.

 Z_i
↑これはそれぞれ独立に同分布 F に従うとする.

 D_i
↑これは右打ち切りの有無を表し,観測が打ち切られた場合に 0, そうでなければ 1 の値をとる.

打ち切りの発生は
 Z_i
↑と独立に生じるものとする.

 d_i=0
↑のとき,観測値
 z_i が示すのは,患者 i の生存時間が観測値以上だったということである.

よってこの場合,尤度は

 \prod^{m}_{i=1} \left\{f(z_i) \right\}^{d_i} \left\{1-F(z_i) \right\}^{1-d_i}

となる.

この辺については,

生存時間解析

生存時間解析

↑を参照した.

 生存時間分析の分野では非常に良く知られた Kaplam-Meier 推定量は,この尤度を最大化するノンパラメトリックな推定量と見ることができる(らしい).

Statistical Models Based on Counting Processes (Springer Series in Statistics)

Statistical Models Based on Counting Processes (Springer Series in Statistics)

 次に,観察の原点 o の前に到着した l 人の患者の,o からの生存時間も観測されているものとする.これを
 (W_j=w_j, E_j=e_j)
↑(j=1, ... ,l)
と表す.

 E_i
↑これは右打ち切りの有無を表し,観測が打ち切られた場合に 0, そうでなければ 1 の値をとる.標本の大きさ(sample size)は N = m + l である.

 再生過程における定理により,時刻 o で生存している患者の余命が w 以上の割合は,o が十分大きいとき,
 G(w) = \int^{\infty}_{w} (1-F(u))/\mu \, du
となる.
この辺については,

確率過程の基礎

確率過程の基礎

↑を参照した.

ここで μ は分布 F の平均を表す.
G は生存関数であるから,密度関数は
 g(w) = - \frac{d}{dw} G(w) = ( 1-F(w) )/\mu
である.

入門・演習 数理統計

入門・演習 数理統計

 Laslett (1982) では,患者の到着が定常ポアソン過程に従うことを仮定し,m がパラメータ N=m+l, μ /(μ + t) の二項分布に従うことを指摘した上で.上記の状況で最大化すべき尤度はこれらをすべて掛け合わせたものになると述べている.
すなわち尤度 L は,

f:id:abrahamcow:20140303182311p:plain

である.

f:id:abrahamcow:20140303183239p:plain
↑「観察期間 (o,o+t) 内のランダムな時間に到着した m 人の患者の生存時間」というのはこういうイメージ.
区間からはみ出した部分は観測されない(打ち切り).

f:id:abrahamcow:20140303183254p:plain
↑「観察の原点 o の前に到着した l 人の患者の,o からの生存時間も観測されている」というのはこういうイメージ.
区間からはみ出した部分は観測されない(打ち切り).

シミュレーション

理解のためにシミュレーションをやる.
そしてノンパラメトリック最尤推定量(NPMLE)という意味がわからないので分布を仮定しちゃうことにする.

 ワイブル分布の場合,密度関数は,
 f(x) = \frac{\alpha}{b} \left( \frac{\alpha}{b} \right)^{\alpha-1} \exp\left(-\left(\frac{x}{b} \right)^{\alpha}\right)

生存関数は,
 1-F(x) = \exp\left(-\left(\frac{x}{b} \right)^{\alpha}\right)

平均は,
 \mu=b\Gamma(1 + 1/\alpha),
 \Gamma(c)=\int^{\infty}_{0}u^{c-1}e^{-u} \, du
である.

そして,
 g(x) = (1-F(x))/\mu
である.G についてはやや煩雑な形になるが,ガンマ関数を用いて
f:id:abrahamcow:20140303183022p:plain
と表すと数値計算上都合がいい.ここで
 \gamma(c,x)=\int^{x}_{0}u^{c-1}e^{-u} \, du
である.

Wikipedia英語版のほうが情報が充実してるな。単純に人口の差なのかな

 ワイブル分布にしたのは単純に,R の survreg のデフォルトがそれだったからだ.

Rによるデータサイエンス - データ解析の基礎から最新手法まで

Rによるデータサイエンス - データ解析の基礎から最新手法まで

Rと生存時間分析(2)


 比較対象として観察の原点 o の前に到着した l 人の患者の,o からの生存時間も右打ち切りの観測とみなして,最初の尤度

 \prod^{m}_{i=1} \left\{f(z_i) \right\}^{d_i} \left\{1-F(z_i) \right\}^{1-d_i}

を最大化するパラメータを推定することもやってみる.

これを簡便法(simple method)と呼ぶことにする.

 ここで,左端点で線分が途切れている状況は,左打ち切りと呼ばれるものとは異なることに注意する.
z が左打ち切りされた観測値であるという場合,z の完全な値が z 未満であることを意味する.
しかし左端点で線分が途切れている状況では,線分の完全な長さは観測値以上である.

 ここで明らかにしたいのは,

  • 簡便法による推定は問題がある
  • それが Laslett の方法によって改善する

ということである.

 最初に述べた状況を乱数を用いたシミュレーションで再現する.シミュレーションの方法を以下に述べる.

ワイブル分布の shape パラメータを 0.5,scale パラメータを 3 とした.

f:id:abrahamcow:20140303184812p:plain

R のコード:

library("ggplot2")
weibull_0.5_3 <- function(x){dweibull(x,shape=0.5, scale=3)}
ggplot(data.frame(x=c(-0.5,3)), aes(x)) +
  stat_function(fun=weibull_0.5_3, geom="line", aes(colour="weibull_0.5_3")) 

観測の幅を t=15 で固定した.

こう設定したパラメータ(真値)と,乱数から推定したパラメータ(推定値)を比較することで推定の精度を見る.

 観測の原点 o は 0 から十分遠い位置とするため o = 1000μと置いてみた.
区間 [o, o+t] の範囲に到達した乱数をデータセットとする.この擬似的なデータセットからパラメータを推定する.

シミュレーションの試行回数は100回とした.

シミュレーションの結果を箱ひげ図で示した.

f:id:abrahamcow:20140303182656p:plain
↑shape(α)の方の推定値.右側が Laslett の方法.


f:id:abrahamcow:20140303182714p:plain
↑scale(b)の方の推定値.右側が Laslett の方法.

 注目して欲しいのは簡便法による推定にバイアスが生じている点である.
対して Laslett の方法は平均的にシミュレーションで仮定した真値に収束している(点線がシミュレーションで仮定した真値).


R のコード:

simLaslett <-function(i=1){
o <- 1000
win<- 15 #t だと t 分布とかぶるので変数名 win にした。
x <- rexp(5000,rate=1) #定常ポアソンなので到着間隔は指数分布に従う。
t1 <- cumsum(x) #到着時間
s <-rweibull(5000,shape=0.5, scale=3) #生存時間
t2 <- t1 + s
      if(max(t1) < o+win){
        print("error")
      }
#
fr <- t2 > o & t1<o
br <- t1>o & t1<o+win
rc <- t2>o+win & t1<o+win
#
#w,e
typ4 =which(fr & rc)
typ3 =which(fr & !rc)
#z,d
typ2 =which(br & rc)
typ1 =which(br & !rc)
#
len1 <- length(typ4)
W2 <- rep(win,len1)
W1 <- t2[typ3]-o
W =c(W1,W2)
E =c(rep(1,length(W1)),rep(0,len1))
#
Z2 <-s[typ2]-(t2[typ2]-(o+win))
Z1 <-s[typ1]
Z =c(Z1,Z2)
D =c(rep(1,length(Z1)),rep(0,length(Z2)))
list(cbind(Z,D),cbind(W,E))
}
#
LLW <- function(z,d,w,e, win, shape, scale){
  mu = scale*gamma(1+1/shape)
  d1 = as.logical(d)
  e1 = as.logical(e)
  lfv =dweibull(z[d1],shape=shape, scale=scale, log=TRUE)
  lfv = replace(lfv, is.nan(lfv), -Inf)
  lFv =pweibull(z[!d1],shape=shape, scale=scale, lower.tail=FALSE, log=TRUE)
  lFv = replace(lFv, is.nan(lFv), -Inf)
  lgv =pweibull(w[e1],shape=shape, scale=scale, lower.tail=FALSE, log=TRUE)- log(mu)
  lgv = replace(lgv, is.nan(lgv), -Inf)
  Gv =   pgamma((w[!e1]/scale)^shape, shape=1+1/shape, lower.tail=FALSE)-
  (w[!e1]*pweibull(w[!e1], shape=shape, scale=scale, lower.tail=FALSE))/mu
  l=length(w)
  m=length(z)
  lbv = dbinom(l, prob=mu/(mu+win),size=l+m,log=TRUE)
  lbv = replace(lbv, is.nan(lbv), -Inf)
  sum(lbv) + sum(lfv) + sum(lFv) + sum(lgv) + sum(log(Gv))
}
library(survival)
#
subfunc<- function(x){
  LLW1 <- function(par){
    LLW(z=x[[1]][,1],d=x[[1]][,2],w=x[[2]][,1], e=x[[2]][,2], win=15, shape=par[1], scale=par[2])
  }
  tmp2 =rbind(x[[1]],x[[2]])
  sr <- survreg(Surv(tmp2[,1], tmp2[,2])~1)
  opts <- optim(c(1/sr$scale, sr$coefficients), fn = LLW1 ,control = list(fnscale=-1))
  return(list(Laslett=opts, simple =sr))
}
simDat <- lapply(1:100, simLaslett)
result <- lapply(simDat, subfunc)
#
library(ggplot2)
resL <- sapply(result, function(X){X$Laslett$par})
resS <- sapply(result, function(X){c(X$simple$scale, X$simple$coefficients)})
#
#shape parametor
tmp1 =  data.frame(method="Simple",value=1/resS[1,])
tmp2 =  data.frame(method="Laslett",value=resL[1,])
resShape <- rbind(tmp1,tmp2)
qp1 <- qplot(method, value,data=resShape, geom = "boxplot")
qp1 + geom_abline(intercept = 0.5, slope = 0, linetype=2)   
#
#scale parametor
tmp3 =  data.frame(method="Simple",value=resS[2,])
tmp4 =  data.frame(method="Laslett",value=resL[2,])
resScale <- rbind(tmp3,tmp4)
qp2 <- qplot(method, value,data=resScale, geom = "boxplot")
qp2+ geom_abline(intercept = 3, slope = 0, linetype=2)   

今後の課題

課題1

 ノンパラメトリック最尤推定量(NPMLE)という意味がわからない.
計数過程とかの勉強をしたら分かるかもしれない.

 あと Cox 回帰(セミパラメトリックモデル)を勉強してみたらいいのかもしれない.

Statistical Models Based on Counting Processes (Springer Series in Statistics)

Statistical Models Based on Counting Processes (Springer Series in Statistics)

Counting Processes and Survival Analysis (Wiley Series in Probability and Statistics)

Counting Processes and Survival Analysis (Wiley Series in Probability and Statistics)

確率過程の基礎

確率過程の基礎

生存時間解析

生存時間解析

課題2

 今回の図とか ggplot2 を使って描いてみた。

 積極的に使っていきたいと思ってる(ぼくの考えたグラフ三原則 - 廿TT)が,ぜんぜんつかいこなせてる気がしない.

統計データの視覚化 (Rで学ぶデータサイエンス 12)

統計データの視覚化 (Rで学ぶデータサイエンス 12)

TeXメモ

Laslett (1982) の尤度関数

\begin{align*}
L=
\left(\begin{array}{c}
m+l \\
l
\end{array}
 \right)
\left(\frac{\mu}{\mu+t} \right)^l \left(\frac{t}{\mu+t} \right)^m
\prod^{m}_{i=1} \left\{f(z_i) \right\}^{d_i} \left\{1-F(z_i) \right\}^{1-d_i}\prod^{l}_{j=1} \left\{g(w_j) \right\}^{e_j} \left\{G(w_j) \right\}^{1-e_j}
\end{align*}

ワイブルの G

\begin{align*}
G(x) &= \int ^{\infty}_{x} g(u) \, du \\
&= 1 - \frac{\gamma \left( 1 + \frac{1}{\alpha}, \left(\frac{x}{b} \right)^\alpha \right) }{\Gamma \left(1 +\frac{1}{\alpha} \right) \exp\left(-\frac{x}{b}\right)} - \frac{x\{1-F(x)\}}{\mu}
\end{align*}