ROC curve analysis에 대한 주위의 요청이 많이 있어 ROC curve에 대한 공부를 해보고 있읍니다. 결국 ROC는 logistic regression 과 관계가 많군요. ROC와 관련된 R 패키지가 많은데 각각 장단점이 있읍니다. Epi 패키지와 pROC 패키지를 분석을 해 보았읍니다. 일단 ROC 그래프 하나만 그리는 것은 매우 쉽습니다. 0과 1로 코딩된 반응변수(response variable), 그리고 이를 설명한 연속변수(explanatory variable) 하나만 있으면 됩니다. moonBook 패키지에 있는 radial 데이타를 사용해보겠읍니다. radial데이타는 115명의 환자 데이타로 이중 male이라는 열이 남자는 1, 여자는 0으로 코딩되어 있고 키(height), 몸무게(weight)로 남자를 구별할 수 있는지 ROC curve analysis를 해보겠읍니다.

require(Epi)
require(pROC)
require(ztable)
require(moonBook)
source("ROC_sub.R")

Epi패키지의 ROC함수를 써서 키와 male에 관한 ROC 그래프를 그려보겠읍니다.

a1=ROC(form=male~height,data=radial,plot="ROC")

비교적 보기 좋은 그래프가 그려집니다. optimal cutoff value 가 표시되며 이때의sensitivity, specificity 등이 표시되고 그래프의 오른쪽 아래 부분에 glm함수의 결과와 model summary, Area under the curve(AUC)가 표시됩니다. ROC 함수의 결과인 a1에는 어떤 값이 들어 있을까요? 다음을 보세요.

a1
$res
                           sens       spec        pvp        pvn
                     1.00000000 0.00000000        NaN 0.49565217
0.0003282059528514   1.00000000 0.03508772 0.00000000 0.48672566
0.000710555991472192 1.00000000 0.05263158 0.00000000 0.48214286
0.00153764716272027  1.00000000 0.07017544 0.00000000 0.47747748
0.00488455998935735  1.00000000 0.08771930 0.00000000 0.47272727
0.00591940529700456  1.00000000 0.10526316 0.00000000 0.46788991
0.00717191333719704  1.00000000 0.17543860 0.00000000 0.44761905
0.0105190713711521   1.00000000 0.22807018 0.00000000 0.43137255
0.0122556995132989   1.00000000 0.24561404 0.00000000 0.42574257
0.0154041231490778   1.00000000 0.35087719 0.00000000 0.38947368
0.0225061899779435   1.00000000 0.36842105 0.00000000 0.38297872
0.0327736695597765   1.00000000 0.42105263 0.00000000 0.36263736
0.0474976383596671   1.00000000 0.54385965 0.00000000 0.30952381
0.0683689492948338   1.00000000 0.59649123 0.00000000 0.28395062
0.0974729628104842   0.98275862 0.61403509 0.02777778 0.27848101
0.11584017881774     0.98275862 0.63157895 0.02702703 0.26923077
0.137142488519708    0.96551724 0.66666667 0.05000000 0.25333333
0.189565696963916    0.96551724 0.68421053 0.04878049 0.24324324
0.256080660565036    0.94827586 0.75438596 0.06521739 0.20289855
0.336251114562426    0.93103448 0.75438596 0.08510638 0.20588235
0.427110239403805    0.86206897 0.85964912 0.14035088 0.13793103
0.523169228894752    0.79310345 0.87719298 0.19354839 0.13207547
0.617544264280943    0.77586207 0.91228070 0.20000000 0.10000000
0.703815216556393    0.72413793 0.94736842 0.22857143 0.06666667
0.777633222738102    0.68965517 0.94736842 0.25000000 0.06976744
0.837306386958018    0.62068966 0.98245614 0.28205128 0.02702703
0.883367752220184    0.60344828 1.00000000 0.28750000 0.00000000
0.917670604017523    0.55172414 1.00000000 0.31325301 0.00000000
0.942540837737788    0.51724138 1.00000000 0.32941176 0.00000000
0.960223862105111    0.39655172 1.00000000 0.38043478 0.00000000
0.972623012123752    0.25862069 1.00000000 0.43000000 0.00000000
0.981232600608132    0.22413793 1.00000000 0.44117647 0.00000000
0.987170337573484    0.17241379 1.00000000 0.45714286 0.00000000
0.99124621893639     0.15517241 1.00000000 0.46226415 0.00000000
0.994035049566317    0.13793103 1.00000000 0.46728972 0.00000000
0.995939038222106    0.08620690 1.00000000 0.48181818 0.00000000
0.997236970503082    0.06896552 1.00000000 0.48648649 0.00000000
0.998120850714131    0.03448276 1.00000000 0.49557522 0.00000000
0.999409670187424    0.00000000 1.00000000 0.50434783        NaN
                          lr.eta
                            -Inf
0.0003282059528514   0.000328206
0.000710555991472192 0.000710556
0.00153764716272027  0.001537647
0.00488455998935735  0.004884560
0.00591940529700456  0.005919405
0.00717191333719704  0.007171913
0.0105190713711521   0.010519071
0.0122556995132989   0.012255700
0.0154041231490778   0.015404123
0.0225061899779435   0.022506190
0.0327736695597765   0.032773670
0.0474976383596671   0.047497638
0.0683689492948338   0.068368949
0.0974729628104842   0.097472963
0.11584017881774     0.115840179
0.137142488519708    0.137142489
0.189565696963916    0.189565697
0.256080660565036    0.256080661
0.336251114562426    0.336251115
0.427110239403805    0.427110239
0.523169228894752    0.523169229
0.617544264280943    0.617544264
0.703815216556393    0.703815217
0.777633222738102    0.777633223
0.837306386958018    0.837306387
0.883367752220184    0.883367752
0.917670604017523    0.917670604
0.942540837737788    0.942540838
0.960223862105111    0.960223862
0.972623012123752    0.972623012
0.981232600608132    0.981232601
0.987170337573484    0.987170338
0.99124621893639     0.991246219
0.994035049566317    0.994035050
0.995939038222106    0.995939038
0.997236970503082    0.997236971
0.998120850714131    0.998120851
0.999409670187424    0.999409670

$AUC
[1] 0.945856

$lr

Call:  glm(formula = form, family = binomial, data = data)

Coefficients:
(Intercept)       height  
   -62.1168       0.3864  

Degrees of Freedom: 114 Total (i.e. Null);  113 Residual
Null Deviance:      159.4 
Residual Deviance: 67.23    AIC: 71.23

ROC함수는 세개의 리스트를 반환하는데 첫번째 res는 데이타프레임으로 각각의 lr.eta값에 대해 민감도, 특이도 등이 들어 있읍니다. 최적의 cutoff point는 이중 민감도+특이도의 합이 제일 큰 것이 되겠읍니다. 즉 다음과 같이 구할 수 있읍니다.

optimal_lr.eta=function(x){
  no=which.max(x$res$sens+x$res$spec)[1]
  result=x$res$lr.eta[no]
  result
}
optimal_lr.eta(a1) 
## [1] 0.4271102

이때의 sensitivity, specificity 등은 ROC 함수의 결과물 중 res 를 보면 알 수 있읍니다. 하지만 키를 얼마로 정하는 것이 남자,여자를 구별할수 있는 최적의 cutoff value일까요? 이것은 일반화선형모델인 glm함수의 결과를 이용해 계산할 수 있읍니다. a1의 리스트 중 세번쨰 lr을 보면 glm함수의 결과가 나와 있는데 절편이 -62.1168, 키에 대한 slope가 0.3864으로 나와 있읍니다. 이것을 회귀공식으로 나타내보면 다음과 같습니다.

\(outcome=\frac{1}{1+e^{(-{\beta}_0+{\beta}_1 Height)}}\)

즉 optimal cutoff value인 height 는 outcome에 lr.eta의 cutoff value인 0.4271을 넣고 \({\beta}_0\)\({\beta}_1\) 에 각각 값을 대입한 후 방정식을 풀면 됩니다. 다음과 같은 함수를 만들 수 있읍니다.

optimal_cutpoint=function(x){
   y=optimal_lr.eta(x)
   b0=unname(x$lr$coeff[1])
   b1=unname(x$lr$coeff[2])
   result=(-log(1/y-1)-b0)/b1
   result
} 
optimal_cutpoint(a1) 
## [1] 160

위의 회귀 공식을 풀어 키 160을 얻었읍니다. 즉 남자, 여자를 키 160을 기준으로 구별하면 민감도 86.2% , 특이도 86.0%로 남,여를 구분할 수 있읍니다. 위와 같은 방법으로 몸무게로 남여를 구별하는 ROC 커브를 그려보겠읍니다.

a2=ROC(form=male~weight,data=radial,plot="ROC")

하지만 불행하게도 두개의 ROC 커브를 하나의 plot으로 그릴 수는 없읍니다. 두개를 하나로 그리려면 a1$res와 a2$res에 있는 sens와 spec를 이용해 다음과 같이 그릴 수 있읍니다.

plot(a1$res$sens~I(1-a1$res$spec),type="l")
par(new=TRUE)
plot(a2$res$sens~I(1-a2$res$spec),type="l",axes=F,xlab="",ylab="",col="red")

이 그래프는 멋이 없읍니다. 또한 두개의 ROC curve를 그리기는 했지만 AUC의 신뢰구간도 알고 싶은 경우에는 Epi패키지로는 알 수 없읍니다. 이 목적으로는 pROC패키지의 roc함수를 쓰면 됩니다.

b1=roc(male~height,radial,ci=T,percent=T)
b2=roc(male~weight,radial,ci=T,percent=T)

plot(b1)

Call:
roc.formula(formula = male ~ height, data = radial, ci = T, percent = T)

Data: height in 57 controls (male 0) < 58 cases (male 1).
Area under the curve: 94.59%
95% CI: 91.07%-98.1% (DeLong)
plot(b2,add=TRUE,col="red")


Call:
roc.formula(formula = male ~ weight, data = radial, ci = T, percent = T)

Data: weight in 57 controls (male 0) < 58 cases (male 1).
Area under the curve: 80.49%
95% CI: 72.53%-88.45% (DeLong)

또한 두개의 커브가 다른지 검정도 하고 싶습니다. 이떄에는 역시 pROC패키지의 roc.test를 하면 됩니다.

roc.test(b1,b2,plot=T)

    DeLong's test for two correlated ROC curves

data:  b1 and b2
Z = 3.8586, p-value = 0.000114
alternative hypothesis: true difference in AUC is not equal to 0
sample estimates:
AUC of roc1 AUC of roc2 
   94.58560    80.49002 

AUC가 0.5보다 큰지 검정도 하고 싶습니다. 이때에는 wilcox.test를 하면 됩니다.

wilcox.test(height~male,data=radial)

    Wilcoxon rank sum test with continuity correction

data:  height by male
W = 179, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0

이상의 여러과정을 한번에 plot하나로 그릴 수 있는 함수를 하나 만들었읍니다. 다음은 사용 예입니다. 내부에서는 두개의 패키지를 모두 이용하고 회귀공식에서 cutoff point를 계산하고 wilcox.test를 실시하고 DeLong’s test까지 실시합니다.

plot_ROC(a1,a2)

하지만 pROC패키지의 roc함수는 하나의 설명변수만 쓸 수 있다는 문제가 있읍니다. 예를 들어 키와 몸무게를 같이 넣고 ROC curve를 그리려면 roc함수로는 불가능합니다. 그런 경우 ROC함수를 써서 만들면 됩니다. 다음은 키와 몸무게를 모두 설명변수로 했을때 ROC 커브를 a3로 하고 plot_ROC함수로 세개의 커브를 모두 그린 예입니다.

a3=ROC(form=male~height+weight,data=radial,plot="")
plot_ROC(a1,a2,a3,show.sens=FALSE)

ROC curve는 glm()함수와 유사합니다. 예를 들어 설명변수를 4개 넣고 glm()함수를 사용하여 모형을 만들고 stepwise backward elimination을 사용해 최종 모형을 선택할 수 있읍니다. 제가 만든 step_ROC함수는 먼저 규정식에 있는 열만을 대상으로 새로운 데이타프레임을 만들고 NA값을 제거 한 후 step함수를 사용하여 AIC값을 기준으로 stepwise backward elimination을 시행합니다. step_ROC 함수는 initial model과 final model, 그리고 table(즉 ANOVA table)을 리스트로 반환합니다. plot_ROC함수는 type=1로 하면 회색 바탕의 그래프를를 그릴 수 있읍니다.

# perform stepwise backward regression
result=step_ROC(male~height+weight+TC,data=radial,plot=FALSE)
Start:  AIC=76.6
male ~ height + weight + TC

         Df Deviance     AIC
- weight  1   11.892  74.665
- TC      1   11.907  74.808
<none>        11.885  76.604
- height  1   20.305 134.589

Step:  AIC=74.67
male ~ height + TC

         Df Deviance     AIC
- TC      1   11.919  72.919
<none>        11.892  74.665
- height  1   27.818 167.845

Step:  AIC=72.92
male ~ height

         Df Deviance     AIC
<none>        11.919  72.919
- height  1   27.991 166.542
plot_ROC(result$initial,result$final,show.lr.eta=FALSE,show.sens=FALSE,type=1)

anova table은 다음과 같이 출력할 수 있읍니다.

print(result$table,type="html")
Analysis of Deviance Table
Model 1: male ~ height + weight + TC
Model 2: male ~ height
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 108 65.08
2 110 65.47 -2.00 -0.39 0.8225

이상 제가 요즘 ROC curve에 대해 공부하면서 만든 내용입니다. 이 내용을 바탕으로 주말을 이용해 제가 만든 shiny app에 ROC curve analysis를 추가하겠읍니다. 제가 통계를 잘 모르니 혹시 이해를 잘못했을 수 있읍니다. 통계에 관한 오류 문제는 언제든지 지적 바랍니다.