우리는 ROC 곡선에서 AUC 값을 구할 수 있다. AUC 값이 1에 가까우면 모델의 성능이 좋다라고 하며, 0.5에 가까워질 수록 성능이 나쁘다- 라고 한다. 그렇다면, ROC에서 신뢰구간은 어떻게 구할까? 특정 AUC값과 비교하여 p value를 구하려면 어떻게 해야 할까?
이를 위해 DeLong의 AUC 표준 오차 구하는 방법을 소개하려 한다. (E. DeLong, 1988)
민감도(Sensitivity, True Positive Rate)
특이도(Specificity, True Negative Rate)
\[ Spec = \frac{True Negatives(TN)}{True\ Negatives(TN) + False \ Positives(FP)} \]
이상적인 ROC 곡선
AUC(Area Under the Curve)
\[ \text{for any real number z,} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \]
\[ Sens(z) = \frac{1}{m}\sum_{i=1}^{m}{I(X_i \geq z)}, \]
\[ Spec(z) = \frac{1}{n}\sum_{j=1}^{n}{I(Y_i < z)}, \]
\[ I(A) = \begin{cases} 1, & A \text{ is true} \\ 0, & \text{otherwise} \end{cases} \]
여기서, Mann-Whitney two sample statistic에 따르면, ROC curve 아래 넓이를 구할때, trapezodial rule로 구한 넓이는 Mann-Whitney two sample statistic으로 구한 넓이로 대체할 수 있다.
\[ \hat{\theta} = \frac{1}{mn}\sum_{j=1}^{n}\sum_{i=1}^{m}{\psi(X_i, Y_j)} \]
\[ \psi(X,Y) = \begin{cases} 1 & \ Y < X \\ 1/2 & \ Y=X \\ 0 & \ Y>X \end{cases} \]
\[ E(\hat{\theta}) = Pr(X>Y) + \frac{1}{2}Pr(X= Y) \]
이제, AUC값의 신뢰성을 측정하기 위해서는 SE(standard Error)를 계산하는 것이 필요하다. 이는 DeLong(1988)이 제시한 방법을 참고하자.
ξ를 각각의 집단 간 공분산이라 하자.
ξ₁₀은 \(C_1\)의 \(X_i\)와 \(C_2\)의 \(Y_j, Y_k\)간 공분산,
ξ₀₁은 \(C_2\)의 \(Y_j\)와 \(C_1\)의 \(X_i, Y_k\)간 공분산,
ξ₁₁은 \(C_1\)의 \(X_i\)와 \(C_2\)의 \(Y_j\)간 자기공분산이다.
ξ₁₀, ξ₀₁, ξ₁₁의 기대값은 아래와 같다.
\[ ξ_{10} = E[\psi(X_i, Y_j) \psi(X_i,Y_k)]-\theta^2, \]
\[ ξ_{01} = E[\psi(X_i, Y_j) \psi(X_k,Y_j)]-\theta^2, \]
\[ ξ_{11} = E[\psi(X_i, Y_j) \psi(X_i,Y_j)]-\theta^2, \]
기댓값들을 이용하여 AUC 추정치의 분산을 계산할 수 있다.
\[ var(\hat{\theta}) = \frac{(n-1)ξ_{10}+(m-1)ξ_{01}}{mn} + \frac{ξ_{11}}{mn} \]
이와 같은 방법으로 단일 표본 집합에서의 AUC의 표준 오차를 구할 수 있다.
\[ ξ_{10}^{rs} = E[\psi(X_i^r, Y_j^r) \psi(X_i^s,Y_k^s)]-\theta^r \theta^s, \]
\[ ξ_{01}^{rs} = E[\psi(X_i^r, Y_j^r) \psi(X_i^s,Y_k^s)]-\theta^r \theta^s, \]
\[ ξ_{11}^{rs} = E[\psi(X_i^r, Y_j^r) \psi(X_i^s,Y_k^s)]-\theta^r \theta^s, \]
이제, 아래 식을 통해 표본 집합의 AUC 값 간의 공분산을 계산할 수 있다.
\[ cov(\hat{\theta^r},\hat{\theta^s}) = \frac{(n-1)ξ_{10}^{rs}+(m-1)ξ_{01}^{rs}}{mn} + \frac{ξ_{11}^{rs}}{mn} \] 이 수식을 통해 여러 표본 집합 간의 공분산을 반영하여 AUC의 표준 오차를 더 정확하게 추정할 수 있다.
\[ V^r_{10}(X_i) = \frac{1}{n} \sum_{j=1}^n\psi(X_i^r,Y_j^r)\ \ \ \ \ (i= 1,2,...,m) \]
\[ V^r_{01}(Y_j) = \frac{1}{m} \sum_{i=1}^m\psi(X_i^r,Y_j^r)\ \ \ \ \ (j= 1,2,...,m) \]
\(\psi(X_i^r,Y_j^r)\)는 \(X_i^r,Y_j^r\)의 관계를 나타내며, 이를 통해 분산을 구한다. \(V^r_{01}(Y_j)\)은 집합 \(Y^r_j\)에서의 분산이다.
두 값을 통해 각 표본 집합에 대한 분산을 따로 계산한 후, 이를 결합하여 최종적으로 표준 오차를 추정할 수 있다. 각 분산 값이 표본 크기에 따라 가중평균되며, AUC 표준 오차 S는 아래와 같다.
\[ \\ S = \frac{1}{m}S_{10} + \frac{1}{n}S_{01} \\ \]
표준 오차를 통해 우리는 AUC 값이 얼마나 신뢰할 수 있는지 평가할 수 있으며, 이는 모델의 성능을 명확히 이해하는데 중요한 역할을 한다.
이를 R에서는 pROC package를 통해 이용할 수 있다.
library(pROC);library(dplyr);library(data.table)
data(aSAH)
SAH <- as.data.table(aSAH)
head(SAH) gos6 outcome gender age wfns s100b ndka
<ord> <fctr> <fctr> <int> <ord> <num> <num>
1: 5 Good Female 42 1 0.13 3.01
2: 5 Good Female 37 1 0.14 8.54
3: 5 Good Female 42 1 0.10 8.09
4: 5 Good Female 27 1 0.04 10.42
5: 1 Poor Female 42 3 0.13 17.40
6: 1 Poor Male 48 2 0.10 12.75데이터는 pROC 패키지 내장 데이터인 aSAH 데이터를 사용하였다. aSAH는 뇌동맥류 파열로 인해 발생하는 SAH에 대한 데이터이다. n=113명이며, 변수는 아래와 같다.
gos6 (Glasgow Outcome scale),
outcome (clinical outcome: Good/Poor)
gender (성별: Male/female)
age (나이)
wfns (신경학적 평가 점수: 1~5, 1= 가장 양호),
s100b (혈액 내 S100B 단백질 농도: 뇌손상 biomarker)
ndka (혈액 내 neuron-specific enoalse 농도: 신경 손상 biomarker)위 변수 중 100b를 사용하여 outcome을 예측하는 ROC curve와 AUC 값을 구해보자.
roc_result <- roc(response = SAH$outcome,
predictor = SAH$s100b,
levels = c("Good", "Poor"), # outcome의 순서 지정
direction = "<")
print(roc_result)
Call:
roc.default(response = SAH$outcome, predictor = SAH$s100b, levels = c("Good", "Poor"), direction = "<")
Data: SAH$s100b in 72 controls (SAH$outcome Good) < 41 cases (SAH$outcome Poor).
Area under the curve: 0.7314위와 같이 지정한 변수 및 AUC 값을 확인할 수 있다.
여기서 ROC curve를 직접 보고싶으면 아래와 같이 확인할 수 있다.
plot(
roc_result,
col = "blue", # curve 색깔 지정
lwd = 2, # curve 두께 지정
main = "ROC Curve for s100b" #그래프 제목 지정
# xlim = c(1, 0), # x축 범위, 일반적으로 1-spec(FPR)
# ylim = c(0, 1) # y축 범위, 일반적으로 Sensitivity
# type = "1" # 그래프의 형태, 기본값은 선그래프
# xlab/ylab = "1-Specificity" # x/y축 레이블 지정.
)
text( # curve와 함께 AUC값 표기.
x = 0.6, y = 0.4,
labels = paste("AUC =", round(auc(roc_result), 3)),
col = "red"
)
이처럼 ROC curve, AUC값을 구할 수 있다. 이제 AUC 값에 대해 추가적으로 알아보자.
auc_result <- auc(roc_result) #이전에 구한 roc_result 사용
se_auc <- sqrt(var(roc_result))
ci_result <- ci.auc(roc_result)
print(auc_result);print(ci_result);print(se_auc)Area under the curve: 0.731495% CI: 0.6301-0.8326 (DeLong)[1] 0.05165929z_value <- (auc_result - 0.7)/se_auc
p_value <- 2*(1- pnorm(abs(z_value)))
cat("p-value for AUC > 0.7:", p_value, "\n")p-value for AUC > 0.7: 0.5437048 위 경우는 유의하지 않다.
library(pROC);library(data.table);library(dplyr)
#data
data(aSAH)
SAH <- as.data.table(aSAH)
#subgroup
subgroups <- list(
data = SAH,
age_over_50 = SAH[age > 50],
age_50_or_below = SAH[age <= 50],
male = SAH[gender == "Male"],
female = SAH[gender == "Female"]
)
# ROC 및 성능 평가 함수(threshold: 0.8)
evaluate_roc <- function(data, predictor_var) {
#predictor_var는 numeric variable이어야.
if (!is.numeric(data[[predictor_var]])) {
data[[predictor_var]] <- as.numeric(data[[predictor_var]])
}
# outcome variable은 factor variable이어야.
if (!is.factor(data$outcome)) {
data$outcome <- factor(data$outcome, levels = c("Poor", "Good"))
}
#ROC
roc_result <- roc(data$outcome, data[[predictor_var]], levels = c("Poor", "Good"))
#AUC, CI
auc_value <- auc(roc_result) %>% round(3)
auc_CI <- ci.auc(roc_result) %>% round(3)
# Handle AUC == 1
if (auc_value == 1) {
warning("AUC is 1. SE and p-value may be undefined.")
se_auc <- NA
p_value <- NA
} else {
# SE
se_auc <- sqrt(var(roc_result))
# P-value
z_value <- (auc_value - 0.8) / se_auc
p_value <- 2 * (1 - pnorm(abs(z_value)))
}
p_value_display <- ifelse(p_value < 0.001, "<0.001", round(p_value, 3))
# 비교하고픈 cutoff가 있다면 비교하여 Sensitivity, Specificity, PPV, NPV, Accuracy등 다양한 값을 구할 수 있다.
#ex. median으로 cutoff를 지정하여 비교할 수 있다. (결과에는 포함하지 않음.)
#median cutoff.
median_cutoff <- median(data[[predictor_var]], na.rm = TRUE)
#
data$predicted <- data[[predictor_var]] >= median_cutoff
table1 <- table(data$outcome, data$predicted)
# 구조 점검.
row_levels <- rownames(table1)
col_levels <- colnames(table1)
if (!("Good" %in% row_levels)) table1 <- rbind(table1, "Good" = c(0, 0))
if (!("Poor" %in% row_levels)) table1 <- rbind(table1, "Poor" = c(0, 0))
if (!("TRUE" %in% col_levels)) table1 <- cbind(table1, "TRUE" = c(0, 0))
if (!("FALSE" %in% col_levels)) table1 <- cbind(table1, "FALSE" = c(0, 0))
TP <- table1["Good", "TRUE"]
FP <- table1["Poor", "TRUE"]
FN <- table1["Good", "FALSE"]
TN <- table1["Poor", "FALSE"]
Sensitivity <- (TP / (TP + FN)) %>% round(3)
Specificity <- (TN / (TN + FP)) %>% round(3)
PPV <- (TP / (TP + FP)) %>% round(3)
NPV <- (TN / (TN + FN)) %>% round(3)
Accuracy <- ((TP + TN) / (TP + FP + TN + FN)) %>% round(3)
list(
AUC = auc_value,
AUC_CI = auc_CI,
`P value` = p_value_display
)
}
# ROC 분석
predictor_vars <- c("gos6", "age", "wfns", "s100b", "ndka")
results <- lapply(subgroups, function(group) {
lapply(predictor_vars, function(var) evaluate_roc(group, var))
})
result_tables <- rbindlist(lapply(names(results), function(group_name) {
#subgroup 가져오기
group_results <- results[[group_name]]
# 예측 변수 가져오기
group_data <- rbindlist(lapply(seq_along(group_results), function(i) {
# 인덱스(i)
data <- group_results[[i]]
if (length(data) > 0) {
as.data.table(data) %>%
mutate(
Group = group_name,
Predictor = paste("Predictor", i),
`AUC(95% CI)` = paste0(min(AUC_CI), "-", max(AUC_CI)),
`P value` = ifelse(`P value` < 0.001, "<0.001", round(`P value`, 3))
) %>%
select(
Group, Predictor, AUC, `AUC(95% CI)`, `P value`
)
}
}))
return(group_data)
}), use.names = TRUE, fill = TRUE)
result_tables[, Predictor := case_when(
Predictor == "Predictor 1" ~ "gos6",
Predictor == "Predictor 2" ~ "age",
Predictor == "Predictor 3" ~ "wfns",
Predictor == "Predictor 4" ~ "s100b",
Predictor == "Predictor 5" ~ "ndka"
)]
# cause CI, 3 row created. rm.
result_tables_filtered <- result_tables[seq(1, nrow(result_tables), by = 3), ]
print(result_tables_filtered) Group Predictor AUC AUC(95% CI) P value
<char> <char> <num> <char> <char>
1: data gos6 1.000 1-1 <NA>
2: data age 0.615 0.508-0.722 <0.001
3: data wfns 0.824 0.749-0.899 0.531
4: data s100b 0.731 0.63-0.833 0.182
5: data ndka 0.612 0.501-0.723 <0.001
6: age_over_50 gos6 1.000 1-1 <NA>
7: age_over_50 age 0.517 0.36-0.675 <0.001
8: age_over_50 wfns 0.738 0.613-0.864 0.334
9: age_over_50 s100b 0.725 0.588-0.861 0.28
10: age_over_50 ndka 0.626 0.473-0.778 0.025
11: age_50_or_below gos6 1.000 1-1 <NA>
12: age_50_or_below age 0.541 0.387-0.695 <0.001
13: age_50_or_below wfns 0.910 0.837-0.983 0.003
14: age_50_or_below s100b 0.702 0.529-0.874 0.266
15: age_50_or_below ndka 0.612 0.453-0.771 0.021
16: male gos6 1.000 1-1 <NA>
17: male age 0.680 0.514-0.845 0.154
18: male wfns 0.876 0.773-0.979 0.148
19: male s100b 0.773 0.632-0.914 0.707
20: male ndka 0.552 0.371-0.734 0.007
21: female gos6 1.000 1-1 <NA>
22: female age 0.635 0.49-0.78 0.026
23: female wfns 0.779 0.67-0.887 0.705
24: female s100b 0.720 0.57-0.87 0.296
25: female ndka 0.667 0.526-0.808 0.064
Group Predictor AUC AUC(95% CI) P value데이터 상 크게 유의미한 변수는 없다.
[1] E. DeLong, D. DeLong, and D. Clarke-Pearson, “Comparing the areasunder two or more correlated receiver operating characteristic curves: anonparametric approach,” Biometrics, pp. 837–845, 1988
@online{lee2024,
author = {LEE, Hojun},
title = {DeLong’s {Method;} for {ROC} {AUC}},
date = {2024-11-20},
url = {https://blog.zarathu.com/posts/2024-11-20-DeLongsMethod/},
langid = {en}
}