Survival Analysis

1. Cox Proportional Hazard Model

생존 분석은 이벤트(사망, 고장 등) 발생까지의 시간을 분석하는 통계 기법입니다. 이때의 '이벤트'는 여러가지가 될 수 있으며, 이를 '사망'이라고 표현하는 것은 통계학에서의 관례입니다. 가장 널리 사용되는 생존 분석 모델 중 하나는 Cox의 비례 위험 모델입니다. 이 모델은 각 변수가 '사망'에 얼마나 영향을 미치는지를 분석합니다.

생존분석에서 '이벤트(event)'는 관심 이벤트가 발생했는지 여부를 나타내는 이진 변수입니다. 즉, 이벤트가 발생하지 않았다면 0, 발생했다면 1의 값을 가집니다. 이에 따라, '연소' 이벤트를 생존분석에 적용하려면, 각 관측치가 연소를 했는지 여부를 나타내는 변수가 필요합니다. 이를 위해 '연소'가 발생한 시점에 해당하는 데이터에는 1을, 그렇지 않은 경우에는 0을 할당해야 합니다.

Cox의 비례 위험 모델은 '비례 위험 가정'이라는 중요한 가정에 의존합니다. 이 가정은 시간에 따라 위험 비율이 일정하다는 것을 의미하며, 다시 말해, 어떤 변수의 효과가 시간에 따라 변하지 않는다는 가정입니다.(공변량의 효과가 시간에 따라 변하지 않는다.) 이 가정이 위반될 경우, 모델의 추정 결과는 신뢰할 수 없게 됩니다.

from lifelines import CoxPHFitter

cph = CoxPHFitter().fit(surv_df, 'TTF', 'event')
cph.check_assumptions(surv_df) # 비례 위험 가정을 검사

d_cols = [col1, col2] # 위 검사 결과에 따른 feature 제외 
surv_df = surv_df.drop(columns=d_cols)

cph = CoxPHFitter().fit(surv_df, 'TTF', 'event')
cph.print_summary()

cph.print_summary()는 Cox의 비례 위험 모델의 결과를 보여줍니다. 다음과 같은 항목들을 확인할 수 있습니다:

• coef: 설명 변수의 회귀 계수입니다. 회귀 계수가 양수인 경우, 그 변수의 증가는 '사망' 위험 증가와 관련이 있다고 해석할 수 있으며, 음수인 경우 '사망' 위험 감소와 관련이 있다고 해석할 수 있습니다.

• exp(coef): 회귀 계수를 지수화한 값입니다. 이것은 '위험비'로 해석될 수 있습니다. 예를 들어, exp(coef)가 2인 경우, 해당 변수가 1 단위 증가할 때마다 위험(사망률)이 두 배가 된다는 것을 의미합니다.

• se(coef): 회귀 계수의 표준 오차입니다. 이 값이 작을수록 추정치의 정확도가 높아집니다.

• z: 회귀 계수를 표준 오차로 나눈 값으로, z-통계량이라고도 합니다. 이 값이 크면 해당 변수의 영향력이 크다는 것을 의미합니다.

• p: 귀무 가설을 기각할 수 있는 p-값입니다. 이 값이 작을수록 해당 변수가 '사망'에 미치는 영향이 통계적으로 유의미하다는 것을 의미합니다.

따라서 이를 통해 각 변수가 연소 발생에 얼마나 영향을 미치는지, 그 영향력이 통계적으로 유의미한지를 알 수 있습니다. 이는 제조 공정에서 연소 발생에 영향을 미치는 주요 요인을 파악하는 데 유용할 것입니다.

그러나 모든 해석은 데이터의 특성과 비즈니스 컨텍스트를 고려하여 이루어져야 하며, p-값 등의 통계적 추정치는 그 자체로 결론을 내리는 것이 아니라 결론을 지원하는 증거로 사용되어야 합니다.

summ = cph.summary
summ[summ['p'] < 0.05].sort_values('se(coef)')

p-value가 0.05 미만인 변수 중 다른 여러 수치들을 고려하여 가장 주요한 영향인자를 도출

2. ML Survival Model

• Gradient Boosting Survival Analysis

from sksurv.ensemble import GradientBoostingSurvivalAnalysis
from sksurv.metrics import brier_score
from sksurv.metrics import integrated_brier_score
from sklearn.preprocessing import MinMaxScaler

# 생존함수 적용을 위한 'event'와 'time' 생성
events = [0]*(df.shape[0]-20) + [1]*20
times = list(np.arange(df.shape[0]+1)
df['event'] = events
df['time'] = times

event에는 실제 event(고장, 사망)가 발생한 시점을 1로, 나머지는 0으로 적용한다.
원래 이게 맞지만, data에 과도하게 0이 많은 경우이면서 model fitting에 알 수 없는 에러가 발생한다면 실제 고장(event)난 시점으로부터 직전 10구간 또는 20 구간을 모두 1로 적용해보면 해결할 수 있다.

# Train-test-split
tr_x = train.drop(columns=['time', 'event'])
tr_y = train[['time', 'event']]
te_x = test.drop(columns=['time', 'event'])
te_y = test[['time', 'event']]

scaler = MinMaxScaler()
tr_sc = scaler.fit_transform(tr_x)
te_sc = scaler.transform(te_x)

tr_x = pd.DataFrame(data = tr_sc, columns=tr_x.columns, index = tr_x.index)
te_x = pd.DataFrame(data = te_sc, columns=te_x.columns, index = te_x.index)
# y값 형식 변경 (모델 입력 구조에 맞추기 위해)
tr_y_arr = np.array(list(zip(tr_y['event'], tr_y['time'])), dtype=[('event','?'), ('Survival_hours', '<f8')])
te_y_arr = np.array(list(zip(te_y['event'], te_y['time'])), dtype=[('event','?'), ('Survival_hours', '<f8')])

#gbs = GradientBoostingSurvivalAnalysis(random_state=42)
gbs = GradientBoostingSurvivalAnalysis(
	n_estimators=300,
    min_sample_split=10,
    min_sample_leaf=5,
    min_depth =5,
    max_features='sqrt',
    random_state=42)
gbs.fit(tr_x, tr_y_arr)
gbs.score(te_x, te_y_arr)  # C-index 

Concordance Index (C-index)는 생존 분석 모델의 성능을 평가하는 지표 중 하나로, 모델이 예측한 생존 시간이 실제 생존 시간과 얼마나 일치하는지를 측정한다. 값의 범위는 0부터 1까지이며, 값이 클수록 모델의 예측 성능이 좋다고 해석된다.

간단히 말해, C-index는 다음과 같은 방식으로 계산된다:

테스트 데이터 세트에 있는 모든 가능한 환자 쌍을 살펴본다.
두 환자 중 하나가 더 빨리 이벤트(예: 사망)를 경험했다면, 모델이 그 환자에게 더 높은 위험을 할당했는지 확인한다.
모든 환자 쌍에 대해 이를 반복하고, 모델이 올바르게 예측한 비율을 계산한다.
예를 들어, C-index가 0.8이라면, 모델이 80%의 환자 쌍에 대해 올바르게 누가 더 빨리 이벤트를 경험할지 예측했다는 것을 의미한다.

C-index가 0.5인 경우, 모델의 예측이 무작위와 같다는 것을 의미하며, C-index가 1인 경우는 모델이 완벽하게 모든 환자 쌍을 분류했다는 것을 의미한다.

# Evaluation
risk = gbs.predict(te_x) # 위험 점수 예측 
te_y['risk'] = risk  

plt.figure(figsize=(10,8))
sns.scatterplot(data=te_y, x=te_y.index, y='event', label='event');
sns.scatterplot(data=te_y, x=te_y.index, y='risk', label='risk');
plt.xlabel("Time(X_test Index)")
plt.title("Event - Predicted Risk Score")
plt.show()

sns.boxplot(data=te_y, y='risk', x='event')
plt.ylabel("Risk score")
plt.show()

from sksurv.metrics import integrated_brier_score
pred_survs = gbs.predict_survival_function(te_x)  # 생존 함수 예측 
times = np.percentile(te_y['time'], [50, 75, 90]) # 특정 분위수에 해당하는 시간
preds = np.asarray([[fn(t) for t in times] for in pred_survs]) # 특정 시간에 대한 생존 확률 
integrated_brier_score(tr_y_arr, te_y_arr, preds, times)

Brier Score는 예측된 확률과 실제 발생 여부 사이의 차이를 측정하는 지표이다. 생존 분석에서는 특정 시간에 대한 생존 확률을 예측하고, 이 확률이 실제로 얼마나 정확한지 평가하는 데 사용된다.
Brier Score가 0에 가까울수록 예측이 정확하다는 것을 의미하며, 1에 가까울수록 예측이 부정확하다는 것을 의미한다.
예를 들어, 모델이 어떤 환자가 1년 안에 생존할 확률을 0.7로 예측했고, 실제로 그 환자가 1년 동안 생존했다면, Brier Score는 $(0.7-1)^2=0.09$ 가 된다.
Integrated brier score는 여러 시간 포인트에서의 brier score를 통합하여 종합적으로 판단하는 지표로 활용된다.