Parameter Estimation for Time-Scaled Inhomogeneous Phase-Type Distributions from Discrete Observations

기술 요약: 이산 관측치를 이용한 시간 척도 변환 불균질 단계형 분포의 매개변수 추정

1. 문제 정의

고전적인 단계형(Phase-Type, PH) 분포는 유한 상태의 균질 마르코프 점프 과정(Markov Jump Process, MJP)의 흡수 시간까지의 시간을 모델링한다. 불균질 단계형(Inhomogeneous Phase-Type, IPH) 분포는 전이 강도가 시간에 따라 변하도록 허용함으로써 헤비 테일(heavy-tailed) 특성, 비상수 위험률(non-constant hazard rates), 그리고 시간 의존적 역학을 포착할 수 있도록 이 프레임워크를 확장한다. 구체적이고 다루기 쉬운 특정 하위 클래스인 시간 척도 변환 IPH 분포는 $\Lambda(t) = h_\beta(t)\Lambda$ 형태의 시간 척도 변환된 부강도 행렬(sub-intensity matrix)을 가정하며, 여기서 $\Lambda$는 상수 기저 행렬이고 $h_\beta(t)$는 $\beta$에 의해 매개변수화된 알려진 비음수 스케일링 함수이다.

이러한 모델에 대한 통계적 추론은 흡수 시간이나 전체 궤적이 정확하게 관찰되는 경우를 위해 개발되었으나, 많은 실제 응용 분야에서는 이산적이고 불규칙하게 간격이 배치된 관측치를 포함한다. 이러한 설정(예: 의료 모니터링에서의 다중 상태 패널 데이터)에서는 정확한 전이 시간과 흡수 시간이 관찰되지 않으므로, 문제는 결측 데이터(missing-data) 문제로 귀결된다. 직접적인 최대 가능도 추정(Maximum Likelihood Estimation, MLE)은 복잡한 행렬 지수 함수를 포함하며, $\Lambda$와 $\beta$의 원소들에 대한 제약 조건이 있는 비선형 최적화를 요구하기 때문에 어렵다. 본 논문은 기초 과정이 이산적인 시간 지점에서만 관찰되는 경우, 시간 척도 변환 IPH 분포의 매개변수 $(\pi, \Lambda, \beta)$를 추정하는 문제를 다룬다.

2. 방법론

저자들은 마르코프 브리지(Markov-bridge) 데이터 증강과 확률적 기대-최대화(Stochastic Expectation-Maximization, SEM) 알고리즘을 결합한 통계적 추론 프레임워크를 제안한다. 핵심 전략은 시간 변환을 통해 불균질 문제를 균질 문제로 변환하고, 잠재적인 연속 시간 궤적을 시뮬레이션한 후, 생성된 완전 데이터 가능도(complete-data likelihood)를 기반으로 매개변수를 업데이트하는 것이다.

2.1 시간 변환

$\Lambda(t) = h_\beta(t)\Lambda$인 하위 클래스의 경우, 불균질 과정 $X$는 시간 변환 $s = g_\beta^{-1}(t) = \int_0^t h_\beta(u) du$를 통해 균질 MJP $Y$로 매핑될 수 있다. 이 변환 하에서, 불균질 과정의 흡수 시간 $\tau$는 $\tau \stackrel{d}{=} g_\beta(\rho)$를 만족하며, 여기서 $\rho$는 매개변수 $(\pi, \Lambda)$를 갖는 표준 균질 PH 분포를 따른다.

2.2 SEM 알고리즘

추정 절차는 누락된 데이터를 시뮬레이션하는 단계와 완전 데이터 가능도를 최대화하는 단계를 교대로 수행한다:

1. 초기화:

   • 초기 분포 $\pi$는 시작 상태의 경험적 빈도를 통해 추정된다.
   • $\beta$에 대한 초기 추측값 $\hat{\beta}_0$가 제공된다.
   • 초기 부강도 행렬 $\hat{\Lambda}_0$는 이산 관측치를 (초기에 시간 척도 변환을 무시하고) 균질 과정으로부터의 연속 관측치로 취급하여 얻는다.
   • 초기 $\hat{\beta}$의 정교화는 (시간 척도 변환 없이) 시뮬레이션된 경로로부터 도출된 완전 데이터 로그 가능도에 대한 경사 상승법(gradient ascent)을 사용하여 수행된다.

2. 반복 절차 (SE-단계 및 M-단계):

   • SE-단계 (Stochastic E-step): 현재의 추정치 $\hat{\beta}$를 사용하여, 관측된 이산 시간 지점들을 균질 시간 척도로 변환한다. 관측 지점들 사이의 잠재적인 연속 시간 궤적은 관측된 상태를 조건으로 하는 마르코프 브리지를 사용하여 시뮬레이션된다. 아직 흡수되지 않은 경로의 경우, 궤적은 흡수 상태까지 확장된다. 이러한 시뮬레이션된 경로는 "완전 데이터" 실현체가 된다.
   • M-단계 (Maximization):
     • $\Lambda$ 업데이트: 시뮬레이션된 균질 경로가 주어졌을 때, 부강도 행렬 $\Lambda$는 균질 MJP에 대한 폐쇄형 최대 가능도 추정량(MLE)(전이 횟수와 점유 시간을 계산)을 사용하여 업데이트된다.
     • $\beta$ 업데이트: 시뮬레이션된 경로는 다시 불균질 시간 척도로 매핑된다. 매개변수 $\beta$는 경사 상승법을 사용하여 $\beta$에 대해 완전 데이터 로그 가능도를 최대화함으로써 정교화된다. 이 단계는 점프 시간과 스케일링된 점유 시간에 의존하는 로그 가능도의 기울기를 평가하는 과정을 포함한다.

3. 종료: 단일 업데이트 후 내부 경사 상승법이 수렴하면 알고리즘이 종료되며, 이는 추정치가 안정되었음을 나타낸다.

2.3 식별 가능성 (Identifiability)

본 논문은 매개변수 식별 가능성을 위한 조건을 확립한다. 구체적으로, 매개변수화 $(\Lambda, \beta)$는 $\beta_1 \neq \beta_2$일 때 스케일링 함수 $h_{\beta_1}(t)$와 $h_{\beta_2}(t)$가 곱셈 상수 이상으로 다를 때만 식별 가능하다. 이 조건은 기저 행렬 $\Lambda$와 시간 스케일링 함수가 서로 혼동되지 않도록 보장한다.

3. 주요 결과

제안된 방법의 성능은 시뮬레이션 연구와 실제 데이터 적용을 통해 평가되었다.

3.1 시뮬레이션 연구

• 행렬-곰페르츠(Matrix-Gompertz) 모델: $h_\beta(t) = e^{\beta t}$를 사용하여, 저자들은 3개의 전이 상태를 가진 1,000개의 경로를 시뮬레이션하였다. 이 방법은 다양한 관측 지평($T$)에 걸쳐 실제 매개변수 $(\pi, \Lambda, \beta)$를 성공적으로 회복하였다. 흡수 시간에 대한 콜모고로프-스미르노프(Kolmogorov-Smirnov, KS) 검정으로 측정된 적합도는 $T$가 클 때 높게 유지되었으나, 충분한 흡수 이벤트를 포착하기에 $T$가 너무 짧을 경우(예: $T=36$은 낮은 p-값을 보임) 유의미하게 저하되었다.
• 행렬-와이불(Matrix-Weibull) 모델: $h_\beta(t) = \beta t^{\beta-1}$를 사용하여, 2개의 전이 상태를 가진 1,000개의 경로에 대해 테스트하였다. 알고리즘은 23번의 반복 후에 정확한 추정치(예: 실제 $\beta=3$, 추정 $\hat{\beta} \approx 2.97$)로 수렴하였다. KS 검정은 p-값이 0.5004로 나타나, 실제 흡수 시간 분포를 재현하는 모델의 능력을 확인하였다.

3.2 실제 데이터 적용

본 방법은 관상동맥 이식편 질환(Coronary Allograft Vasculopathy, CAV) 진행을 모니터링하는 622명의 심장 이식 수혜자 데이터셋에 적용되었다. 데이터는 환자의 상태(CAV 미발생, 경증, 중등도/중증)와 사망(흡수 상태)에 대한 불규칙하게 간격이 배치된 관측치로 구성되었다.

• 모델 적합: 곰페르츠 스케일링 함수($h_\beta(t) = e^{\beta t}$)를 갖는 시간 척도 변환 IPH 모델이 피팅되었다. 추정된 매개변수는 전이 강도가 시간에 따라 지수적으로 증가함을 나타냈다.
• 비교: 피팅된 불균질 모델은 균질 PH 모델과 비교되었다. 불균질 모델은 균질 모델(KS p-값 = 0.01066)보다 유의미하게 더 나은 적합도(KS p-값 = 0.5966)를 제공하였으며, 이는 이 맥락에서 시간 가변적 강도를 고려하는 것의 필요성을 입증하였다.
• 해석: 결과는 중등도/중증 상태에 있는 환자들이 사망에 대한 순간 위험도가 가장 높음을 강조하였으며, 모델은 질병 진행과 관련된 가속화되는 위험 역학을 포착하였다.

4. 의의 및 주장

본 논문은 불규칙한 다중 상태 데이터에 시간 척도 변환 IPH 모델을 피팅하기 위한 유연하고 계산적으로 다루기 쉬운 도구를 제공한다고 주장한다. 주요 기여는 다음과 같다:

• 결측 데이터 처리: 이 프레임워크는 마르코프 브리지를 통해 전이 및 흡수 시간을 재구성함으로써 unobserved(관찰되지 않은) 전이 및 흡수 시간의 문제를 효과적으로 해결한다.
• 제약 조건이 있는 최적화 회피: 완전 데이터 업데이트를 활용함으로써, 관측 데이터 로그 가능도의 어려운 제약 조건이 있는 비선형 최적화를 피한다. $\Lambda$에 대한 업데이트는 폐쇄형(closed-form)이며, $\beta$에 대한 업데이트는 1차원 경사 상승법이다.
• 이론적 토대: 저자들은 시간 척도 하위 클래스에 대한 식별 가능성 조건을 확립하여, 기저 행렬과 스케일링 함수 사이의 매개변수 혼동을 배제하였다.
• 실용적 유용성: CAV 적용을 통해, 제안된 SEM 기반 추론 방법이 신뢰할 수 있는 매개변수 추정치를 제공하며, 불규칙한 관측이 존재하는 실제 환경에서 균질 모델보다 우수한 성능을 보이며 시간적 진화 양상을 정확하게 포착함을 입증하였다.

저자들은 보고된 추정량들이 시뮬레이션된 완전 데이터 실현치에서 평가된 최대 가능도 추정량임을 강조한다. 표준 SEM 조건 하에서, 추정량 시퀀스는 관측 데이터 극대값 근처에 집중하는 에르고딕 마르코프 체인을 형성하며, 몬테카를로 변동성(Monte Carlo variability)을 제외하면 최대 가능도 추정량으로서의 해석을 유지한다.