절단자료의 최대가능도¶

앞 절들에서 여러 모수적 생존 모형 --- 지수, 와이불, 로그정규, 로그로지스틱 --- 을 소개하고 지수 모형의 MLE를 유도했다. 이 절에서는 모든 모수적 생존 모형의 모수 추정을 떠받치는 절단자료의 일반적인 가능도 틀을 전개한다.

핵심 착상은 관측된 사건과 절단된 관측치가 가능도에 서로 다른 인자를 기여한다는 것이다. 시점 \(t_i\)의 사건은 밀도 \(f(t_i)\)를, 시점 \(t_i\)의 절단은 생존확률 \(S(t_i)\)를 기여한다.

절단자료의 가능도¶

\((t_1, \delta_1), \ldots, (t_n, \delta_n)\)을 관측 자료라 하자. \(t_i\)는 관측된 시간이고 \(\delta_i \in \{0, 1\}\)은 사건 지시자다(사건이면 \(\delta_i = 1\), 절단이면 \(\delta_i = 0\)).

독립 절단 가정 아래에서 모수벡터 \(\boldsymbol{\theta}\)의 가능도는

이다.

해석:

• \(\delta_i = 1\)이면 그 대상이 정확히 \(t_i\)에 사건을 겪었으므로 기여는 밀도 \(f(t_i)\)다.
• \(\delta_i = 0\)이면 그 대상이 \(t_i\)에 절단되어 우리가 아는 것은 \(T_i > t_i\)뿐이므로 기여는 \(P(T_i > t_i) = S(t_i)\)다.

로그가능도¶

\(f(t) = h(t) S(t)\)와 \(S(t) = \exp(-H(t))\)를 쓰면 가능도는

이 되고, 로그가능도는

이다. 이 분리된 형태는 계산에 편리하다. 위험함수와 누적위험함수가 어떤 모수 모형에서든 기본 구성요소이기 때문이다.

왼쪽은 지수 모형에서 \(t = 10\)에 관측된 대상 한 명의 가능도 기여를 \(\lambda\)의 함수로 그린 것이다. 사건으로 관측되었으면 기여가 \(f(10) = \lambda e^{-10\lambda}\)이고, 이 함수는 \(\hat\lambda = 1/10 = 0.10\)에서 뚜렷한 봉우리를 갖는다. "10개월에 사건이 일어났다"는 진술이 \(\lambda\)가 대략 얼마여야 하는지를 위아래 양쪽에서 좁혀 주기 때문이다.

절단으로 관측되었으면 기여가 \(S(10) = e^{-10\lambda}\)이고, 이 함수에는 봉우리가 없다. \(\lambda\)가 작을수록 단조 증가할 뿐이다. "적어도 10개월은 버텼다"는 진술은 \(\lambda\)가 너무 크지 않다는 상한만 말해 주고 하한은 전혀 말해 주지 않는다. 절단 관측이 정보를 주기는 하되 한쪽으로만 준다는 것이 이 두 곡선의 차이다. 앞의 상자에서 "절단된 관측치를 버리면 편향된다"고 한 것과 "절단만으로는 모수를 정할 수 없다"는 것이 동시에 참인 이유다.

오른쪽은 그 결과다. 와이불 \(k = 1.7\), \(\lambda = 50\)에서 \(n = 300\)을 뽑고 절단의 강도만 바꿔 가며 형상모수의 프로파일 로그가능도 \(\ell(k) - \ell(\hat k)\)를 그렸다. 절단율 \(13\%\)(사건 262건)에서는 곡선이 뾰족해 프로파일 95% 구간이 \((1.41,\ 1.69)\)로 폭 \(0.28\)이다. 절단율 \(70\%\)(사건 90건)에서는 폭이 \(0.54\)로 두 배, 절단율 \(95\%\)(사건 15건)에서는 \((1.17,\ 2.70)\)으로 폭이 \(1.53\)까지 부푼다. 마지막 경우는 "위험이 감소한다"(\(k<1\))만 겨우 배제할 뿐 그 이상 아무 말도 하지 못한다.

표본 크기 \(n = 300\)은 세 경우 모두 같았다는 점이 핵심이다. 달라진 것은 사건 수뿐이다. 왼쪽 그림이 그 이유를 이미 말했다. 봉우리를 만드는 것은 \(\delta_i \ln h(t_i)\) 항이고, 절단 관측은 그 항을 기여하지 않는다. 그래서 생존분석에서는 표본 크기가 아니라 사건 수를 세야 하며, 임상시험 설계가 "환자 몇 명"이 아니라 "사건 몇 건"을 목표로 잡는 것도 이 때문이다.

점수방정식¶

점수함수는 로그가능도의 기울기다.

MLE \(\hat{\boldsymbol{\theta}}\)는 \(U(\hat{\boldsymbol{\theta}}) = \mathbf{0}\)을 푼다.

정보행렬¶

관측 피셔 정보행렬은

를 \(\hat{\boldsymbol{\theta}}\)에서 평가한 것이다. MLE의 점근 공분산은 \(\text{Var}(\hat{\boldsymbol{\theta}}) \approx \mathcal{I}(\hat{\boldsymbol{\theta}})^{-1}\) 이다.

개별 모수의 신뢰구간은

이다.

수치 최적화¶

지수 모형을 제외하면 생존 모형의 점수방정식에는 닫힌 형태의 해가 없다. 표준적인 수치 방법은 다음과 같다.

1. 뉴턴-랩슨. 헤세행렬(관측정보)을 써서 매 반복 모수 추정치를 갱신한다.

1. 피셔 점수법. 관측정보를 기대정보로 대체한다. 모형이 옳게 지정되어 있으면 뉴턴-랩슨과 동등하다.

2. 준뉴턴법(BFGS, L-BFGS). 이차도함수 계산을 피하려고 헤세행렬을 근사한다. 소프트웨어 패키지에서 널리 쓰인다.

와이불 가능도¶

모수가 \((k, \lambda)\)인 와이불 모형에서

이므로 로그가능도는

이다.

가능도를 이용한 모형 비교¶

모수 모형은 가능도 기반 기준으로 비교할 수 있다.

가능도비 검정(내포 모형용). 지수 모형(\(k=1\))을 와이불 모형(\(k\) 자유)과 견주려면

를 쓴다.

아카이케 정보기준(비내포 모형용):

여기서 \(p\)는 모수의 개수다. AIC가 낮을수록 적합도와 복잡도의 절충이 낫다. 21.5절에서 AIC와 다른 선택 기준을 자세히 다룬다.

좌측절단과 구간절단 다루기¶

가능도는 다른 절단 유형으로도 확장된다.

• 좌측절단(\(T_i < t_i\)): \(F(t_i) = 1 - S(t_i)\)를 기여한다.
• 구간절단(\(L_i < T_i \leq R_i\)): \(S(L_i) - S(R_i) = F(R_i) - F(L_i)\)를 기여한다.

여러 절단 유형이 섞인 경우의 일반적인 가능도는

이다. 여기서 \(\mathcal{E}\), \(\mathcal{R}\), \(\mathcal{L}\), \(\mathcal{I}\)는 각각 정확 관측, 우측절단, 좌측절단, 구간절단 관측치의 집합이다.

연습문제¶

import numpy as np
from scipy.optimize import minimize
from scipy import stats

rng = np.random.default_rng(31)
n = 400
k_true, lam_true = 1.7, 50.0
T = lam_true * rng.weibull(k_true, n)
C = rng.exponential(70.0, n)
t = np.minimum(T, C); d = (T <= C).astype(int)

def neg_loglik(p):
    k, lam = np.exp(p)                 # optimize on the log scale
    h = (k / lam) * (t / lam) ** (k - 1)
    H = (t / lam) ** k
    return -np.sum(d * np.log(h) - H)

res = minimize(neg_loglik, [0.0, np.log(t.mean())], method='BFGS')
k_hat, lam_hat = np.exp(res.x)
ll_weib = -res.fun

lam_exp = t.sum() / d.sum()            # exponential MLE, closed form
ll_exp = np.sum(d * np.log(1 / lam_exp) - t / lam_exp)

LR = 2 * (ll_weib - ll_exp)
print(f"Weibull : k={k_hat:.4f}, lambda={lam_hat:.3f}, loglik={ll_weib:.4f}")
print(f"Exp     : lambda={lam_exp:.3f}, loglik={ll_exp:.4f}")
print(f"LR = {LR:.4f}, p = {stats.chi2.sf(LR, 1):.3e}")

Weibull : k=1.5960, lambda=49.300, loglik=-1088.6557
Exp     : lambda=52.155, loglik=-1124.6066
LR = 71.9018, p = 2.262e-17

정리하며¶

절단자료의 가능도는 두 종류의 기여로 이루어진다.

• 사건이면 밀도, 절단이면 생존함수를 곱한다. 사건 관측은 "정확히 \(t_i\) 에 일어났다"이고, 절단 관측은 "\(t_i\) 를 넘겼다" 는 정보만 준다.
• 로그를 취하면 \(\sum[\delta_i\log h(t_i)-H(t_i)]\) 가 된다. 위험함수와 누적위험으로 적히며, 모든 모수 모형에 공통인 형태다.
• 절단 관측이 버려지지 않는다. 정보의 양이 적을 뿐 기여는 한다. 이것이 생존분석이 절단을 다루는 방식이다.
• 무작위 절단 가정이 여기에 들어 있다. 절단 시점이 사건시간과 독립이어야 이 가능도가 옳다.
• 6장의 최대가능도 이론이 그대로 적용된다. 수치 최적화, 피셔 정보량, 표준오차 모두 같은 방식이다.

다음 절 모수적 생존 모형 (코드) 로 넘어간다.