생존 분석 - 카플란 마이어 추정 [파이썬 예제]
카플란 마이어 하드코딩
이번 예제에서는 강아지 보호소를 예시로 사용한다. 강아지 보호소에는 새로운 강아지가 입소하고, 입양되어 퇴소하는 사건이 발생한다. 따라서 강아지가 보호소에 남아있는 기간을 생존기간으로, 입양되어 퇴소하는 사건을 위험으로 간주한다.
데이터 셋 구성
가상의 강아지 보호소 데이터 구성은 아래와 같다.
- start : 보호소 입소 시간
- end : 입양 사건 발생 시간. 9는 관측 종료 시점이므로 해당 컬럼이 9일 경우 cencored
- status : 1일 경우 사건 발생, 0일 경우 사건 미발생
1
2
3
4
5
6
7
8
9
10
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
obs = pd.DataFrame()
obs['start'] = 0,1,2,2,4,6,7
obs['end'] = 5,2,6,9,9,8,9
obs['status'] = 1,1,1,0,0,1,0
obs
| start | end | status | |
|---|---|---|---|
| 0 | 0 | 5 | 1 |
| 1 | 1 | 2 | 1 |
| 2 | 2 | 6 | 1 |
| 3 | 2 | 9 | 0 |
| 4 | 4 | 9 | 0 |
| 5 | 6 | 8 | 1 |
| 6 | 7 | 9 | 0 |
A. 레코드 plot
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def plot_lifelines(obs):
for y,row in obs.iterrows():
start = row['start']
end = row['end']
status = row['status']
if status == 0:
plt.hlines(y,start,end,color='C0')
else :
plt.hlines(y,start,end,color='C1')
plt.plot(end,y,marker='o',color='C1')
plt.xlabel('Time (weeks)')
plt.ylabel('Dog Index')
plt.gca().invert_yaxis()
plot_lifelines(obs)
카플란 마이어 추정
카플란 마이어 추정은 두가지 아이디어에서 시작한다.
- 레코드의 시작지점을 무시하고 duration time만을 고려한다.
- 위험함수를 추정한다. 위험함수는 각 기간 내 발생한 사건을 위험에 처한 총 레코드 수로 나누어 도출한다.
A. Duration time 을 기준으로 데이터세트 생성
1
2
3
4
5
6
7
duration = obs['end'] - obs['start']
shifted = obs.copy()
shifted['start'] = 0
shifted['end'] = duration
plot_lifelines(shifted)
plt.xlabel('Duration (week)');
B. 위험함수 도출을 위한 At Risk 도출
At Risk는 특정 기간동안 사건이 발생하지 않았거나 절단되지 않은 레코드의 수를 뜻한다. 아래 코드는 numpy의 meshgrid 함수를 이용하여 사건이 발생한 Duration Time과 각 레코드의 end 시점을 각각 비교한다. 사건이 발생한 $Duration Time_i$ 보다 레코드의 $end_i$ 시점이 같거나 크다면, $DurationTime_i$ 내에 생존한 레코드이다. 즉, 사건이 발생한 각각의 Duration Time 내 위험에 노출된 레코드들의 수를 계산한다.
1
2
3
ts = duration.unique()
ts.sort()
ts
1
⇒ array([1, 2, 4, 5, 7])
1
2
E,T = np.meshgrid(shifted['end'],ts)
E
1
2
3
4
5
array([[5, 1, 4, 7, 5, 2, 2],
[5, 1, 4, 7, 5, 2, 2],
[5, 1, 4, 7, 5, 2, 2],
[5, 1, 4, 7, 5, 2, 2],
[5, 1, 4, 7, 5, 2, 2]])
1
2
at_risk = (T <= E).sum(axis=1)
at_risk
1
array([7, 6, 4, 3, 1])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
T<br
[<br>
[1, 1, 1, 1, 1, 1, 1],<br>
[2, 2, 2, 2, 2, 2, 2],<br>
[4, 4, 4, 4, 4, 4, 4],<br>
[5, 5, 5, 5, 5, 5, 5],<br>
[7, 7, 7, 7, 7, 7, 7]<br>
]<br>
<=
E<br>
[<br>
[5, 1, 4, 7, 5, 2, 2],<br>
[5, 1, 4, 7, 5, 2, 2],<br>
[5, 1, 4, 7, 5, 2, 2],<br>
[5, 1, 4, 7, 5, 2, 2],<br>
[5, 1, 4, 7, 5, 2, 2]<br>
]<br>
This results in a boolean array:
[<br>
[True, True, True, True, True, True, True],<br>
[True, False, True, True, True, True, True],<br>
[True, False, True, True, True, False, False],<br>
[True, False, False, True, True, False, False],<br>
[False, False, False, True, False, False, False]<br>
]<br>
[7, 6, 4, 3, 1]
This gives us the number of individuals who are considered "at risk" for each time point.
C. 위험함수 도출을 위한 사건 발생 수 도출
각 t(Duration Time)마다 입양이 발생한 수를 계산한다.
1
2
3
4
5
6
7
adopted = pd.Series(0,index=ts)
for t in ts:
k = (shifted['status']==1) & (t == shifted['end'])
adopted[t] = k.sum()
adopted
1
2
3
4
5
6
1 1
2 1
4 1
5 1
7 0
dtype: int64
1
2
3
4
d = dict(adopted=adopted,
at_risk=at_risk)
df= pd.DataFrame(d,index=ts)
df
| adopted | at_risk | |
|---|---|---|
| 1 | 1 | 7 |
| 2 | 1 | 6 |
| 4 | 1 | 4 |
| 5 | 1 | 3 |
| 7 | 0 | 1 |
D. 위험함수 계산. 사건발생 수 / 위험에 노출된 레코드 수
1
2
df['hazard'] = df['adopted']/df['at_risk']
df
| adopted | at_risk | hazard | |
|---|---|---|---|
| 1 | 1 | 7 | 0.142857 |
| 2 | 1 | 6 | 0.166667 |
| 4 | 1 | 4 | 0.250000 |
| 5 | 1 | 3 | 0.333333 |
| 7 | 0 | 1 | 0.000000 |
생존함수 계산
위에서 계산한 위험함수는 특정 기간에 사건이 발생할 확률이다. 이 케이스에서 사건은 입양이다. 따라서 반대의 케이스는 입양이 되지 않은 생존이며, 이 생존의 누적곱이 생존함수가 된다.
1
2
df['surv'] = (1 - df['hazard']).cumprod()
df
| adopted | at_risk | hazard | surv | |
|---|---|---|---|---|
| 1 | 1 | 7 | 0.142857 | 0.857143 |
| 2 | 1 | 6 | 0.166667 | 0.714286 |
| 4 | 1 | 4 | 0.250000 | 0.535714 |
| 5 | 1 | 3 | 0.333333 | 0.357143 |
| 7 | 0 | 1 | 0.000000 | 0.357143 |
누적 분포함수는 1 - 생존함수이다.
1
2
df['cdf'] = 1 - df['surv']
df
| adopted | at_risk | hazard | surv | cdf | |
|---|---|---|---|---|---|
| 1 | 1 | 7 | 0.142857 | 0.857143 | 0.142857 |
| 2 | 1 | 6 | 0.166667 | 0.714286 | 0.285714 |
| 4 | 1 | 4 | 0.250000 | 0.535714 | 0.464286 |
| 5 | 1 | 3 | 0.333333 | 0.357143 | 0.642857 |
| 7 | 0 | 1 | 0.000000 | 0.357143 | 0.642857 |
확률 질량함수는 cdf의 변화량이다
1
2
df['pdf'] = np.diff(df['cdf'],prepend=0)
df
| adopted | at_risk | hazard | surv | cdf | ||
|---|---|---|---|---|---|---|
| 1 | 1 | 7 | 0.142857 | 0.857143 | 0.142857 | 0.142857 |
| 2 | 1 | 6 | 0.166667 | 0.714286 | 0.285714 | 0.142857 |
| 4 | 1 | 4 | 0.250000 | 0.535714 | 0.464286 | 0.178571 |
| 5 | 1 | 3 | 0.333333 | 0.357143 | 0.642857 | 0.178571 |
| 7 | 0 | 1 | 0.000000 | 0.357143 | 0.642857 | 0.000000 |
라이브러리를 이용한 카플란마이어 추정
위의 하드코딩으로 구현한 카플란 마이어 추정은 lifelines 라이브러리를 활용하면 간단히 구현할 수 있다.
1
2
3
4
5
6
7
from lifelines import KaplanMeierFitter
kmf = KaplanMeierFitter()
T = shifted['end']
E = shifted['status']
kmf.fit(T,E)
1
<lifelines.KaplanMeierFitter:"KM_estimate", fitted with 7 total observations, 3 right-censored observations>
1
kmf.survival_function_
| KM_estimate | |
|---|---|
| timeline | |
| 0.0 | 1.000000 |
| 1.0 | 0.857143 |
| 2.0 | 0.714286 |
| 4.0 | 0.535714 |
| 5.0 | 0.357143 |
| 7.0 | 0.357143 |
생존함수의 confidence interval을 구할 수 있다. 현재 예제의 자료가 매우 작기에 CI는 상당히 넓다.
1
2
ci = kmf.confidence_interval_survival_function_
ci
| KM_estimate_lower_0.95 | KM_estimate_upper_0.95 | |
|---|---|---|
| 0.0 | 1.000000 | 1.000000 |
| 1.0 | 0.334054 | 0.978561 |
| 2.0 | 0.258154 | 0.919797 |
| 4.0 | 0.131988 | 0.824997 |
| 5.0 | 0.051977 | 0.698713 |
| 7.0 | 0.051977 | 0.698713 |
1
2
3
4
5
ts = ci.index
low, high = np.transpose(ci.values)
plt.fill_between(ts,low,high,color='grey',alpha=0.3)
kmf.survival_function_.plot(ax=plt.gca())
plt.ylabel('Survival function')
1
kmf.median_survival_time_
5.0
카플란 마이어 곡선 해석
생존 곡선을 통해 다양한 의사결정을 내릴 수 있다.
- 해당 강아지 보호소의 중위 생존값은 5. 따라서 입소 후 5주째 50%의 입양이 발생한다.
- 공간 활용 최적화를 위해 어느 주기로 신규 입소견을 받아야 할지 결정할 수 있다.
- 또한, 입양 사건 발생시 매출이 발생한다면, 생존기간과 매출을 곱하여 LTV와 예상 매출을 계산할 수 있다.
This post is licensed under CC BY 4.0 by the author.