비례위험모형

3. 콕스 모델

콕스 모델은 비례위험의 가정이 성립하고 생존함수와 변수 사이에 로그-선형관계가 있다는 가정 하에, 회귀분석 방법을 이용하여 중도 절단된 자료를 처리한다.

피험자 $i$에 대한 공변량 벡터를 $X\_i\; =\; (X\_\{i1\},\; \backslash dots,\; X\_\{ip\})$라고 할 때, 콕스 비례 위험 모형의 위험 함수 $\backslash lambda(t)$는 다음과 같은 형태를 가진다.

$$$$

\begin{align}

\lambda(t|X_i) &= \lambda_0(t)\exp(\beta_1X_{i1} + \cdots + \beta_pX_{ip}) \\

&= \lambda_0(t)\exp(X_i \cdot \beta)

\end{align}



이 식은 공변량 벡터 $X\_i$를 가진 피험자 $i$의 시간 $t$에서의 위험 함수를 나타낸다. 여기서 $\backslash lambda\_0(t)$는 기준선 위험(baseline hazard)으로, 모든 피험자에게 동일하게 적용되는 시간 $t$에 따른 기본적인 위험 수준을 의미한다. 피험자 간의 위험률 차이는 오직 공변량의 효과를 나타내는 $\backslash exp(X\_i\; \backslash cdot\; \backslash beta)$ 항에 의해서만 결정된다.

데이비드 콕스 경은 비례 위험 가정이 성립한다면(또는 성립한다고 가정하면), 전체 위험 함수($\backslash lambda\_0(t)$ 포함)를 구체적으로 알지 못하더라도 효과 매개변수 $\backslash beta$를 추정할 수 있다는 점을 발견했다. 이러한 생존 데이터 분석 접근 방식을 콕스 비례 위험 모형,^[2] 또는 줄여서 콕스 모형, 비례 위험 모형이라고 부른다.^[3] 그러나 콕스는 비례 위험 가정을 생물학적으로 해석하는 것이 매우 까다로울 수 있다고 지적하기도 했다.^[4][5]

콕스 모델은 부분 우도(partial likelihood)라는 개념을 사용하여 매개변수 $\backslash beta$를 추정한다. 이는 시간 경과에 따른 위험의 변화를 직접 모델링하지 않고도 공변량의 영향을 추정할 수 있게 해준다. 이 모델을 통해 추정된 공변량의 효과는 종종 위험비(hazard ratio)로 해석된다.

3. 1. 비례성

3. 2. 절편 항의 부재

3. 3. 유일한 시간을 위한 가능도

3. 4. 동률 시간이 존재하는 경우의 가능도

4. 예시

다음은 콕스 비례위험모형의 실제 사용에 대한 몇 가지 예시이다.

4. 1. 단일 이진 공변량

단일 공변량을 사용하는 콕스 비례위험모형은 다음과 같이 표현할 수 있다. 여기서 $\backslash beta\_1$은 병원의 영향을 나타내는 계수이고, ''i''는 각 환자를 의미한다.

$$

\overbrace{\lambda(t|X_{i})}^{\text{환자 i의 위험률}} = \underbrace{\lambda_0(t)}_{\text{기준선 위험률}}\cdot\overbrace{\exp(\beta_1 X_{i})}^{\text{환자 i의 위험 스케일링 팩터}}



통계 소프트웨어를 사용하여 $\backslash beta\_1$ 값을 추정한 결과 약 2.12를 얻었다고 가정하자. 위험비(Hazard Ratio, HR)는 이 값의 지수인 $\backslash exp(\backslash beta\_1)\; =\; \backslash exp(2.12)$로 계산된다. 두 병원 간의 위험률 비율을 비교하면 다음과 같다.

$$

\frac{\lambda(t|X=1)}{\lambda(t|X=0)} = \frac{\cancel{\lambda_0(t)}\exp(\beta_1 \cdot 1)}{\cancel{\lambda_0(t)}\exp(\beta_1 \cdot 0)} = \exp(\beta_1)



따라서 병원 A(X=1)와 병원 B(X=0)의 위험비는 $\backslash exp(2.12)\; \backslash approx\; 8.32$이다. 통계적 유의성을 잠시 고려하지 않는다면, 이 결과는 병원 A에서 수술받은 환자가 병원 B에서 수술받은 환자에 비해 특정 시점에서의 순간적인 사망 위험이 약 8.3배 더 높다는 것을 의미한다.

이 해석에는 몇 가지 중요한 주의사항이 따른다.

# 사망 위험이 8.3배 더 높다는 것이 병원 A에서 8.3배 더 많은 환자가 사망한다는 의미는 아니다. 생존 분석은 사건 발생 여부뿐만 아니라 사건 발생 속도(얼마나 빨리 발생하는지)를 함께 고려한다.

# "사망 위험"은 비율(rate)의 개념으로, 단위 시간당 사건 발생률을 의미한다 (예: 속도). 하지만 상대적 위험률인 위험비는 단위가 없는 비교값이다. 즉, 병원 A의 사망 위험률이 병원 B(기준 집단)의 사망 위험률보다 8.3배 더 높다는 의미이다.

# 반대로 병원 B의 위험률은 병원 A에 비해 $1/8.32\; =\; \backslash exp(-2.12)\; \backslash approx\; 0.12$배이다.

# 이 모형만으로는 각 병원에서의 생존 확률 자체를 직접 추론하기는 어렵다. 생존 확률을 계산하려면 $\backslash beta\_1$ 추정치뿐만 아니라 기준선 위험률 $\backslash lambda\_0(t)$의 추정치가 필요하지만, 일반적인 콕스 모형 추정 과정에서는 기준선 위험률을 직접 추정하지 않는 경우가 많다.

# 콕스 모형에서 시간에 따라 변하는 유일한 부분은 기준선 위험률 $\backslash lambda\_0(t)$이다. 위험비 $\backslash exp(\backslash beta\_1)$는 시간에 따라 변하지 않는다고 가정하므로, 시간 단위를 월 대신 년으로 측정했더라도 $\backslash beta\_1$ 추정치와 위험비는 동일하게 계산된다.

# 이 분석 결과는 병원 위치와 사망 위험 간의 연관성을 보여줄 뿐, 병원 자체가 위험률 차이를 유발했다고 단정하기는 어렵다. 데이터가 관찰 연구에서 얻어진 경우, 인과관계를 추론하는 데에는 한계가 있다.

4. 2. 단일 연속형 공변량

'''*''' 종료 날짜는 파산, 매각 등 이벤트 발생일 또는 관찰 종료일(2022-01-01)을 의미한다.

이전 예시들과 달리 이 데이터에는 연속형 변수인 P/E가 포함되어 있다. 콕스 모델은 다음과 같이 설정할 수 있다.

$$\backslash lambda(t|P\_\{i\})\; =\; \backslash lambda\_0(t)\backslash cdot\backslash exp(\backslash beta\_1\; P\_\{i\})$$

여기서 $P\_i$는 회사 ''i''의 P/E 비율을 나타낸다. 이 데이터를 콕스 모델에 적용하면, 모델의 계수 $\backslash beta\_1$의 추정값을 얻을 수 있으며, 이 예시에서는 -0.34라고 가정하자. 따라서 특정 회사의 위험 함수 추정치는 다음과 같다.

$$\backslash lambda(t|P\_\{i\})\; =\; \backslash lambda\_0(t)\backslash cdot\backslash exp(-0.34\; P\_\{i\})$$

콕스 모델의 특성상 기준 위험 $\backslash lambda\_0(t)$ 자체는 추정되지 않으므로, 특정 시점의 절대적인 위험 수준을 계산할 수는 없다. 하지만 두 회사 ''i''와 ''j'' 사이의 상대적인 위험 비율은 계산할 수 있다.

$$\backslash begin\{align\}$$

\frac{\lambda(t|P_{i})}{\lambda(t|P_{j})}

&= \frac{ \cancel{\lambda_0(t)}\cdot\exp(-0.34 P_{i})}{\cancel{\lambda_0(t)}\cdot \exp(-0.34 P_{j})} \\

&= \exp(-0.34 (P_{i} - P_{j}))

\end{align}

이 식의 우변에 있는 값들은 모두 알려진 값(또는 추정된 값 $\backslash beta\_1$과 데이터 $P\_i,\; P\_j$)이므로, 두 회사 간의 위험 비율 계산이 가능하다. 또한, 계산된 위험 비율 식에는 시간에 따라 변하는 항($t$)이 없으므로, 두 회사의 위험 비율은 시간이 지나도 일정하게 유지된다. 이것이 바로 '비례 위험' 가정이다. 예를 들어, P/E가 각각 6.3과 3.0인 회사 5와 회사 2의 위험 비율은 다음과 같이 계산된다: $\backslash exp(-0.34\; (6.3\; -\; 3.0))\; =\; 0.33$. 이는 주어진 기간 내에서 회사 5가 파산이나 매각 등으로 사라질 위험이 회사 2의 약 1/3 수준임을 의미한다. P/E가 높은 회사 5의 위험이 더 낮게 나타난 것이다($\backslash beta\_1$이 음수이므로).

결과 해석 시 다음과 같은 점에 유의해야 한다.

# 위험 비율 ($\backslash exp(\backslash beta\_1)$) 해석: 흔히 "위험 비율"이라고 부르는 값은 $\backslash exp(\backslash beta\_1)$이며, 이 예시에서는 $\backslash exp(-0.34)\; \backslash approx\; 0.71$이다. 이 값은 공변량(여기서는 P/E) 값이 1단위 차이 나는 두 대상 간의 위험 비율로 해석될 수 있다. 즉, 회사 ''i''의 P/E가 회사 ''j''보다 정확히 1만큼 높다면($P\_\{i\}\; =\; P\_\{j\}\; +\; 1$), 회사 ''i''의 위험은 회사 ''j''의 위험의 $\backslash exp(-0.34\; \backslash times\; 1)\; =\; 0.71$배가 된다. P/E가 1단위 높을수록 종료 위험이 약 29% 낮아지는 것이다.

# 기준 위험 ($\backslash lambda\_0(t)$) 해석의 한계: 기준 위험 $\backslash lambda\_0(t)$는 모든 공변량의 값이 0일 때의 위험 함수, 즉 $\backslash lambda(t|P\_\{i\}=0)\; =\; \backslash lambda\_0(t)\backslash cdot\backslash exp(-0.34\; \backslash cdot\; 0)\; =\; \backslash lambda\_0(t)$로 생각할 수 있다. 따라서 기준 위험을 P/E가 0인 "기준 회사"의 위험으로 해석할 수도 있다. 하지만 이 예시에서 P/E가 0이라는 것은 현실적으로 불가능하며(주가가 0이거나 이익이 무한대여야 함), 의미를 부여하기 어렵다. 따라서 "모든 공변량이 0일 때의 위험"이라는 해석이 더 적절하지만, 이 예시처럼 공변량 0이 불가능한 경우에는 기준 위험 자체의 직접적인 해석은 제한적이다.

# 개별 위험 인자 ($\backslash exp(\backslash beta\_1\; P\_\{i\})$) 해석의 문제점: $\backslash exp(\backslash beta\_1\; P\_\{i\})$ 값 자체를 회사 ''i''의 위험 수준을 나타내는 지표로 해석하고 싶을 수 있다. 하지만 이 값은 $\backslash exp(\backslash beta\_1\; P\_\{i\})\; =\; \backslash frac\{\backslash exp(\backslash beta\_1\; P\_\{i\})\}\{\backslash exp(\backslash beta\_1\; \backslash cdot\; 0)\}\; =\; \backslash frac\{\backslash lambda(t|P\_\{i\})\}\{\backslash lambda(t|P=0)\}$와 같이, P/E가 0인 가상의 기준 회사 대비 회사 ''i''의 위험 비율을 의미한다. 위에서 언급했듯이 P/E=0인 기준 회사는 현실성이 없으므로, $\backslash exp(\backslash beta\_1\; P\_\{i\})$ 값 자체를 독립적으로 해석하는 것은 이 예시에서는 의미가 없다. 중요한 것은 P/E 값이 다른 두 실제 회사 간의 위험 비율이다.

5. 시간 종속 예측 변수 및 계수

안데르센과 길(Andersen and Gill)의 계수 과정 공식을 사용하면 콕스 모형을 확장하여 시간에 따라 변하는 변수, 시간에 따라 달라지는 층화, 그리고 각 대상자별로 여러 사건이 발생하는 경우까지 다룰 수 있다.^[8][29] 예를 들어, 시간에 따라 변하는 회귀 변수를 사용하는 위험 모형은 실업 보험이 실업 기간에 미치는 영향을 추정하는 데 사용될 수 있다.^{[9][10][30][31]}

콕스 모형은 시간 종속 공변량(즉, 예측 변수)을 허용하도록 일반화될 수 있다. 즉, 예측 변수의 값이 시간에 따라 변할 수 있다는 의미이다.

나아가, 콕스 모형은 시간 종속 계수를 허용하도록 더 일반화될 수 있다. 이는 어떤 처치(예: 약물 투여)의 효과가 시간에 따라 달라질 수 있음을 의미한다. 예를 들어, 특정 약물이 이환율 발생 직후에는 매우 효과적이지만 시간이 지남에 따라 그 효과가 점차 감소하는 경우가 있을 수 있다.^[11][32] 이러한 경우, 계수가 시간에 따라 변하지 않는다는 가설(정상성 가설)을 검증해 볼 수 있다. 더 자세한 내용과 관련 소프트웨어(R 패키지) 정보는 Martinussen과 Scheike (2006)의 연구에서 찾아볼 수 있다.^{[11][12][32][33]} 신뢰성 계산 분야에서도 시간 종속 공변량을 사용한 콕스 모형의 적용이 고려되고 있다.^[34]

이와 관련하여, 공변량의 효과를 곱셈적이 아닌 덧셈 형태로 지정하는 부가적 위험 모형도 이론적으로 가능하다.^[13][26] 이 모형은 다음과 같이 표현된다.

$$$$

\lambda(t|X_i) = \lambda_0(t) + \beta_1X_{i1} + \cdots + \beta_pX_{ip} = \lambda_0(t) + X_i \cdot \beta.



하지만 로그 가능도를 최대화하는 방식으로 모형을 추정할 경우, 위험 함수 $\backslash lambda(t\backslash mid\; X\_i)$의 값이 음수가 되지 않도록 제약을 가해야 하는 복잡성이 있다. 이러한 이유로 부가적 위험 모형은 실제 분석에서는 자주 사용되지 않는다. 반면, 최소 제곱법을 목표로 할 경우에는 음수 제한이 엄격하게 요구되지는 않는다.

6. 기준 위험 함수 지정

콕스 모형은 특정 상황에서 기준 위험 함수 $\backslash lambda\_0(t)$를 구체적인 함수 형태로 지정하여 특화할 수 있다. 예를 들어, 기준 위험 함수가 와이블 분포의 위험 함수 형태를 따른다고 가정하면, 이를 '와이블 비례 위험 모형'이라고 부른다. 와이블 분포를 기준 위험 함수로 사용하는 경우는 해당 모형이 비례 위험 모형이면서 동시에 가속 고장 시간 모형의 성질도 만족하는 유일한 경우이다.

이처럼 기준 위험 함수를 특정 분포로 지정하는 모형을 '모수적 비례 위험 모형'이라고 하며, 기준 위험 함수를 특정하지 않는 일반적인 콕스 비례 위험 모형은 '반모수 모형'으로 분류된다.

일부 연구자들은 기준 위험 함수를 특정 모수적 형태로 지정하는 경우에도, 이 분야에 대한 데이비드 콕스의 기여를 인정하여 '콕스 비례 위험 모형'이라는 용어를 사용하기도 한다.^[14][35]

'콕스 회귀 모형'이라는 용어는 시간 의존적 요인을 포함하도록 확장된 모형을 지칭하는 데 쓰이기도 하지만, 콕스 비례 위험 모형 자체가 회귀분석의 일종으로 볼 수 있으므로 용어 사용에 혼란이 있을 수 있다.

7. 푸아송 모델과의 관계

비례 위험 모형과 푸아송 회귀 모형 사이에는 관계가 있다. 이 관계를 이용하여, 푸아송 회귀 분석용 소프트웨어를 통해 비례 위험 모형을 근사적으로 적합시키기도 한다. 이렇게 하는 주된 이유는 계산 속도가 훨씬 빠르다는 장점 때문이다. 과거 컴퓨터 성능이 낮았던 시절에는 이러한 방식이 더욱 중요하게 여겨졌으나, 오늘날에도 대규모 데이터 집합이나 복잡한 문제를 다룰 때 여전히 유용하게 사용될 수 있다. Laird와 Olivier (1981)는 이 관계에 대한 수학적 세부 사항을 제시했으며^[15], 그들은 "우리는 [푸아송 모형]이 참이라고 가정하는 것이 아니라, 단지 우도(likelihood)를 유도하기 위한 장치로 사용할 뿐이다"라고 언급했다. McCullagh와 Nelder의 저서 『일반화 선형 모형』에도 비례 위험 모형을 일반화 선형 모형으로 변환하는 방법에 대한 장이 있다.^[16]

8. 고차원 설정

고차원 설정에서는 공변량의 수 p가 표본 크기 n에 비해 클 수 있다. 이러한 경우 Lasso 방법은 모형을 선택하는 데 사용되는 고전적인 전략 중 하나이다. Tibshirani(1997)는 비례 위험 회귀 모수 추정을 위한 Lasso 절차를 제안했다.^[17] 회귀 모수 β의 Lasso 추정량은 L¹-노름 제약 조건 하에서 Cox 부분 로그-우도 함수의 음의 값을 최소화하는 값으로 정의된다.

:$$

\ell(\beta) = \sum_j \left(\sum_{i\in H_j} X_i \cdot \beta -\sum_{\ell=0}^{m_j-1}\log\left(\sum_{i:Y_i\ge t_j}\theta_i - \frac{\ell}{m_j}\sum_{i\in H_j}\theta_i\right)\right)

+ \lambda \|\beta\|_1



최근 이 주제에 대한 이론적인 발전이 있었다.^{[18][19][20][21]}

9. 소프트웨어 구현

• '''매스매티카(Wolfram Mathematica)''': ''CoxModelFit'' 함수.^[22]
• '''R''': '''survival''' 패키지 내의 ''coxph()'' 함수.
• '''SAS''': ''phreg'' 프로시저.
• '''Stata''': ''stcox'' 명령어.
• '''파이썬(Python)''': '''lifelines''' 라이브러리 내의 ''CoxPHFitter'', '''statsmodels''' 라이브러리 내의 ''phreg''.
• '''SPSS''': '''Cox 회귀''' 분석 기능에서 사용 가능.
• '''MATLAB''': ''fitcox'' 또는 ''coxphfit'' 함수.
• '''줄리아(Julia)''': '''Survival.jl''' 라이브러리에서 사용 가능.
• '''JMP''': '''비례 위험 적합''' 플랫폼에서 사용 가능.
• '''프리즘(Prism)''': 생존 분석 및 다변수 분석 기능에서 사용 가능.

참조

_[1] 논문 Analysis of Survival Data under the Proportional Hazards Model
_[2] 논문 Regression Models and Life-Tables
_[3] 논문 Fifty Years of the Cox Model 2023-03-10
_[4] 논문 A Conversation with Sir David Cox
_[5] 간행물 Some remarks on the analysis of survival data
_[6] 문서 Each failure contributes to the likelihood function
_[7] 논문 The Efficiency of Cox's Likelihood Function for Censored Data
_[8] 논문 Cox's regression model for counting processes, a large sample study.
_[9] 논문 Unemployment Insurance and Unemployment Spells http://www.nber.org/[...]
_[10] 논문 Unemployment Duration, Benefit Duration, and the Business Cycle http://www.bde.es/f/[...]
_[11] 서적 Dynamic Regression Models for Survival Data Springer
_[12] 웹사이트 timereg: Flexible Regression Models for Survival Data https://cran.r-proje[...]
_[13] 간행물 Some remarks on the analysis of survival data
_[14] 논문 Generating survival times to simulate Cox proportional hazards models
_[15] 논문 Covariance Analysis of Censored Survival Data Using Log-Linear Analysis Techniques
_[16] 서적 Generalized Linear Models Chapman & Hall/CRC
_[17] 논문 The Lasso method for variable selection in the Cox model
_[18] 논문 Regularization for Cox's proportional hazards model with NP-dimensionality
_[19] 논문 Structured Estimation in Nonparametric Cox Model
_[20] 논문 Non-asymptotic oracle inequalities for the high-dimensional Cox regression via Lasso
_[21] 논문 Oracle inequalities for the lasso in the Cox model
_[22] 웹사이트 CoxModelFit https://reference.wo[...]
_[23] 논문 Analysis of Survival Data under the Proportional Hazards Model
_[24] 논문 Regression Models and Life-Tables
_[25] 논문 A Conversation with Sir David Cox
_[26] 간행물 Some remarks on the analysis of survival data
_[27] 문서 Each failure contributes to the likelihood function
_[28] 논문 The Efficiency of Cox's Likelihood Function for Censored Data
_[29] 논문 Cox's regression model for counting processes, a large sample study.
_[30] 논문 Unemployment Insurance and Unemployment Spells http://www.nber.org/[...]
_[31] 논문 Unemployment Duration, Benefit Duration, and the Business Cycle http://www.bde.es/f/[...]
_[32] 서적 Dynamic Regression Models for Survival Data Springer
_[33] 웹사이트 timereg: Flexible Regression Models for Survival Data https://cran.r-proje[...] 2021-10-17
_[34] 논문 Decline and repair, and covariate effects http://usir.salford.[...]
_[35] 논문 Generating survival times to simulate Cox proportional hazards models
_[36] 논문 Covariance Analysis of Censored Survival Data Using Log-Linear Analysis Techniques
_[37] 서적 Generalized Linear Models Chapman & Hall/CRC
_[38] 논문 The Lasso method for variable selection in the Cox model
_[39] 논문 Regularization for Cox's proportional hazards model with NP-dimensionality
_[40] 논문 Structured Estimation in Nonparametric Cox Model
_[41] 논문 Non-asymptotic oracle inequalities for the high-dimensional Cox regression via Lasso
_[42] 논문 Oracle inequalities for the lasso in the Cox model
_[43] 저널 Survival Analysis Part II: Multivariate data analysis – an introduction to concepts and methods 2003-08-04