딥러닝에서 클래스 불균형을 다루는 방법

현실 데이터에는 클래스 불균형 (class imbalance) 문제가 자주 있다. 어떤 데이터에서 각 클래스 (주로 범주형 반응 변수) 가 갖고 있는 데이터의 양에 차이가 큰 경우, 클래스 불균형이 있다고 말한다. 예를 들어, 병원에서 질병이 있는 사람과 질병이 없는 사람의 데이터를 수집했다고 하자. 일반적으로 질병이 있는 사람이 질병이 없는 사람에 비해 적다. 비단 병원 데이터뿐 아니라 대부분의 "현실 데이터" 에 클래스 불균형 문제가 있다. 

왜 클래스 균형이 필요한가?

왜 데이터가 클래스 균형을 이루어야할까? 그리고 언제 클래스 균형이 필요할까? 핵심은 다음과 같다. 클래스 균형 클래스 균형은 소수의 클래스에 특별히 더 큰 관심이 있는 경우에 필요하다. 

예를 들어 현재 재정 상황 및 집의 특성 등을 토대로 집을 사야할지 말아야할지를 예측하는 모델을 만들고 싶다고 하자. 사지말라고 예측하는 것과 사라고 예측하는 것은 그 무게가 다르다. 집을 사라고 예측하는 것은 훨씬 더 큰 리스크를 수반한다. 잘못된 투자는 큰 손실로 이루어질 수 있기 때문이다. 따라서 '집을 사라' 라고 예측하는 것에 대해서는 더 큰 정확도를 가져야한다. 하지만 데이터가 '집을 사지마라' 클래스에 몰려있는 경우, '집을 사지마라' 예측에 있어서는 높은 정확도를 가질 수 있어도 '집을 사라' 라고 예측하는 것에 관해서는 예측 성능이 좋지 않게 된다. 따라서 클래스 불균형이 있는 경우, 클래스에 따라 정확도가 달라지게 된다. 이를 해결하기 위해서는 따라서 '집을 사라' 클래스에는 더욱 큰 비중 (weight) 를 두고 정확한 예측을 할 수 있도록 만들어야한다.   

만약 소수 클래스에 관심이 없다면 어떻게 할까? 예를 들어, 이미지 분류 문제를 예로 들어보자. 그리고 오직, 전체 예측의 정확도 (accuracy) 에만 관심이 있다고 하자. 이 경우에는 굳이 클래스 균형을 맞출 필요가 없다. 왜냐하면 트레이닝 데이터에 만에 데이터를 위주로 학습하면, 모델의 정확도가 높아질 것이기 때문이다. 따라서 이런 경우에는 소수 클래스를 무시하더라도 전체 성능에 큰 영향을 주지 않기 때문에, 클래스 균형을 맞추는 것이 굳이 필요하지 않다고 할 수 있다. 

클래스 균형이 필요한 상황과 불필요한 상황을 예로 들어 설명했다. 다음으로는 딥러닝에서 클래스 균형을 맞추기 위한 두 가지 테크닉을 소개한다. 

(1) Weight balancing

Weight balancing 은 training set 의 각 데이터에서 loss 를 계산할 때 특정 클래스의 데이터에 더 큰 loss 값을 갖도록 하는 방법이다. 예를 들어, 이전 예에서 집을 사라는 클래스에 관해서는 더 큰 정확도가 필요하므로, 트레이닝 할 때, 집을 사라는 클래스의 데이터에 관해서는 loss 가 더 크도록 만드는 것이다. 이를 구현하는 한 가지 간단한 방법은 원하는 클래스의 데이터의 loss 에는 특정 값을 곱하고, 이를 통해 딥러닝 모델을 트레이닝하는 것이다.  

예를 들어, "집을 사라" 클래스에는 75 %의 가중치를 두고, "집을 사지마라" 클래스에는 25 %의 가중치를 둘 수 있다. 이를 python keras 를 통해 구현하면 아래와 같다. class_weight 라는 dictionary 를 만들고, keas model 의 class_weight parameter 로 넣어주면 된다. 

import keras

class_weight = {"buy": 0.75,
                "don't buy": 0.25}

model.fit(X_train, Y_train, epochs=10, batch_size=32, class_weight=class_weight)

물론 이 값을 예를 든 값이며, 분야와 최종 성능을 고려해 가중치 비율의 최적 세팅을 찾으면 된다. 다른 한 가지 방법은 클래스의 비율에 따라 가중치를 두는 방법인데, 예를 들어, 클래스의 비율이 1:9 라면 가중치를 9:1로 줌으로써 적은 샘플 수를 가진 클래스를 전체 loss 에 동일하게 기여하도록 할 수 있다. 

Weight balancing 에 사용할 수 있는 다른 방법은 Focal loss 를 사용하는 것이다. Focal loss 의 메인 아이디어는 다음과 같다. 다중 클래스 분류 문제에서, A, B, C 3개의 클래스가 존재한다고 하자. A 클래스는 상대적으로 분류하기 쉽고, B, C 클래스는 쉽다고 하자. 총 100번의 epoch 에서 단지 10번의 epoch 만에 validation set 에 대해 99 % 의 정확도를 얻었다. 그럼에도 불구하고 나머지 90 epoch 에 대해 A 클래스는 계속 loss 의 계산에 기여한다. 만약 상대적으로 분류하기 쉬운 A 클래스의 데이터 대신, B, C 클래스의 데이터에 더욱 집중을 해서 loss 를 계산을 하면 전체적인 정확도를 더 높일 수 있지 않을까? 예를 들어 batch size 가 64 라고 하면, 64 개의 sample 을 본 후, loss 를 계산해서 backpropagation 을 통해 weight 를 업데이트 하게 되는데 이 때, 이 loss 의 계산에 현재까지의 클래스 별 정확도를 고려한 weight 를 줌으로서 전반적인 모델의 정확도를 높이고자 하는 것이다. 

Focal loss 는 어떤 batch 의 트레이닝 데이터에 같은 weight 를 주지 않고, 분류 성능이 높은 클래스에 대해서는 down-weighting 을 한다. 이 때, gamma (위 그림) 를 주어, 이  down-weighting 의 정도를 결정한다. 이 방법은 분류가 힘든 데이터에 대한 트레닝을 강조하는 효과가 있다. Focal loss 는 Keras 에서 아래와 같은 custom loss function 을 정의하고 loss parameter 에 넣어줌으로써 구현할 수 있다. 

import keras
from keras import backend as K
import tensorflow as tf

# Define our custom loss function
def focal_loss(y_true, y_pred):
    gamma = 2.0, alpha = 0.25
    pt_1 = tf.where(tf.equal(y_true, 1), y_pred, tf.ones_like(y_pred))
    pt_0 = tf.where(tf.equal(y_true, 0), y_pred, tf.zeros_like(y_pred))
    return -K.sum(alpha * K.pow(1. - pt_1, gamma) * K.log(pt_1))-K.sum((1-alpha) * K.pow( pt_0, gamma) * K.log(1. - pt_0))

# Compile our model
adam = Adam(lr=0.0001)
model.compile(loss=[focal_loss], metrics=["accuracy"], optimizer=adam) 

(2) Over and under sampling

클래스 불균형을 다루기 위한 다른 방법은 바로 샘플링을 이용하는 것이다.

Under and and Over Sampling

예를 들어, 위 그림에서 파란색 데이터가 주황색 데이터에비해 양이 현저히 적다. 이 경우 두 가지 방법 - Undersampling, Oversampling 으로 샘플링을 할 수 있다. 

Undersampling 은 Majority class (파란색 데이터) 의 일부만을 선택하고, Minority class (주황색 데이터) 는 최대한 많은 데이터를 사용하는 방법이다. 이 때 Undersampling 된 파란색 데이터가 원본 데이터와 비교해 대표성이 있어야한다. Oversampling 은 Minority class 의 복사본을 만들어, Majority class 의 수만큼 데이터를 만들어주는 것이다. 똑같은 데이터를 그대로 복사하는 것이기 때문에 새로운 데이터는 기존 데이터와 같은 성질을 갖게된다.  

Reference

https://towardsdatascience.com/handling-imbalanced-datasets-in-deep-learning-f48407a0e758

https://towardsdatascience.com/methods-for-dealing-with-imbalanced-data-5b761be45a18

Adaptive Histogram Equalization 이란 무엇인가? 

Adaptive Histogram Equalization (AHE) 이란 이미지 전처리 기법으로 이미지의 contrast 를 늘리는 방법이다. Histogram Equalization (HE) 방법을 조금 더 개선해서 만든 방법이라고 할 수 있는데, 그 차이점은 HE 는 하나의 histogram 을 통해 uniform 분포를 가정한 equalization 을 수행하는 반면, AHE 는 여러개의 histogram 을 만든다는 차이점이 있다. HE 의 문제점은 하나의 histogram 으로 pixel intensity 의 redistribution 을 수행하므로, pixel intensity 의 분포가 이미지의 전 범위에 걸쳐서 동일해야 적절히 수행될 수 있는 방법이다. 만약 이미지의 일부분이 다른 지역과 비교해 다른 분포를 갖고 있을 경우, HE 를 적용한 이미지는 왜곡된다. 

AHE 방법은 이미지를 grid 를 이용해서 여러 개로 쪼개고 이 쪼갠 이미지에 대해서 각각 HE 를 적용한다. 따라서 그 grid 안에서 이미지의 contrast 가 좋아지기 때문에 이미지의 local contrast 를 조정하기에 좋은 방법이라고 할 수 있다. 위 그림처럼 어떤 pixel 의 AHE 를 적용한 값은 각 box 안에서 HE 를 수행한 값이다. 경계에 위치한 pixel 에 대한 transformation 은 특별하게 다루어지는데, 그 반대편에 위치한 pixel 값들을 mirroring 함으로써 가상의 픽셀 값을 만든 후, AHE 를 적용하게 된다. 

AHE 의 단점은 어떤 grid 안에 있는 pixel intensity 들이 매우 좁은 지역에 몰려있는 경우에 발생한다. 이 경우 AHE 는 매우 좁은 지역에 위치해 있는 pixel intensity 들을 넓은 지역으로 퍼트리기 때문에 이 안에 다른 지역과 매우 작은 차이를 보이는 noise 가 있다고 하더라도 상당한 크기의 peak 를 만들어내게 된다. 예를 들어, 8x8 grid 안에서 하나의 pixel 이 1이고 나머지는 0인 경우, 1이 255로 변환이 되게 된다. 즉, 의미 없는 pixel 의 값이 크게 변화되는 문제가 생긴다. 이를 noise amplification 이라고 부르기도 한다. 또한 의미없는 지역의 pixel intensity 를 enhance 시키기도 한다. 즉, 이로 인해 원래 이미지에 왜곡이 생기게된다. 

Contrast Limited Adaptive Histogram Equalization 

Adaptive Histogram Equalization 과 늘 따라다니는 것이 Contrast limited adaptive histogram equalization (CLAHE) 방법이다. 이 방법은 AHE 의 변형으로 AHE 의 중대한 문제점인 noise amplification 을 해결하기 위해 contrast limit 을 활용한다. 만약 pixel intensity 가 좁은 구역에 몰려있다면, 아래와 같은 두 가지 현상이 생긴다. 

  1. 그 부근에서 slope of cdf 가 급격하게 증가한다.

  2. histogram의 높이가 매우 높다.

1, 2 는 결국 같다. slope of cdf 가 급격하게 증가하는 부분이 histogram 의 높이가 높은 지점이다. 따라서 CLAHE 방법은 위 그림과 같이 cdf 를 계산하기 전에 histogram 에 높이에 제한을 걸어 특정 높이 이상에 있는 pixel 값들을 redistribution 하는 과정을 먼저 거치게 된다. 이 때 histogram 의 높이를 제한하는 값을 clip limit 이라고 부른다. redistribution 이후에, 다시 clip limit 의 값을 초과할 수 있는데, 이를 그대로 놔두고 할 수 있지만, 제거하고 싶은 경우, 반복적으로 redistribution 을 실시한 후, 초과하는 pixel 값들이 무시할만 한 정도까지 반복할 수 있다. 결국 이러한 과정은 cdf function 의 기울기가 너무 높지 않게 제한한다. 최종적으로 가상의 cdf 함수를 구할 수 있고, 이를 통해 HE 를 적용하면 CLAHE 를 적용한 이미지를 얻을 수 있다. 

그러면 이 이미지는 좁은 범위에 pixel intensity 가 몰려 있는 경우 생기는 noise 에 robust 하도록 만든다. 왜냐하면 몰려 있는 픽셀의 값을 redistribution 을 하기 때문에 급격한 기울기가 없는 cdf를 통해 pixel intensity 를 변환하기 때문에 noise 의 pixel intensity 를 너무 큰 값으로 변환시키지 않도록 만든다. 그렇기 때문에 CLAHE 를 실시한 후의 local region 의 pixel intensity 가 0~255 사이에 위치해 있는 것이 아니라, 더 좁은 범위에 몰려 있게 만들 수 있다. 즉 CLAHE 를 사용했을 때 장점은, contrast 가 낮은 지역에 위치한 noise 에 robust 하도록 변환함으로써, 변환된 이미지가 실제 이미지와 비슷한 특성을 갖도록 만드는 것이다. 

참고

https://en.wikipedia.org/wiki/Adaptive_histogram_equalization

Histogram Equalization 이란 무엇인가? 

IDEA

Histogram equalization (HE) 은 이미지를 전처리하는 방법이다. 위 그림만 보고도 HE 의 기본적인 아이디어를 알 수 있는데, HE 는 이미지의 contrast 가 적을 때 매우 유용하게 사용할 수 있는 방법이다. "Before Histogram Equalization" 이미지는 contrast 가 매우 떨어진다는 볼 수 있다. 오른쪽 histogram 을 보면 이미지 픽셀이 0~255에 고르게 퍼져있는 게 아니라, 일부분에 몰려있다는 것을 확인할 수 있다. "After Histogram Equalization" 을 보면, 이미지의 contrast 가 증가한 것을 확인할 수 있다. 이렇게 이미지의 밝기를 조정하는 것이 HE 로, HE 이후에는 밝기가 0~255 사이에 고르게 분포함을 알 수 있다. 이는 Cumulative density function (cdf) 을 통해 볼 수 있는데, 0~255 사이에서 cdf 가 일직선 형태를 나타내게 된다. 

Formula

L = 266 (일반적으로 0~255 사이에서 pixel intensity 를 표현하기 때문에)    

Example

HE 를 실제로 하는 방법은 해당 pixel intensity 의 cdf function 값을 구한 후, 0~255 사이의 uniform 분포에 맞게 값을 변환을 시켜서 넣어주면 된다. 예를 들어, 총 픽셀 수가 64 인 이미지가 있다고 하자 이 이미지에서 가장 낮은 pixel intensity 는 52 이다. 이를 어떻게 변환하는지를 보면 cdf(v) = 1, cdf_min = 1 MxN=64, cdf_min = 1, L-1 = 255 이고 분자가 0이므로 값은 0이 된다. pixel intensity 가 55 일 때, cdf(v) = 4, cdf_min = 1, MxN = 64 이므로, h(v) = 12 이다. cdf_min 은 cdf(v) 중 가장 작은 값으로 이 값을 분자 분모에서 빼주는 이유는 h(v) 를 0부터 시작하도록 하기 위함이다. 

참고

https://en.wikipedia.org/wiki/Histogram_equalization

Linearly weighted kappa 란 무엇인가?

Formula

$$ k_w = (P_0 - P_e) / (1 - P_e) $$

사용하는 이유

두 명의 측정자 (rater) 의 일치도를 평가 (inter-rater agreement) 할 때 쓰이며, 측정치가 순서가 있는 정성적 지표 (ordinal categorical variable) 일 때 쓰인다.  그냥 kappa coefficient 대신에 linearly weighted kappa 를 사용하는 이유는, 이러한 순서형 척도 변수에서 근접하게 측정된 값을 그냥 버리기 아깝기 때문에 이를 고려해서 일치도를 평가하기 위함이다. 

방법

P_0 =  proportion weighted observed agreement

P_e =  proportion weighted chance agreement

위의 두 값을 구하면 된다. 이를 위해 우선 아래와 같은 weight matrix 를 생성한다.  weight matrix는 측정치가 5개의 범주로 나누어지는 경우 아래와 같이 구할 수 있다. 

두 명의 측정자의 측정치가 아래와 같다고 한다. 예를 들어, 아래 테이블은 측정자 A와 측정자 B가 동시에 1이라고 판단한 경우는 8건이고 이는 전체의 0.07 이며, 측정자 A가 1, 측정자 B 가 2라고 판단한 경우는 2건으로 전체의 0.02 임을 의미한다.  

이 때, $$ P_0 = 0.07 * 1 + 0.02 * 0.75 + 0.009 * 0.5 + 0.026 * 0.75 + 0.09 * 1 ... + 0.02 * 1 = 0.93 $$ 로 구할 수 있다.

또한, $$ P_e =1 *(0.094 * 0.094) + 0.75 * (0.094 * 0.1709) + 0.50 * (0.094 * 0.530) ... + 0.02 * (0.0171*0.0171) = 0.78 $$ 로 구할 수 있다. 

결국 $$ k_w = (P_0 - P_e) / (1 - P_e) = (0.93 - 0.78) / (1 - 0.78) = 0.68 $$

Reference

http://www.anestesiarianimazione.com/DWLDocuments/The%20Linearly%20Weighted%20Kappa.pdf

로지스틱 회귀분석의 최대우도추정법 직관적 이해

Logistic regression의 Random component 

$$ Y \sim B(n, \pi) $$ 

Logistic regression은 반응변수 y가 위와 같은 이항 분포를 따른다고 가정한다.  

Model 

$$ \pi = \frac{1}{1+e^{-{(X\beta+\alpha)}}} $$ 

이 때, X를 통해 모수 pi 를 위와 같이 추정하는 것이다. 또는 일반적으로 아래와같이 표현되기도 한다. 

$$ logit(\pi) = X\beta + \alpha $$ 

Logistic regression에서 하고자하는건 결국, Random component에 주어진 pi를 추정하는 것인데, 이것이 0-1 사이의 값을 갖기 때문에 logit 변환을 수행하는 것이다. logit link function은 0~1사이의 값을 음의 무한대에서 양의 무한대로 바꾸어준다. 

Random component가 이항 분포이므로 likelihood는 아래와 같이 계산된다.

이곳에 Model을 집어넣음으로써, beta와 alpha 를 maximum likelihood 추정법으로 계산할 수 있게 된다. 근데 여기서 n과 y는 n번 시도했을 때 성공횟수 y를 나타낸다. 근데, 실제 각각의 데이터를 보면 시행횟수 n=1이며, y=0 or 1이다. 따라서 위 likelihood 는 베르누이 likelihood 와 같게된다. 

결국, 위의 likelihood를 최대화하면 계수 추정값을 얻을 수 있는데, log likelihood를 최대화하는 방식으로 쉽게 계수를 추정할 수 있다. 로지스틱 회귀분석에서는 보통 iteratively reweighted least squares (IRLS) 라는 방법을 통해 계수를 추정한다. 

$$ \sum_{i=1}^{N} y_i log(\pi_i) + (1-y_i) log (1-\pi_i) $$

Logistic regression의 log likelihood이다. 이는 negative binary cross entropy 이다. 결국, log likelihood를 최대화 하는 것은 negative binary cross entropy를 최대화 하는 것과 같고 이는 binary cross entropy를 최소화화는 것과 같다. 

고급통계 - Net reclassification improvements

Net reclassification index란 예측 모형에서 새로운 모델이 과거 모델보다 얼마나 얼마나 더 나아졌는지를 측정하는 지표이다. AUC 값을 기준으로 이를 판단하는게 과거에는 일반적이었으나, 그 한계점이 지적되외서 NRI가 등장하게 되었다. Framingham Heart Study에서 HDL 관련 연구를 수행할 때, AUC의 한계를 볼 수 있었는데, HDL을 예측 모형에 추가했을 때, AUC값이 크게 증가하지 않는 것이었다. 하지만 AUC 측면에서가 아니라, 모형의 계수를 통해 HDL의 effect를 파악할 때, HDL은 심장질환에 중요한 예측 인자였다. 

AUC의 한계는 다음과 같이 제시되었다. 1) 임상적인 의미와 연결시키기 힘들다. 2) 작은 변화는 잡아내기가 힘들다.

하지만 NRI는 이전 모델과 비교해 새로운 모델에서 얼마나 많은 사람들이 제대로 재분류되었는가를 파악함으로써 모델의 개선을 측정한다. 

출처 - https://en.wikipedia.org/wiki/Net_reclassification_improvement

만약 환자 30명, 비환자 70명에 대해서 환자-비환자 여부를 판단하는 모델 Test1, Test2를 만들었다고 하자. Test 2가 Test1에 비해 얼마나 우월한지를 판단하고자 한다. 위 테이블에서 여기서 검은색 글씨는 Test 1, Test 2 둘 다 맞춘거고, 하얀색 글씨는 Test1, Test 2 둘 다 틀린 경우, 초록색 글씨는 Test 2에서만 맞추고, Test 1에서는 틀린 것, 빨간색 글씨는 Test 1에서 맞추고, Test 2 에서 틑린 것이다. 

이 때, NRI의 계산은 매우 간단한데, Test 2를 통해서 얻은 이득이 얼마나 더 많은지를 계산해주면 된다. 이것은 Test 2를 도입했을 때, 추가적으로 맞춘것-오히려 틀린것 을 해줌으로써 계산되며, Event, Non-event 두 그룹의 대상자에 대해서 계산한 다음에 더함으로써 최종 NRI를 계산하게 된다. 

즉, Event가 발생한 사람들에 대한 NRI는 아래와 같이 계산된다.

$$ NRI_e = (4-2)/30 = 0.067 $$ 

다음으로 Non-event에 대한 NRI는 아래와 같이 계산된다.

$$ NRI_{ne} = (8-4)/70 = 0.057 $$ 

따라서 NRI는 아래와 같이 계산된다.

$$ NRI = NRI_e + NRI_{ne}= 0.124 $$ 

Gradient Boosting Algorithm의 직관적인 이해

실패를 통해 성공을 발전시켜라. 낙담과 실패는 성공으로 가는 가장 확실한 두 개의 디딤돌이다. -데일 카네기

Gradient Boosting Algorithm (GBM)은 회귀분석 또는 분류 분석을 수행할 수 있는 예측모형이며 예측모형의 앙상블 방법론 중 부스팅 계열에 속하는 알고리즘입니다. Gradient Boosting Algorithm은 Tabular format 데이터 (엑셀형태와 같이 X-Y Grid로 되어있는 데이터)에 대한 예측에서 엄청난 성능을 보여주고, 머신러닝 알고리즘 중에서도 가장 예측 성능이 높다고 알려진 알고리즘입니다. 그렇기 때문에 Gradient Boosting Algorithm을 구현한 패키지들이 많습니다. LightGBM, CatBoost, XGBoost 같은 파이썬 패키지들이 모두 Gradient Boosting Algorithm을 구현한 패키지들입니다. GBM은 계산량이 상당히 많이 필요한 알고리즘이기 때문에, 이를 하드웨어 효율적으로 구현하는 것이 필요한데, 위 패키지들은 모두 GBM을 효율적으로 구현하려고한 패키지들이라고 볼 수 있습니다.

Boosting 이란?

주요 앙상블 알고리즘은 bagging과 boosting으로 나눌 수 있고, Gradient boosting은 boosting 계열의 앙상블 알고리즘입니다. Gradient Boosting 알고리즘은 Gradient 를 이용하여 Boosting하는 알고리즘입니다. Gradient boosting 외의 boosting 계열의 알고리즘의 예로는 Adaptive boosting을 들 수 있습니다. 

먼저 Boosting의 개념부터 살펴봅시다. Boosting이란 약한 분류기를 결합하여 강한 분류기를 만드는 과정입니다. 분류기 A, B, C 가 있고, 각각의 0.3 정도의 accuracy를 보여준다고 합시다. A, B, C를 결합하여 더 높은 정확도, 예를 들어 0.7 정도의 accuracy를 얻는 게 앙상블 알고리즘의 기본 원리입니다. Boosting은 이 과정을 순차적으로 실행합니다. A 분류기를 만든 후, 그 정보를 바탕으로 B 분류기를 만들고, 다시 그 정보를 바탕으로 C 분류기를 만듭니다. 그리고 최종적으로 만들어진 분류기들을 모두 결합하여 최종 모델을 만드는 것이 Boosting의 원리입니다. 

GBM의 직관적인 이해

GBM을 이해하는 가장 쉬운 방법은 Residual fitting으로 이해하는 것입니다. 아주 간단한 모델 A를 통해 y를 예측하고 남은 잔차 (residual)을 다시 B라는 모델을 통해 예측하고 A+B 모델을 통해 y를 예측한다면 A보다 나은 B 모델을 만들 수 있게 되죠. 이러한 방법을 계속하면 잔차는 계속해서 줄어들게되고, training set을 잘 설명하는 예측 모형을 만들 수 있게 됩니다. 하지만 이러한 방식은 bias는 상당히 줄일 수 있어도, 과적합이 일어날 수도 있다는 단점이 있습니다. 따라서 실제로 GBM을 사용할 때는 sampling, penalizing 등의 regularization 테크닉을 이용하여 더 advanced 된 모델을 이용하는 것이 보편적입니다. 하지만 본 포스팅의 목적은 GBM 에 대하여 직관적인 이해를 해보는 것이기 때문에 이 부분은 다음에 기회가 있을 때 정리해보도록 하겠습니다. 

위 그림을 보시면 tree 1을 통해 예측하고 남은 잔차를 tree2를 통해 예측하고, 이를 반복함으로서 점점 잔차를 줄여나가는 것을 볼 수 있습니다. 이 때, 각각의 모델 tree1,2,3 을약한 분류기 (weak learner), 이를 결합한 분류기를 강한 분류기 (strong learner)라고도 합니다. 보통 약한 분류기로는 간단한 의사결정나무 (decision tree)를 많이 사용합니다. 이를 Gradient boosting tree라고도 하는데, 구현한 대표적인 라이브러리로 XGboost를 들 수 있습니다. XGBoost는 python, R로 이용할 수 있습니다. 

Gradient란 무엇인가?

그렇다면 residual fitting과 gradient가 무슨 관계인지 궁금하실텐데요. residual은 loss function을 squared error로 설정하였을 때, negative gradient 입니다. 따라서 residual에 fitting해서 다음 모델을 순차적으로 만들어 나가는 것은 negative gradient를 이용해 다음 모델을 순차적으로 만들어 나가는 것으로 볼 수 있습니다. 우선 아래 식을 통해 이를 확인해봅시다. 

loss function을 아래와 같이 정의하면,

$$ j(y_i, f(x_i)) = \frac{1}{2}(y_i - f(x_i))^2 $$ 

이 때의 negative gradient는 residual이 됩니다.

residual이 $$ y_i - f(x_i) $$ 이기 때문에 negative gradient가 residual이 된다는 것을 알 수 있습니다. 

따라서 gradient boosting은 다음 모델을 만들 때, negative gradient를 이용해 만들기 때문에 gradient boosting 이라고 할 수 있습니다. 또한 residual fitting model은 gradient boosting 모델의 한 가지 종류입니다. 왜냐하면 다른 loss function을 이용해서 gradient boosting을 할 수 있기 때문입니다. 예를 들어 회귀 문제가 아닌 분류 문제라면 loss function으로 squared error 를 이용할 수 없기 때문에 새로운 loss function을 만들어야하는데, 이 경우에도 negative gradient를 이용해서 새로운 model 을 fitting하고 이를 합산해 나가시는 식으로 최종 모델을 만들게 됩니다. 

직관적으로 이것은 무엇을 뜻할까요? 왜 negative gradient를 이용해서 새로운 모델을 만드는 걸까요? negative gradient는 pseudo-residual이라고도 불리며, 이것은 어떤 데이터 포인트에서 loss function이 줄어들기 위해 f(x)가 가려고하는 방향입니다. 이 방향에 새로운 모델을 fitting해서 이것을 이전 모델과 결합하면, f(x) 는 loss function이 줄어드는 방향으로 업데이트가 되겠죠. 이것이 gradient boosting의 아이디어입니다. Gradient boosting을 gradient descent + boosting 이라고 하기도 합니다. 어쨌든, loss function을 줄이는 방향의 negative gradient를 얻고, 이를 활용해 boosting을 하는 것이기 때문에 gradient descent와 boosting이 결합된 방법이다. 라고 이해하셔도 괜찮습니다. 

위 그림은 위키피디아에서 gradient boosting에 대하여 설명한 알고리즘도입니다. 이제 위 내용을 이해할 수 있습니다. gradient boosting에서는 learning rate를 통하여 pseudo-residual에 fitting된 모델을 어느정도로 update할지에 대해서 정하게 되는데, 이 값은 위 알고리즘도의 3번처럼, loss를 최소화하는 값을 optimization 과정을 통해 얻을 수도 있겠지만, 일반적으로 0.001~0.01 정도로 매우 작은 값으로 설정하는 것이 보편적입니다. 왜냐하면 작은 값으로 설정하여야 굉장히 세밀한 분류기를 얻을 수 있습니다. learning rate가 높을 수록 빠르게 모델의 bias를 줄여나가지만, learning rate가 적으면 디테일한 부분을 놓칠 수 있다고 이해하시면 됩니다. 이와 관련해서는 gradient boosting 과정을 simulation 해볼 수 있는 사이트를 소개합니다. 이 곳에서 직접 gradient boosting 방법이 어떻게 함수를 fitting 해나가는지를 시각적으로 확인할 수 있습니다. 

참고

https://machinelearningmastery.com/gentle-introduction-gradient-boosting-algorithm-machine-learning/

생존분석 - 모수적 모형

Cox의 비례위험 모형

일반적으로 생존분석에 사용되는 모델인 Cox의 비례위험모형의 경우 준모수적 (Semi-parametric) 방법입니다. 이 방법은 생존함수 (Survival function)을 추정하기 위해 Hazard function에 대한 모델링을 합니다. 시간에 따라 변하는 Baseline hazard가 존재하고, Covariate 들의 effect에 따라 hazard가 달라진다고 모델링합니다. 그리고 hazard의 effect인 Hazard ratio가 시간에 따라 동일하다는 비례위험을 가정하여 Hazard function과 Survival function 을 추정하게 됩니다. 이 방법의 장점은 Survival function에 분포에 대해서 가정을하지 않기 때문에 강건하며, Covariate를 보정할 수 있다는 것입니다. 또한 Baseline hazard function이 어떤 함수인지를 알 필요가 없이 계수 추정을 통해 Hazard ratio만 구할 수 있다면 Covariate의 effect를 파악할 수 있습니다. 

모수적 모형

반면 모수적 모형은 Hazard function과 Survival function의 모양에 대해 Cox 보다 더 강한 '가정'을 합니다. 이 가정은 Hazard function과 Survival function이 어떤 분포를 따를 것이라는 믿음에서 출발합니다. 가장 많이 사용되는 모수적 모델은 Weibull 모형입니다. 

Weibull 모형의 Hazard function과 Survival function 은 다음과 같이 정의됩니다.

$$ h(t) = \lambda p t^{p-1} $$

$$ S(t) = exp(-\lambda t^p) $$ 

h(t)와 S(t)에 정의에 따라 h(t)가 저렇게 정의되면 S(t)가 저렇게 계산되는 것을 확인할 수 있습니다.

$$ \lambda = exp(\beta_0 + \beta_1 TRT) $$ 

만약 Treatment 변수에 대해 Weibull model을 만들 경우, 이 때 위와 같이 식을 놓고, beta0, beta1, p 3개의 계수를 추정하게 됩니다. 

또는 Hazard를 기준으로 계수를 추정하는 것이 아니라 Survival time을 기준으로 추정할 수 있는데, 이 경우가 모수적 방법에서 많이 사용되는 방법입니다. 이를 Acceleration Failure Time Model이라고 합니다. 

$$ S(t) = exp(-\lambda t^p) $$

$$ t = [-lnS(t)]^{1/p} \times \frac{1}{\lambda^{\frac{1}{p}}} $$

$$ \frac{1}{\lambda^{\frac{1}{p}}} = exp(\alpha_0 + \alpha_1 TRT) $$

AFT 모델에서는 위와 같이 놓고 alpha0, alpha1, p를 추정하게 됩니다. 

이렇게 두 가지 방법으로 구한 계수 beta와 alpha의 의미는 다른데, 만약 p가 고정되있다고 한다면, beta는 exponential을 취하면 Hazard ratio를 뜻하며, alpha는 exponential을 취하면 Accerleration factor가 됩니다. Acceleration factor에 대해서는 본 포스팅에서는 다루지 않을 예정이므로 이 자료를 참고 바랍니다. 

Strata 별로 p가 동일하다고 가정한 경우, Weibull model은 Proportional Hazard 가정과 Acceleration Failure Time 가정을 만족합니다. (https://en.wikipedia.org/wiki/Accelerated_failure_time_model)

PBC 데이터를 통한 Weibull model의 구축

pbc.dat

library(survival)

data <- read.table("pbc.dat")

# Column rename
colnames(data) <- c('id', 'futime', 'status', 'drug', 'age', 'sex', 'ascites', 'hepato', 'spiders', 'edema', 'bili', 'chol', 'albumin', 'copper', 'alk_phos', 'sgot', 'trig', 
                   'platelet', 'protime', 'stage')

head(data)

data <- data[data$drug != '.', ]
data$drug <- as.character(data$drug)
data$drug[(data$drug == 2)] = 0
data$drug <- factor(data$drug)

Drug = 0 플라시보, Drug = 1 페니실린에 대해 Death에 대한 Time to event 분석을 수행합니다.

Weibull Model Fitting

library(survival)
library(eha)

# AFT  
swei = survreg(s~drug, dist = "weibull",data = data) 

# PH  
sweiph = phreg(s~drug, dist = "weibull",data = data)

summary(swei) 
summary(sweiph)

Call:
survreg(formula = s ~ drug, data = data, dist = "weibull")
              Value Std. Error     z      p
(Intercept)  8.3167     0.1094 75.99 <2e-16
drug1       -0.0438     0.1431 -0.31  0.759
Log(scale)  -0.1532     0.0729 -2.10  0.036

Scale= 0.858 

Weibull distribution
Loglik(model)= -1348.2   Loglik(intercept only)= -1348.3
	Chisq= 0.09 on 1 degrees of freedom, p= 0.76 
Number of Newton-Raphson Iterations: 5 
n= 312 

Call:
phreg(formula = s ~ drug, data = data, dist = "weibull")

Covariate          W.mean      Coef Exp(Coef)  se(Coef)    Wald p
drug 
               0    0.491     0         1           (reference)
               1    0.509     0.051     1.052     0.167     0.759 

log(scale)                    8.317               0.109     0.000 
log(shape)                    0.153               0.073     0.036 

Events                    144 
Total time at risk        625985 
Max. log. likelihood      -1348.2 
LR test statistic         0.09 
Degrees of freedom        1 
Overall p-value           0.759387

Median Survival Time Comparison

predict(swei,newdata=data.frame(drug='0'),type=c("quantile"),se.fit=FALSE,p=0.5)
predict(swei,newdata=data.frame(drug='1'),type=c("quantile"),se.fit=FALSE,p=0.5)

Weibull Model을 통해 구한 Survival Function과 Kaplen-Meier estimation의 비교

data$event_all=ifelse(data$status==0,0,1)
data$futime = ifelse(data$futime==0,0.1,data$futime)
s = Surv(data$futime,data$event_all)

sKM <- survfit(s ~ drug, data=data)
plot(sKM)
lines(predict(swei,newdata=list(drug='0'),type="quantile", p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col='red')
lines(predict(swei,newdata=list(drug='1'),type="quantile", p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col='red')

위 그림에서 검은색 직선이 Kaplen-Meier Survival function estimation이며, 붉은색과 파란색이 Weibull 모델로 추정한 Survival function 입니다. Weibull 모델의 S(t)의 모양에 대한 가정이 맞다면 두 직선은 가까울 것입니다. 그리고 위 직선들은 p가 같다는 가정하에 구해진 직선들입니다. 즉, PH, AFT 가정을 한 후에 그려진 직선입니다. 그렇기 때문에 Weibull 모형의 Survival function에 대한 가정, 그리고, PH, AFT 가정을 해도 될지에 대해서 체크를 해보아야합니다.

Weibull 가정의 체크 

(https://stat.ethz.ch/education/semesters/ss2011/seminar/contents/handout_9.pdf)

Weibull 모형에서 가정을 체크하는 방법은 X 축을 log(t), Y 축을 log(-log(S(t)) 로 그래프를 그려보는 것입니다. 이 때, S(t)는 Kaplen-Meier 방법으로 구한 Survival function 입니다. 이 때, 그룹간의 그래프의 모양을 보고 Weibull 생존모형 가정이 맞는지, 또한 두 그룹에서 PH, AFT 가정이 맞는지를 체크할 수 있습니다. 

먼저, 두 그룹의 직선이 평행하게 나오면 Weibull 가정을 만족합니다. 그렇지 않으면 Weibull 가정이 만족되지 못하며 AFT 가정도 만족되지 못합니다. 또한 두 그룹의 기울기의 비가 일정하면 PH 가정을 만족합니다. 

sKM <- survfit(s~ drug, data=data)
kmfit.TRT1 = summary(sKM)
kmfit.TRT2 = data.frame(kmfit.TRT1$strata, kmfit.TRT1$time, kmfit.TRT1$surv)
names(kmfit.TRT2) = c("TRT","time","survival")
TRT0=kmfit.TRT2[kmfit.TRT2$TRT=="drug=0", ]
TRT1=kmfit.TRT2[kmfit.TRT2$TRT=="drug=1", ]
plot(log(TRT0$time), log(-log(TRT0$survival)), pch=16, xlab="time in days using log-arithmic scale", ylab="log-log survival",type="p", xlim=c(3,10), ylim=c(-6,1), font=2, lwd=2, font.lab=2, col="Black",main="log-log curves by treatment")
par(new=T)
plot(log(TRT1$time), log(-log(TRT1$survival)), pch=3, col="grey50", ylab="",xlab="",  xlim=c(3,10), ylim=c(-6,1), font=2, lwd=2, font.lab=2)

이 그래프의 경우 우선 그룹에서 거의 직선이라는 것을 확인할 수 있지만, 두 그래프가 평행하지는 않습니다. (교차가 일어납니다.) 그렇기 때문에 Weibull 가정을 만족하지만 PH 가정과 AFT 가정을 만족하지 않으며, 두 그룹에서 p가 다르다는 알 수 있습니다.

Time-dependent variable의 정의 

Definition : 생존 분석에 있어서 시간에 따라 변수의 값이 바뀌는 변수를 time-dependent variable이라고한다. 

예 : 흡연 여부, 직장에서의 지위 등

Extended cox model : Time-dependent variable을 생존분석에서 사용하는 방법 

개념 : 만약 시간에 따라 변수의 값이 변화하는 변수가 있다면 extended cox model을 수행하여야한다. 

- 기존 데이터의 형태

- Extended cox model을 위한 데이터의 형태

모든 time에 대해 데이터를 쪼개서 start와 time으로 정의한다. 하나의 row를 시간에 따라 쪼개므로 한 사람당 여러개의 row가 생기게 된다. 만약 해당 시간 내에서 event가 발생한 경우, event는 1이된다.

library(data.table)
library(survival)
data <- read.csv("addict.csv")

- Extended cox model을 위한 데이터의 형태 만드는 법

data.cp <- survSplit(data, cut=data$time[data$event==1], end="time", event="event", start="start", id="ID2")

만약 아래와 같이 logdose라는 time-dependent variable이 있다고 하면 

data.cp$logdose <- data.cp$dose * log(data.cp$time)

Extended cox model 은 아래와 같이 수행할 수 있다.

Y = Surv(data.cp$start, data.cp$time, data.cp$event)
coxph(Y ~ clinic + logdose +prison+ cluster(ID2), data=data.cp)

Time-dependent variable의 활용 1 : PH (proportional hazard) Assumption check

개념 : PH Assumption에 위반이 되는 것이 의심되는 변수가 있을 때, 이 변수에 대한 time-dependent variable을 정의함으로서 PH 가정을 체크할 수 있다. 

만약에 위와 같은 데이터에 대해서 dose에 대해 PH Assumption 가정이 위반된다고 해보자. 즉, 이말은 시간에 따라 dose의 hazard ratio가 달라진다는 것이다. 이를 테스트하기위해 dose*g(t)라는 새로운 변수를 정의하고 이 변수의 계수가 유의한지 확인하는 방법을 통해 해당 변수에 대한 PH Assumption을 체크할 수가 있다. 

Y = Surv(data.cp$start, data.cp$time, data.cp$event)
coxph(Y ~ clinic + dose + logdose +prison+ cluster(ID2), data=data.cp)

Call:
coxph(formula = Y ~ clinic + dose + logdose + prison + cluster(ID2), 
    data = data.cp)

            coef exp(coef) se(coef) robust se     z       p
clinic  -1.01988   0.36064  0.21542   0.23637 -4.31 1.6e-05
dose    -0.08262   0.92070  0.03598   0.02960 -2.79  0.0053
logdose  0.00862   1.00865  0.00645   0.00525  1.64  0.1007
prison   0.34063   1.40584  0.16747   0.15972  2.13  0.0329

Likelihood ratio test=66.35  on 4 df, p=1e-13
n= 18708, number of events= 150 
   (1379 observations deleted due to missingness)

위 결과에서 logdose의 p-value > 0.05이므로, dose가 PH Assumption을 만족한다고 결론내릴 수 있다. 위 방법은 Wald test를 한 것이며, 다른 방법으로는 Likelihood ratio test를 통해 p-value를 구할 수도 있다. 

Likelihood ratio test

Y = Surv(data.cp$start, data.cp$time, data.cp$event)
mod1 <- coxph(Y ~ clinic + dose + logdose +prison+ cluster(ID2), data=data.cp)
mod2 <- coxph(Y ~ clinic + dose +prison+ cluster(ID2), data=data.cp)

print(mod1$loglik[2])
print(mod2$loglik[2])

LRT <- (-2) * (mod2$loglik[2]-mod1$loglik[2])
pvalue <- 1 - pchisq(LRT, 2) 

이는 logdose를 포함한 full model과 logdose를 포함하지 않은 restricted model의 likelihood ratio를 구하고 이것이 카이제곱 분포를 따른다는 것을 이용한 검정이다. p-value를 계산한 결과 0.4 정도로 위의 wald test와 마찬가지로 logdose는 의 계수는 유의하지 않았다. 

따라서 dose에 대해 PH Assumption이 어긋나지 않는다고 결론 내릴 수 있다.

Time-dependent variable의 활용 2 : PH Assumption 만족 안하는 경우, 시간에 따라 달라지는 Hazard ratio 모델링

개념 : 어떤 변수에 대하여 heavy side function을 이용한 time-dependent variable을 정의함으로서 시간에 따라 달라지는 hazard ratio를 모델링할 수 있다. 

heavy side function은 어떤 값을 기준으로 그 이상일 때는 1, 미만일 때는 0을 돌려주는 함수를 말한다. 

예제 : 만약 clinic에 대해 시간에 따라 달라지는 HR을 모델링 하고 싶다고 해보자. 

fit <- survfit(Surv(data$time, data$event)~data$clinic)
plot(fit)

이것이 clinic=1과 clinic=2일 때 각각 그린 카플란 마이어 생존함수 추정값이다. 위 그림에서 hazard ratio가 시간에 따라 바뀐다는 것이 의심스러울 수 있다. 이 때, 어떤 시점을 기준으로 heavy side function을 이용한 time-dependent variable을 정의함으로써 이를 모델링할 수 있다. 

365를 기준으로 데이터를 쪼갠다. 

data.cp365 <- survSplit(data, cut=365, end="time", event="event", start="start", id="ID2")

data.cp365$hv1 <- data.cp365$clinic * (data.cp365$start < 365)
data.cp365$hv2 <- data.cp365$clinic * (data.cp365$start >= 365)

위와 같이 두 개의 heavy side function을 원래 clinic의 값에 곱해준 것을 두 개의 변수로 정의한다.  

Extended cox model

Y365 <- Surv(data.cp365$start, data.cp365$time, data.cp365$event)
coxph(Y365 ~ dose + hv2 + hv1 + prison + cluster(ID2), data=data.cp365)

Call:
coxph(formula = Y365 ~ dose + hv2 + hv1 + prison + cluster(ID2), 
    data = data.cp365)

           coef exp(coef) se(coef) robust se     z       p
dose   -0.03548   0.96514  0.00643   0.00652 -5.44 5.3e-08
hv2    -1.83052   0.16033  0.38595   0.39838 -4.59 4.3e-06
hv1    -0.45937   0.63168  0.25529   0.25998 -1.77   0.077
prison  0.37795   1.45929  0.16841   0.16765  2.25   0.024

Likelihood ratio test=74.25  on 4 df, p=3e-15
n= 360, number of events= 150 

이 때 결과를 보면 hv1의 HR은 0.63, hv2의 HR은 0.16이다. 이 말은 clinic 1과 clinic2의 hazard ratio가 365미만에서는 0.63이었는데, 365 이상에서는 0.16이라고 해석할 수 있다. 즉, 두 clinic에서 보이는 생존률에 대한 경향성이 시간이 지날 수록 강해지고 있는 것이다. 이는 위의 카플란 마이서 생존곡선에서도 확인할 수 있다. 

Cox 비례위험 모형의 가정 체크

개념 : Cox 비례위험 모형은 "비례위험 가정" 을한다. 

비례위험가정은 시간에 따라 어떤 변수에서의 Hazard ratio가 일정해야한다는 것이다.

비례위험가정이 어긋나면, Stratified cox regression 또는 time-dependant variable 정의를 통한 extended cox regression을 수행할 수 있다. 

비례위험 가정의 체크는 아래방식으로 할 수 있다. 

1. Log-log plot

2. Observed-expected plot

3. Goodness of fit test

본 포스팅에서는 R의 survival, survminer 패키지를 통한 실습을 통해 cox 비례 위험 모형의 가정을 체크하는 방법을 알아본다. 

addict.csv

library(data.table)
library(survival)
data <- read.csv("addict.csv")

head(data)

먼저 데이터의 형태는 위와 같이 생겼다. 

관심변수 : clinic 에 대하여 Kaplan-Meier 플롯 그리기

# Treatment 별로 그룹지어서 KM Estimation
fit = survfit(y~clinic, data=data)
par(mai=c(1,1,1,1))
plot(fit,main="KM Curve", xlab="Time(Week)", ylab=expression(paste(hat(S),"(t)")), col=c("blue", "red"))

우리는 이것이 비례위험 가정을 만족하는지를 확인하고 싶다. 비례위험 가정을 만족하려면 모든 timepoints에서 Hazard ratio가 같아야한다. hazard란 f(t)/S(t) (참고 : https://3months.tistory.com/349?category=743476) 로 위 그림만 봐서는 가정을 만족하는지 한 눈에 알 수 없다.  

로그로그 플롯

개념 : 로그로그 플롯을 그리고 그룹간에 평행하게 나오면 비례위험가정을 만족한다. 

plot(fit, fun="cloglog", lty=1:2, col=c("Black", "Grey50"), lwd=2, font=2, ylim=c(-3,2), font.lab=2, main="log-log KM curves by Rx", ylab="log-log survival", xlab="time in weeks using logarithmic scale")
legend("bottomright",lty=1:2,legend=c("Group1","Group2"),bty="n", text.font=2, lwd=2, col=c("Black", "Grey50"))

관찰값-기댓갓 플롯

개념 : Kaplen-Meier 생존 함수와 Cox regression을 통해 구한 생존함수를 동시에 그려본다.

이 때, Expected plot은 Cox 비례위험 모형에서 나오는 생존함수로 관찰한 생존함수와 모양의 차이가 많이나면 비례위험 가정을 만족하지 않는 것이다. 

kmfit <- survfit(y~clinic, data=data)
#observed
plot(kmfit, lty="dashed", col=c("Black", "Grey50"), lwd=2, font=2, font.lab=2, main="Observed Versus Expected Plots by Clinic", ylab="Survival probability", xlab="Survival time")
par(new=T)
#expected
exp <- coxph(y~clinic, data=data)
new_df <- data.frame(clinic=c(1,2))
kmfit.exp <- survfit(exp, newdata=new_df)
plot(kmfit.exp, lty="solid", col=c("Black", "Grey50"), lwd=2, font=2, font.lab=2)

Goodness of fit test

개념 : Goodness of fit test(GoF Test) 는 관찰값-기댓값의 차이를 통계적으로 검정하는 방법이다. (카이스퀘어 분포 이용)

# Gooodness of fit test
mod1 <- coxph(y~clinic, data=data)
cox.zph(mod1, transform=rank)

plot(cox.zph(mod1, transform=rank), se=F, var="clinic")
abline(h=0, lty=3)
plot(cox.zph(mod1, transform=rank), se=F, var="prison")
abline(h=0, lty=3)

        rho chisq       p
clinic -0.255  10.5 0.00118

p-value가 0.001이므로 clinic 변수는 비례위험 가정을 만족하지 않는다. 이는 그림을 통해서도 확인할 수 있다. 위 그림에서 직선이 평행하지를 확인하면되는데, 직선은 추정된 beta값을 의미하는데 직선이 평행하지 않고 어떠한 경향성이 보이면 beta 값이 시간에 따라 바뀐다는 의미로 비례위험 가정에 어긋난다. 

참고 : https://www.r-bloggers.com/cox-model-assumptions/