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

廿TT

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

R によるターンブルのアルゴリズム

生存時間分析 R

勉強のためターンブル(Turnbull)のアルゴリズムを愚直に書いてみる。

生存時間解析

生存時間解析

クライン&メシュベルガー『生存時間解析』 p.17 に以下のような表が与えられている。

age event right left
10 4 0 0
11 12 0 0
12 19 2 0
13 24 15 1
14 20 24 2
15 13 18 3
16 3 14 2
17 1 6 3
18 0 0 1
>18 4 0 0

これはカリフォルニア州の高校生のマリファナ使用に関するデータだそうだ。「はじめてマリファナを使用したのはいつですか?」age は年齢(これを t_i とする)、>18 は 18 より上という意味。event はイベント発生数(d_i とする)、right は右打ち切り数(r_i)、left は左打ち切り数( c_i )。

例えばこの表の 4 行目は 13 歳ではじめてマリファナを使用したと答えた少年が 24 人、マリファナを使用したことがないと答えた 13 歳の少年が 15 人、マリファナを使用したことがあるがいくつのときだったか覚えていないと答えた少年が 1 人いることを表す。

このデータに対して、少年がはじめてマリファナを使用するまでの生存関数を推定してみる。普通、時間 t での生存関数 S(t) の値は、その時間まで生存している(故障していない/死亡していない)モノの割合を表すが、この場合はマリファナを使用していない少年の割合を表す。以下コードはすべて R。

#まずは表のデータ
#「18 より上」にはダミーで 999 を入れておく

tab117 <- data.frame(age =c(10:18,999), 
event =c(4,12,19,24,20,13,3,1,0,4),
right =c(0, 0, 2,15,24,18,14,6,0,0),
left =c(0, 0, 0, 1, 2, 3, 2,3,1,0)
)

ステップ 0:各時間(t_j とする)に対して生存関数の初期推定値 S_0(t_j) を定める。(上記の本にはここで、「いかなる論理的な推論もうまく進むようにすることに注意すべきである」と書かれているが意味不明)初期値には左打ち切りを無視したカプラン・マイヤー推定量を用いることが推奨されているらしいのでそれに従う。

カプラン・マイヤー推定量は
\displaystyle \hat S_0(t_j) = \prod^{m}_{j=1} (1-d_j/Y_j)

ここで、
\displaystyle Y_i = \sum^m_{j=1}(d_j+r_j)

第 K+1 回目の更新式は


ステップ 1j \le i に対し、条件付き確率

 p_{ij} = P( t_{j-1} < X \le t_j | X \le t_i )

\hat p_{ij} = \frac{S_K(t_{j-1})-S_K(t_j)}{1-S_K(t_i)}

で推定する。

ステップ 2:時間 t_j におけるイベント発生数を

\displaystyle \hat d_j = d_j + \sum^{m}_{i=j} c_i p_{ij}

で推定する。

ステップ 3:左打ち切りを無視し、通常のカプラン・マイヤー推定量を用いて生存関数

\displaystyle \hat S_{K+1}(t_i) = \prod^{i}_{j=1} (1-\hat d_i/\hat Y_i)

を推定する。ここで、

\displaystyle \hat Y_i = \sum^m_{j=1}(\hat d_j+r_j)

第一ステップは、こう。

num <- which(tab117$left>0)
lnum <- length(num)
Y = numeric()
Y2 =Y
S0 = numeric(10+1)
S0[1] = 1
S1 =S0
dhat <- numeric(10)
dhat[10] <- tab117$event[10]
dhat2 =dhat
pij <- matrix(rep(NA, 9*9),c(9,9))
pij2 <- pij
for(i in 1:10){ Y[i] =  sum(tab117$event[i:10] + tab117$right[i:10]) }
for(i in 1:10){ S0[i+1] = prod(1-tab117$event[1:i]/Y[1:i]) }
for(j in num){ for(i in 1:j) {pij[i,j] =(S0[i]-S0[i+1])/(1- S0[j+1]) } }
for(i in 1:9) dhat[i] <- tab117$event[i] + sum(tab117$left[num] *  pij[,num][i,], na.rm = NA)

第一ステップでの結果を見てみる。

pij の値(i:行、j:列)

NA NA NA 0.067 0.049 0.039 0.036 0.034 0.034
NA NA NA 0.202 0.146 0.116 0.108 0.102 0.102
NA NA NA 0.32 0.23 0.184 0.17 0.161 0.161
NA NA NA 0.41 0.295 0.235 0.218 0.207 0.207
NA NA NA NA 0.28 0.223 0.207 0.196 0.196
NA NA NA NA NA 0.204 0.189 0.179 0.179
NA NA NA NA NA NA 0.071 0.068 0.068
NA NA NA NA NA NA NA 0.054 0.054
NA NA NA NA NA NA NA NA 0


dhat、Y、S0 の値

dhat Y S0
4.488 179 0.977
13.464 175 0.906
21.318 163 0.794
26.971 142 0.651
22.427 103 0.516
14.704 59 0.392
3.414 28 0.345
1.215 11 0.308
0 4 0.308
4 4 0

#ここからがメインのループ

#shokichi =S0
for(k in 1:100){
for(i in 1:10){ Y2[i] =  sum(dhat[i:10] + tab117$right[i:10]) }
for(i in 1:10){ S1[i+1] = prod(1-dhat[1:i]/Y2[1:i]) }
for(j in num){ for(i in 1:j) {pij2[i,j] =(S1[i]-S1[i+1])/(1- S1[j+1]) } }
for(i in 1:9){ dhat2[i] = tab117$event[i] + sum(tab117$left[num] *  pij2[,num][i,], na.rm = NA)}
if(all(abs(S1-S0) < 0.001)) break
S0=S1
pij =pij2
dhat =dhat2
}

結果、

> i
[1] 9
> S1
 [1] 1.0000000 0.9765031 0.9060124 0.7944021 0.6513060 0.5157529 0.3921066 0.3453536
 [9] 0.3079407 0.3079407 0.0000000

となり、9 回のイテレーションで収束し S1 が生存関数の最終的な推定値である。

グラフはこう。
f:id:abrahamcow:20121209235212p:plain

#png()
plot(S1,type="S", xlab="", ylab="", xaxt="n")
points(shokichi, type="S", xlab="", ylab="", xaxt="n", lty=2)
axis(side=1,at=1:10 , labels=tab117[,1])
par(family="HiraKakuPro-W3")
legend(2,0.4, lty=1:2, legend=c("S1", "初期値"))
#dev.off()

考察
カリフォルニア州の少年は大麻吸いすぎですね。