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

廿TT

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

2重打ち切りデータについてのターンブルのアルゴリズム

以前のこれです。関数化した。

2重打ち切りは区間打ち切りとはまた別。解説はいずれ。

surv.doub <- function(time, exact, right, left){
	n0 <- length(time)
	num <- which(left>0)   
	lnum <- length(num)
	Y = numeric()
	Y2 =Y
	S0 = numeric(n0+1)
	S0[1] = 1
	S1 =S0
	dhat <- numeric(n0)
	dhat[n0] <- exact[n0]
	dhat2 =dhat
	pij <- matrix(rep(NA, n0*n0),c(n0,n0))
	pij2 <- pij
	for(i in 1:n0){ Y[i] =  sum( exact[i:n0] +  right[i:n0]) }
	for(i in 1:n0){ S0[i+1] = prod(1- exact[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:(n0-1)) dhat[i] <-  exact[i] + sum( left[num] *  pij[,num][i,], na.rm = NA)
	shokichi =S0
	#
	for(k in 1:100){
		for(i in 1:n0){ Y2[i] =  sum(dhat[i:n0] +  right[i:n0]) }
		for(i in 1:n0){ 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:n0){ dhat2[i] =  exact[i] + sum( left[num] *  pij2[,num][i,], na.rm = NA)}
		if(all(abs(S1-S0) < 0.001)) break
		S0=S1
		pij =pij2
		dhat =dhat2
	}
	return(list(result =data.frame(time=c(0, time), surv =S1), it=k))
}


生存時間解析

生存時間解析

この本では「18歳以上」の行を一行別に作っているが、それって18歳での右打ち切りと同じことじゃないの? と思った。なので「〇〇歳以上」も一つの行にまとめることにした。

tab117_2 <- data.frame(age =c(10:18), 
event =c(4,12,19,24,20,13,3,1,0),
right =c(0, 0, 2,15,24,18,14,6,4),
left =c(0, 0, 0, 1, 2, 3, 2,3,1)
)
res1 <- surv.doub(time=tab117_2[,1], exact =tab117_2[,2], right=tab117_2[,3], left =tab117_2[,4])
plot(res1$result,type="S")

結果、

> res1
$result
   time      surv
1     0 1.0000000
2    10 0.9765031
3    11 0.9060124
4    12 0.7944021
5    13 0.6513060
6    14 0.5157529
7    15 0.3921066
8    16 0.3453536
9    17 0.3079407
10   18 0.3079407

$it
[1] 2

数字はだいたい一緒なんだけどなぜかたった2回で収束してる。

f:id:abrahamcow:20130303190803p:plain