§ 9.1–9.2 — Time-Dependent Covariates

1 도입 — 왜 공변량이 “움직여야” 하는가

Ch.8의 표준 Cox 모형은 연구 시작 시점에 관측된 고정 공변량(fixed-time covariates) 만 허용한다.

\[ h(t \mid \mathbf{Z}) = h_0(t) \exp(\boldsymbol{\beta}^\top \mathbf{Z}) \]

여기서 \(\mathbf{Z}\) 는 시점 0에서 측정한 후 변하지 않는다. 성별, 질병 유형, 이식 시 연령 등이 대표적이다.

그러나 실제 생존 연구에서는 추적 기간 동안 값이 바뀌는 변수가 흔하다:

변수 유형	예시	왜 고정 공변량으로 불충분한가
이진 중간 사건	혈소판 회복, aGVHD 발생	사건 발생 전후로 위험이 근본적으로 달라진다
연속 반복 측정	혈압, 종양 크기, 바이오마커	현재 값이 미래 위험을 결정한다
실험자 통제 처치	치료 스케줄 변경	계획된 개입이 위험을 변동시킨다

이러한 변수를 시간의존 공변량(time-dependent covariates) 이라 한다. 직관적으로, “환자의 상태가 시간에 따라 진화하므로, 위험률 역시 그 진화에 반응해야 한다”는 것이 핵심 동기이다 (Klein & Moeschberger, 2003, Ch.9).

2 확장 모형 — 식 (9.1.1)

고정 공변량 \(\mathbf{Z}\) 를 시간의 함수 \(\mathbf{Z}(t)\) 로 바꾸면 다음과 같다:

\[ h[t \mid \mathbf{Z}(t)] = h_0(t) \exp\!\left[\sum_{k=1}^{p} \beta_k Z_k(t)\right] \]

직관

“시점 \(t\) 에서의 순간 위험은 바로 그 시점의 공변량 값에 의해 결정된다.”

혈소판이 아직 회복되지 않은 환자( \(Z_P(t)=0\) )와 이미 회복된 환자( \(Z_P(t)=1\) )는 같은 시점이라도 서로 다른 위험을 경험한다. 고정 공변량 모형에서는 이를 표현할 방법이 없다.

이 모형은 Ch.8의 식 (8.1.2)를 포함하는 상위 모형이다 — \(Z_k(t)\) 가 상수이면 표준 Cox 모형으로 축소된다.

3 편우도(Partial Likelihood) 확장 — 식 (9.2.1)

3.1 데이터 구조

표본 크기 \(n\) 에서 \(j\) 번째 개인의 데이터는 삼중체(triple)로 구성된다:

\[ \bigl[T_j,\; \delta_j,\; \{\mathbf{Z}_j(t),\; 0 \le t \le T_j\}\bigr], \quad j = 1, \ldots, n \]

\(T_j\) : 관찰 종료 시점 (사건 또는 중도절단)
\(\delta_j\) : 사건 지시자 ( \(\delta_j = 1\) 이면 사건 발생)
\(\mathbf{Z}_j(t)\) : 관찰 기간 동안의 공변량 경로(covariate path)

고정 공변량 모형과의 핵심 차이는, 각 개인의 정보가 하나의 벡터가 아닌 함수(프로세스) 라는 점이다.

3.2 확장 편우도

사건 시점을 \(t_1 < t_2 < \cdots < t_D\) 로 정렬하면, 편우도는 다음과 같다:

\[ L(\boldsymbol{\beta}) = \prod_{i=1}^{D} \frac{\exp\!\left[\boldsymbol{\beta}^\top \mathbf{Z}_{(i)}(t_i)\right]}{\sum_{j \in R(t_i)} \exp\!\left[\boldsymbol{\beta}^\top \mathbf{Z}_j(t_i)\right]} \]

직관 — 각 분수의 의미

분자: 시점 \(t_i\) 에서 실제로 사건이 발생한 개인의 그 시점 공변량 값으로 계산한 상대 위험.

분모: 시점 \(t_i\) 직전에 아직 위험에 노출된(at risk) 모든 개인의 그 시점 공변량 값으로 계산한 상대 위험의 합.

고정 공변량 편우도(식 8.2.1)에서는 분모의 \(\mathbf{Z}_j\) 가 시점과 무관했지만, 여기서는 매 사건 시점마다 위험집합 구성원의 현재 공변량 값을 재평가한다. 이것이 “시간이 흐르면서 환자 상태가 바뀐다”는 현실을 반영하는 메커니즘이다.

3.3 추정과 검정

편우도를 \(\boldsymbol{\beta}\) 에 대해 최대화하는 절차는 Ch.8과 동일하다 — Newton-Raphson 반복법으로 점수 방정식 \(U(\boldsymbol{\beta}) = 0\) 을 풀고, 관측 정보 행렬로 표준 오차를 구한다. 동점(ties)이 존재하면 Breslow 또는 Efron 근사를 동일하게 적용한다.

4 내부 vs 외부 공변량 — 해석의 경계

Kalbfleisch and Prentice (1980)는 시간의존 공변량을 두 유형으로 구분한다.

유형	정의	예시	생존함수 추정
외부(external)	공변량의 경로가 실패 과정과 독립적이다	일일 기온, 계획된 투약 스케줄, 고정 공변량	가능 — \(S(t \mid \mathbf{Z}(t))\) 직접 추정
내부(internal)	공변량이 관측되려면 개인이 살아 있어야 한다	혈소판 회복 시점, aGVHD 발생, 종양 크기	불가능 — \(P[X \ge t \mid \mathbf{Z}(t)] = 1\) (자명)

왜 구분이 중요한가

내부 공변량의 경우, 조건부 생존함수 \(S(t \mid \mathbf{Z}(t))\) 를 추정하려면 “시점 \(t\) 에서 \(\mathbf{Z}(t)\) 가 관측되었다”는 사실 자체가 “개인이 시점 \(t\) 까지 생존했다”를 내포하므로, 결과가 항상 1이 되어 무의미하다.

따라서 내부 공변량이 포함된 모형에서는 위험비(hazard ratio) 해석에 집중하고, 생존 확률 예측은 별도의 다상태 프레임워크(§9.5)를 필요로 한다.

외부 공변량은 개인의 생사와 무관하게 값이 결정되므로, Ch.8의 생존함수 추정법(§8.8)을 그대로 적용할 수 있다.

5 PH 가정의 형식적 검정 — 시간의존 공변량 활용

5.1 아이디어

시간의존 공변량은 PH 가정이 맞는지 형식적으로(formally) 검정하는 도구로도 쓰인다.

고정 공변량 \(Z_1\) 에 대해 PH 가정이 성립한다면, 위험비 \(h(t \mid Z_1) / h(t \mid Z_1^*)\) 는 시간에 의존하지 않는다. 이를 검정하기 위해 인위적(artificial) 시간의존 공변량을 만든다:

\[ Z_2(t) = Z_1 \times g(t) \]

여기서 \(g(t)\) 는 알려진 시간 함수이다. 실무에서 가장 많이 쓰는 선택은 \(g(t) = \ln t\) 이다.

5.2 검정 모형과 해석

\(Z_1\) 과 \(Z_2(t) = Z_1 \cdot \ln t\) 를 함께 모형에 넣으면:

\[ h(t \mid Z_1) = h_0(t) \exp\!\bigl[\beta_1 Z_1 + \beta_2 (Z_1 \cdot \ln t)\bigr] \]

두 개인( \(Z_1\) vs \(Z_1^*\) )의 위험비를 구하면:

\[ \frac{h(t \mid Z_1)}{h(t \mid Z_1^*)} = \exp\!\bigl\{(\beta_1 + \beta_2 \ln t)(Z_1 - Z_1^*)\bigr\} \]

직관

\(\beta_2 = 0\) 이면 위험비가 시간과 무관하게 일정하다 — PH 가정이 성립한다.

\(\beta_2 \ne 0\) 이면 위험비가 \(t\) 에 따라 변한다 — PH 가정이 깨진다.

따라서 \(H_0: \beta_2 = 0\) 을 검정하는 것이 곧 PH 가정의 검정이다.

이 검정의 검정력(power)은 \(g(t)\) 의 선택에 의존한다. \(\ln t\) 는 대부분의 단조 이탈(departure)을 잡아내지만, 비단조적 패턴에는 다른 선택이 필요할 수 있다.

5.3 예제: 신장투석 카테터 데이터 (Example 9.2)

신장 부전 환자에게 경피적(percutaneous, \(Z_1=1\) ) 또는 수술적(surgical, \(Z_1=0\) ) 카테터를 배치한 후, 첫 감염까지의 시간을 추적한다.

\(Z_2(t) = Z_1 \times \ln t\) 로 검정한 결과:

우도비 검정: \(\chi^2 = 12.22\) ( \(p = 0.0005\) )
Wald 검정: \(\chi^2 = 6.19\) ( \(p = 0.013\) )

결론: PH 가정이 강하게 기각된다 — 경피적 카테터의 상대 위험이 시간이 지남에 따라 변한다.

6 PH 위반 시 대안 — 변환점(Change-Point) 모형

6.1 문제

\(g(t) = \ln t\) 처럼 특정 함수 형태를 가정하면, 실제 이탈 패턴을 제대로 포착하지 못할 수 있다. 대안은 위험비가 특정 시점 \(\tau\) 를 기준으로 구간별(piecewise) 상수라고 가정하는 것이다.

6.2 모형 구성 — 식 (9.2.3)–(9.2.4)

코딩 1 (원본 \(Z_1\) + 지시 공변량):

\[ Z_2(t) = \begin{cases} Z_1 & \text{if } t > \tau \\ 0 & \text{if } t \le \tau \end{cases} \]

이 모형의 위험률:

\[ h[t \mid Z(t)] = \begin{cases} h_0(t)\, e^{\beta_1 Z_1} & \text{if } t \le \tau \\ h_0(t)\, e^{(\beta_1 + \beta_2) Z_1} & \text{if } t > \tau \end{cases} \]

직관 — 변환점의 의미

\(e^{\beta_1}\) : 시점 \(\tau\) 이전의 상대 위험
\(e^{\beta_1 + \beta_2}\) : 시점 \(\tau\) 이후의 상대 위험
\(e^{\beta_2}\) : \(\tau\) 전후의 상대 위험 변화량

PH 가정은 “한 조건에서 모든 시간에 걸쳐 위험비가 동일”을 요구하지만, 변환점 모형은 이를 “두 구간에서 각각 PH가 성립”으로 완화한다.

코딩 2 (구간별 직접 추정):

\[ Z_2(t) = \begin{cases} Z_1 & t > \tau \\ 0 & t \le \tau \end{cases}, \quad Z_3(t) = \begin{cases} Z_1 & t \le \tau \\ 0 & t > \tau \end{cases} \]

이 코딩에서 \(e^{\theta_3}\) 은 \(\tau\) 이전의 상대 위험, \(e^{\theta_2}\) 는 \(\tau\) 이후의 상대 위험을 직접 추정한다. 두 코딩은 동일한 로그우도를 산출하며 \(\beta_1 = \theta_3\) , \(\beta_1 + \beta_2 = \theta_2\) 관계를 갖는다.

6.3 최적 변환점 \(\tau\) 의 결정

우도는 사건 시점에서만 값이 바뀌므로, 모든 사건 시점 \(t_1, t_2, \ldots, t_D\) 에 대해 모형을 적합하고 로그 편우도를 최대화하는 \(\tau\) 를 선택한다.

절차	설명
후보 격자	모든 사건 시점을 \(\tau\) 후보로 설정한다
모형 적합	각 \(\tau\) 에서 확장 편우도를 최대화한다
선택	최대 로그 편우도를 주는 \(\tau\) 를 최적 변환점으로 채택한다
재검증	선택된 구간 내에서 다시 PH를 검정하여, 추가 분할이 필요한지 확인한다

6.4 예제 계속: 카테터 데이터의 변환점

모든 사건 시점에 대해 로그 편우도를 계산한 결과, \(\tau = 3.5\) 개월에서 최대값을 기록했다.

구간	\(b\)	SE\((b)\)	Wald \(\chi^2\)	\(p\)	해석
\(t \le 3.5\) 개월	\(-2.089\)	0.760	7.56	0.006	경피적 카테터가 유의하게 우수 ( \(e^{-2.089} = 0.124\) )
\(t > 3.5\) 개월	\(1.081\)	0.783	1.91	0.167	두 군 간 차이 없음

임상적 직관

이식 초기(3.5개월 이내)에는 경피적 카테터의 감염 위험이 수술적 카테터의 약 1/8이다. 그러나 3.5개월 이후에는 카테터 유형에 따른 차이가 사라진다. 이는 초기 상처 치유 과정에서의 차이가 장기적으로는 소멸함을 시사한다.

구간 내 PH 검정 결과(추가 \(Z_1 \cdot \ln t\) 항의 Wald p-value: 0.82, 0.28), 각 구간 내에서 PH 가정이 유지됨을 확인했다.

7 BMT 데이터에서의 종합 활용 (Example 9.1)

7.1 시간의존 공변량 정의

골수이식(BMT) 데이터에서 세 가지 중간 사건을 지시 공변량으로 코딩한다:

\[ Z_A(t) = \begin{cases} 1 & t \ge \text{aGVHD 발생 시점} \\ 0 & \text{otherwise} \end{cases} \]

\[ Z_P(t) = \begin{cases} 1 & t \ge \text{혈소판 회복 시점} \\ 0 & \text{otherwise} \end{cases} \]

\[ Z_C(t) = \begin{cases} 1 & t \ge \text{cGVHD 발생 시점} \\ 0 & \text{otherwise} \end{cases} \]

각 공변량은 해당 사건이 발생하기 전에는 0, 발생 후에는 1로 전환되는 계단 함수(step function) 이다.

7.2 단변량 분석 결과

질병 유형( \(Z_1, Z_2\) )을 보정한 단변량 분석에서:

시간의존 공변량	\(\chi^2\) (LR)	\(p\)	해석
\(Z_A(t)\) : aGVHD	1.17	0.28	무의미
\(Z_C(t)\) : cGVHD	0.46	0.50	무의미
\(Z_P(t)\) : 혈소판 회복	9.64	0.002	유의

혈소판 회복만이 무병생존(disease-free survival)에 유의한 영향을 미친다. \(b_P = -1.130\) 이므로 \(e^{-1.130} = 0.323\) — 혈소판 회복 후의 사망/재발 위험은 회복 전의 약 1/3로 감소한다.

7.3 고정 + 시간의존 결합 모형

고정 인자(FAB, 연령)와 시간의존 인자( \(Z_P(t)\) )를 동시에 포함한 결합 모형에서:

시간의존 효과의 추가 검정: \(-2[-356.99 - (-353.31)] = 7.36\) (1 df, \(p = 0.0067\) )
고정 효과의 추가 검정: \(-2[-361.82 - (-353.31)] = 17.02\) (4 df, \(p = 0.0019\) )

두 유형 모두 유의하므로, 질병 유형 비교 시 고정·시간의존 인자를 함께 보정해야 한다.

7.4 교호작용 — 혈소판 회복 전후의 효과 차이

고정 인자와 혈소판 회복의 교호작용을 전진 선택법으로 추가한 결과, 세 교호작용(Group \(\times Z_P\) , FAB \(\times Z_P\) , Age \(\times Z_P\) ) 모두 유의했다.

핵심 해석 (Table 9.4의 \(Z_1\) 과 \(Z_7(t) = Z_1 \times Z_P(t)\) ):

혈소판 회복 전: AML low-risk vs ALL의 상대 위험 = \(e^{1.3073} = 3.70\) (유의하지 않음, CI: 0.74–18.39)
혈소판 회복 후: 상대 위험 = \(e^{1.3073 + (-3.0374)} = e^{-1.73} = 0.18\) (유의, CI: 0.08–0.41)

직관

혈소판 회복 이전에는 AML low-risk과 ALL 환자의 예후 차이가 명확하지 않다. 그러나 혈소판이 회복된 후에는 AML low-risk 환자의 재발/사망 위험이 ALL의 18%에 불과하다. 즉, 위험 유형 간 예후 차이는 혈소판 회복이라는 이정표를 지나야 드러난다.

7.5 BMT 데이터에서의 PH 가정 검정

각 고정 인자에 대해 \(Z_i \cdot \ln t\) 형태의 인위적 시간의존 공변량을 만들어 PH를 검정한다 (Table 9.5):

인자	Wald \(\chi^2\)	df	\(p\)	결론
Group (질병 유형)	1.74	2	0.42	PH 유지
Waiting time	0.01	1	0.94	PH 유지
FAB status	0.44	1	0.51	PH 유지
MTX	4.32	1	0.038	PH 위반
Sex	0.22	3	0.97	PH 유지
CMV status	1.69	3	0.64	PH 유지
Age	4.76	3	0.19	PH 유지

MTX(methotrexate) 사용 여부만 PH 가정을 위반한다. 이 결과는 §9.3에서 MTX를 층화 변수로 활용하는 근거가 된다.

8 데이터 구조 — 계수과정(Counting Process) 형식

시간의존 공변량을 소프트웨어에서 다루려면, 각 개인의 관찰 기간을 공변량이 일정한 (시작, 끝] 구간들로 분할한다.

예시: 혈소판이 20일에 회복되고 45일에 사건(사망)이 발생한 환자

구간 시작	구간 끝	사건	\(Z_P\)
0	20	0	0
20	45	1	1

직관

하나의 환자가 여러 행(row)으로 표현된다. 각 행은 “이 구간 동안 공변량이 이 값으로 일정했다”를 나타낸다. 사건 지시자는 마지막 구간의 끝점에서만 1이 된다.

이 형식은 R의 Surv(start, stop, event) 또는 Python의 lifelines 에서 episode-splitting으로 구현한다.

9 실무 주의사항 (Practical Notes)

항목	내용
동점 처리	SAS PHREG: Breslow (기본) / Efron / Exact; R `coxph`: Efron (기본)
고정 vs 시간의존	연구 시작 시점에 알 수 없는 변수(중간 사건)는 반드시 시간의존 코딩한다
생존함수 추정	시간의존 공변량이 있으면 \(\int h_0(t) e^{\beta^\top Z(t)} dt\) 가 확률 과정에 의존하므로, 외부 공변량이 아닌 한 단순 Breslow 추정이 불가하다
불멸 시간 편향(immortal time bias)	중간 사건 발생 전 시간을 처리군에 잘못 포함하면 결과가 왜곡된다 — 시간의존 코딩이 이를 방지한다

10 코드 예시

10.1 Step 1: 순수 Python으로 확장 편우도 계산

import numpy as np

# 간소화 데이터: (사건시점, 사건여부, 공변량경로)
# 구간 형식: [(start, stop, event, Z_P)]
episodes = [
    # 환자1: 혈소판 20일 회복, 45일 사망
    (0, 20, 0, 0),
    (20, 45, 1, 1),
    # 환자2: 혈소판 미회복, 30일 중도절단
    (0, 30, 0, 0),
    # 환자3: 혈소판 15일 회복, 60일 중도절단
    (0, 15, 0, 0),
    (15, 60, 0, 1),
]

def extended_partial_likelihood(beta, episodes):
    """확장 편우도의 로그값을 계산한다."""
    # 사건 시점 추출
    event_times = sorted(set(
        stop for start, stop, event, z in episodes if event == 1
    ))
    
    log_lik = 0.0
    for t in event_times:
        # 사건 발생자의 공변량
        z_event = None
        for start, stop, event, z in episodes:
            if stop == t and event == 1:
                z_event = z
                break
        
        # 위험집합: 시점 t 직전에 아직 관찰 중인 개인들
        risk_set_z = []
        for start, stop, event, z in episodes:
            if start < t <= stop:
                risk_set_z.append(z)
        
        # 편우도 기여분
        numerator = np.exp(beta * z_event)
        denominator = sum(np.exp(beta * z) for z in risk_set_z)
        log_lik += np.log(numerator / denominator)
    
    return log_lik

# 격자 탐색으로 beta 추정
betas = np.linspace(-3, 1, 100)
log_liks = [extended_partial_likelihood(b, episodes) for b in betas]
beta_hat = betas[np.argmax(log_liks)]
print(f"beta_hat = {beta_hat:.3f}")
print(f"hazard ratio = exp(beta) = {np.exp(beta_hat):.3f}")

10.2 Step 2: lifelines / R survival 패키지 활용

# Python: lifelines (episode-splitting 형식)
import pandas as pd
from lifelines import CoxPHFitter

# 계수과정 형식 데이터프레임
df = pd.DataFrame({
    'id': [1, 1, 2, 3, 3],
    'start': [0, 20, 0, 0, 15],
    'stop': [20, 45, 30, 15, 60],
    'event': [0, 1, 0, 0, 0],
    'platelet_recovery': [0, 1, 0, 0, 1],
})

cph = CoxPHFitter()
cph.fit(
    df,
    duration_col='stop',
    event_col='event',
    entry_col='start',  # 좌절단 시작점
    formula='platelet_recovery',
)
cph.print_summary()

# R: survival 패키지
library(survival)

# 계수과정 형식
df <- data.frame(
  id     = c(1, 1, 2, 3, 3),
  start  = c(0, 20, 0, 0, 15),
  stop   = c(20, 45, 30, 15, 60),
  event  = c(0, 1, 0, 0, 0),
  zp     = c(0, 1, 0, 0, 1)
)

fit <- coxph(
  Surv(start, stop, event) ~ zp,
  data = df
)
summary(fit)

# PH 가정 검정 (Schoenfeld 잔차 기반)
cox.zph(fit)

11 요약 — §9.1–9.2 핵심 흐름

고정 공변량 Cox 모형 (Ch.8)
    │
    ▼ 한계: 시간에 따라 값이 변하는 변수를 포함 불가
    │
확장: Z → Z(t) — 시간의존 공변량 모형 (식 9.1.1)
    │
    ├─ 편우도 확장 (식 9.2.1): 매 사건 시점의 현재 공변량 값 사용
    │
    ├─ 내부 vs 외부 공변량 구분 → 해석 범위 결정
    │
    ├─ PH 가정 검정: Z₂(t) = Z₁ × g(t) 추가 후 β₂ = 0 검정
    │
    └─ 변환점 모형: τ 전후 구간별 PH 적합
         └─ 모든 사건 시점에 대해 우도 최대화 → 최적 τ

12 관련 주제

선행 지식

후속 주제

관련 개념