이 챕터를 마치면 다음을 할 수 있습니다:
ggsurvfit 패키지로 제작metafor)이 챕터에서는 R 내장 데이터셋과 시뮬레이션 데이터를 사용합니다.
R 패키지의 데이터를 사용하려면 반드시 패키지를 먼저 불러와야 합니다:
library(tidyverse)
library(survival) # ⚠️ 먼저 패키지 로드!
library(ggsurvfit)
library(metafor)
# 생존 분석: survival 패키지 내장 데이터
data(lung) # 폐암 환자 생존 데이터 (N=228)
data(colon) # 대장암 생존 데이터
# 메타분석: metafor 패키지 내장 데이터
data(dat.bcg) # BCG 백신 효과 메타분석 (13개 연구)
# 임상시험 데이터 (시뮬레이션)
clinical <- read_csv(here("data", "processed", "clinical_trial.csv"))출처: North Central Cancer Treatment Group 환자 수: 228명 (진행성 폐암 환자) 주요 변수:
| 변수명 | 설명 | 타입 |
|---|---|---|
time |
생존 시간 (일) | numeric |
status |
생존 상태 (1=사망, 2=생존) | numeric |
age |
나이 (세) | numeric |
sex |
성별 (1=남성, 2=여성) | numeric |
ph.ecog |
ECOG 수행 점수 (0-3) | numeric |
ph.karno |
Karnofsky 점수 (의사 평가) | numeric |
pat.karno |
Karnofsky 점수 (환자 평가) | numeric |
meal.cal |
식사 칼로리 | numeric |
wt.loss |
최근 6개월 체중 감소 (파운드) | numeric |
사용 예시: Kaplan-Meier 생존 곡선, Cox 비례위험 모델
참고: ?lung 명령어로 자세한 설명을 볼 수 있습니다.
생존 분석(Survival Analysis)은 사건(event)이 발생할 때까지의 시간을 분석하는 통계 기법입니다.
주요 용어:
예시:
환자 A: 진단 → 5년 추적 → 사망 (event=1, time=5)
환자 B: 진단 → 3년 추적 → 탈락 (event=0, time=3, censored)
환자 C: 진단 → 10년 추적 → 생존 (event=0, time=10, censored)중도절단(censoring) 때문입니다!
문제: - 환자 B는 3년 후 탈락했지만, 실제로는 5년째 사망했을 수 있음 - 이를 무시하고 “3년 생존”으로 분석하면 과소추정
해결책: - 생존 분석은 중도절단을 정보(information)로 활용 - “최소한 3년은 생존했다”는 정보를 반영
Kaplan-Meier 추정량은 생존 확률을 시간에 따라 계산합니다.
S(t) = Π [1 - (사건 발생 수 / 위험 인구)]R 내장 데이터로 실습:
# ⚠️ 중요: survival 패키지를 먼저 로드해야 lung 데이터 사용 가능!
library(survival) # 1단계: 패키지 로드
library(ggsurvfit)
# 2단계: lung 데이터 불러오기
data(lung) # survival 패키지의 내장 데이터
# 데이터 구조 확인
head(lung, 3)#> inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
#> 1 3 306 2 74 1 1 90 100 1175 NA
#> 2 3 455 2 68 1 0 90 90 1225 15
#> 3 3 1010 1 56 1 0 90 90 NA 15예시: - 환자 1: time=306일, status=2 → 306일 후 사망 - 환자 2: time=455일, status=1 → 455일까지 추적 후 중도절단 (생존)
TIP: library(survival) 후에는 data(lung) 없이 바로 lung을 사용해도 됩니다!
# 생존 객체 생성
# Surv(time, status): time=생존 시간(일), status=1(사망), 0(censored)
surv_obj <- Surv(lung$time, lung$status)
# 전체 환자의 생존 곡선
fit_all <- survfit(surv_obj ~ 1, data = lung)
# 기본 생존 곡선
ggsurvfit(fit_all, linewidth = 1) +
labs(title = "전체 폐암 환자 생존 곡선",
x = "시간 (일)",
y = "생존 확률") +
theme_minimal()
# 성별(sex)로 구분: 1=남성, 2=여성
fit_sex <- survfit2(Surv(time, status) ~ sex, data = lung)
# ggsurvfit로 시각화
ggsurvfit(fit_sex, linewidth = 1) +
add_confidence_interval(alpha = 0.2) + # 95% CI 추가
labs(title = "폐암 환자 생존 곡선 (성별)",
x = "시간 (일)",
y = "생존 확률",
color = "성별") +
scale_color_manual(values = c("1" = "#2196F3", "2" = "#E91E63"),
labels = c("1" = "남성", "2" = "여성")) +
scale_fill_manual(values = c("1" = "#2196F3", "2" = "#E91E63"),
labels = c("1" = "남성", "2" = "여성")) +
theme_minimal()
ggsurvfit(fit_sex, linewidth = 1) +
add_confidence_interval(alpha = 0.2) +
add_risktable(
risktable_stats = c("n.risk", "cum.event"), # 위험 인구, 누적 사건
size = 3.5
) +
labs(title = "폐암 환자 생존 곡선 (위험 테이블 포함)",
x = "시간 (일)",
y = "생존 확률") +
scale_color_manual(values = c("1" = "#2196F3", "2" = "#E91E63"),
labels = c("1" = "남성", "2" = "여성")) +
scale_fill_manual(values = c("1" = "#2196F3", "2" = "#E91E63"),
labels = c("1" = "남성", "2" = "여성")) +
scale_ggsurvfit() + # 축 스케일 조정
theme_minimal()
두 생존 곡선이 통계적으로 다른지 검정합니다.
가설: - H0: 두 그룹의 생존 곡선이 동일하다 - H1: 두 그룹의 생존 곡선이 다르다
#> Call:
#> survdiff(formula = Surv(time, status) ~ sex, data = lung)
#>
#> N Observed Expected (O-E)^2/E (O-E)^2/V
#> sex=1 138 112 91.6 4.55 10.3
#> sex=2 90 53 73.4 5.68 10.3
#>
#> Chisq= 10.3 on 1 degrees of freedom, p= 0.001결과 해석: - p-value < 0.05: 두 그룹의 생존에 유의한 차이가 있음 - p-value ≥ 0.05: 두 그룹 간 차이가 없음
# p-value를 그래프에 추가
ggsurvfit(fit_sex, linewidth = 1) +
add_confidence_interval(alpha = 0.2) +
add_pvalue(location = "annotation", size = 5) + # p-value 추가
labs(title = "폐암 환자 생존 곡선 (Log-rank p-value 포함)",
x = "시간 (일)",
y = "생존 확률") +
scale_color_manual(values = c("1" = "#2196F3", "2" = "#E91E63"),
labels = c("1" = "남성", "2" = "여성")) +
scale_fill_manual(values = c("1" = "#2196F3", "2" = "#E91E63"),
labels = c("1" = "남성", "2" = "여성")) +
theme_minimal()
생존 확률이 50%가 되는 시점입니다.
# 중앙 생존 기간 추출
summary(fit_sex)$table#> records n.max n.start events rmean se(rmean) median 0.95LCL 0.95UCL
#> sex=1 138 138 138 112 326.0841 22.91156 270 212 310
#> sex=2 90 90 90 53 460.6473 34.68985 426 348 550ggsurvfit(fit_sex, linewidth = 1) +
add_confidence_interval(alpha = 0.2) +
add_quantile(y_value = 0.5, linewidth = 0.8, linetype = "dashed") + # 50% 라인
labs(title = "폐암 환자 생존 곡선 (중앙 생존 기간)",
x = "시간 (일)",
y = "생존 확률") +
scale_color_manual(values = c("1" = "#2196F3", "2" = "#E91E63"),
labels = c("1" = "남성", "2" = "여성")) +
scale_fill_manual(values = c("1" = "#2196F3", "2" = "#E91E63"),
labels = c("1" = "남성", "2" = "여성")) +
theme_minimal()
다변량 분석으로 위험비(Hazard Ratio)를 계산합니다.
# Cox 모델: 성별, 나이, ECOG 수행 상태 보정
cox_model <- coxph(Surv(time, status) ~ sex + age + ph.ecog, data = lung)
summary(cox_model)#> Call:
#> coxph(formula = Surv(time, status) ~ sex + age + ph.ecog, data = lung)
#>
#> n= 227, number of events= 164
#> (1 observation deleted due to missingness)
#>
#> coef exp(coef) se(coef) z Pr(>|z|)
#> sex -0.552612 0.575445 0.167739 -3.294 0.000986 ***
#> age 0.011067 1.011128 0.009267 1.194 0.232416
#> ph.ecog 0.463728 1.589991 0.113577 4.083 4.45e-05 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> exp(coef) exp(-coef) lower .95 upper .95
#> sex 0.5754 1.7378 0.4142 0.7994
#> age 1.0111 0.9890 0.9929 1.0297
#> ph.ecog 1.5900 0.6289 1.2727 1.9864
#>
#> Concordance= 0.637 (se = 0.025 )
#> Likelihood ratio test= 30.5 on 3 df, p=1e-06
#> Wald test = 29.93 on 3 df, p=1e-06
#> Score (logrank) test = 30.5 on 3 df, p=1e-06Hazard Ratio (HR) 해석: - HR = 1: 위험 동일 - HR > 1: 위험 증가 (나쁜 예후) - HR < 1: 위험 감소 (좋은 예후)
예시: - sex HR = 0.59: 여성이 남성보다 사망 위험 41% 낮음 - age HR = 1.02: 나이 1세 증가마다 사망 위험 2% 증가
메타분석(Meta-Analysis)은 여러 연구 결과를 통합하여 하나의 추정값을 도출하는 통계 기법입니다.
목적: 1. 검정력 증가: 개별 연구보다 큰 표본 크기 2. 정밀도 향상: 더 좁은 신뢰구간 3. 일관성 평가: 연구 간 이질성(heterogeneity) 탐색 4. 근거 종합: 의사결정을 위한 최선의 추정
언제 사용하나요? - 여러 RCT 결과 통합 - 체계적 문헌고찰 (Systematic Review) - 임상 진료 가이드라인 개발 - 정책 의사결정
고정 효과 모델 (Fixed-Effect Model) - 가정: 모든 연구가 동일한 참값(true effect)을 추정 - 사용: 연구 간 이질성이 낮을 때 - 가중치: 큰 연구에 더 많은 weight
랜덤 효과 모델 (Random-Effects Model) - 가정: 각 연구의 참값이 정규분포를 따름 - 사용: 연구 간 이질성이 있을 때 - 가중치: 작은 연구에도 더 많은 weight
권장: 대부분의 경우 랜덤 효과 모델 사용
메타분석 결과를 시각화하는 표준 방법입니다.
library(metafor)
# BCG 백신 메타분석 데이터
data(dat.bcg)
# 효과 크기 계산 (Risk Ratio)
dat <- escalc(measure = "RR",
ai = tpos, bi = tneg, # 치료군
ci = cpos, di = cneg, # 대조군
data = dat.bcg,
slab = paste(author, year, sep = ", "))
# 랜덤 효과 모델
res <- rma(yi, vi, data = dat)
# Forest plot
forest(res,
atransf = exp, # log RR → RR로 변환
at = log(c(0.05, 0.25, 1, 4)),
xlim = c(-16, 6),
header = "Author(s) and Year",
xlab = "Risk Ratio (95% CI)",
mlab = "Random-Effects Model",
refline = 1) # RR=1 기준선
Forest Plot 읽는 법:
연구 간 차이가 얼마나 큰지 평가합니다.
# 이질성 통계량
print(res)#>
#> Random-Effects Model (k = 13; tau^2 estimator: REML)
#>
#> tau^2 (estimated amount of total heterogeneity): 0.3132 (SE = 0.1664)
#> tau (square root of estimated tau^2 value): 0.5597
#> I^2 (total heterogeneity / total variability): 92.22%
#> H^2 (total variability / sampling variability): 12.86
#>
#> Test for Heterogeneity:
#> Q(df = 12) = 152.2330, p-val < .0001
#>
#> Model Results:
#>
#> estimate se zval pval ci.lb ci.ub
#> -0.7145 0.1798 -3.9744 <.0001 -1.0669 -0.3622 ***
#>
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1주요 지표:
| 지표 | 의미 | 해석 |
|---|---|---|
| Q 통계량 | 연구 간 변동 | p < 0.05: 이질성 존재 |
| I² | 전체 변동 중 이질성 비율 | <25%: 낮음, 25-75%: 중간, >75%: 높음 |
| τ² | 연구 간 분산 | 클수록 이질성 큼 |
예시 해석:
I² = 85% → 전체 변동의 85%가 이질성 때문
→ 연구 간 차이가 큼 → 하위그룹 분석 필요출판 편향: 통계적으로 유의한 결과만 출판되어 메타분석이 왜곡되는 현상
탐지 방법: Funnel Plot
# Funnel plot
funnel(res,
xlab = "Log Risk Ratio",
ylab = "Standard Error",
main = "Funnel Plot (BCG 백신 메타분석)")
Funnel Plot 해석:
Egger’s Test:
# Egger's regression test
regtest(res)#>
#> Regression Test for Funnel Plot Asymmetry
#>
#> Model: mixed-effects meta-regression model
#> Predictor: standard error
#>
#> Test for Funnel Plot Asymmetry: z = -0.8033, p = 0.4218
#> Limit Estimate (as sei -> 0): b = -0.5104 (CI: -1.1182, 0.0974)metafor의 기본 forest plot은 base R 그래픽입니다. ggplot2로 더 예쁘게 만들 수 있습니다.
library(tidyverse)
# 메타분석 결과 데이터프레임 변환
forest_data <- tibble(
study = dat$slab,
estimate = exp(dat$yi), # RR
lower_ci = exp(dat$yi - 1.96 * sqrt(dat$vi)),
upper_ci = exp(dat$yi + 1.96 * sqrt(dat$vi)),
weight = weights(res)
)
# 통합 추정값 추가
pooled_estimate <- tibble(
study = "Pooled (Random-Effects)",
estimate = exp(res$beta[1]),
lower_ci = exp(res$ci.lb),
upper_ci = exp(res$ci.ub),
weight = 100
)
forest_data <- bind_rows(forest_data, pooled_estimate)
# ggplot2 Forest Plot
ggplot(forest_data, aes(y = reorder(study, estimate), x = estimate,
xmin = lower_ci, xmax = upper_ci)) +
geom_pointrange(aes(size = weight)) +
geom_vline(xintercept = 1, linetype = "dashed", color = "red") +
scale_x_log10(breaks = c(0.1, 0.25, 0.5, 1, 2, 4)) +
labs(title = "BCG 백신 효과 메타분석",
subtitle = "Risk Ratio (95% CI)",
x = "Risk Ratio (log scale)",
y = "Study") +
theme_minimal() +
theme(legend.position = "none")
ROC 곡선은 진단 검사의 성능을 평가하는 도구입니다.
주요 개념:
| 용어 | 정의 | 예시 |
|---|---|---|
| True Positive (TP) | 질병 O, 검사 양성 | 암 환자를 양성으로 진단 |
| False Positive (FP) | 질병 X, 검사 양성 | 정상인을 양성으로 오진 |
| True Negative (TN) | 질병 X, 검사 음성 | 정상인을 음성으로 진단 |
| False Negative (FN) | 질병 O, 검사 음성 | 암 환자를 음성으로 오진 |
성능 지표:
민감도 (Sensitivity) = TP / (TP + FN) # 실제 양성 중 양성으로 진단
특이도 (Specificity) = TN / (TN + FP) # 실제 음성 중 음성으로 진단ROC 곡선: - x축: 1 - 특이도 (False Positive Rate) - y축: 민감도 (True Positive Rate) - 모든 임계값(threshold)에서의 성능을 시각화
AUC: ROC 곡선 아래 면적
| AUC 값 | 해석 |
|---|---|
| 1.0 | 완벽한 진단 |
| 0.9 - 1.0 | 우수 (Excellent) |
| 0.8 - 0.9 | 좋음 (Good) |
| 0.7 - 0.8 | 보통 (Fair) |
| 0.5 - 0.7 | 나쁨 (Poor) |
| 0.5 | 무작위 추측 |
library(pROC)
library(tidyverse)
# 시뮬레이션 데이터 생성
set.seed(123)
n <- 200
disease_status <- rbinom(n, 1, 0.3) # 30% 유병률
test_score <- ifelse(disease_status == 1,
rnorm(n, mean = 70, sd = 15), # 환자군
rnorm(n, mean = 50, sd = 15)) # 정상군
# ROC 객체 생성
roc_obj <- roc(disease_status, test_score)
# ROC 곡선
ggroc(roc_obj, linewidth = 1.5) +
geom_abline(intercept = 1, slope = 1,
linetype = "dashed", color = "gray") + # 무작위 추측선
annotate("text", x = 0.3, y = 0.3,
label = paste0("AUC = ", round(auc(roc_obj), 3)),
size = 5, color = "darkblue") +
labs(title = "ROC 곡선 (진단 검사 성능)",
x = "1 - 특이도 (False Positive Rate)",
y = "민감도 (True Positive Rate)") +
theme_minimal()
| 통계량 | 의미 | 해석 | 시각화 용도 |
|---|---|---|---|
| SE (Standard Error) | 평균의 정밀도 | 표본 평균이 모평균과 얼마나 가까운지 | 평균의 불확실성 |
| 95% CI | 모평균에 대한 신뢰 | 95% 확률로 모평균이 이 구간에 포함 | 통계적 유의성 검정 |
| SD (Standard Deviation) | 데이터의 산포 | 개별 관측값이 평균에서 얼마나 퍼져있는지 | 데이터 분포 |
문제: SE를 CI로 착각하거나 SD를 SE로 혼동
올바른 사용: - 데이터 분포 표현: SD 사용 - 평균 비교: 95% CI 사용 - 평균의 정밀도: SE 사용
임상시험/메타분석: 거의 항상 95% CI 사용!
library(tidyverse)
# 시뮬레이션 데이터
set.seed(456)
treatment_data <- tibble(
group = rep(c("A", "B", "C"), each = 30),
value = c(rnorm(30, 70, 15),
rnorm(30, 80, 15),
rnorm(30, 75, 15))
)
# 요약 통계
summary_stats <- treatment_data %>%
group_by(group) %>%
summarize(
n = n(),
mean = mean(value),
sd = sd(value),
se = sd / sqrt(n),
ci_lower = mean - 1.96 * se,
ci_upper = mean + 1.96 * se
)
# 시각화
p1 <- ggplot(summary_stats, aes(x = group, y = mean)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.2) +
labs(title = "평균 ± SE", y = "값") +
theme_minimal()
p2 <- ggplot(summary_stats, aes(x = group, y = mean)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.2) +
labs(title = "평균 ± 95% CI", y = "값") +
theme_minimal()
p3 <- ggplot(summary_stats, aes(x = group, y = mean)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width = 0.2) +
labs(title = "평균 ± SD", y = "값") +
theme_minimal()
library(patchwork)
p1 + p2 + p3 +
plot_annotation(title = "SE vs CI vs SD 비교")
관찰: - SE: 가장 좁음 (평균의 정밀도) - 95% CI: SE의 1.96배 - SD: 가장 넓음 (데이터 분포)
✅ 생존 분석 - Kaplan-Meier 생존 곡선 - 위험 테이블 (Risk Table) - Log-rank 검정 - 중앙 생존 기간 - Cox 비례위험 모델
✅ 메타분석 - 고정 효과 vs. 랜덤 효과 - Forest Plot 시각화 - 이질성 평가 (I², τ²) - Funnel Plot (출판 편향) - ggplot2 커스터마이징
✅ ROC 곡선 - 민감도 vs. 특이도 - AUC 해석 - 최적 임계값
✅ 통계 지표 - SE vs. CI vs. SD 차이 - 올바른 사용법
생존 곡선:
fit <- survfit2(Surv(time, status) ~ group, data)
ggsurvfit(fit) +
add_confidence_interval() +
add_risktable() +
add_pvalue()Forest Plot:
ROC 곡선:
roc_obj <- roc(disease, test_score)
ggroc(roc_obj) +
annotate("text", label = paste0("AUC = ", round(auc(roc_obj), 3)))Chapter 6에서는 학술지 출판을 위한 전문적인 그래프 제작 기법을 배웁니다:
Nature, Science, NEJM 수준의 그래프를 만들 수 있습니다!