까먹을까봐 적어놓는 블로그

백혈병 데이터

중도절단 시점을 제외한 그룹1과 그룹2의 관측 시점은 다음과 같이 17개이다.

관측 시점 => 1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 13, 15, 16, 17, 22, 23 (17개)

예) 7시점에서 관측 도수 분할표

 예를 들어, 7 시점에서 그룹1의 절단(failed) 사람의 수 $d_{17}$는 1명이고, 그룹2에서는 7시점에서 절단되지 않았으므로 $d_{27}$은 0이다. 또한, 그룹1의 7 시점에서 위험에 처한(at risk) 사람의 수 $Y_{17}$은 7시점을 포함하여 17명이고, 7 시점에서 그룹 2의 위험 사람의 수 $Y_{27}$는 12명이다.

위와 같은 분할표를 각 관측 시점에 대하여 만들면 다음과 같은 표를 작성하는 것이 편리하다.

  $d_{1i}$과 $d_{i2}$는 각각 관측시점 $i$에서 그룹1과 그룹2의 절단 수이고, $Y_{i1}$과 $Y_{i2}$는 각각 관측시점 $i$에서 그룹1과 그룹2의 위험(at risk) 수이다.

로그-순위 검정통계량이 $ \chi ^2 _{0.05}(1)=3.84 $보다 크므로 두 집단의 생존 곡선의 차이가 있다.

Python의 lifelines를 이용한 로그 순위검정

from lifelines.statistics import logrank_test
import numpy as np

time1 = np.array([6, 6, 6, 6, 7, 9, 10, 10, 11, 13, 16, 17, 19, 20, 22, 23, 25, 25, 32, 32, 34])
status1 = np.array([1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0])

time2 = np.array([1, 1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 8, 8, 11, 11, 12, 12, 15, 17, 22, 23])
status2 = np.repeat(1, 21)

result = logrank_test(time1, time2,
                      status1, status2)
print(result.test_statistic)
print(result.p_value)

> 16.792940989216547

> 4.168809109334511e-05

from lifelines import KaplanMeierFitter
import matplotlib.pyplot as plt

kmf = KaplanMeierFitter(alpha = 0.05)
fig = plt.figure(figsize = (10, 10))

kmf.fit(time1, status1, label = 'Group 1')
kmf.plot()
kmf.fit(time2, status2, label = 'Group 2')
kmf.plot()
plt.xlabel('time')
plt.ylabel('surv prop')
plt.show()

R의 survival 패키지를 이용한 로그순위 검정

library(survival)
time1 <- c(6, 6, 6, 6, 7, 9, 10, 10, 11, 13, 
           16, 17, 19, 20, 22, 23, 25, 25, 32, 32, 34)
status1 <- c(1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 
             1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0)

time2 <- c(1, 1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 8, 8, 
           11, 11, 12, 12, 15, 17, 22, 23)
status2 <- rep(1, 21)

time <- c(time1, time2)
status <- c(status1, status2)
group <- c(rep(1, 21), rep(2, 21)) # 1: trt, 2: placebo

fit = survdiff(Surv(time, status) ~ group)
fit

plot(survfit(Surv(time, status) ~ group), 
     lty=c(1,2), xlab = "time", ylab = "surv prob")
legend(20, 1, c('group1', 'group2'), lty = c(1,2))

$i$번째 개체의 생존시간 $T_i$는 서로 독립이고, 중도절단시간 $C_i$도 서로 독립이라고 하자. 그러면 관측시간 $Y_i=min(T_i, C_i)$이고 각각의 관측값은 $(Y_i, \delta _i)$로 대응하게 된다. 여기서 $\delta _i=0$이면 중도절단된 경우이고 $\delta _i=1$이면 중도절단 되지 않은 경우다. 시점 t를 넘어 생존할 확률을 나타내는 생존함수(survival function) $S(t)$는

이고, 시점 t까지는 생존했다가 시점 t 바로 직후에 사망하게 되는 순간 위험률을 나타내는 위험함수 $h(t)$는

로 정의된다.

 T가 확률밀도함수 $f(x)$와 분포함수 $F(x)$를 갖는다고 하자. $S(t)$는 $S(0)=1, \lim _{t→∞}S(t)=0$이므로 생존함수 $S(t)$는

따라서 양변을 미분하면 $F(t)$는

또한, 생존함수 $f(x)$ 다음과 같이도 표현할 수 있다.

한편, 위험함수 $h(t)$는

$F(0)=1$을 이용하여 t에 대하여 적분하면 누적위험함수(culmulative harzard function) $H(t)$는

이므로 생존함수 $S(t)$는 누적위험함수 $H(t)$를 이용해 표현할 수 있다.

또한 위 식을 t에 대해 미분하면 생존함수 $f(t)$는 다음과 같이 표현할 수 있다.