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

廿TT

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

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

生存時間分析 R

生存時間解析

生存時間解析


この本の pp.147-150 の内容をなるべく愚直に書いてみた。
ところどころよくわかってないんだけどこれであってるんだろうか。

データセットこちらから入手できるが、R 用に加工した。

Data1_18 <-
structure(list(lower = c(0L, 0L, 0L, 0L, 0L, 4L, 4L, 4L, 5L, 
5L, 5L, 6L, 7L, 7L, 8L, 8L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 
11L, 11L, 12L, 13L, 13L, 13L, 13L, 14L, 14L, 15L, 15L, 16L, 16L, 
16L, 16L, 17L, 17L, 17L, 17L, 17L, 17L, 18L, 18L, 18L, 18L, 19L, 
19L, 21L, 22L, 22L, 23L, 24L, 24L, 24L, 24L, 25L, 26L, 27L, 30L, 
30L, 31L, 32L, 32L, 33L, 33L, 34L, 34L, 34L, 34L, 35L, 35L, 36L, 
36L, 36L, 36L, 37L, 37L, 37L, 37L, 38L, 40L, 44L, 45L, 46L, 46L, 
46L, 46L, 46L, 46L, 46L, 46L, 48L), upper = c(5L, 7L, 8L, 5L, 
22L, 11L, 8L, 9L, 11L, 12L, 8L, 10L, 14L, 16L, 12L, 21L, 17L, 
35L, 15L, 18L, 13L, 17L, 20L, NA, NA, 20L, 39L, NA, NA, NA, 17L, 
19L, 22L, NA, 20L, 24L, 24L, 60L, 25L, 25L, 23L, 26L, 27L, NA, 
26L, 24L, 25L, NA, 35L, 32L, NA, 32L, NA, NA, 30L, 31L, NA, NA, 
37L, 40L, 34L, 34L, 36L, NA, NA, NA, 40L, NA, 34L, NA, NA, NA, 
39L, NA, 44L, 48L, NA, NA, 44L, NA, NA, NA, NA, NA, 48L, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, 48L), treatment = c(0L, 0L, 0L, 
1L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 
0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 
1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 
1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 
0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L)), .Names = c("lower", 
"upper", "treatment"), row.names = c(NA, -95L), class = "data.frame")


アルゴリズム。

surv.int <- function(lower, upper){
dat <- cbind(lower, upper)
#dat <- dat_1
tau <- unique(c(dat[,1],dat[,2]))
tau <- tau[!is.na(tau)]
tau <- sort(tau)
dat <- replace(dat, is.na(dat), max(tau))
n0 <- dim(dat)[1]
m0 <- length(tau)-1
#
#initial val.
subp <- numeric()
for(i in 1:n0){subp[i] <- sum( tau > dat[i,1] & tau <= dat[i,2] )}
subp <- replace(subp, subp==0, 1)
subt <- numeric()
subd <- numeric()
for(i in 1:n0){
	wch <- which(dat[i,2]==tau)
	subt <- c(subt, tau[(wch-subp[i]+1):wch])
	subd <- c(subd, rep(1/subp[i],subp[i]))
}
#
S0 <- numeric()
for(i in 1:(length(tau)-1)){ S0[i] <- (n0 - sum(subd[subt <= tau[1+i]]))/n0 }
S0 <- c(1,S0)
#
alpha <- function(x, y){as.integer(x[,1] >= y[,1] & x[,2] <= y[,2])}# "<=" であってる
#
for(it in 1:500){
	p0 <- rev(diff(rev(S0)))
	#
	d1 <- numeric(m0)
	#
	tmp <- dat[,1]==dat[,2] #exact
	rp1 <- (1:n0)
	if(any(tmp)){ 
		wch <- which(tmp)
		exact <- table(dat[wch,1])
		d1[tau[-1] %in% names(exact)] <- exact
		rp1 <- rp1[-wch] #skip
	}
	for(i in rp1){
		bunbo = 0
		bunbo <- bunbo+sum(alpha(x=cbind(tau[1:m0], tau[(1:m0)+1]), y=cbind(dat[i,1], dat[i,2]))*p0)
		d1 <- d1 + (p0*alpha(x=cbind(tau[1:m0], tau[(1:m0)+1]), y=cbind(dat[i,1], dat[i,2])))/bunbo
	}
	d1 <- c(0,d1)
	Y1 <- rev(cumsum(rev(d1)))
	S1 <- numeric(m0+1)
	for(i in 1:(m0+1)){ S1[i] = prod(1-d1[1:i]/Y1[1:i]) }
	if(max(abs(S1 - S0)) < 1e-07) break
	S0 =S1
}
dup1 <- duplicated(round(S1,8))
#dup1 <- duplicated(S1)
return(list(result=data.frame(time=tau[!dup1], survival =S1[!dup1]), it =it))
}


例。

#source("Data1_18.R")
dat_0 <- Data1_18[Data1_18$treatment == 0, 1:2]
res0 <- surv.int(dat_0[,1],dat_0[,2])
dat_1 <- Data1_18[Data1_18$treatment == 1, 1:2]
res1 <- surv.int(dat_1[,1],dat_1[,2])

survival パッケージの survfit 関数は区間打ち切りデータにも対応しているらしい。
使いかた、これであってるんだろうか。

library(survival)
sf1 <- survfit(Surv(time=lower,time2=upper, type="interval2")~treatment, data=Data1_18)

プロット。

#png("jisaku.png")
par(family="HiraKakuPro-W3")
plot(res0$result,type="S", main="自作")
points(res1$result,type="S", lty=2)
#dev.off()
#png("survfit.png")
plot(sf1,lty=1:2, main="survfit")
#dev.off()

f:id:abrahamcow:20130303221732p:plain
f:id:abrahamcow:20130303221747p:plain

解説はいずれ。