Logistic Regression
Seo Young Park
2024-01-26
\(~\)
R을 이용해서 logistic regression을 수행하는 방법을 알아보도록 하겠습니다. 이 강의에서는 logistic regression 자체에 대한 내용은 소개하지 않습니다. 회귀분석, logistic regression, odds ratio에 대한 배경지식이 없는 분은 무작정 분석하지 마시고 반드시 logistic regression이 무슨 분석인지 공부를 하시기 바랍니다. 여기저기에 잘 소개되어있지만, 예를 들어 이런 책을 보시면 잘 설명되어있습니다: 바이오통계학
그런데, R에서 기본이 되는 data structure가 벡터(vector)와 행렬(matrix)입니다. 벡터와 행렬을 잘 이해하고 이용하면 아주 번거로울 수 있는 작업을 아주 간단하게 수행 할 수도 있습니다. 꼭 필요한 개념만 잠시 살펴보고 넘어가겠습니다.
\(~\)
Playing around with the vectors and matrices in R
지난 시간에 살펴봤듯이, 벡터는 여러개의 elements들이 일렬로 저장되어있습니다. 간단하게 세 개의 벡터를 만들어봅시다.
aa<-1:4
bb<-5:8
cc<-c(10, 20, 30, 40)저장된 벡터를 확인해봅시다.
aa## [1] 1 2 3 4bb## [1] 5 6 7 8cc## [1] 10 20 30 40벡터 여러개를 각 행으로 합치거나, 각 열로 합쳐서 행렬을 만들 수 있습니다.
mat1<-rbind(aa, bb, cc)
mat1## [,1] [,2] [,3] [,4]
## aa 1 2 3 4
## bb 5 6 7 8
## cc 10 20 30 40mat2<-cbind(aa, bb, cc)
mat2## aa bb cc
## [1,] 1 5 10
## [2,] 2 6 20
## [3,] 3 7 30
## [4,] 4 8 40꼭 rbind()나 cbind()를 이용해야하만 행렬을
만들 수 있는 것은 아닙니다. matrix()안에 각 원소의 값을
벡터로 넣어주고, 행의 갯수 또는 열의 갯수를 지정해서도 만들 수
있습니다.
mat3<-matrix(11:18, nrow=4)
mat3## [,1] [,2]
## [1,] 11 15
## [2,] 12 16
## [3,] 13 17
## [4,] 14 18mat4<-matrix(11:18, ncol=4)
mat4## [,1] [,2] [,3] [,4]
## [1,] 11 13 15 17
## [2,] 12 14 16 18행렬이 몇행 몇열을 가지고 있는지 알고 싶으면,
dim()함수를 쓰면 됩니다.
dim(mat1)## [1] 3 4dim(mat2)## [1] 4 3dim(mat3)## [1] 4 2dim(mat4)## [1] 2 4벡터 뿐만 아니라, 행렬도 rbind()나
cbind()를 이용해서 합칠 수 있습니다.
mat5<-rbind(mat1, mat4)
mat5## [,1] [,2] [,3] [,4]
## aa 1 2 3 4
## bb 5 6 7 8
## cc 10 20 30 40
## 11 13 15 17
## 12 14 16 18mat6<-cbind(mat2, mat3)
mat6## aa bb cc
## [1,] 1 5 10 11 15
## [2,] 2 6 20 12 16
## [3,] 3 7 30 13 17
## [4,] 4 8 40 14 18x행, y열의 원소를 뽑아내려면 행렬 이름 옆에 [x, y]를
붙여주면 됩니다.
mat6[1, 2]## bb
## 5mat6[1, ]## aa bb cc
## 1 5 10 11 15mat6[-1, ]## aa bb cc
## [1,] 2 6 20 12 16
## [2,] 3 7 30 13 17
## [3,] 4 8 40 14 18mat6[ , 2]## [1] 5 6 7 8mat6[ , 2:3]## bb cc
## [1,] 5 10
## [2,] 6 20
## [3,] 7 30
## [4,] 8 40자 이제 벡터와 행렬에 대한 아주 기본적인 내용은 보여드렸습니다.
그럼 본격적으로 logistic regression을 해봅시다.
\(~\)
Logistic Regression
지난주까지 정리한 데이터 dat3가 필요합니다. 지난주에
보여드린 코드를 다시 돌려도 되지만, 분석을 할 때마다 데이터 클리닝 하는
코드를 매번 다시 돌리기 귀찮을 수도 있습니다. 그럴 때는
save()함수를 이용해서, 저장하고 싶은 object를 .RData
확장자를 가진 파일로 저장하는 것이 좋습니다. 아래는 지난주까지 돌렸던
코드입니다. 맨 아랫줄의 코드가 dat3를 clean_data.RData라는
파일에 저장하는 명령어입니다.
setwd("C:/Users/biostat81/Dropbox/Teaching/AMC Radiology R 2024")
library(plyr)
library(dplyr)
library(tableone)
dat0<-read.csv("ex_data.csv", nrow=52)
dat1<-dat0 %>% distinct()
dat2<-dat1[-45, ]
dat3<-dat2 %>% rename(Recur=Recur..1..Recur..0..Censored.) %>%
mutate_at(vars(OP_date, Recur_date),
list(~as.Date(., format="%Y-%m-%d"))) %>%
mutate(CA19.9=as.double(replace(CA19.9, CA19.9=="<1.0", 1)),
adc_decrease=as.double(replace(adc_decrease, adc_decrease=="na" | adc_decrease==".", NA)),
log.CA19.9=log(CA19.9),
log.CEA=log(CEA),
CA19.9.group=ifelse(CA19.9>37, 1, 0),
TNM.b=mapvalues(TNM, from=c(1, 2, 3, 4, 5), to=c("1or2", "1or2", "3", "4", "5")),
RFS=as.double(Recur_date - OP_date)
) %>%
mutate_at(vars(sex, Recur, local_6m, TNM, CCI, CA19.9.group, TNM.b),
list(~as.factor(.)))
save(dat3, file="clean_data.RData")따라서, 위의 코드를 실행해서 얻은 clean_data.RData 파일이 working
directory에 저장되어있다면, RStudio를 종료했다가 다시 열었을 때, 위의
코드를 실행하지 않아도 간단히 load() 함수를 이용해서 이미
만들어놓은 dat3를 불러올 수 있습니다.
load("clean_data.RData")분석에 필요한 패키지 두가지를 불러옵니다.
library(lmtest)
library(pROC)먼저, local_6m을 결과변수로, age를
설명변수로 하는 univariate logistic regression 모델을 fitting하고 결과를
보는 방법은 다음과 같습니다.
obj<-glm(local_6m~age, data=dat3, family="binomial")
summary(obj)##
## Call:
## glm(formula = local_6m ~ age, family = "binomial", data = dat3)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.44512 1.36291 -1.060 0.289
## age 0.02269 0.02345 0.968 0.333
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 68.994 on 49 degrees of freedom
## Residual deviance: 68.031 on 48 degrees of freedom
## AIC: 72.031
##
## Number of Fisher Scoring iterations: 4summary() 함수는 회귀분석 결과를 보는 가장 기본적인
방법이지만, 아쉬운 점이 많습니다. OR을 보여주지도 않고, OR의 95%
confidence interval도 보여주지 않죠. 이것들을 뽑아내서 보기좋게
정렬해봅시다. 일단 summary() 외에도 결과를 뽑아내는 세가지
함수를 알아두시면 좋습니다.
coef(obj) #### beta coefficients## (Intercept) age
## -1.445116 0.022693coefci(obj) #### 95% CI## 2.5 % 97.5 %
## (Intercept) -4.11636956 1.22613782
## age -0.02326803 0.06865403coeftest(obj) #### p-value##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.445116 1.362910 -1.0603 0.2890
## age 0.022693 0.023450 0.9677 0.3332일단 beta coefficients에 지수함수를 취해서 OR (odds ratio)을 계산해봅시다.
exp(coef(obj))## (Intercept) age
## 0.2357188 1.0229524OR의 95% CI 역시, beta coefficient의 95% CI에 지수함수를 취해주면 됩니다.
exp(coefci(obj))## 2.5 % 97.5 %
## (Intercept) 0.0163036 3.408042
## age 0.9770006 1.071066앞에서 배운 행렬 합치기를 이용하여, OR과 OR의 95% CI를 나란히 붙여볼까요?
cbind(exp(coef(obj)), exp(coefci(obj)))## 2.5 % 97.5 %
## (Intercept) 0.2357188 0.0163036 3.408042
## age 1.0229524 0.9770006 1.071066이제 p-value만 얻어내서 오른쪽 끝내 붙여주면 완벽합니다.
coeftest()의 output이 행렬이고, p-value는 네번째 열에
있다는 점을 이용하면 p-value만 뽑아낼 수 있겠죠.
coeftest(obj)[, 4]## (Intercept) age
## 0.2890005 0.3331837그러면 이 p-value도 cbind() 안에 넣어서 OR, 95% CI,
p-value 순서로 보이게 정리해줍니다. 열의 이름(column name)이
정의되어있지 않은 열은 이름을 붙여줄 수도 있습니다.
cbind(OR=exp(coef(obj)), exp(coefci(obj)), p.value=coeftest(obj)[, 4])## OR 2.5 % 97.5 % p.value
## (Intercept) 0.2357188 0.0163036 3.408042 0.2890005
## age 1.0229524 0.9770006 1.071066 0.3331837사실 univariate logistic regression에서는 intercept는 관심없을 경우도 많죠. intercept 텀이 첫번째 행에 있으니까, 첫번째 행을 지워줄 수도 있습니다.
cbind(OR=exp(coef(obj)), exp(coefci(obj)), p.value=coeftest(obj)[, 4])[-1, ]## OR 2.5 % 97.5 % p.value
## 1.0229524 0.9770006 1.0710656 0.3331837따라서, 정리하자면, 아래의 두 줄을 반복하되 설명변수 이름만 바꿔가면서 넣어주면, 여러 설명변수에 대한 univariate logistic regression 결과를 얻어낼 수 있습니다.
obj<-glm(local_6m~age, data=dat3, family="binomial")
cbind(OR=exp(coef(obj)), exp(coefci(obj)), p.value=coeftest(obj)[, 4])[-1, ]## OR 2.5 % 97.5 % p.value
## 1.0229524 0.9770006 1.0710656 0.3331837그래서 관심있는 변수들에 대해서 각각 위의 두 줄을 실행시켜서 결과를
얻어내면 되는데, 후보 설명변수가 많은 경우, 결과를 한줄 씩 일일이
메모장에 붙였다가 엑셀에 붙이는 게 좀 번거로울 수 있죠. 그래서 아래와
같이 각 모델의 결과를 벡터로 저장한 후, rbind()를 이용해서
행렬로 묶어버리면, 한번에 모든 결과를 복/붙 할 수 있습니다.
obj<-glm(local_6m~age, data=dat3, family="binomial")
result1<-cbind(OR=exp(coef(obj)), exp(coefci(obj)), p.value=coeftest(obj)[, 4])[-1, ]
obj<-glm(local_6m~TNM.b, data=dat3, family="binomial")
result2<-cbind(OR=exp(coef(obj)), exp(coefci(obj)), p.value=coeftest(obj)[, 4])[-1, ]
obj<-glm(local_6m~log.CA19.9, data=dat3, family="binomial")
result3<-cbind(OR=exp(coef(obj)), exp(coefci(obj)), p.value=coeftest(obj)[, 4])[-1, ]
obj<-glm(local_6m~adc_decrease, data=dat3, family="binomial")
result4<-cbind(OR=exp(coef(obj)), exp(coefci(obj)), p.value=coeftest(obj)[, 4])[-1, ]
obj<-glm(local_6m~apt_1_mean, data=dat3, family="binomial")
result5<-cbind(OR=exp(coef(obj)), exp(coefci(obj)), p.value=coeftest(obj)[, 4])[-1, ]
rbind(age=result1, result2, log.CA19.9=result3, adc_decrease=result4, apt_1_mean=result5)## OR 2.5 % 97.5 % p.value
## age 1.0229524 0.9770006 1.0710656 0.333183691
## TNM.b3 2.5000000 0.3409325 18.3320730 0.367380415
## TNM.b4 1.9642857 0.3182447 12.1240634 0.467211915
## TNM.b5 5.0000000 0.4720500 52.9604881 0.181364351
## log.CA19.9 1.4997746 1.0295015 2.1848671 0.034735885
## adc_decrease 0.9387088 0.8984004 0.9808256 0.004734741
## apt_1_mean 2.1528407 1.1435493 4.0529280 0.017524140TNM.b가 설명변수인 모델은, TNM.b가 총 네 개의 범주가 있기 때문에 4-1 = 3개의 dummy variable에 대한 결과가 나옵니다. 그럼 각 범주를 pairwise로 비교한 결과가 아니라, 모든 범주간의 차이를 종합적으로 고려한 p-value를 계산하려면 어떻게 해야할까요? 설명변수에서 TNM.b를 제외한 모델, 즉 상수항만 있는 모델을 먼저 fitting 한 다음, 원 모델과 상수항만 있는 모델을 비교하는 likelihood ratio test를 수행하면 됩니다.
obj<-glm(local_6m~TNM.b, data=dat3, family="binomial")
obj0<-update(obj, .~.-TNM.b)
lrtest(obj, obj0)## Likelihood ratio test
##
## Model 1: local_6m ~ TNM.b
## Model 2: local_6m ~ 1
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 4 -33.473
## 2 1 -34.497 -3 2.0484 0.5624Pr(>Chisq) 밑에 나온 값이 두 모델을 비교한 p-value입니다. 즉, TNM.b에 대한 overall p-value가 됩니다.
이번에는 여러개의 설명변수가 포함된 multivariable logistic regression
모델을 fitting해보겠습니다. univariate logistic regression 과 같은
방법으로 하되 설명변수를 +로 연결해주면 되고, output를 뽑아내는 방법과
범주가 세 개 이상인 범주형 변수의 p-value 구하는 법도 동일합니다. 한가지
주의할 점은, 결측이 있는 경우에는 na.action=na.exclude라는
옵션을 포함해주는 것이 좋습니다. 이 옵션의 역할은 뒤에서 다시
설명하겠습니다.
obj.m<-glm(local_6m~age+TNM.b+log.CA19.9+adc_decrease+apt_1_mean, family="binomial", data=dat3, na.action=na.exclude)
summary(obj.m)##
## Call:
## glm(formula = local_6m ~ age + TNM.b + log.CA19.9 + adc_decrease +
## apt_1_mean, family = "binomial", data = dat3, na.action = na.exclude)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.24747 2.53194 -1.283 0.1996
## age 0.01811 0.03241 0.559 0.5763
## TNM.b3 0.50435 1.32151 0.382 0.7027
## TNM.b4 0.05256 1.22987 0.043 0.9659
## TNM.b5 0.85461 1.63482 0.523 0.6011
## log.CA19.9 0.41150 0.26393 1.559 0.1190
## adc_decrease -0.06830 0.02761 -2.473 0.0134 *
## apt_1_mean 0.70900 0.38260 1.853 0.0639 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 66.208 on 47 degrees of freedom
## Residual deviance: 45.289 on 40 degrees of freedom
## (2 observations deleted due to missingness)
## AIC: 61.289
##
## Number of Fisher Scoring iterations: 6cbind(OR=exp(coef(obj.m)), exp(coefci(obj.m)), p.value=coeftest(obj.m)[, 4])[-1, ]## OR 2.5 % 97.5 % p.value
## age 1.0182766 0.95559911 1.0850651 0.57631545
## TNM.b3 1.6559086 0.12421183 22.0754597 0.70272400
## TNM.b4 1.0539697 0.09461546 11.7407039 0.96590939
## TNM.b5 2.3504489 0.09540920 57.9043754 0.60114626
## log.CA19.9 1.5090839 0.89961411 2.5314569 0.11896336
## adc_decrease 0.9339835 0.88478096 0.9859222 0.01338241
## apt_1_mean 2.0319593 0.95993625 4.3011801 0.06386873multivariable logistic regression 모델을 fit할 때에도 TNM.b의 overall p-value를 구하기 위해선 likelihood ratio test를 해야합니다.
obj.m0<-update(obj.m, .~.-TNM.b)
lrtest(obj.m, obj.m0)## Likelihood ratio test
##
## Model 1: local_6m ~ age + TNM.b + log.CA19.9 + adc_decrease + apt_1_mean
## Model 2: local_6m ~ age + log.CA19.9 + adc_decrease + apt_1_mean
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 8 -22.645
## 2 5 -22.946 -3 0.6028 0.8958logistic regression model은 결과변수 Y값이 1일 확률(우리 예제에서는
local_6m=1일 확률)을 독립변수들을 이용해서 추정하는 수식이라고
생각하시면 됩니다. 그러면 데이터를 이용하여 fit한 logistic regression
model을 이용하여, 우리 데이터의 환자들 각각의 6개월내 국소 재발 확률을
계산하려면 어떻게 해야할까요? predict() 함수를 이용하면
쉽게 계산할 수 있습니다.
predict(obj.m, type="response")## 1 2 3 4 5 6
## 0.996159313 0.824968098 0.434773172 0.167485627 0.515421686 NA
## 7 8 9 10 11 12
## 0.787463979 0.266181849 0.684303949 0.865112205 0.918082899 0.526567315
## 13 14 15 16 17 18
## 0.364978753 0.080087874 0.046303323 0.455567860 0.125818166 0.246216348
## 19 20 21 22 23 24
## 0.743047539 0.752886233 0.610323932 0.256621921 0.728282855 0.158503103
## 25 26 27 28 29 30
## 0.039707425 0.123367250 0.290681592 0.050641847 0.065728466 0.191942046
## 31 32 33 34 35 36
## 0.774829899 0.803690718 0.317085678 0.623887711 0.346718205 NA
## 37 38 39 40 41 42
## 0.310859089 0.441498353 0.782050282 0.087717242 0.007887333 0.502961916
## 43 44 46 47 48 49
## 0.995700002 0.496823280 0.789630164 0.999999995 0.034392855 0.741003375
## 50 51
## 0.345620997 0.280416282환자 두 명은 6개월 내 국소 재발 확률의 추정값이 NA, 즉 결측으로
표시된 것이 보이시나요? 이 환자들은 독립변수 중에 결측이 있어서, 6개월
내 국소 재발 확률을 추정할 수 없는 환자들입니다. 이렇게 추정 확률이
결측인 환자가 그냥 제외되지 않고 결측이라고 NA가 찍혀서 나오는 이유는,
앞에서 logistic regression을 할 때 glm() 함수 내에서
na.action=na.exclude 옵션을 사용했기 때문입니다. 만약 이
옵션을 쓰지 않으면 어떻게 되는지 봅시다.
obj.x<-glm(local_6m~age+TNM.b+log.CA19.9+adc_decrease+apt_1_mean, family="binomial", data=dat3)
predict(obj.x, type="response")## 1 2 3 4 5 7
## 0.996159313 0.824968098 0.434773172 0.167485627 0.515421686 0.787463979
## 8 9 10 11 12 13
## 0.266181849 0.684303949 0.865112205 0.918082899 0.526567315 0.364978753
## 14 15 16 17 18 19
## 0.080087874 0.046303323 0.455567860 0.125818166 0.246216348 0.743047539
## 20 21 22 23 24 25
## 0.752886233 0.610323932 0.256621921 0.728282855 0.158503103 0.039707425
## 26 27 28 29 30 31
## 0.123367250 0.290681592 0.050641847 0.065728466 0.191942046 0.774829899
## 32 33 34 35 37 38
## 0.803690718 0.317085678 0.623887711 0.346718205 0.310859089 0.441498353
## 39 40 41 42 43 44
## 0.782050282 0.087717242 0.007887333 0.502961916 0.995700002 0.496823280
## 46 47 48 49 50 51
## 0.789630164 0.999999995 0.034392855 0.741003375 0.345620997 0.280416282국소 재발 확률의 추정치가 결측인 환자가 그냥 사라진 채 48명 환자의
결과값만 나왔습니다. 결과가 이렇게 나와도 상관없는 경우는
na.action=na.exclude을 쓰지 않아도 되지만, 이런 결과가
나오면 곤란한 경우가 있습니다. 그 경우는 뒤에서 다시 살펴보겠습니다.
\(~\)
회귀분석에서 multivariable model에 들어갈 변수는 어떻게 정하느냐, 하는 것은 간단한 문제가 아니고, 각 연구의 목적과 가설, 자료 특성에 따라 달라져야 하는 문제입니다. 이 주제로 두꺼운 책 한권이 나올만큼 깊고 어려운 문제라서 이 강의에서 짧게 말씀드리기는 어렵습니다. 그렇지만 강조하고 싶은 것은, 임상과학에서 공식처럼 쓰이고 있는 ‘univariate 모델에서 유의한 것 골라낸 후 backward elimination하기’ 같은 방식을 무턱대로 사용하는 것은 지양해야합니다. 그래도 경우에 따라 backward elimination을 해야만 하는 경우도 있긴 합니다. 그래서 R을 이용해서 backward elimination을 수행하는 방법을 보여드리겠습니다.
먼저 후보 변수가 모두 들어간 모델을 먼저 fitting 합니다. 우리는
앞에서 obj.m에 저장한 모델을 사용하겠습니다. 그런 다음에
step()함수 안에 결과를 넣어줍니다. 그러면 backward
elimination 적용한 이후 최종 선택된 모델에 대한 정보를 output으로 내놓게
됩니다. 그것을 일반적인 multivariable model 결과 뽑아내는 방식으로
뽑아내면 됩니다.
obj.b<-step(obj.m, direction = "backward")## Start: AIC=61.29
## local_6m ~ age + TNM.b + log.CA19.9 + adc_decrease + apt_1_mean
##
## Df Deviance AIC
## - TNM.b 3 45.892 55.892
## - age 1 45.604 59.604
## <none> 45.289 61.289
## - log.CA19.9 1 48.242 62.242
## - apt_1_mean 1 49.577 63.577
## - adc_decrease 1 53.200 67.200
##
## Step: AIC=55.89
## local_6m ~ age + log.CA19.9 + adc_decrease + apt_1_mean
##
## Df Deviance AIC
## - age 1 46.201 54.201
## <none> 45.892 55.892
## - log.CA19.9 1 49.508 57.508
## - apt_1_mean 1 50.003 58.003
## - adc_decrease 1 54.014 62.014
##
## Step: AIC=54.2
## local_6m ~ log.CA19.9 + adc_decrease + apt_1_mean
##
## Df Deviance AIC
## <none> 46.201 54.201
## - log.CA19.9 1 49.995 55.995
## - apt_1_mean 1 50.286 56.286
## - adc_decrease 1 54.545 60.545summary(obj.b)##
## Call:
## glm(formula = local_6m ~ log.CA19.9 + adc_decrease + apt_1_mean,
## family = "binomial", data = dat3, na.action = na.exclude)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.96453 1.22920 -1.598 0.1100
## log.CA19.9 0.43295 0.25000 1.732 0.0833 .
## adc_decrease -0.06715 0.02657 -2.527 0.0115 *
## apt_1_mean 0.65424 0.37152 1.761 0.0782 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 66.208 on 47 degrees of freedom
## Residual deviance: 46.201 on 44 degrees of freedom
## (2 observations deleted due to missingness)
## AIC: 54.201
##
## Number of Fisher Scoring iterations: 6cbind(OR=exp(coef(obj.b)), exp(coefci(obj.b)), p.value=coeftest(obj.b)[, 4])[-1, ]## OR 2.5 % 97.5 % p.value
## log.CA19.9 1.5417943 0.9445499 2.5166797 0.08331307
## adc_decrease 0.9350562 0.8876002 0.9850496 0.01151092
## apt_1_mean 1.9236854 0.9287393 3.9845038 0.07824189predict(obj.b, type="response")## 1 2 3 4 5 6 7
## 0.98860727 0.84752382 0.39543883 0.24927445 0.40983357 NA 0.62963339
## 8 9 10 11 12 13 14
## 0.18596725 0.66693796 0.90306187 0.93276152 0.51252422 0.47943068 0.08903005
## 15 16 17 18 19 20 21
## 0.06471112 0.50940013 0.11743565 0.22766609 0.83507470 0.76907455 0.63959317
## 22 23 24 25 26 27 28
## 0.20409517 0.74595116 0.11883699 0.04628036 0.09424343 0.39604274 0.09772262
## 29 30 31 32 33 34 35
## 0.08817852 0.30804632 0.64120798 0.81424297 0.32314417 0.53403821 0.35983084
## 36 37 38 39 40 41 42
## NA 0.40520016 0.35662241 0.65324177 0.08010164 0.01060379 0.51142674
## 43 44 46 47 48 49 50
## 0.99295449 0.44451878 0.79929045 0.99999998 0.03799788 0.75543055 0.40959670
## 51
## 0.31817290\(~\)
ROC analysis
logistic regression model의 성능을 평가하는 여러가지 방법 중 아마도 가장 많이 쓰이는 것은 c-index일 것입니다. 이것은, logistic regression model을 이용하여 추정한 P(Y=1)이, 실제로 Y=1인 그룹과 Y=0인 그룹을 얼마나 잘 분류하느냐를 나타내는 지표입니다. 구체적으로는, logistic regression model을 이용하여 추정한 P(Y=1)을 가지고 ROC (receiver operating characteristics) curve라는 것을 그려서 그 curve의 아래쪽 면적(AUC=area under the curve)을 계산하면 그게 바로 c-index가 됩니다. ROC 분석은 알고보면 그렇게 어려운 개념이 아니지만 제대로 모르는 채로 사용하는 경우가 매우 많습니다. 그러면 안됩니다! 반드시 ROC 분석이 무엇인지, logistic regression model의 성능을 평가하는데 어떻게 쓰이는 것인지 공부하시기 바랍니다. 이 역시 앞에서 소개드린 바이오통계학 책에 매우 자세히 나와있습니다.
자, 그럼 c-index를 계산하기 위해서 ROC 분석을 해봅시다. 이것은
pROC라는 패키지의 roc() 함수를 이용하는 것이 가장
간단합니다. 이 함수 안에 예측하려는 결과변수의 실제 값~ 예측에 쓰는
마커의 값을 넣어주면 AUC를 계산해주는데요, AUC의 95% CI까지 얻고 싶으면
ci=T 옵션을 넣어주면 됩니다. 이 AUC가 곧 logistic regression model의
c-index 입니다.
roc(dat3$local_6m~predict(obj.m, type="response"), ci=T)## Setting levels: control = 0, case = 1## Setting direction: controls < cases##
## Call:
## roc.formula(formula = dat3$local_6m ~ predict(obj.m, type = "response"), ci = T)
##
## Data: predict(obj.m, type = "response") in 26 controls (dat3$local_6m 0) < 22 cases (dat3$local_6m 1).
## Area under the curve: 0.8462
## 95% CI: 0.7386-0.9537 (DeLong)roc(dat3$local_6m~predict(obj.b, type="response"), ci=T)## Setting levels: control = 0, case = 1
## Setting direction: controls < cases##
## Call:
## roc.formula(formula = dat3$local_6m ~ predict(obj.b, type = "response"), ci = T)
##
## Data: predict(obj.b, type = "response") in 26 controls (dat3$local_6m 0) < 22 cases (dat3$local_6m 1).
## Area under the curve: 0.8357
## 95% CI: 0.7248-0.9466 (DeLong)그런데 아까 na.action=na.exclude 옵션을 쓰지 않고 모델을
돌린 후 저장했던 obj.x라는 객체를 이용해서 같은 계산을 하면
어떻게 될까요?
roc(dat3$local_6m~predict(obj.x, type="response"), ci=T)## Error in model.frame.default(formula = dat3$local_6m ~ predict(obj.x, : variable lengths differ (found for 'predict(obj.x, type = "response")')obj.x를 만들 때 na.action=na.exclude 옵션을
포함하지 않았기 때문에, predicted()함수는 결측이 있는
환자는 무시하고 결측이 없는 환자 48명의 predicted probability만으로
이루어진 벡터를 만들어낸다는 걸 앞에서 확인했었죠. 그런데 outcome인
dat3$local_6m은 총 50개의 값이 있는 벡터이니,
roc()함수는 이 중 어떤 48개의 값을 써야할지 알 수 없게
되어버리는 겁니다. 즉, 실제 학생은 50명인데 그 중 48명의 시험점수만 주고
이들을 성적 순으로 줄 세우라고 명령하는 것과 같아요. 그래서 오류가 나는
것입니다. 이런 종류의 오류를 예방하기 위해서 결측이 있는 데이터를 다룰
때는 na.action=na.exclude 옵션을 꼭 써주는 것이
좋습니다.
ROC 분석의 결과를 저장해놓으면 더 다양한 분석도 할 수 있습니다. 일단
저장한 이후 결과를 좀더 깔끔하게 뽑아낼 수도 있고요. plot()
함수 안에 넣어 ROC curve를 그릴 수도 있습니다.
roc.m<-roc(dat3$local_6m~predict(obj.m, type="response"), ci=T)## Setting levels: control = 0, case = 1## Setting direction: controls < casesroc.m$auc## Area under the curve: 0.8462roc.m$ci## 95% CI: 0.7386-0.9537 (DeLong)plot(roc.m)
두개의 ROC curve를 하나의 그래프에 겹쳐 그리고 싶을 때도 있죠. 그럴
때는 add=T 옵션이 유용합니다. 그리고 두 직선을 구분하기
위해 lty 옵션의 값을 바꿔줍니다. (디폴트가
lty=1 입니다.)
roc.b<-roc(dat3$local_6m~predict(obj.b, type="response"), ci=T)## Setting levels: control = 0, case = 1## Setting direction: controls < casesplot(roc.m)
plot(roc.b, add=T, lty=2)
legend(0.4, 0.5, lty=1:2, legend=c("Model M", "Model B"))
두 모델의 AUC가 다른지 테스트하려면, roc.test안에, roc
함수 돌린 결과물을 넣어주면 됩니다.
roc.test(roc.m, roc.b)##
## DeLong's test for two correlated ROC curves
##
## data: roc.m and roc.b
## Z = 0.65203, p-value = 0.5144
## alternative hypothesis: true difference in AUC is not equal to 0
## 95 percent confidence interval:
## -0.02104148 0.04202050
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.8461538 0.8356643이렇게 해서 logistic regression을 R을 수행하는 방법과 ROC 분석을 이용하여 c-index를 계산하는 방법을 알아보았습니다. 다음 시간엔 R을 이용하여 생존분석을 수행하는 방법을 알아보겠습니다.