우리는 ROC 곡선에서 AUC 값을 구할 수 있다. AUC 값이 1에 가까우면 모델의 성능이 좋다라고 하며, 0.5에 가까워질 수록 성능이 나쁘다- 라고 한다. 그렇다면, ROC에서 신뢰구간은 어떻게 구할까? 특정 AUC값과 비교하여 p value를 구하려면 어떻게 해야 할까?
이를 위해 DeLong의 AUC 표준 오차 구하는 방법을 소개하려 한다. (E. DeLong, 1988)
다양한 임계값에서의 민감도와 False Positive Rate의 관계를 시각화한 것으로,
민감도(Sensitivity)를 Y축으로, 1 - 특이도(False Positive Rate)을 X축으로 하여 관계를 그린 그래프이다.
민감도(Sensitivity, True Positive Rate)
실제 양성 사례를 얼마나 잘 분류하는지/ 실제 양성인 샘플을 양성으로 올바르게 분류한 비율
특이도(Specificity, True Negative Rate)
실제 음성 사례를 얼마나 잘 분류하는지/ 실제 음성인 샘플을 음성으로 올바르게 분류한 비율
False Positive Rate(FPR) : 1-Spec, 실제 음성인 샘플을 잘못 양성으로 분류한 비율
이상적인 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이라 할 때, 실제로 이벤트가 발생한 집단을 (,n=m), 발생하지 않은 집단을 (,n=n)이라 하자.
이때 실수 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도 선 위에 놓일 것이다.
Figure 2: AUC
AUC
ROC curve는 위와 같이 구할 수 있다. 그럼 이의 넓이: AUC 값은 어떻게 구할까? 보통, 곡선 아래 넓이는 trapezoidal rule을 통해 구한다. 고등학교 때 배운 적분을 떠올리면 된다. 수많은 사다리꼴로 쪼개어 넓이를 근사하던 기억을 되살려 보자.
Figure 3: Trapezes
여기서, Mann-Whitney two sample statistic에 따르면, ROC curve 아래 넓이를 구할때, trapezodial rule로 구한 넓이는 Mann-Whitney two sample statistic으로 구한 넓이로 대체할 수 있다.
Mann-Whitney statistic는 확률를 예측한다. 에서 무작위 추출한 값이 의 값보다 같거나 작을 확률을 추정한다. (test를 시행한 총 집단의 수를 N이라 할 때, 실제로 이벤트가 발생한 집단을 (,n=m), 발생하지 않은 집단을 (,n=n)이라 하자)
모든 (X, Y) 쌍에 대해 X > Y이면 1을, X = Y이면 , X < Y 이면 0을 부여하여 확률을 구한다.
직관적으로, ROC curve와 AUC 값은 곧 모델의 추정이 옳을 확률이며, 이의 성능의 최고값은 1, 최저값은 1/2이라는 점을 고민하면 위의 추정은 그럴듯 하다.
그럼, 확률(AUC)은 위와 같이 정리할 수 있다.
Standard Error
이제, AUC값의 신뢰성을 측정하기 위해서는 SE(standard Error)를 계산하는 것이 필요하다. 이는 DeLong(1988)이 제시한 방법을 참고하자.
ξ를 각각의 집단 간 공분산이라 하자.
ξ₁₀은 의 와 의 간 공분산,
ξ₀₁은 의 와 의 간 공분산,
ξ₁₁은 의 와 의 간 자기공분산이다.
ξ₁₀, ξ₀₁, ξ₁₁의 기대값은 아래와 같다.
기댓값들을 이용하여 AUC 추정치의 분산을 계산할 수 있다.
이와 같은 방법으로 단일 표본 집합에서의 AUC의 표준 오차를 구할 수 있다.
이제, 단일 표본 집합이 아닌 다른 표본집합 r과 s에 대해 다뤄보자. 여러 표본 집합이 있을 경우, 각 표본 간의 상호 공분산 또한 고려되어야한다. 두 표본 집합 r과 s에 대해 AUC의 공분산 계산은 아래와 같다.
이제, 아래 식을 통해 표본 집합의 AUC 값 간의 공분산을 계산할 수 있다.
이 수식을 통해 여러 표본 집합 간의 공분산을 반영하여 AUC의 표준 오차를 더 정확하게 추정할 수 있다.
이를 바탕으로 우리가 궁금한 값 “표준 오차”에 접근해 보자. (Hoeffding_1948, Bamber_1975, Sen_1960)과 같은 분들 덕분에, 우리는 AUC 표준 오차를 보다 정확히 추정할 수 있다.
은 집합r에서 값을 기반으로 한 분산이다.
는 의 관계를 나타내며, 이를 통해 분산을 구한다. 은 집합 에서의 분산이다.
두 값을 통해 각 표본 집합에 대한 분산을 따로 계산한 후, 이를 결합하여 최종적으로 표준 오차를 추정할 수 있다. 각 분산 값이 표본 크기에 따라 가중평균되며, AUC 표준 오차 S는 아래와 같다.
표준 오차를 통해 우리는 AUC 값이 얼마나 신뢰할 수 있는지 평가할 수 있으며, 이는 모델의 성능을 명확히 이해하는데 중요한 역할을 한다.
데이터는 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_aucp_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)#datadata(aSAH)SAH<-as.data.table(aSAH)#subgroupsubgroups<-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"))}#ROCroc_result<-roc(data$outcome, data[[predictor_var]], levels =c("Poor", "Good"))#AUC, CIauc_value<-auc(roc_result)%>%round(3)auc_CI<-ci.auc(roc_result)%>%round(3)# Handle AUC == 1if(auc_value==1){warning("AUC is 1. SE and p-value may be undefined.")se_auc<-NAp_value<-NA}else{# SE se_auc<-sqrt(var(roc_result))# P-valuez_value<-(auc_value-0.8)/se_aucp_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_cutofftable1<-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