콘텐츠로 이동
Data Prep
상세

Survival Analysis (생존 분석)

개요

항목 내용
분류 Statistics / Biostatistics / Machine Learning
핵심 개념 사건 발생까지의 시간(time-to-event) 분석
기원 의학 연구 (생존 기간 분석), 보험 수리학
적용 분야 의료 (환자 생존), 고객 이탈, 장비 고장, 금융 (부도 예측)
핵심 난제 중도절단(Censoring) -- 관측 기간 내 사건 미발생

수학적 기초

핵심 함수 정의

생존 함수 (Survival Function):

\[S(t) = P(T > t) = 1 - F(t)\]
  • \(T\): 사건 발생까지의 시간 (확률 변수)
  • \(F(t) = P(T \leq t)\): 누적 분포 함수
  • \(S(0) = 1\), \(\lim_{t \to \infty} S(t) = 0\)

위험 함수 (Hazard Function):

\[h(t) = \lim_{\Delta t \to 0} \frac{P(t \leq T < t + \Delta t \mid T \geq t)}{\Delta t} = \frac{f(t)}{S(t)}\]
  • 시각 \(t\)까지 생존한 개체가 그 직후 사건을 경험할 순간 확률
  • \(f(t)\): 확률 밀도 함수

누적 위험 함수 (Cumulative Hazard Function):

\[H(t) = \int_0^t h(u) \, du = -\ln S(t)\]

함수 간 관계

f(t) ---> S(t) = 1 - integral(f) ---> h(t) = f(t)/S(t)
  ^                                        |
  |                                        v
  +---------- f(t) = h(t) * S(t) <--- H(t) = integral(h)
                                            |
                                            v
                                    S(t) = exp(-H(t))
\[S(t) = \exp\left(-\int_0^t h(u) \, du\right) = \exp(-H(t))\]

주요 모수적 분포

분포 위험 함수 h(t) 생존 함수 S(t) 특징
Exponential \(\lambda\) (상수) \(e^{-\lambda t}\) 메모리리스, 일정한 위험
Weibull \(\lambda p (\lambda t)^{p-1}\) \(e^{-(\lambda t)^p}\) p>1 증가, p<1 감소, p=1 지수
Log-normal 복잡 \(1 - \Phi\left(\frac{\ln t - \mu}{\sigma}\right)\) 비단조 위험
Log-logistic \(\frac{(\lambda p)(\lambda t)^{p-1}}{1+(\lambda t)^p}\) \(\frac{1}{1+(\lambda t)^p}\) 비단조, 닫힌 형태
Gompertz \(\lambda e^{\alpha t}\) \(e^{-\frac{\lambda}{\alpha}(e^{\alpha t}-1)}\) 노화 모델링

중도절단 (Censoring)

중도절단 유형

시간 --->

대상 1: |====X              (사건 관측)
대상 2: |========|-->        (우측 절단: 연구 종료)
대상 3:     |====|-->        (우측 절단: 추적 소실)
대상 4:          |====X      (사건 관측)
대상 5: <--|=====X           (좌측 절단: 이전 사건 미관측)
대상 6:   |-----|            (구간 절단: 검진 사이에 발생)
          t1    t2

X = 사건 발생,  --> = 절단 (사건 미관측)
유형 정의 예시
우측 절단 (Right) 관측 종료 시점까지 사건 미발생 연구 종료 시 생존 환자
좌측 절단 (Left) 관측 시작 이전에 사건 발생 가능 HIV 검진 전 감염
구간 절단 (Interval) 두 관측 시점 사이에 사건 발생 정기 검진 사이 질병 진단

우측 절단 데이터 표현

각 관측치 \(i\)에 대해: - \(t_i\): 관측된 시간 (사건 시간 또는 절단 시간) - \(\delta_i\): 사건 지시자 (\(\delta_i = 1\): 사건 발생, \(\delta_i = 0\): 절단)

우도 함수 (Likelihood):

\[L = \prod_{i=1}^{n} \left[ f(t_i) \right]^{\delta_i} \left[ S(t_i) \right]^{1-\delta_i} = \prod_{i=1}^{n} \left[ h(t_i) \right]^{\delta_i} S(t_i)\]

비모수적 방법

Kaplan-Meier 추정량

\[\hat{S}(t) = \prod_{t_i \leq t} \left(1 - \frac{d_i}{n_i}\right)\]
  • \(t_i\): 사건이 발생한 시점 (순서대로)
  • \(d_i\): 시점 \(t_i\)에서의 사건 수
  • \(n_i\): 시점 \(t_i\) 직전의 위험 집합(risk set) 크기

Greenwood 공식 (분산 추정):

\[\widehat{\text{Var}}\left[\hat{S}(t)\right] = \hat{S}(t)^2 \sum_{t_i \leq t} \frac{d_i}{n_i(n_i - d_i)}\]

Nelson-Aalen 추정량

\[\hat{H}(t) = \sum_{t_i \leq t} \frac{d_i}{n_i}\]
  • 누적 위험 함수의 비모수 추정
  • Kaplan-Meier보다 작은 표본에서 안정적

Log-Rank 검정

두 그룹의 생존 곡선 비교:

\[\chi^2 = \frac{\left(\sum_{j} (d_{1j} - e_{1j})\right)^2}{\sum_{j} v_j}\]
  • \(e_{1j} = n_{1j} \cdot d_j / n_j\): 그룹 1의 기대 사건 수
  • \(v_j\): 분산 항
  • 귀무가설: 두 그룹의 생존 함수가 동일

반모수적 방법: Cox 비례 위험 모델

모델 정의

\[h(t \mid \mathbf{x}) = h_0(t) \cdot \exp(\boldsymbol{\beta}^T \mathbf{x})\]
  • \(h_0(t)\): 기저 위험 함수 (비모수적, 형태 가정 없음)
  • \(\exp(\boldsymbol{\beta}^T \mathbf{x})\): 공변량에 의한 위험 비율
  • 비례 위험 가정 (Proportional Hazards): 위험 비율이 시간에 무관

Partial Likelihood

\[L(\boldsymbol{\beta}) = \prod_{i: \delta_i=1} \frac{\exp(\boldsymbol{\beta}^T \mathbf{x}_i)}{\sum_{j \in R(t_i)} \exp(\boldsymbol{\beta}^T \mathbf{x}_j)}\]
  • \(R(t_i)\): 시점 \(t_i\)의 위험 집합 (risk set)
  • 기저 위험 \(h_0(t)\)가 소거됨 (핵심 장점)

위험 비율 (Hazard Ratio) 해석

\[\text{HR} = \frac{h(t \mid x_k + 1)}{h(t \mid x_k)} = \exp(\beta_k)\]
HR 값 해석
HR = 1 효과 없음
HR > 1 위험 증가 (예: HR=2 -> 위험 2배)
HR < 1 보호 효과 (예: HR=0.5 -> 위험 절반)

비례 위험 가정 검정

  1. Schoenfeld 잔차 검정: 잔차와 시간의 상관관계 검정
  2. Log(-log(S(t))) 그래프: 두 그룹의 곡선이 평행해야 함
  3. 시간 의존 공변량 추가: \(\beta(t) = \beta + \gamma \cdot g(t)\) 에서 \(\gamma = 0\) 검정

모수적 방법: Accelerated Failure Time (AFT)

모델 정의

\[\ln T = \boldsymbol{\alpha}^T \mathbf{x} + \sigma \epsilon\]
  • 사건 시간의 로그를 선형 모델로 설명
  • \(\epsilon\): 오차 분포 (표준 극값, 정규, 로지스틱)
  • Cox 모델과 달리 사건 시간에 직접 효과를 모델링

가속 인자 (Acceleration Factor):

\[\frac{S_1(t)}{S_0(t)} = S_0\left(\frac{t}{\exp(\boldsymbol{\alpha}^T \mathbf{x})}\right)\]
  • \(\exp(\alpha_k) > 1\): 사건 시간 연장 (보호 효과)
  • \(\exp(\alpha_k) < 1\): 사건 시간 단축 (위험 증가)

머신러닝 확장

Random Survival Forest (RSF)

항목 내용
논문 Ishwaran et al. (2008), "Random Survival Forests"
기반 Random Forest + 생존 분석
분할 기준 Log-rank 통계량 최대화

알고리즘:

RSF 알고리즘
===========================================
1. B개의 부트스트랩 샘플 생성
2. 각 트리에서:
   a. 무작위로 p개 중 m개 변수 선택
   b. Log-rank 기준으로 최적 분할
   c. 터미널 노드에서 Nelson-Aalen 추정
3. 앙상블:
   S_hat(t|x) = (1/B) * sum(S_b(t|x))
===========================================

DeepSurv

항목 내용
논문 Katzman et al. (2018), "DeepSurv: Personalized Treatment Recommender System Using A Cox Proportional Hazards Deep Neural Network"
핵심 Cox PH의 선형 항을 신경망으로 대체
\[h(t \mid \mathbf{x}) = h_0(t) \cdot \exp\left(\text{NN}(\mathbf{x})\right)\]

손실 함수 (Negative Log Partial Likelihood):

\[\mathcal{L} = -\sum_{i:\delta_i=1} \left[ \text{NN}(\mathbf{x}_i) - \log \sum_{j \in R(t_i)} \exp\left(\text{NN}(\mathbf{x}_j)\right) \right]\]

Cox-nnet

  • Cox 모델의 단일 은닉층 신경망 확장
  • Ridge penalty 적용
  • DeepSurv보다 단순하지만 효과적

DeepHit

항목 내용
논문 Lee et al. (2018), "DeepHit: A Deep Learning Approach to Survival Analysis with Competing Risks"
특징 비례 위험 가정 불필요, 경쟁 위험 직접 모델링

손실 함수:

\[\mathcal{L} = \mathcal{L}_1 + \alpha \cdot \mathcal{L}_2\]
  • \(\mathcal{L}_1\): Log-likelihood (사건 발생 확률)
  • \(\mathcal{L}_2\): Ranking loss (Concordance 최적화)

경쟁 위험 모델 (Competing Risks)

정의

여러 종류의 사건이 경쟁하는 상황: - 환자: 질병 사망 vs 다른 원인 사망 - 고객: 자발적 이탈 vs 계약 만료 vs 경쟁사 전환

원인별 위험 함수

\[h_k(t) = \lim_{\Delta t \to 0} \frac{P(t \leq T < t+\Delta t, \ D=k \mid T \geq t)}{\Delta t}\]
  • \(k\): 사건 유형 (1, 2, ..., K)
  • \(D\): 사건 유형 변수

누적 발생 함수 (Cumulative Incidence Function, CIF)

\[F_k(t) = P(T \leq t, D = k) = \int_0^t h_k(u) \cdot S(u) \, du\]

여기서 \(S(u) = \exp\left(-\sum_{k=1}^{K} H_k(u)\right)\)

주의: 경쟁 위험에서 \(1 - \text{KM}\)은 CIF와 다름 (편향 발생)

Fine-Gray 모델

\[h_k^*(t \mid \mathbf{x}) = h_{k,0}^*(t) \cdot \exp(\boldsymbol{\beta}_k^T \mathbf{x})\]
  • 부분분포 위험(subdistribution hazard)에 대한 Cox 모형
  • CIF에 직접 공변량 효과를 모델링

시간 의존 공변량

확장 Cox 모델

\[h(t \mid \mathbf{x}(t)) = h_0(t) \cdot \exp\left(\boldsymbol{\beta}_1^T \mathbf{x}_{\text{fixed}} + \boldsymbol{\beta}_2^T \mathbf{x}(t)\right)\]
  • \(\mathbf{x}(t)\): 시간에 따라 변하는 공변량 (예: 혈압, 약물 투여량)
  • 데이터 구조: 시작-종료 형식 (counting process)

데이터 구조 예:

환자 start stop 사건 치료 혈압
1 0 30 0 A 120
1 30 60 0 B 135
1 60 85 1 B 142
2 0 50 0 A 118
2 50 120 1 A 128

평가 지표

Concordance Index (C-index)

\[C = P(\hat{T}_i > \hat{T}_j \mid T_i > T_j)\]
  • 비교 가능한 쌍(comparable pairs)에서 순서가 맞는 비율
  • C = 0.5: 무작위, C = 1.0: 완벽한 순서 예측
  • Harrell's C-index: 절단 고려

계산:

\[C = \frac{\sum_{i,j} \mathbf{1}[\hat{\eta}_i > \hat{\eta}_j] \cdot \mathbf{1}[T_i < T_j] \cdot \delta_i}{\sum_{i,j} \mathbf{1}[T_i < T_j] \cdot \delta_i}\]

Brier Score

\[\text{BS}(t) = \frac{1}{n} \sum_{i=1}^{n} \left[\hat{S}(t \mid \mathbf{x}_i) - \mathbf{1}(T_i > t)\right]^2 \cdot W_i(t)\]
  • \(W_i(t)\): IPCW (Inverse Probability of Censoring Weighting) 가중치
  • 시점 \(t\)에서의 예측 정확도
  • 0에 가까울수록 좋음

Integrated Brier Score (IBS)

\[\text{IBS} = \frac{1}{t_{\max} - t_{\min}} \int_{t_{\min}}^{t_{\max}} \text{BS}(t) \, dt\]
  • 전체 시간 구간에 걸친 평균 Brier Score
  • 모델 간 종합 비교에 사용

지표 비교

지표 측정 대상 범위 절단 고려 시점 의존
C-index 판별력 (순위) 0.5~1.0 O (IPCW) X (전체)
Brier Score 교정(calibration) + 판별력 0~1 O (IPCW) O (특정 t)
IBS 종합 예측 성능 0~1 O X (통합)
AUC(t) 시점별 판별력 0.5~1.0 O O

실무 적용 사례

고객 이탈 분석 (Customer Churn)

요소 매핑
시간 가입 후 경과 일수
사건 이탈 (서비스 해지)
절단 관측 시점에 아직 활성 고객
공변량 이용 빈도, 결제 금액, 고객 등급
경쟁 위험 자발적 이탈 vs 계약 만료 vs 비자발적 해지

장비 고장 예측 (Predictive Maintenance)

요소 매핑
시간 가동 시간 / 사이클 수
사건 장비 고장
절단 정기 교체, 관측 종료
공변량 센서값, 운영 환경, 유지보수 이력
시간 의존 공변량 실시간 센서 데이터

의료 생존 분석

요소 매핑
시간 진단 후 생존 기간
사건 사망 / 재발
절단 추적 소실, 연구 종료
공변량 나이, 성별, 병기, 치료법
경쟁 위험 질병 사망 vs 기타 원인 사망

실무 코드 예시

lifelines를 이용한 Kaplan-Meier & Cox PH

import pandas as pd
import numpy as np
from lifelines import KaplanMeierFitter, CoxPHFitter
from lifelines.statistics import logrank_test

# -------------------------------------------------------
# 데이터 준비
# -------------------------------------------------------
np.random.seed(42)
n = 500
data = pd.DataFrame({
    "duration": np.random.weibull(1.5, n) * 100,
    "event": np.random.binomial(1, 0.7, n),
    "age": np.random.normal(50, 12, n),
    "treatment": np.random.choice(["A", "B"], n),
    "biomarker": np.random.exponential(1, n),
})
data["treatment_B"] = (data["treatment"] == "B").astype(int)

# -------------------------------------------------------
# Kaplan-Meier 추정
# -------------------------------------------------------
kmf = KaplanMeierFitter()

# 그룹별 추정
for name, group in data.groupby("treatment"):
    kmf.fit(
        durations=group["duration"],
        event_observed=group["event"],
        label=f"Treatment {name}",
    )
    print(f"Treatment {name}:")
    print(f"  Median survival: {kmf.median_survival_time_:.1f}")
    print(f"  S(50): {kmf.predict(50):.3f}")

# -------------------------------------------------------
# Log-Rank 검정
# -------------------------------------------------------
group_a = data[data["treatment"] == "A"]
group_b = data[data["treatment"] == "B"]

result = logrank_test(
    group_a["duration"], group_b["duration"],
    event_observed_A=group_a["event"],
    event_observed_B=group_b["event"],
)
print(f"\nLog-Rank Test p-value: {result.p_value:.4f}")

# -------------------------------------------------------
# Cox Proportional Hazards
# -------------------------------------------------------
cph = CoxPHFitter()
cph.fit(
    data[["duration", "event", "age", "treatment_B", "biomarker"]],
    duration_col="duration",
    event_col="event",
)

cph.print_summary()
# 출력 예:
#              coef  exp(coef)  se(coef)      z      p
# age          0.02    1.02      0.005    4.12  <0.001
# treatment_B -0.35    0.70      0.12    -2.92   0.004
# biomarker    0.18    1.20      0.06     3.01   0.003

# 비례 위험 가정 검정
cph.check_assumptions(data[["duration", "event", "age", "treatment_B", "biomarker"]])

scikit-survival을 이용한 Random Survival Forest

import numpy as np
import pandas as pd
from sksurv.datasets import load_veterans_lung_cancer
from sksurv.ensemble import RandomSurvivalForest
from sksurv.preprocessing import OneHotEncoder
from sksurv.metrics import concordance_index_censored, integrated_brier_score
from sklearn.model_selection import train_test_split

# -------------------------------------------------------
# 데이터 로드
# -------------------------------------------------------
data_x, data_y = load_veterans_lung_cancer()
# data_y: structured array with (event, time) fields

# 범주형 인코딩
encoder = OneHotEncoder()
data_x_enc = encoder.fit_transform(data_x)

X_train, X_test, y_train, y_test = train_test_split(
    data_x_enc, data_y, test_size=0.25, random_state=42
)

# -------------------------------------------------------
# Random Survival Forest
# -------------------------------------------------------
rsf = RandomSurvivalForest(
    n_estimators=100,
    min_samples_split=10,
    min_samples_leaf=5,
    max_features="sqrt",
    n_jobs=-1,
    random_state=42,
)
rsf.fit(X_train, y_train)

# -------------------------------------------------------
# 평가
# -------------------------------------------------------
# C-index
pred_risk = rsf.predict(X_test)
c_index = concordance_index_censored(
    y_test["Status"], y_test["Survival_in_days"], pred_risk
)
print(f"C-index: {c_index[0]:.4f}")

# 생존 함수 예측
surv_funcs = rsf.predict_survival_function(X_test)
print(f"Patient 0, S(100): {surv_funcs[0](100):.3f}")
print(f"Patient 0, S(200): {surv_funcs[0](200):.3f}")

# Integrated Brier Score
times = np.linspace(10, 500, 50)
surv_prob = np.row_stack([fn(times) for fn in surv_funcs])
ibs = integrated_brier_score(y_train, y_test, surv_prob, times)
print(f"IBS: {ibs:.4f}")

# 변수 중요도
importances = pd.Series(
    rsf.feature_importances_, index=data_x_enc.columns
).sort_values(ascending=False)
print("\nFeature Importances:")
print(importances.head(5))

PyCox를 이용한 DeepSurv / DeepHit

import numpy as np
import pandas as pd
import torch
import torchtuples as tt
from pycox.models import CoxPH, DeepHitSingle
from pycox.evaluation import EvalSurv
from sklearn.preprocessing import StandardScaler
from sklearn_pandas import DataFrameMapper

# -------------------------------------------------------
# 데이터 준비 (SUPPORT 데이터셋)
# -------------------------------------------------------
from pycox.datasets import support

df = support.read_df()
# columns: duration, event, x0, x1, ..., x13

# 훈련/검증/테스트 분할
df_train = df.sample(frac=0.7, random_state=42)
df_remain = df.drop(df_train.index)
df_val = df_remain.sample(frac=0.5, random_state=42)
df_test = df_remain.drop(df_val.index)

# 전처리
cols_numeric = [f"x{i}" for i in range(14)]
standardize = [(col, StandardScaler()) for col in cols_numeric]
mapper = DataFrameMapper(standardize, df_out=True)

x_train = mapper.fit_transform(df_train).values.astype("float32")
x_val = mapper.transform(df_val).values.astype("float32")
x_test = mapper.transform(df_test).values.astype("float32")

get_target = lambda df: (df["duration"].values, df["event"].values)
y_train = get_target(df_train)
y_val = get_target(df_val)
durations_test, events_test = get_target(df_test)

# -------------------------------------------------------
# DeepSurv (CoxPH with Neural Network)
# -------------------------------------------------------
in_features = x_train.shape[1]
num_nodes = [64, 64, 64]
out_features = 1
batch_norm = True
dropout = 0.2

net = tt.practical.MLPVanilla(
    in_features, num_nodes, out_features,
    batch_norm, dropout
)

model = CoxPH(net, optimizer=tt.optim.Adam)
model.optimizer.set_lr(1e-3)

# 학습
epochs = 100
callbacks = [tt.callbacks.EarlyStopping()]
log = model.fit(
    x_train, y_train, batch_size=256,
    epochs=epochs, callbacks=callbacks,
    val_data=(x_val, y_val),
)
print(f"Best epoch: {log.to_pandas().val_loss.argmin()}")

# Breslow 추정으로 기저 위험 계산
_ = model.compute_baseline_hazards()

# 예측 및 평가
surv = model.predict_surv_df(x_test)

ev = EvalSurv(surv, durations_test, events_test, censor_surv="km")
c_index = ev.concordance_td()
ibs = ev.integrated_brier_score(np.linspace(0, 365, 100))
print(f"C-index (td): {c_index:.4f}")
print(f"IBS: {ibs:.4f}")

# -------------------------------------------------------
# DeepHit
# -------------------------------------------------------
# 시간 이산화
num_durations = 50
labtrans = DeepHitSingle.label_transform(num_durations)
y_train_dh = labtrans.fit_transform(*y_train)
y_val_dh = labtrans.transform(*y_val)

net_dh = tt.practical.MLPVanilla(
    in_features, [64, 64, 64], labtrans.out_features,
    batch_norm=True, dropout=0.2,
)

model_dh = DeepHitSingle(net_dh, optimizer=tt.optim.Adam, alpha=0.2, sigma=0.1)
model_dh.optimizer.set_lr(1e-3)

log_dh = model_dh.fit(
    x_train, y_train_dh, batch_size=256,
    epochs=100, callbacks=[tt.callbacks.EarlyStopping()],
    val_data=(x_val, y_val_dh),
)

surv_dh = model_dh.predict_surv_df(x_test)
ev_dh = EvalSurv(surv_dh, durations_test, events_test, censor_surv="km")
print(f"DeepHit C-index: {ev_dh.concordance_td():.4f}")

lifelines를 이용한 경쟁 위험 분석

from lifelines import AalenJohansenFitter

# 경쟁 위험 데이터
data_cr = pd.DataFrame({
    "duration": [5, 10, 15, 20, 7, 12, 3, 25, 8, 18],
    "event_type": [1, 2, 0, 1, 2, 1, 0, 2, 1, 0],
    # 0: 절단, 1: 사건 A, 2: 사건 B
})

aj = AalenJohansenFitter()

# 사건 유형 1에 대한 CIF
aj.fit(
    durations=data_cr["duration"],
    event_observed=data_cr["event_type"],
    event_of_interest=1,
)
print("CIF for event 1:")
print(aj.cumulative_density_.head(10))

모델 비교표

모델 유형 PH 가정 비선형성 경쟁 위험 해석 가능성 확장성
Kaplan-Meier 비모수 -- -- X 높음 낮음
Cox PH 반모수 필요 X X 높음 중간
AFT 모수 X X X 높음 중간
Fine-Gray 반모수 필요 X O 높음 중간
RSF ML X O O (확장) 중간 높음
DeepSurv DL 필요 O X 낮음 높음
DeepHit DL X O O 낮음 높음
Cox-nnet DL 필요 O X 낮음 높음

한계와 발전 방향

현재 한계

  1. 중도절단 가정: 대부분 독립적 중도절단(independent censoring) 가정 -- 실제로는 정보적 중도절단이 흔함
  2. 고차원 데이터: 유전체/영상 데이터 등 초고차원 공변량 처리의 어려움
  3. 시간 의존 효과: 비례 위험 가정 위반 시 모델링 복잡
  4. 인과 추론: 관측 데이터에서 치료 효과 추정의 편향

발전 방향

  1. 동적 예측 (Dynamic Prediction): 새로운 정보 유입 시 실시간 예측 갱신
  2. 인과적 생존 분석: Causal Survival Analysis (IPTW, G-computation)
  3. 딥러닝 확장: Transformer 기반 생존 모델 (SurvTRACE, TransSurv)
  4. 연합학습: 의료 기관 간 데이터 공유 없이 협력 학습 (Federated Survival Analysis)
  5. 해석 가능한 DL: SHAP for survival models, SurvSHAP(t)