R 생존분석 완벽 실전 가이드: 카플란-마이어(Kaplan-Meier) 곡선, 로그순위검정, 콕스 비례위험 회귀분석(Cox PH) 총정리

0. 생존분석(Survival Analysis)이란 무엇인가?

데이터 분석에서 “특정 사건(Event)이 일어났는가?”만을 다루는 로지스틱 회귀분석이나 분류 모델과 달리, “그 사건이 일어날 때까지 시간이 얼마나 걸렸는가?”를 분석해야 하는 상황이 매우 빈번하게 발생합니다.

이처럼 관심 있는 사건이 발생할 때까지 경과된 시간(Time-to-Event)을 통계적으로 모델링하고 분석하는 방법론을 생존분석(Survival Analysis) 또는 신뢰성 분석(Reliability Analysis)이라고 합니다.

0-1. 생존분석이 광범위하게 응용되는 도메인

의료 및 임상 연구: 암 환자의 수술/항암 치료 후 생존 기간, 특정 질환의 재발(Recurrence)까지 걸리는 시간.
산업 공학 및 신뢰성 공학: 기계 장비 및 자동차 부품이 마모되어 고장(Failure) 날 때까지의 수명, 시스템 장애 후 복구 시간(MTTR).
비즈니스 및 고객 마케팅: 신규 고객이 가입한 후 첫 구매 또는 서비스 해지(Churn)까지 걸리는 기간.
사회 과학 및 경제학: 실직자가 구직 활동 후 새로운 일자리를 얻을 때까지의 실업 기간.
재난 및 안전 방재: 건축물 준공 후 화재 발생까지의 주기, 특정 단층대에서 지진이 재발생하기까지의 시간 간격.

1. 생존분석의 핵심 개념: 중도절단과 두 가지 기본 함수

일반적인 선형 회귀분석(Linear Regression)을 생존 시간 데이터에 직접 적용할 수 없는 가장 결정적인 이유는 오른쪽 왜도(Right-skewed) 분포와 생존 데이터 특유의 결측값인 중도절단(Censoring) 때문입니다.

1-1. 중도절단(Censoring)의 원인과 유형

중도절단 데이터(Censored Data)는 관찰 기간 동안 대상자에게서 관심 사건이 끝내 발생하지 않은 경우를 의미합니다.

추적 실패 (Loss to Follow-up): 환자가 이사를 가거나 연락이 두절되어 관찰 불가.
연구 조기 중단 (Drop out): 환자가 부작용이나 개인 변심으로 임상시험 참여 중단.
연구 종료 (Termination of Study): 정해진 연구 기간이 끝났으나 여전히 생존 중인 경우.
무관한 원인으로 인한 사건 발생: 분석 대상 질환이 아닌 교통사고 등 무관한 사고로 관찰 종료.

⚠️ 중도절단 데이터 처리의 절대 원칙: 중도절단 데이터를 “미완료된 결측치”로 보고 단순히 제거(Drop)하면, 상대적으로 오래 생존한 개체들이 데이터셋에서 대거 누락되어 생존 기간이 심각하게 과소평가(Underestimation)되는 치명적인 편향(Bias)이 발생합니다.

1-2. R에서 Surv(time, event) 객체 생성

R의 표준 생존분석 라이브러리인 survival 패키지는 Surv() 함수를 통해 생존 시간(time)과 사건 발생 여부(event / status)를 결합한 특수 객체를 생성합니다.

library(survival)
library(dplyr)

# 8명의 가상 환자 데이터셋 생성
df_sample <- tibble(
  id     = 1:8,
  time   = c(1, 1, 4, 5, 6, 9, 9, 22),
  status = c(1, 1, 1, 1, 0, 1, 0, 1) # 1 = 사망(사건), 0 = 중도절단
)

# Surv 객체 생성
surv_obj <- Surv(time = df_sample$time, event = df_sample$status)
print(surv_obj)
# [1]  1   1   4   5   6+  9   9+ 22

출력 결과에서 6+, 9+와 같이 숫자 뒤에 붙은 + 기호가 바로 중도절단(Censored)된 관측치를 나타냅니다. 이는 “최소 6일(혹은 9일)까지는 사건 없이 확실히 생존했다”는 정보를 담고 있습니다.

1-3. 생존 함수 S(t) 와 위험 함수 h(t)

생존분석은 상호 보완적인 두 가지 핵심 수학적 함수로 현상을 설명합니다:

생존 함수 (Survival Function, $S(t)$):
- 특정 시점 $t$를 넘어서까지 생존(사건이 발생하지 않음)할 확률입니다. \[\hat{S}(t) = P(T > t)\]
- $S(0) = 1$ (시작 시점에는 100% 생존)이며, $t o \infty$ 일 때 $S(t) o 0$으로 단조 감소합니다.
위험 함수 (Hazard Function, $h(t)$):
- 시점 $t$까지 생존한 개체가 바로 그 순간($\Delta t o 0$)에 사건을 겪을 조건부 순간 위험률(Instantaneous Rate)입니다. \[h(t) = \lim_{\Delta t o 0} rac{P(t < T \le t + \Delta t \mid T > t)}{\Delta t}\]
- 생존 함수와 위험 함수는 누적 위험 함수 $H(t) = \int_0^t h(u) du$를 매개로 $S(t) = \exp(-H(t))$의 관계를 가집니다.

2. 카플란-마이어(Kaplan-Meier) 분석 및 로그순위 검정

카플란-마이어(Kaplan-Meier, KM) 추정법은 중도절단 정보를 엄밀히 반영하여 생존 함수 $S(t)$를 비모수적(Non-parametric)으로 추정하는 표준 기법입니다.

2-1. 카플란-마이어 누적 생존율 계산 공식

사건이 발생한 각 시점 $t_i$에서, 해당 시점 직전까지 생존해 있던 위험 집단(Number at Risk) $n_i$ 중 실제 사건(사망)이 발생한 수 $d_i$를 이용해 조건부 생존 비율 $p_i$를 구하고, 이를 순차적으로 곱합니다:

\[p_i = 1 - rac{d_i}{n_i}\] \[\hat{S}(t) = \prod_{t_i \le t} \left(1 - rac{d_i}{n_i} ight) = \hat{S}(t-1) imes \left(1 - rac{d_i}{n_i} ight)\]

2-2. R 실습: lung 데이터셋을 활용한 KM 모델링

survival 패키지에 내장된 228명의 진행성 폐암 환자 임상 데이터인 lung 데이터셋을 활용합니다.

library(survival)
library(survminer)
library(dplyr)

# 1. 데이터 로드 및 팩터 변환
data("lung", package = "survival")

lung <- lung %>%
  mutate(
    sex = factor(sex, levels = c(1, 2), labels = c("Male", "Female")),
    event = (status == 2) # status: 1=중도절단, 2=사망
  )

# 2. 전체 환자 대상 KM 모델 적합 (~ 1)
km_all <- survfit(Surv(time, event) ~ 1, data = lung)
print(km_all)

출력 결과 해석: * n = 228: 전체 환자 수 * events = 165: 관찰 기간 내 사망 환자 수 * median = 310: 중앙 생존 기간(Median Survival Time). 환자의 50%가 생존해 있는 시점이 310일임을 의미합니다. * 0.95LCL / 0.95UCL: 중앙 생존 기간의 95% 신뢰구간 (284일 ~ 363일).

2-3. 집단 간 생존 곡선 비교 및 로그순위 검정 (Log-rank Test)

성별(Male vs Female)에 따라 생존 기간에 통계적으로 유의한 차이가 존재하는지 검정합니다.

# 성별에 따른 KM 모델 적합
km_sex <- survfit(Surv(time, event) ~ sex, data = lung)

# 로그순위 검정 실행
logrank_res <- survdiff(Surv(time, event) ~ sex, data = lung)
print(logrank_res)

Call:
survdiff(formula = Surv(time, event) ~ sex, data = lung)

            N Observed Expected (O-E)^2/E (O-E)^2/V
sex=Male   138      112     91.6      4.55      10.3
sex=Female  90       53     73.4      5.67      10.3

 Chisq= 10.3  on 1 degrees of freedom, p= 0.00131

통계적 해석: 카이제곱 통계량 $\chi^2 = 10.3$ ($df=1$), $p = 0.00131 < 0.05$로 귀무가설($H_0$: 두 집단의 생존 곡선은 동일하다)을 기각합니다. 즉, 남성과 여성 환자 간 생존 기간에 통계적으로 매우 유의한 차이가 존재합니다.

2-4. survminer::ggsurvplot() 기반 출판용 고품질 시각화

survminer 패키지의 ggsurvplot() 함수를 사용하면 신뢰구간, 위험 테이블(Risk Table), 생존 중앙값 가이드라인, 로그순위 검정 p-value가 통합된 저널 게재용 시각화를 완성할 수 있습니다.

# 출판용 Kaplan-Meier 생존 곡선 플롯
ggsurvplot(
  km_sex,
  data = lung,
  pval = TRUE,                 # Log-rank 검정 p-value 자동 표시
  conf.int = TRUE,             # 95% 신뢰구간 리본 표시
  risk.table = "abs_pct",      # 하단 위험 집단 수 및 백분율 테이블
  risk.table.col = "strata",   # 집단별 색상 일치
  surv.median.line = "hv",     # 생존 중앙값(50%) 수직/수평 가이드라인
  legend.title = "성별",
  legend.labs = c("남성 (Male)", "여성 (Female)"),
  palette = c("#2563eb", "#e11d48"),
  xlab = "추적 관찰 기간 (일)",
  ylab = "누적 생존 확률 S(t)",
  ggtheme = theme_bw()
)

2-5. 재난 방재 분야 활용: 이벤트 누적 발생 확률 곡선

ggsurvplot(..., fun = "event") 옵션을 지정하면 생존 확률 대신 누적 사건 발생 확률($1 - S(t)$)을 시각화할 수 있어, 지역별 화재·지진 등 재난 발생 위험도 분석에 탁월하게 활용됩니다.

3. 콕스 비례위험 회귀분석 (Cox Proportional Hazards Regression)

카플란-마이어 분석은 1개의 범주형 변수에 대해서만 그룹별 생존 곡선을 비교할 수 있으며, 환자의 나이, 혈압, 종양 크기 등 연속형 변수를 통제할 수 없습니다.

콕스 비례위험 회귀모형(Cox Proportional Hazards Model)은 여러 개의 예측 변수(연속형 및 범주형)가 생존 시간에 미치는 영향을 동시에 통제하고 평가하는 준모수적(Semi-parametric) 다변량 생존분석 모델입니다.

3-1. 콕스 회귀모형의 수학적 구조

\[h_i(t) = h_0(t) \exp(eta_1 X_{i1} + eta_2 X_{i2} + \cdots + eta_p X_{ip})\]

양변에 자연로그를 취하면: \[\ln\left(rac{h_i(t)}{h_0(t)} ight) = eta_1 X_{i1} + eta_2 X_{i2} + \cdots + eta_p X_{ip}\]

$h_i(t)$: $i$번째 개체의 시점 $t$에서의 위험률(Hazard Rate).
$h_0(t)$: 기저 위험 함수(Baseline Hazard). 모든 공변량 $X$의 값이 0일 때의 위험률로, 특정 확률분포(예: 와이블, 지수분포)를 가정하지 않는 비모수적 성격을 가집니다.
$\exp(\sum eta_j X_j)$: 공변량들의 상대적 위험 기여도를 나타내는 모수적 성격.

3-2. 위험비 (Hazard Ratio, HR)의 완벽 해석

공변량 $X_1$이 1단위 증가할 때의 위험비(HR)는 다음과 같이 정의됩니다:

\[ ext{HR} = rac{h(t \mid X_1 = k + 1)}{h(t \mid X_1 = k)} = rac{h_0(t) e^{eta_1(k+1) + \cdots}}{h_0(t) e^{eta_1 k + \cdots}} = e^{eta_1}\]

$ ext{HR} = 1$ ($eta = 0$): 해당 변수는 사건 발생 위험에 아무런 영향을 미치지 않음.
$ ext{HR} > 1$ ($eta > 0$): 위험 요인 (Risk Factor). 변수값이 증가할수록 위험률이 증가하고 생존 기간이 단축됨. (예: $ ext{HR}=1.58 o$ 위험 58% 증가)
$ ext{HR} < 1$ ($eta < 0$): 보호 요인 (Protective Factor). 변수값이 증가할수록 위험률이 감소하고 생존 기간이 연장됨. (예: $ ext{HR}=0.57 o$ 위험 43% 감소)

3-3. R 실습: coxph() 다변량 모델 적합 및 summary 해석

나이(age), 성별(sex), 신체활동수행점수(ph.ecog, 0=양호 ~ 3=불량)를 공변량으로 투입합니다.

# 1. 결측치 정제 및 콕스 모델 적합
lung_clean <- lung %>% filter(!is.na(ph.ecog))

cox_fit <- coxph(Surv(time, event) ~ age + sex + ph.ecog, data = lung_clean)
summary(cox_fit)

Call:
coxph(formula = Surv(time, event) ~ age + sex + ph.ecog, data = lung_clean)

  n= 227, number of events= 164 

              coef exp(coef) se(coef)      z Pr(>|z|)    
age       0.01107   1.01113  0.00927  1.194  0.23242    
sexFemale -0.55861  0.57199  0.16723 -3.340  0.00084 ***
ph.ecog   0.46309   1.58897  0.11347  4.081 4.48e-05 ***
---
            exp(coef) exp(-coef) lower .95 upper .95
age            1.0111     0.9890    0.9929    1.0297
sexFemale      0.5720     1.7483    0.4121    0.7938
ph.ecog        1.5890     0.6293    1.2721    1.9848

Concordance= 0.637  (se = 0.025 )
Likelihood ratio test= 30.5  on 3 df,   p=1.08e-06
Wald test            = 29.93 on 3 df,   p=1.43e-06
Score (logrank) test = 30.5  on 3 df,   p=1.08e-06

핵심 해석: 1. sexFemale ($ ext{HR} = 0.572, p = 0.00084$): 나이와 ECOG 점수를 통제했을 때, 여성 환자는 남성 환자에 비해 사망 위험이 약 42.8% 유의하게 낮습니다 ($1 - 0.572 = 0.428$). 2. ph.ecog ($ ext{HR} = 1.589, p < 0.0001$)**: ECOG 점수가 1점 나빠질 때마다 사망 위험이 약 **58.9% 유의하게 증가합니다**. 3. **`age` ($ ext{HR} = 1.011, p = 0.232$): 다른 요인을 통제한 상태에서 나이의 단독 영향은 $p > 0.05$로 통계적으로 유의하지 않았습니다.

3-4. survminer::ggforest() 포레스트 플롯 시각화

# 콕스 회귀모형 위험비 포레스트 플롯
ggforest(cox_fit, data = lung_clean)

3-5. 콕스 모델의 핵심 전제: 비례위험 가정 검정 (cox.zph)

콕스 회귀모형이 타당하려면 “시간이 흘러도 집단 간의 위험비(HR) 비율이 일정하게 유지되어야 한다”는 비례위험 가정(Proportional Hazards Assumption)을 반드시 만족해야 합니다.

# 비례위험 가정 검정
ph_check <- cox.zph(cox_fit)
print(ph_check)

# Schoenfeld 잔차 플롯
ggcoxzph(ph_check)

cox.zph() 결과 모든 변수의 $p > 0.05$이면 비례위험 가정을 만족하는 것으로 판단합니다.
Schoenfeld 잔차 플롯에서 평활선이 수평($y=0$)을 유지하면 시간 의존적 효과가 없음을 시각적으로 입증합니다.

4. [심화] 다요인 생존분석 및 반복 사건·프레일티 모델

4-1. 다요인 생존 곡선 분할 시각화 (colon 데이터셋)

치료법(rx), 성별(sex), 종양 분화도(differ) 등 3개 이상의 범주형 요인이 복합적으로 작용할 때는 ggsurvplot 객체에 facet_grid를 결합하여 다차원 패널 분할 그래프를 구성합니다.

data("colon", package = "survival")

colon_death <- colon %>%
  filter(etype == 2) %>% # 사망 사건
  mutate(
    rx     = factor(rx, labels = c("Obs", "Lev", "Lev+5FU")),
    sex    = factor(sex, labels = c("Male", "Female")),
    differ = factor(differ, labels = c("Well", "Mod", "Poor")),
    event  = (status == 1)
  )

fit_multi <- survfit(Surv(time, event) ~ rx + sex + differ, data = colon_death)
gg_plot <- ggsurvplot(fit_multi, data = colon_death, conf.int = FALSE)

# facet_grid로 분화도 x 치료법 패널 분할
gg_plot$plot +
  facet_grid(differ ~ rx) +
  theme_bw() +
  labs(title = "Colon Cancer Survival: Stratified by Differentiation and Treatment")

4-2. 반복 사건(Recurrent Events) 및 프레일티 모델 (Frailty Model)

환자 1명에게 질병이 여러 번 재발하거나, 동일 시설에서 화재가 반복 발생하는 경우 데이터 내 개체 간 상관성(Intra-cluster Correlation)이 발생합니다.

# 1. 클러스터 강건 표준오차 콕스 모델
cox_recurrent <- coxph(
  Surv(start, stop, event) ~ treatment + cluster(id),
  data = recurrent_data
)

# 2. 프레일티 모델 (랜덤 효과를 통한 개체 간 이질성 모델링)
# library(frailtypack)
# frailty_fit <- frailtyPenal(Surv(time, event) ~ treatment + cluster(id), data = recurrent_data)

5. 카플란-마이어 분석 vs 콕스 회귀분석 비교 요약

비교 항목	카플란-마이어 (Kaplan-Meier) 분석	콕스 비례위험 회귀분석 (Cox PH Model)
모델링 유형	비모수적 (Non-parametric)	준모수적 (Semi-parametric)
분석 목적	단일 집단/요인별 생존 곡선 추정 및 비교	다변량 공변량의 상대적 위험 기여도 평가
설명변수(X) 형태	범주형 요인(Group)만 가능	연속형 및 범주형 변수 동시 투입 가능
핵심 출력 지표	누적 생존율 $\hat{S}(t)$, 중앙 생존 기간	위험비(Hazard Ratio, $ ext{HR} = e^eta$)
통계적 가설검정	로그순위 검정 (Log-rank Test)	Wald Test, 우도비 검정 (Likelihood Ratio Test)
핵심 통계적 가정	독립적 중도절단(Independent Censoring)	독립적 중도절단 + 비례위험 가정(PH Assumption)

6. 실무자를 위한 핵심 체크리스트 요약

중도절단 보존: 관찰 기간 내 사건이 발생하지 않은 데이터를 제거하지 않고 Surv(time, event)의 우측 중도절단으로 정확히 코딩했는가?
생존 대표값 선정: 비대칭 오른쪽 왜도 분포를 고려하여 산술 평균 대신 중앙 생존 기간(Median Survival Time)을 보고했는가?
로그순위 검정 적용: 2개 이상 집단 간 생존 곡선의 통계적 차이를 $p < 0.05$ 기준으로 올바르게 판정했는가?
위험비(HR) 직관적 해석: $ ext{HR} > 1$(위험 요인)과 $ ext{HR} < 1$(보호 요인)의 백분율 위험 증감 수치를 정확히 설명했는가?
비례위험 가정 검토: 콕스 회귀분석 적합 후 반드시 cox.zph()를 실행하여 비례위험 가정이 충족되었는지 검증했는가?

📚 참고 문헌 및 공식 출처

구분	자료명 및 공식 사이트	주요 내용	바로가기
R 공식 재단	CRAN: survival - Survival Analysis in R	Terry Therneau 박사의 R 표준 생존분석 패키지 공식 도큐멘테이션	CRAN 공식 포털
R 시각화 패키지	survminer: Drawing Survival Curves using ggplot2	Alboukadel Kassambara의 출판용 생존 곡선 및 포레스트 플롯 공식 가이드	공식 문서
학술 교재	Survival Analysis: A Self-Learning Text (3rd ed.)	David G. Kleinbaum & Mitchel Klein 저, 생존분석 및 콕스 모형 바이블	Springer 도서 정보
고급 통계 서적	Regression Modeling Strategies (with R)	Frank E. Harrell Jr. 교수의 다변량 생존분석 및 모형 검증 전략	공식 연구 사이트
임상 통계 저널	Statistical Analysis of Survival Data in Clinical Trials	생존분석의 중도절단 처리와 비례위험 가정 검정 실무 가이드라인	NCBI 논문 전문