Concept of Survival Analysis

고객 이탈예측 모형 개발을 위해 생존분석(Survival Analysis)관련 자료를 정리하고자 포스팅을 작성한다. 생존분석은 어떠한 사건이 발생할 확률을 시간이라는 변수와 같이 고려하는 분석 방법으로 의료 통계에서 많이 사용되는데 이 방법을 고객 이탈을 예측하는 모델을 개발하는 과정에서도 사용한다.

생존분석을 공부하기 위해 평소 통계학 문헌에서 언급되지 않았던 몇가지 생소한 개념에 대한 사전 정의가 필요하고 이번 포스팅은 그러한 개념들을 정리하기 위한 용도로 작성할 것이다. 생존분석 기법을 고객 이탈예측 시 적용하기 위해 공부하는 것이므로 몇가지 단어 선택도 변경하여 정리하기로 한다.

시간(time)

생존분석은 고객이 서비스를 가입(또는 구매)한 시점으로부터 이탈(사건이 발생)할 때까지의 시간 구간(time interval) 데이터에 관심이 있다. 그래서 우리의 관심 대상이 될 반응 변수는 고객이 이탈하는 사건이 발생할 때까지 걸리는 시간이다.

생존분석 실생 시, 시간 경과에 따른 고객의 잔존 확률을 구하게 된다. 이 때 시간은 하나의 독립 변수로 사용된다. 대신 이 때 시간은 상대적 시간 개념으로 분석 대상이 관측되기 시작한 시점(즉, 서비스를 가입한 시점)을 0으로 하여 시간을 계산하게 된다.

사건(event)

잔존의 반대인 이탈을 의미한다. 고객이 이탈하는 것은 한번 발생하며 0과 1로 구분 가능하다.

생존함수(survival function)

생존함수는 관찰대상이 특정 기준 시간보다 더 늦게 이탈하거나 이탈 자체가 일어나지 않을 확률을 계산하는 함수다. 즉 사건이 발생하기까지의 시간에 대한 함수로 다음과 같이 함수 \(S(t)\)로 표현 한다.

\[S(t) = p(T > t) = \int_{t}^{infty}f(x)dx = 1- p(T \leq t) = 1- F(t)\]

생존함수는 확률을 계산하는 함수이므로 \(0 \leq S(t) \leq 1\) 조건을 만족한다. 그리고 이탈이 관측되는 시간 t는 관측 즉시 발생할 수도 있고 또는 영영 발생하지 않을 수도 있다. 따라서 이 함수에서 t의 범위는 \(t \in [0,\infty)\) 이며, \(S(0) = 1, S(\infty) = 0\) 인 것도 쉽게 받아들일 수 있다. 더불어 생존함수가 non-increasing function인 것도 받아들일 수 있다.

그래서 \(p(T>t)\) 는 고객의 이탈이 t시점 이후 발생할 확률이며 이는 적어도 t시점까지는 고객이 이탈하지 않고 서비스를 이용한다고 해석할 수 있다.

위험함수 (hazard function)

조건부 확률에 대한 함수 \(h(t)\) 이며, t시점에 잔존했다는 조건 하, t시점에 이탈이 발생할 확률을 계산한다. 즉 t시점 이전까지는 이탈이 발생하지 않았는데 특정한 시점 t에 이탈이 발생할 확률을 의미한다.

\[\begin{align} h(t) &= \text{lim}_{\Delta \to 0} \frac{ p(t \leq T \leq T + \Delta t\vert T \geq t) }{ \Delta t} = \frac{f(t)}{S(t)} = \frac{ \frac{d}{dt}F(t) }{S(t)} \nonumber \\ &= \frac{ \frac{d}{dt}(1-S(t))}{S(t)} = \frac{-S'(t) }{S(t)} = -d\log{S(t)} \end{align}\]

위험함수는 증가/감소/상수형태 모두 가능하다.

만약 서비스 가입 기간 후 오랜 시간이 지나 평균적인 서비스 이용 기간을 넘어섰다면 이탈 위험은 증가할 것이다.

반대로 서비스 가입 초기는 한번 사용해보고 즉시 이탈할 수 있는 가능성이 있다. 그러나 일정기간 후 서비스 이용에 대한 안정기에 접어들면 위험함수는 (일시적으로라도) 감소가 가능하다.

누적위험함수(cumulative hazard function)

시간 0에서부터 t시점까지의 위험함수를 적분한 값이며 이는 t시점까지 이탈이 발생할 확률을 모두 더한 것이다.

\[H(t) = \int_{0}^{t}h(u)du = \int_{0}^{t}\frac{ \frac{d}{du}(1-S(u)) }{ S(u) } = -\log{S(t)}\]

수식을 통해 알 수 있는 것처럼 \(S(t) = \text{exp}\left(- H(t)\right)\) 이다.

중도절단(censoring / censored)

중도절단은 생존 데이터의 가장 큰 특징 중 하나라고 볼 수 있다. 중도 절단 데이터는 한정된 분석 기간 안에서 사건이 발생하지 않고 종료될 경우에 발생할 수 있고 분석의 대상이 되는 고객이 추적 불가능한 경우에도 발생할 수 있다.

Survival

그림처럼 20주 전부터 관찰을 시작해 현재까지의 사용자 생존시간을 추적할 수 있다. 사용자가 서비스를 사용하기 시작하면, 그 때부터 사용자의 수명이 시작된다. 관측 기간 내에 사용자가 서비스를 이탈한 경우, x자로 표기한다. 반면에 분석종료 시점까지 이탈하지 않거나(사용자 C) 도중에 사용자 추적이 불가능한 경우(사용자 A)가 중도 절단 케이스로 간주된다.

생존 데이터에 대한 표본에 대해 중도 절단된 대상자 \(i\)에 대한 관측 중단 시점을 \(c_{i}\) 로 표기한다. 사건 발생시간은 \(y_{i}\) 로 표기하기로 하자. 만약 중도 절단 이후 사건이 발생한다면 \(c_{i} = y_{i}\) 라고 볼 수 있다.

따라서 생존 시간 \(t_{i}\) 는 다음과 같이 표현할 수 있을 것이다. \(\delta_{i}\)는 i번째 데이터가 사건발생이 관측된 것을 나타내는 identity function이다. 즉, \(\delta_{i} = 0\)인 경우 중도절단 된 것을 의미한다.

\[t_{i} = \text{min}(y_{i},c_{i}) \quad \delta_{i} = \mathcal{I}(y_{i} < c_{i})\]

Kaplan-Meier 추정량

생존함수와 위험함수는 데이터를 통해 추정해야할 대상이다. 생존시간에 대한 분포가정을 할 수 있다면 parametric한 추정이 가능하다. 그러나 실제 데이터에 대해서는 parametric 가정을 할 수 없는 경우가 많다. 이러한 경우는 non-parametric 한 추정 방법을 고려해야 하는데 가장 대표적인 생존함수 추정 방법이 Kaplan-Meier 추정이다.

그림과 같이 Kaplan Meier 추정량에 대한 Plot을 확인할 수 있다. 시간 경과에 따라 그룹별 누적생존확률이 보여진다. 관측 시작부터 점차 누적 생존확률이 떨어지는 것을 확인할 수 있다. 이 Plot을 바탕으로 그룹별 생존 확률이 유의미한 차이가 있는 것으로 해석된다. 따라서 그룹을 표기하는 이 변수는 각 집단별 생존확률을 구분하는데 유의미한 변수일 것이라고 추측할 수 있게 된다.

생존 데이터는 i번째 사용자에 대한 생존시간 \(y_{i}\)와 중도절단여부 \(\delta_{i}\)를 포함해 \(\{(y_{i}, \delta_{i}), i=1,2,...,n \}\) 을 표기한다. 중도절단 대상자의 관측중단 시점은 \(c_{i}\) 로 표현한다.

\(y_{i} < c_{i}\) 인 경우는 \(y_{i}=t_{i}\) 인데 그렇지 못한 경우는 \(c_{i} = t_{i}\) 이다.

\[t_{i} = \text{min}(y_{i},c_{i}) \quad \delta_{i} = \mathcal{I}(y_{i}< c_{i})\]

사건이 발생한 unique한 시간이 총 D개가 있다면 \(t_{1} \leq t_{2} \leq \cdot\cdot\cdot \leq t_{D}\) 로 표현할 수 있을 것이고 \(t_{i}\) 시점에 발생한 사건의 수를 d_{i}로 표현한다. 만약 동시간에 사건의 중복발생이 없다면, \(d_{i}=1\) 이다.

\(Y_{i}\) 는 \(t_{i}\) 시점에서 사건발생 가능성이 있는 개체수를 의미하며, \(\frac{d_{i}}{Y_{i}}\) 는 \(t_{i}\) 시점 이전까지는 잔존하나 \(t_{i}\) 시점에 이탈할 조건부 확률이다.

이렇다면, \(1-\frac{d_{i}}{Y_{i}}\) 는 \(t_{i}\) 시점 이후로도 잔존할 확률을 의미한다.

시점 \(t_{i}\)에서의 생존함수 \(S(t_{i})\)는 \(t_{i-1}\) 시점의 생존함수 \(S(t_{i-1})\) 과 \(t_{i}\) 시점에서의 생존율 \(p_{i} = 1-\frac{d_{i}}{Y_{i}}\) 을 곱한 값이 된다.

생존확률을 구해보면 다음과 같다. 그리고 이 값이 생존함수 \(S(t_{i})\) 의 추정량이 된다.

\[\begin{align} S(t_{i}) &= p(T>t_{i}) = S(t_{i-1})\times p_{i} \nonumber \\ &= S(t_{i-2})\times p_{i-1} \times p_{i} = p_{1} \times p_{2} \times ... \times p_{i} \nonumber \\ &= \prod_{j=1}^{t_{j} \leq t_{i}} (1-\frac{d_{j}}{Y_{j}}) \nonumber \end{align}\]