Introduction

우리는 ROC 곡선에서 AUC 값을 구할 수 있다. AUC 값이 1에 가까우면 모델의 성능이 좋다라고 하며, 0.5에 가까워질 수록 성능이 나쁘다- 라고 한다. 그렇다면, ROC에서 신뢰구간은 어떻게 구할까? 특정 AUC값과 비교하여 p value를 구하려면 어떻게 해야 할까?

이를 위해 DeLong의 AUC 표준 오차 구하는 방법을 소개하려 한다. (E. DeLong, 1988)

ROC curve (Recevier Operating Characteristic Curve)

ROC 곡선은 이진 분류 문제에서 모델의 성능을 평가하는데 사용되는 시각적 도구이다.
다양한 임계값에서의 민감도와 False Positive Rate의 관계를 시각화한 것으로,
민감도(Sensitivity)를 Y축으로, 1 - 특이도(False Positive Rate)을 X축으로 하여 관계를 그린 그래프이다.

민감도(Sensitivity, True Positive Rate)

실제 양성 사례를 얼마나 잘 분류하는지/ 실제 양성인 샘플을 양성으로 올바르게 분류한 비율 $S e n s = \frac{T r u e P o s i t i v e s (T P)}{T r u e P o s i t i v e s (T P) + F a l s e N e g a t i v e s (F N)}$

특이도(Specificity, True Negative Rate)

실제 음성 사례를 얼마나 잘 분류하는지/ 실제 음성인 샘플을 음성으로 올바르게 분류한 비율
False Positive Rate(FPR) : 1-Spec, 실제 음성인 샘플을 잘못 양성으로 분류한 비율

$S p e c = \frac{T r u e N e g a t i v e s (T N)}{T r u e N e g a t i v e s (T N) + F a l s e P o s i t i v e s (F P)}$

이상적인 ROC 곡선

이상적인 분류는 ROC가 (0,1)을 지날 때, 즉 FPR이 0이고, TPR이 1인 경우.
무작위 분류는 ROC가 y=x일 때, TPR = FPR인 경우.

AUC(Area Under the Curve)

AUC는 ROC 곡선 아래 면적을 의미한다.
AUC가 1에 가까울수록 성능이 좋은 모델. 0.5에 가까울수록 성능이 안좋은 모델.

Figure 1: ROC

Empirical ROC Curve

아래와 같이 ROC Curve를 추정할 수 있다.
test를 시행한 총 집단의 수를 N이라 할 때, 실제로 이벤트가 발생한 집단을 $C_{1}$ ( $X_{i}$ ,n=m), 발생하지 않은 집단을 $C 2$ ( $Y_{i}$ ,n=n)이라 하자.

$for any real number z,$

$S e n s (z) = \frac{1}{m} \sum_{i = 1}^{m} I (X_{i} \geq z),$

$S p e c (z) = \frac{1}{n} \sum_{j = 1}^{n} I (Y_{i} < z),$

$I (A) = {\begin{cases} 1, & A is true \\ 0, & otherwise \end{cases}$

이때 실수 z가 variable 내 가능한 값들 내에서 움직인다면, ROC curve는 [1 - spec(z)]를 X로, Sens(z)를 Y로 갖는 plot이라 할 수 있다. 만약 z가 가능한 최대값보다 크다면 curve는 (0,0)을, 최솟값보다 작다면(1,1)을 지날 것이다. Sens(z) = 1 - Spec(z)라면 y=x 위, 45도 선 위에 놓일 것이다.

AUC

ROC curve는 위와 같이 구할 수 있다. 그럼 이의 넓이: AUC 값은 어떻게 구할까? 보통, 곡선 아래 넓이는 trapezoidal rule을 통해 구한다. 고등학교 때 배운 적분을 떠올리면 된다. 수많은 사다리꼴로 쪼개어 넓이를 근사하던 기억을 되살려 보자.

여기서, Mann-Whitney two sample statistic에 따르면, ROC curve 아래 넓이를 구할때, trapezodial rule로 구한 넓이는 Mann-Whitney two sample statistic으로 구한 넓이로 대체할 수 있다.

Mann-Whitney statistic는 확률 $θ$ 를 예측한다. $C_{2}$ 에서 무작위 추출한 값이 $C_{1}$ 의 값보다 같거나 작을 확률을 추정한다. (test를 시행한 총 집단의 수를 N이라 할 때, 실제로 이벤트가 발생한 집단을 $C_{1}$ ( $X_{i}$ ,n=m), 발생하지 않은 집단을 $C 2$ ( $Y_{i}$ ,n=n)이라 하자)

$\hat{θ} = \frac{1}{m n} \sum_{j = 1}^{n} \sum_{i = 1}^{m} ψ (X_{i}, Y_{j})$

$ψ (X, Y) = {\begin{cases} 1 & Y < X \\ 1 / 2 & Y = X \\ 0 & Y > X \end{cases}$

모든 (X, Y) 쌍에 대해 X > Y이면 1을, X = Y이면 $\frac{1}{2}$ , X < Y 이면 0을 부여하여 확률을 구한다.
직관적으로, ROC curve와 AUC 값은 곧 모델의 추정이 옳을 확률이며, 이의 성능의 최고값은 1, 최저값은 1/2이라는 점을 고민하면 위의 추정은 그럴듯 하다.

$E (\hat{θ}) = P r (X > Y) + \frac{1}{2} P r (X = Y)$

그럼, 확률(AUC)은 위와 같이 정리할 수 있다.

Standard Error

이제, 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 [ψ (X_{i}, Y_{j}) ψ (X_{i}, Y_{k})] - θ^{2},$

$ξ_{01} = E [ψ (X_{i}, Y_{j}) ψ (X_{k}, Y_{j})] - θ^{2},$

$ξ_{11} = E [ψ (X_{i}, Y_{j}) ψ (X_{i}, Y_{j})] - θ^{2},$

기댓값들을 이용하여 AUC 추정치의 분산을 계산할 수 있다.

$v a r (\hat{θ}) = \frac{(n - 1) ξ_{10} + (m - 1) ξ_{01}}{m n} + \frac{ξ_{11}}{m n}$

이와 같은 방법으로 단일 표본 집합에서의 AUC의 표준 오차를 구할 수 있다.

이제, 단일 표본 집합이 아닌 다른 표본집합 r과 s에 대해 다뤄보자. 여러 표본 집합이 있을 경우, 각 표본 간의 상호 공분산 또한 고려되어야한다. 두 표본 집합 r과 s에 대해 AUC의 공분산 계산은 아래와 같다.

$ξ_{10}^{r s} = E [ψ (X_{i}^{r}, Y_{j}^{r}) ψ (X_{i}^{s}, Y_{k}^{s})] - θ^{r} θ^{s},$

$ξ_{01}^{r s} = E [ψ (X_{i}^{r}, Y_{j}^{r}) ψ (X_{i}^{s}, Y_{k}^{s})] - θ^{r} θ^{s},$

$ξ_{11}^{r s} = E [ψ (X_{i}^{r}, Y_{j}^{r}) ψ (X_{i}^{s}, Y_{k}^{s})] - θ^{r} θ^{s},$

이제, 아래 식을 통해 표본 집합의 AUC 값 간의 공분산을 계산할 수 있다.

$c o v (\hat{θ^{r}}, \hat{θ^{s}}) = \frac{(n - 1) ξ_{10}^{r s} + (m - 1) ξ_{01}^{r s}}{m n} + \frac{ξ_{11}^{r s}}{m n}$ 이 수식을 통해 여러 표본 집합 간의 공분산을 반영하여 AUC의 표준 오차를 더 정확하게 추정할 수 있다.

이를 바탕으로 우리가 궁금한 값 “표준 오차”에 접근해 보자. (Hoeffding_1948, Bamber_1975, Sen_1960)과 같은 분들 덕분에, 우리는 AUC 표준 오차를 보다 정확히 추정할 수 있다.

$V_{10}^{r} (X_{i}) = \frac{1}{n} \sum_{j = 1}^{n} ψ (X_{i}^{r}, Y_{j}^{r}) (i = 1, 2, . . ., m)$

$V_{01}^{r} (Y_{j}) = \frac{1}{m} \sum_{i = 1}^{m} ψ (X_{i}^{r}, Y_{j}^{r}) (j = 1, 2, . . ., m)$

$V_{10}^{r} (X_{i})$ 은 집합r에서 값을 기반으로 한 분산이다.

$ψ (X_{i}^{r}, Y_{j}^{r})$ 는 $X_{i}^{r}, Y_{j}^{r}$ 의 관계를 나타내며, 이를 통해 분산을 구한다. $V_{01}^{r} (Y_{j})$ 은 집합 $Y_{j}^{r}$ 에서의 분산이다.

두 값을 통해 각 표본 집합에 대한 분산을 따로 계산한 후, 이를 결합하여 최종적으로 표준 오차를 추정할 수 있다. 각 분산 값이 표본 크기에 따라 가중평균되며, AUC 표준 오차 S는 아래와 같다.

$S = \frac{1}{m} S_{10} + \frac{1}{n} S_{01}$

표준 오차를 통해 우리는 AUC 값이 얼마나 신뢰할 수 있는지 평가할 수 있으며, 이는 모델의 성능을 명확히 이해하는데 중요한 역할을 한다.

Code

이를 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 값에 대해 추가적으로 알아보자.

pROC package을 이용해 AUC의 표준오차(SE)와 신뢰구간(CI)를 구할 수 있다.

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.7314

95% CI: 0.6301-0.8326 (DeLong)

[1] 0.05165929

위처럼 SE/CI를 구하는 것에 더해, ROC curve에서 특정 threshold(0.7)와 비교한 p-value를 구할 수 있다.

z_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

위 경우는 유의하지 않다.

Example

위와 같은 방법으로 aSAH의 다양한 변수들 중, outcome 예측에 가장 좋은 변수를 확인하자. 또한, age/gender를 기준으로 subgroup analysis도 진행하자.

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

데이터 상 크게 유의미한 변수는 없다.

Conclusion

ROC curve를 그리고, AUC 값을 구할 수 있다.
AUC 값의 SE, CI를 구할 수 있다.
threshold와 비교하여 p-value를 구할 수 있다.
이를 실제 데이터에 적용/연습해 보았다.

Reference

[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

Reuse

https://creativecommons.org/licenses/by-nc/4.0/

Citation

BibTeX citation:

@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}
}

For attribution, please cite this work as:

LEE, Hojun. 2024. “DeLong’s Method; for ROC AUC.” November 20, 2024. https://blog.zarathu.com/posts/2024-11-20-DeLongsMethod.