ROC Curve and AUC

[안재형]

> library(MASS)
> data(cats)
> cats[1:10,]
   Sex Bwt Hwt
1    F 2.0 7.0
2    F 2.0 7.4
3    F 2.0 9.5
4    F 2.1 7.2
5    F 2.1 7.3
6    F 2.1 7.6
7    F 2.1 8.1
8    F 2.1 8.2
9    F 2.1 8.3
10   F 2.1 8.5

### Sex ~ Bwt and Sex ~ Hwt

# AUC

> install.packages("caTools")

> library(caTools)

> colAUC(cats[,2:3],cats[,1],plotROC=TRUE)
              Bwt      Hwt
F vs. M 0.8338451 0.759048

# ROC Curves

### Epi 패키지 이용

> library(Epi)
> ROC(form=cats$Sex~cats$Bwt, plot="ROC")

# ROC()는 Bwt, Hwt를 predictor로 넣는게 가능하네요(Sex ~ Bwt + Hwt)

> ROC(form=cats$Sex~cats$Bwt+cats$Hwt, plot="ROC")

### ROC 결과물

> data(cats, package="MASS")

> library(Epi)

> out = ROC(form=cats$Sex~cats$Bwt+cats$Hwt)

> out

$res

                        sens       spec       PV+       PV-    lr.eta

0.191243960703729 1.00000000 0.00000000 0.6736111       NaN        NA

0.193025369106466 0.97938144 0.00000000 0.6690141 0.0000000 0.1912440

0.194459564487476 0.97938144 0.02127660 0.6737589 0.3333333 0.1930254

0.202121627955293 0.97938144 0.04255319 0.6785714 0.5000000 0.1944596

0.255310074155313 0.97938144 0.06382979 0.6834532 0.6000000 0.2021216

0.255746729288057 0.97938144 0.08510638 0.6884058 0.6666667 0.2553101

0.257059632742871 0.97938144 0.10638298 0.6934307 0.7142857 0.2557467

.

.

.

0.994208761189084 0.04123711 1.00000000 1.0000000 0.3357143 0.9910133

0.994467173239425 0.03092784 1.00000000 1.0000000 0.3333333 0.9942088

0.995896529418839 0.02061856 1.00000000 1.0000000 0.3309859 0.9944672

0.99643076153607  0.01030928 1.00000000 1.0000000 0.3286713 0.9958965

                          0.00000000 1.00000000 1.0000000 0.3286713 0.9964308

$AUC

[1] 0.8329678

$lr

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

Coefficients:

(Intercept)     cats$Bwt     cats$Hwt  

   -8.69896      3.55391      0.02295  

Degrees of Freedom: 143 Total (i.e. Null);  141 Residual

Null Deviance:      181.9 

Residual Deviance: 132.2        AIC: 138.2 

### logistic regression 결과와 동일

> out2=glm(factor(Sex)~Bwt+Hwt,data=cats,family=binomial)

### lr.eta는 p의 추정치

> fitted(out2)

### 추정치를 데이터 cats에 pred라는 변수로 추가. 이 pred가 lr.eta입니다. lr이 Logistic Regression의 약자인듯합니다.

> cats$pred=fitted(out2)

> head(cats)

  Sex Bwt Hwt      pred

1   F 2.0 7.0 0.1930254

2   F 2.0 7.4 0.1944596

3   F 2.0 9.5 0.2021216

4   F 2.1 7.2 0.2553101

5   F 2.1 7.3 0.2557467

6   F 2.1 7.6 0.2570596

[강성찬]

그래프 예쁘네요^^

[안재형]

R이 그래프는 죽이죠. 그런데 두번째 그래프가 Bwt, Hwt 두 변수를 predictor로 사용한 ROC인데 AUC가 0.7809827네요.
위에 Bwt하나만 사용한 AUC가 0.8338451로 더 높은데 좀 이상하네요. 제가 뭘 잘못 이해하고 있는건지... 시간날때 조사해봐야겠습니다.

[안재형]

음... 알겠습니다. Bwt+Hwt 하면 R에서 보통 두 변수를 predictor로 이용하는데, 이 경우는 그냥 합해서 한변수로 사용한듯하네요. 이런 낭패가^^

[potassium]

안녕하세요, 또 질문을 남기는군요. 여기에 댓글로 남겨도 보시는지 모르겠네요. ROC curve 를 그리고 난후 결국은 cutoff 값을 설정을 하려 합니다. 그런데, ROC curve 에서 직선식과 만나는점의 원래 변수 값을 모르겠네요. 여기서는 Bwt 나 Hwt 값이 되겠죠. 이것을 볼수 있는 방법이 있나요? 이와 함께 가능하면 Bwt cutoff 에 따라, Sens/Spec 등이 변하는 table 을 볼수 있는 함수가 있나요?

[안재형]

ROC는 결국은 logistic regression입니다.
out = ROC(form=cats$Sex~cats$Bwt+cats$Hwt, plot="ROC")
하면 out에 sens, spec, PV+, PV-, lr.eta의 결과물이 있습니다. Btw+Hwt의 조합이 나열되어 있지는 않고 Btw+Hwt에 의해서 추정되는 lr.eta가 있습니다. lr.eta는 p의 추정치 같습니다. 이 lr.eta와 원래 데이터에 붙이면 Btw와 Hwt가 무슨값인지 알수있습니다. 위에 내용을 추가했습니다.

[potassium]

감사합니다. 많은 도움이 되었습니다.

[potassium]

한가지 더 질문을 드리자면, glm 함수의 독립변수는 직접 이름을 넣어야 하는지요?
변수가 많아서 각각 glm 을 할때 반복문으로 변수를 계속 교체 하려 하는데, 예를 들어
for (i in 1: length(val.name)) 과 함께 glm(factor(sex)~ val.name[i], data=aa, bi....) 형식으로 하면 val.name[i] 이 변수 인식이 되지 않는군요. glm 함수 자체의 문제인것 같습니다만, 해결할 수 있는 방법이 있을지요?

[안재형]

R에 내장되어 있는 데이터 attitude를 예를 들면,
lm(attitude[,1] ~ attitude[,2]) 이런식으로 하면 첫번째 변수를 y, 두번째 변수를 x로 회귀분석을 합니다. y, x자리에 넣어야하는건 vector이지 "이름"이 아닙니다. glm(y~x, data=aa) 이런식이면 data=aa라고 했기에 y는 aa$y, x는 aa$x로 vector를 넣은것이지 변수이름을 넣은건 아니라는걸 염두하시기 바랍니다.

[potassium]

아하. 그렇군요. 매번 감사합니다.

[김은주]

Epi 패키지 install한 후에 ROC 함수를 돌렸는데 (out1=ROC(form=na$group~na$sex+na$age+na$BMI+na$SBP+na$DBP+na$AST+na$ALT+na$Glucose+na$Uricacid+na$TG+na$HDL+na$LDL) ROC function을 찾을 수 없다는 에러가 뜨네요.. 뭐가 잘못된걸까요?

[안재형]

library(Epi) 하셨나요?

[kejjjj]

너무 기초적인 실수를햇네요.. R을 처음 다뤄봐서 ^^ 감사합니다!!

[안재형]

저두 가끔 빼먹어요^^