Logistic Regression and ROC analysis
Seo Young Park, PhD
\(~\)
R을 이용해서 logistic regression 과 ROC analysis를 수행하는 방법을 알아보도록 하겠습니다.
그런데, 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)벡터에 들어간 element의 갯수가 그 벡터의 길이(length)가 됩니다.
length(aa)## [1] 4length(bb)## [1] 4length(cc)## [1] 4벡터 여러개를 각 행으로 합치거나, 각 열로 합쳐서 행렬을 만들 수 있습니다.
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\\biost\\Dropbox\\Teaching\\AMC Radiology R 2020")
library(plyr)
library(dplyr)
library(tableone)
dat0<-read.csv("example_data.csv", stringsAsFactors = F, nrow=52)
dat1<-dat0 %>% distinct()
length(unique(dat1$id))
dat1$id[which(duplicated(dat1$id))]
dat1[dat1$id==41, ]
dat2<-dat1[-45, ]
length(unique(dat2$id))
summary(dat2)
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")),
CCI.b=factor(mapvalues(CCI, from=c(0, 1, 2, 3, 4, 6), to=c("0or1", "0or1", "2", "3", ">=4", ">=4")),
levels=c("0or1", "2", "3", ">=4")),
RFS=as.double(Recur_date - OP_date)
) %>%
mutate_at(vars(volume_change, adc_1_mean),
list(~as.double(replace(., .=="n/a" | .=="", NA)))) %>%
mutate_at(vars(sex, Recur, local_6m, TNM, CCI, CA19.9.group, TNM.b, CCI.b),
list(~as.factor(.)))
summary(dat3)
save(dat3, file="clean_data.RData")따라서, 위의 코드를 실행해서 얻은 clean_data.RData 파일이 working directory에 저장되어있다면, 간단히 load() 함수를 이용해서 dat3를 불러올 수 있습니다.
load("clean_data.RData")분석에 필요한 패키지 세가지를 불러옵니다.
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(pROC)## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, var먼저, 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)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2883 -1.1422 -0.8995 1.2107 1.4328
##
## 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)
anova(obj, obj0, test="LR")## Analysis of Deviance Table
##
## Model 1: local_6m ~ TNM.b
## Model 2: local_6m ~ 1
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 46 66.946
## 2 49 68.994 -3 -2.0484 0.5624Pr(>Chi) 밑에 나온 값이 두 모델을 비교한 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)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8045 -0.7606 -0.2746 0.7716 1.6493
##
## 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.06386873obj.m<-glm(local_6m~age+TNM.b+log.CA19.9+adc_decrease+apt_1_mean, family="binomial", data=dat3, na.action=na.exclude)
obj.m0<-update(obj.m, .~.-TNM.b)
anova(obj.m, obj.m0, test="LR")## Analysis of Deviance Table
##
## 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
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 40 45.289
## 2 43 45.892 -3 -0.60279 0.8958na.action=na.exclude라는 옵션을 포함하지 않으면 무슨 문제가 생길까요? 사실 회귀계수나 OR, p-value등을 계산하는데에는 아무 영향을 끼치지 않습니다. 그런데 그 모델을 이용해 계산한 predictied probability를 이용해야하는 분석 (예를 들어 ROC analysis)에는 지장을 초래할 수 있습니다. 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)
obj.x<-glm(local_6m~age+TNM.b+log.CA19.9+adc_decrease+apt_1_mean, family="binomial", data=dat3)
predict(obj.m, type="response")## 1 2 3 4 5 6 7 8
## 0.996159313 0.824968098 0.434773172 0.167485627 0.515421686 NA 0.787463979 0.266181849
## 9 10 11 12 13 14 15 16
## 0.684303949 0.865112205 0.918082899 0.526567315 0.364978753 0.080087874 0.046303323 0.455567860
## 17 18 19 20 21 22 23 24
## 0.125818166 0.246216348 0.743047539 0.752886233 0.610323932 0.256621921 0.728282855 0.158503103
## 25 26 27 28 29 30 31 32
## 0.039707425 0.123367250 0.290681592 0.050641847 0.065728466 0.191942046 0.774829899 0.803690718
## 33 34 35 36 37 38 39 40
## 0.317085678 0.623887711 0.346718205 NA 0.310859089 0.441498353 0.782050282 0.087717242
## 41 42 43 44 45 46 47 48
## 0.007887333 0.502961916 0.995700002 0.496823280 0.789630164 0.999999995 0.034392855 0.741003375
## 49 50
## 0.345620997 0.280416282predict(obj.x, type="response")## 1 2 3 4 5 7 8 9
## 0.996159313 0.824968098 0.434773172 0.167485627 0.515421686 0.787463979 0.266181849 0.684303949
## 10 11 12 13 14 15 16 17
## 0.865112205 0.918082899 0.526567315 0.364978753 0.080087874 0.046303323 0.455567860 0.125818166
## 18 19 20 21 22 23 24 25
## 0.246216348 0.743047539 0.752886233 0.610323932 0.256621921 0.728282855 0.158503103 0.039707425
## 26 27 28 29 30 31 32 33
## 0.123367250 0.290681592 0.050641847 0.065728466 0.191942046 0.774829899 0.803690718 0.317085678
## 34 35 37 38 39 40 41 42
## 0.623887711 0.346718205 0.310859089 0.441498353 0.782050282 0.087717242 0.007887333 0.502961916
## 43 44 45 46 47 48 49 50
## 0.995700002 0.496823280 0.789630164 0.999999995 0.034392855 0.741003375 0.345620997 0.280416282roc(dat3$local_6m~predict(obj.m, type="response"))## Setting levels: control = 0, case = 1## Setting direction: controls < cases##
## Call:
## roc.formula(formula = dat3$local_6m ~ predict(obj.m, type = "response"))
##
## Data: predict(obj.m, type = "response") in 26 controls (dat3$local_6m 0) < 22 cases (dat3$local_6m 1).
## Area under the curve: 0.8462roc(dat3$local_6m~predict(obj.x, type="response"))## Error in model.frame.default(formula = dat3$local_6m ~ predict(obj.x, : variable lengths differ (found for 'predict(obj.x, type = "response")')na.action=na.exclude 옵션을 포함하지 않으면, predicted()함수는 결측이 있는 환자는 무시하고 결측이 없는 환자 50명의 predicted probability만으로 이루어진 벡터를 만들어냅니다. 그런데 outcome인 dat3$local_6m은 총 52개의 값이 있는 벡터이니, roc()함수는 이 중 어떤 50개의 값을 써야할지 알 수 없게 되어버리는 거죠. 그렇지만 na.action=na.exclude 옵션을 포함하면, predicted()함수가 결측이 있는 환자 자리에는 NA값을 부여하기 때문에, NA를 포함해서 총 52개의 predicted probability 값이 있는 벡터를 만들어냅니다. 따라서 roc()함수가 어떤 환자가 결측인지 알아차리고 분석에서 제외할 수 있게 되는 것이죠.
\(~\)
회귀분석에서 multivariable model에 들어갈 변수는 어떻게 정하느냐, 하는 것은 간단한 문제가 아니고, 각 연구의 목적과 가설, 자료 특성에 따라 달라져야 하는 문제입니다. 이 주제로 두꺼운 책 한권이 나올만큼 깊고 어려운 문제라서 이 강의에서 짧게 말씀드리기는 어렵습니다. 그렇지만 강조하고 싶은 것은, 임상과학에서 공식처럼 쓰이고 있는 ‘univariate 모델에서 유의한 것 골라낸 후 backward elimination하기’ 같은 방식을 무턱대로 사용하는 것은 지양해야합니다. 그래도 경우에 따라 backward elimination을 해야만 하는 경우도 있긴 합니다. 그래서 R을 이용해서 backward elimination을 수행하는 방법을 보여드리겠습니다.
먼저 후보 변수가 모두 들어간 모델을 먼저 fitting 합니다.
obj.m<-glm(local_6m~age+TNM.b+log.CA19.9+adc_decrease+apt_1_mean, family="binomial", data=dat3, na.action=na.exclude)그런 다음에 step()함수 안에 결과를 넣어줍니다. 그러면 backward elimination 적용한 이후 최종 선택된 모델에 대한 정보를 output으로 내놓게 됩니다. 그것을 일반적인 multivariable model 결과 뽑아내는 방식으로 뽑아내면 됩니다.
obj.b<-step(obj.m, direction = "backward", trace=0)
summary(obj.b)##
## Call:
## glm(formula = local_6m ~ log.CA19.9 + adc_decrease + apt_1_mean,
## family = "binomial", data = dat3, na.action = na.exclude)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8348 -0.7825 -0.2931 0.9107 1.7828
##
## 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 8 9
## 0.98860727 0.84752382 0.39543883 0.24927445 0.40983357 NA 0.62963339 0.18596725 0.66693796
## 10 11 12 13 14 15 16 17 18
## 0.90306187 0.93276152 0.51252422 0.47943068 0.08903005 0.06471112 0.50940013 0.11743565 0.22766609
## 19 20 21 22 23 24 25 26 27
## 0.83507470 0.76907455 0.63959317 0.20409517 0.74595116 0.11883699 0.04628036 0.09424343 0.39604274
## 28 29 30 31 32 33 34 35 36
## 0.09772262 0.08817852 0.30804632 0.64120798 0.81424297 0.32314417 0.53403821 0.35983084 NA
## 37 38 39 40 41 42 43 44 45
## 0.40520016 0.35662241 0.65324177 0.08010164 0.01060379 0.51142674 0.99295449 0.44451878 0.79929045
## 46 47 48 49 50
## 0.99999998 0.03799788 0.75543055 0.40959670 0.31817290\(~\)
ROC analysis
앞에서 na.action 옵션을 설명하느라 roc()함수를 사용하는 방법을 살짝 보여드렸죠. 이 함수 안에 예측하려는 결과변수의 실제 값~ 예측에 쓰는 마커의 값을 넣어주면 AUC를 계산해주는데요, AUC의 95% CI까지 얻고 싶으면 ci=T 옵션을 넣어주면 됩니다.
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)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.b<-roc(dat3$local_6m~predict(obj.b, 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)roc.b$auc## Area under the curve: 0.8357roc.b$ci## 95% CI: 0.7248-0.9466 (DeLong)plot(roc.m)
plot(roc.b)
두개의 ROC curve를 하나의 그래프에 겹쳐 그리고 싶을 때도 있죠. 그럴 때는 add=T 옵션이 유용합니다. 그리고 두 직선을 구분하기 위해 lty 옵션의 값을 바꿔줍니다. (디폴트가 lty=1 입니다.)
plot(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
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.8461538 0.8356643\(~\)
Added value of a new marker by NRI and IDI
NRI (Net Reclassification Index)와 IDI(Integrated Discrimination Index)는 예측 모델에 새로 추가된 설명변수의 유용성을 평가하기 위한 통계량입니다. 성능이 괜찮은 예측 모델의 경우, 어떤 새로운 설명변수를 추가해도 AUC 값이 드라마틱하게 좋아지기는 힘들어서 통계적으로 유의한 AUC 증가를 보기는 어렵다는 문제의식에서 제안된 통계량이죠. 그래서 한동안 엄청난 인기를 누렸으나, 요즘은 NRI와 IDI가 너무 민감해서 별로 예측력이 없는 설명변수도 통계적으로 유의한 것처럼 보이게 한다는 비판을 받아서 점차 쓰지 않는 추세입니다. (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5837334/, https://journals.lww.com/epidem/FullText/2014/03000/Commentary___On_NRI,_IDI,_and__Good_Looking_.17.aspx) 그래도 저널이나 리뷰어가 요구하는 경우에는 어쩔 수 없이 계산을 할 수 밖에 없겠죠. 그래서 R을 이용해서 NRI와 IDI를 계산하는 방법을 보여드리도록 하겠습니다. NRI와 IDI를 계산해주는 패키지가 여럿 있지만 그 중에서 저는 Hmisc 패키지를 이용하는 방법을 보여드리려고 합니다.
참고로, 원래 처음에 제안된 NRI는 예측 모델이 알려주는 예측 확률의 컷오프를 정해서 계산하는 것이었는데, 이것은 컷오프를 어떻게 정하느냐에 따라 NRI 값이 바뀌고, 컷오프 정하는 가장 좋은 방법이 무엇이냐가 불분명해서, 어쩔 수 없이 arbitrary하게 컷오프를 정해서 써야한다는 단점이 있습니다. 그래서 컷오프 없이 계산하는 Continuous NRI라는 방법이 제안되었고, 제가 보여드리는 NRI는 이 Continuous NRI 입니다.
\(~\)
우리의 예제 데이터에서, local_6m을 예측하는 기존의 모델이 age, TNM.b, log.CA19.9 세가지 설명변수를 사용한 모델이라고 하고, 우리가 새로 제시하는 설명변수는 영상에서 나온 adc_decrease라고 합시다. 그래서 이 adc_decrease가 예측을 얼마나 향상시키는지를 NRI, IDI를 통해서 보여주려고 한다고 하겠습니다. 그럼 먼저, 기존의 모델과, 기존의 모델에 adc_decrease를 추가한 새로운 모델 두가지를 fitting 합니다.
obj.a1<-glm(local_6m~age+TNM.b+log.CA19.9, family="binomial", data=dat3, na.action=na.exclude)
summary(obj.a1)##
## Call:
## glm(formula = local_6m ~ age + TNM.b + log.CA19.9, family = "binomial",
## data = dat3, na.action = na.exclude)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4672 -1.0283 -0.5706 1.1552 1.7308
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.03200 1.89650 -1.599 0.1099
## age 0.02125 0.02639 0.805 0.4207
## TNM.b3 0.40074 1.08307 0.370 0.7114
## TNM.b4 0.34424 1.01439 0.339 0.7343
## TNM.b5 1.17459 1.29031 0.910 0.3627
## log.CA19.9 0.37131 0.20409 1.819 0.0689 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 68.994 on 49 degrees of freedom
## Residual deviance: 61.859 on 44 degrees of freedom
## AIC: 73.859
##
## Number of Fisher Scoring iterations: 4obj.a2<-glm(local_6m~age+TNM.b+log.CA19.9+adc_decrease, family="binomial", data=dat3, na.action=na.exclude)
summary(obj.a2)##
## Call:
## glm(formula = local_6m ~ age + TNM.b + log.CA19.9 + adc_decrease,
## family = "binomial", data = dat3, na.action = na.exclude)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7750 -0.8626 -0.3147 0.9449 1.7687
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.94537 2.10789 -0.448 0.65380
## age 0.01149 0.03012 0.381 0.70289
## TNM.b3 0.04200 1.19295 0.035 0.97192
## TNM.b4 -0.44414 1.14407 -0.388 0.69786
## TNM.b5 0.02314 1.52754 0.015 0.98791
## log.CA19.9 0.51813 0.26507 1.955 0.05061 .
## adc_decrease -0.06889 0.02563 -2.688 0.00718 **
## ---
## 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: 49.577 on 41 degrees of freedom
## (2 observations deleted due to missingness)
## AIC: 63.577
##
## Number of Fisher Scoring iterations: 5그리고 두 모델이 예측해주는, 각 환자가 local_6m=1일 확률을 저장합니다.
pred.a1<-predict(obj.a1, type="response")
pred.a2<-predict(obj.a2, type="response")이제 Hmisc 패키지를 불러오고, improveProb()라는 함수를 써서 NRI와 IDI를 계산합니다. 기존의 모델에서 나온 확률을 먼저 쓰고, 콤마 뒤에 새로운 모델에서 나온 확률을 씁니다. 그 다음 실제 결과변수를 써야하는데, 주의할 점은 꼭 0과 1로 코딩된 numeric type 벡터여야 한다는 점입니다.
library(Hmisc)improveProb(pred.a1, pred.a2, y=as.double(as.character(dat3$local_6m)))##
## Analysis of Proportions of Subjects with Improvement in Predicted Probability
##
## Number of events: 22 Number of non-events: 26
##
## Proportions of Positive and Negative Changes in Probabilities
##
## Proportion
## Increase for events (1) 0.818
## Increase for non-events (2) 0.385
## Decrease for events (3) 0.182
## Decrease for non-events (4) 0.615
##
##
## Net Reclassification Improvement
##
## Index SE Z 2P Lower 0.95 Upper 0.95
## NRI (1-3+4-2) 0.867 0.252 3.44 0.000577 0.373 1.361
## NRI for events (1-3) 0.636 0.164 3.87 0.000109 0.314 0.959
## NRI for non-events (4-2) 0.231 0.191 1.21 0.226533 -0.143 0.605
##
##
## Analysis of Changes in Predicted Probabilities
##
## Mean Change in Probability
## Increase for events (sensitivity) 0.0881
## Decrease for non-events (specificity) 0.0722
##
##
## Integrated Discrimination Improvement
## (average of sensitivity and 1-specificity over [0,1];
## also is difference in Yates' discrimination slope)
##
## IDI SE Z 2P Lower 0.95 Upper 0.95
## 0.16026 0.05393 2.97145 0.00296 0.05455 0.26596output에서 NRI Index 0.867이 Continuous NRI 값이고, 2P 밑에 있는 0.000577이 이 NRI에 대한 (이 NRI가 0인지 아닌지 테스트한) p-value입니다. IDI는 0.16026, IDI의 p-value는 0.00296입니다.
복습할 겸, 이 두 모델의 ROC 커브를 그려보고, AUC값의 차이를 테스트해볼까요?
roc.a1<-roc(dat3$local_6m~pred.a1)## Setting levels: control = 0, case = 1## Setting direction: controls < casesroc.a2<-roc(dat3$local_6m~pred.a2)## Setting levels: control = 0, case = 1
## Setting direction: controls < casesplot(roc.a1)
plot(roc.a2, add=T, lty=2)
roc.test(roc.a1, roc.a2)##
## DeLong's test for two correlated ROC curves
##
## data: roc.a1 and roc.a2
## Z = -1.5921, p-value = 0.1114
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.6818182 0.7919580AUC 값의 차이는 유의하지 않게 나왔습니다. 확실히 NRI와 IDI가 더 민감한 것을 볼 수 있습니다.
\(~\)
Optimism-corrected AUC
예측 모델을 만들 때 사용한 데이터를 그대로 이용해서 그 모델을 평가하면, 모델의 성능이 좋은 것처럼 보일 수 밖에 없습니다. 애초에 그 데이터의 signal과 noise를 모두 반영해서 만든 모델이니까요. 모델을 만들 때 사용한 internal data를 그대로 이용해서 모델을 평가하는 것을 apparent validation이라고 하고, apparent validation에서 모델의 성능을 실제보다 좋게 평가하게 되는 경향/정도를 optimism이라고 합니다.
어떤 예측 모델을 객관적으로 평가하려면, 그 데이터와 같은 signal을 가지고 있되 noise는 다른 external data를 가지고 평가를 해야합니다. 그런데 external data가 available하지 않은 경우가 많이 있죠. 그런 경우에는 external data에 모델을 평가했을 때 성능을 bootstrap을 이용해서 추정하는 optimism-corrected estimate을 계산해볼 수 있습니다. 로지스틱 회귀분석을 이용해서 만든 예측모델에 대한 optimism-corrected AUC를 계산해주는 믿을만한 패키지를 아직까지 발견하지 못해서, 제가 만든 코드를 공유해드립니다. 이 코드가 어떻게 구성되어있는지 이해하신 후 각자의 연구에서 필요한 부분을 고쳐서 쓰셔야 합니다. 아래 코드는 예제 데이터에서 local_6m을 예측하는 모델을 age, log.CA19.9, adc_decrease, apt_1_mean 네 변수를 이용한 로지스틱 회귀분석으로 만들었을 경우 그 모델의 optimism-corrected AUC를 계산하는 코드입니다.
obj.a3<-glm(local_6m~age+log.CA19.9+adc_decrease+apt_1_mean, family="binomial", data=dat3, na.action=na.exclude)
roc.a3<-roc(dat3$local_6m~predict(obj.a3, type="response"))
roc.a3$auc## Area under the curve: 0.8444n.iter<-500
bootstrap.performance<-rep(NA, n.iter)
test.performance<-rep(NA, n.iter)
for(i in 1:n.iter){
# print(i)
set.seed(2021+i)
boot.dat<-sample_n(dat3, nrow(dat3), replace=TRUE)
### build model on boot.dat
obj.b<-glm(local_6m~age+log.CA19.9+adc_decrease+apt_1_mean, data=boot.dat, family="binomial", na.action="na.exclude")
#### calculate bootstrap performance
phat.b<-predict(obj.b, newdata=boot.dat, type="response")
bootstrap.performance[i]<-roc(boot.dat$local_6m~phat.b)$auc
#### apply the model to the original data
phat.o<-predict(obj.b, newdata=dat3, type="response")
test.performance[i]<-roc(dat3$local_6m~phat.o)$auc
}#### calculate optimism = bootstrap performance - test performance
optimism<-mean(bootstrap.performance-test.performance)
optimism## [1] 0.04689721### optimism corrected AUC
roc.a3$auc - optimism## [1] 0.7975084모델에 설명변수로 TNM.b를 넣으면 위의 코드가 잘 돌아가지 않게 됩니다. TNM.b에 환자 수가 너무 적은 카테고리가 있어서 생기는 문제입니다. 구체적인 것은 강의에서 말씀드리겠습니다.
\(~\)
위의 코드는 모델에 들어가는 변수의 종류가 데이터와 관계없이 미리 고정되어있는 경우에만 유효합니다. 데이터를 보고 변수를 선택해서 만든 모델이라면, 변수 선택 과정에서 생기는 optimism이 있을 수 밖에 없는데 그런 optimism은 위의 코드에 반영되어있지 않습니다. 따라서, 데이터를 이용해서 변수를 선택해서 모델을 만들었을 경우, 위의 코드는 optimism을 과소추정할 가능성이 있습니다. 그럼 데이터를 이용해서 변수를 선택해서 모델을 만들었을 경우, 변수 선택 과정에서 생기는 optimism까지 포함한 optimism은 어떻게 추정할 수 있을까요? 매뉴얼하게 변수를 선택한 경우는 사실상 optimism을 제대로 추정하는 것이 불가능합니다. 하지만, backward elimination 같은 automatic procedure를 이용한 경우에는 가능합니다. 아래는 backward elimination을 이용해서 변수를 선택한 모델의 optimism-corrected AUC를 추정하는 코드입니다. 아래의 코드는 결측이 있는 자료에서는 작동하지 않기 때문에 결측을 제거한 dat4를 만들어서 모델을 만든 후 optimism-corrected AUC를 계산하는 코드입니다.
dat4<-dat3 %>% filter(!is.na(adc_decrease))
obj.a40<-glm(local_6m~age+log.CA19.9+adc_decrease+apt_1_mean, family="binomial", data=dat4, na.action=na.exclude)
obj.a4<-step(obj.a40, direction="backward", trace=0)
roc.a4<-roc(dat4$local_6m~predict(obj.a4, type="response"))
roc.a4$auc## Area under the curve: 0.8357n.iter<-500
bootstrap.performance<-rep(NA, n.iter)
test.performance<-rep(NA, n.iter)
for(i in 1:n.iter){
# print(i)
set.seed(2021+i)
boot.dat<-sample_n(dat4, nrow(dat4), replace=TRUE)
### build model on boot.dat using backward elimination
obj.b0<-glm(local_6m~age+log.CA19.9+adc_decrease+apt_1_mean, data=boot.dat, family="binomial", na.action="na.exclude")
obj.b<-step(obj.b0, direction="backward", trace=0)
#### calculate bootstrap performance
phat.b<-predict(obj.b, newdata=boot.dat, type="response")
bootstrap.performance[i]<-roc(boot.dat$local_6m~phat.b)$auc
#### apply the model to the original data
phat.o<-predict(obj.b, newdata=dat4, type="response")
test.performance[i]<-roc(dat4$local_6m~phat.o)$auc
}#### calculate optimism = bootstrap performance - test performance
optimism<-mean(bootstrap.performance-test.performance)
optimism## [1] 0.05126268### optimism corrected AUC
roc.a4$auc - optimism## [1] 0.7844017