adjustedCurves: Estimating Confounder-Adjusted Survival Curves in R

기술 요약: adjustedCurves R 패키지

문제 정의
우측 검열(right-censored)된 사건 발생 시간 종점(time-to-event endpoints)을 다루는 연구에서, 표준적인 접근 방식은 처치 할당에 따라 층화된 카플란-마이어(Kaplan-Meier, KM) 추정량을 시각화하는 것이다. 이러한 방법은 공변량이 균형 잡힌 무작위 대조 시험(RCT)에서는 편향되지 않은 추정치를 제공하지만, 관찰 연구나 순응도(adherence) 문제가 있는 RCT에서는 실패한다. 이러한 시나리오에서 교란(confounding)은 묘사된 생존 곡선을 체계적으로 왜곡하며, 이는 종종 조정된 회귀 분석(예: Cox 모델)과 조정되지 않은 그래픽 표현 사이의 불일치를 초래한다. 더욱이, 회귀 모델은 일반적으로 상대적 효과 척도인 위험비(hazard ratios)를 제공하는데, 이는 인과적 해석력이 부족하며 CONSORT 및 STROBE 가이드라인에서 권장하는 절대적 효과 척도(예: 생존 확률 또는 평균 생존 시간)로 쉽게 변환되지 않는다. 기존의 교란 조정 생존 곡선 소프트웨어 구현체들은 단 하나 또는 두 개의 방법만을 지원하거나, 일관된 프레임워크가 부족하거나, 신뢰 구간 및 출판 가능한 수준의 플로팅 루틴을 제공하지 못하는 등 파편화되어 있다.

방법론
본 논문은 다양한 문헌 기반 방법론을 하나의 일관된 프레임워크 아래 통합한 adjustedCurves R 패키지를 소개한다. 이 패키지는 모든 인구가 처치 $z$를 받았을 경우의 반사실적 생존 확률인 한계 생존 확률(marginal survival probability, $S_z(t)$)을 추정하는 것을 목표로 한다.

패키지는 조건부 교환 가능성(unmeasured confounding이 없다는 가정) 및 기타 표준적인 인과적 식별 가능성 가정을 전제로 다음과 같은 추정 전략을 구현한다:

1. 층화(Stratification): 교란 층 내에서의 가중치 적용 KM 추정치 풀링 (방법: strat_cupples, strat_nieto, strat_amato).
2. G-Computation (직접 표준화): 사건 발생 시간 결과(예: Cox 모델)를 모델링하고 공변량 분포에 대해 예측값을 평균화함 (방법: direct).
3. 역확률 가중치(Inverse Probability of Treatment Weighting, IPTW): 성향 점수의 역수를 사용하여 교란 요인을 균형 있게 조정함 (방법: iptw_km, iptw_cox).
4. 성향 점수 매칭(Propensity Score Matching): 성향 점수에 기반하여 개별 사례를 매칭함으로써 KM 추정을 위한 균형 잡힌 표본을 생성함 (방법: matching).
5. 경험적 가능도(Empirical Likelihood, EL): 모수적 가정 없이 교란 요인을 균형 있게 맞추기 위해 제약된 가능도 극대화를 통해 가중치를 도출함 (방법: emp_lik).
6. 증강 역확률 가중치(Augmented Inverse Probability of Treatment Weighting, AIPTW): 결과 모델링과 처치 가중법을 결합한 이중 강건(doubly robust) 접근법 (방법: aiptw).
7. 의사값(Pseudo-Values, PV): 생존 데이터를 변환하여 일반적인 회귀 기법(예: GEE)이 검열을 처리할 수 있도록 함. IPTW, G-Computation, AIPTW에 적용 가능 (방법: iptw_pseudo, direct_pseudo, aiptw_pseudo).
8. 도구 변수(Instrumental Variable, IV) 추정: 특정 조건 하에서 측정되지 않은 교란을 해결하기 위한 2단계 잔차 포함법(2SRI-F) (방법: iv_2SRIF).
9. 근사 인과 추론(Proximal Causal Inference): 측정되지 않은 교란이 존재하는 상황에서 처치 유도 및 결과 유도 대리 변수를 활용하여 인과 효과를 복구하는 방법 (방법: prox_iptw, prox_aiptw).

또한, 이 패키지는 이러한 곡수로부터 도출되는 한계 처치 효과의 추정을 지원한다:

• 특정 시점에서의 생존 확률 차이.
• 생존 시간 분위수 (예: 중앙 생존 시간).
• 제한 평균 생존 시간 (Restricted Mean Survival Times, RMST).

주요 기여 및 특징

• 통합 프레임워크: 서로 다른 구문을 가진 고립된 방법들을 구현하는 기존 R 패키지들(예: survminer, riskRegression, RISCA)과 달리, adjustedCurves는 모든 주요 조정 방법에 대해 고수준의 사용자 친화적 인터페이스(adjustedsurv())를 제공한다.
• 통계적 추론: 점별 신뢰 구간(point-wise confidence intervals)을 점근적 근사(사용 가능한 경우) 또는 계산 집약적인 부트스트래핑을 통해 계산한다.
• 시각화: 신뢰 밴드, 위험 표(risk tables), 중앙 생존 지표를 포함하여 출판 가능한 그래픽을 생성할 수 있는 ggplot2 기반의 고도로 커스터마이징 가능한 plot() 함수를 포함한다.
• 결측 데이터 처리: 다중 대체 데이터(mids 객체)에 대한 네이티브 지원을 통해 결측 공변량이 존재하는 상황에서도 강건한 분석을 가능하게 한다.
• 사후 추정 보정: AIPTW 및 PV 기반 방법에서 발생하는 알려진 문제인, 생존 곡선이 비감소(non-increasing) 상태를 유지하고 [0, 1] 확률 경계 내에 머물도록 하기 위해 이소토닉 회귀(isotioninc regression) 및 경계 절단(boundary truncation)을 구현한다.
• 경쟁 사건(Competing Events): 교란이 조정된 원인 특이 누적 발생 함수(cause-specific cumulative incidence functions)를 추정하기 위한 adjustedcif() 함수를 포함한다.

결과 및 예시 사례
저자들은 호르몬 요법의 효과를 평가하기 위해 로테르담 종양 은행 데이터(유방암 환자 2,982명)를 사용하여 패키지를 시연한다.

• 미조정 분석: 표준 KM 곡선은 호르몬 요법을 받는 환자가 그렇지 않은 환자에 비해 생존 확률이 현저히 낮음을 시사했다.
• 조정 분석: IPTW, G-Computation, AIPTW를 사용하여 교란 요인(연령, 림프절, 종양 크기, 등급, 폐경 여부)을 조정한 후, 두 그룹의 생존 곡선은 거의 동일해졌으며 신뢰 구간이 겹치는 것으로 나타났다.
• 효과 척도: 교란이 해결된 후, 5, 6, 7년 시점에서의 생존 확률 차이, 중앙 생존 시간, RMST 계산 결과 처치 그룹 간에 통계적으로 유의미한 차이가 없었다.
  이 예시는 미조정 곡선이 교란(예: 더 아픈 환자가 더 많은 처치를 받는 경우)으로 인해 어떻게 오해를 불러일으킬 수 있는지, 그리고 이 패키지가 어떻게 이러한 편향을 교정하여 진정한 인과 효과를 드러내는지 보여준다.

의의 및 주장
본 논문은 adjustedCurves가 고급 사용자를 위한 유연성을 유지하면서도 R 경험이 적은 연구자들이 교란 조정 생존 곡량 추정을 쉽게 사용할 수 있도록 함으로써 통계 소프트웨어의 중요한 공백을 메운다고 주장한다. 다양한 방법론에 대한 일관된 프레임워크를 제공함으로써, 이 패키지는 관찰 연구에서 조정된 회귀 결과와 조정되지 않은 그래픽 표현 사이의 불일치를 줄이는 것을 목표로 한다. 저자들은 이 패키지가 위험비보다 더 해석하기 쉽고 인과적으로 타당한 절대적 효과 척도(RMST 및 생존 차이와 같은)의 보고를 가능하게 함으로써 현재의 보고 가이드라인과 일치한다고 강조한다.

저자들은 한계점도 언급하였다: 현재 이 패키지는 기저(baseline) 교란 변수만을 처리하며(시간 가변적 교변량은 불가), 연속형 변수는 범주형으로 분류해야 하며, CRAN에서 의존성 패키지를 사용할 수 없어 TMLE(Targeted Maximum Likelihood Estimation)는 제외된다는 점이다. 그러나 이 패키지는 관찰 환경에서 인과적 생존 효과를 추정하기 위한 포괄적인 도구 역할을 한다.