[생존분석] 두 그룹에 대한 로그-순위 검정(log-rank test)

백혈병 데이터

	재발까지 걸린 시간 (단위: 주)
그룹1	6, 6, 6, 6+, 7, 9+, 10, 10+, 11+, 13 16, 17+, 19+, 20+, 22, 23, 25+, 25+, 32+, 32+, 34+
그룹2	1, 1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 8, 8, 11, 11, 12, 12, 15, 17, 22, 23

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

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

 

 

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

	failed	unfailed	합
그룹1	1	16	17
그룹2	0	12	12
합	1	28	29

 예를 들어, 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))