Chuẩn bị data

Data tải từ trang web Kaggle về máy tính. Bộ dữ liệu diabetes có nhiều biến số. Trong đó, tôi sử dụng biến số Glucose (mg/dL) và Outcome (chẩn đoán diabetes) để xây dựng ROC và tính các thông số liên quan như AUC, sensitivity, specificity, cutoff points …

library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages -------------------------------- tidyverse 1.3.1 --
v ggplot2 3.3.5     v purrr   0.3.4
v tibble  3.1.4     v dplyr   1.0.7
v tidyr   1.1.4     v stringr 1.4.0
v readr   2.0.2     v forcats 0.5.1
-- Conflicts ----------------------------------- tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
Warning message:
package ‘ROCR’ was built under R version 4.1.2 
library(readr)
library(readxl)
library(ROCR)
library(OptimalCutpoints)
df <- read.csv("E:/Statistics/R/Data/diabetes.csv")
head(df)

Kiểm tra các biến số của dataset xem có những NA (missing) hay không

sapply(df,function(x) sum(is.na(x)))
             Pregnancies                  Glucose 
                       0                        0 
           BloodPressure            SkinThickness 
                       0                        0 
                 Insulin                      BMI 
                       0                        0 
DiabetesPedigreeFunction                      Age 
                       0                        0 
                 Outcome 
                       0 

Như vậy là dataset đã hoàn hảo để phân tích. Tiếp tục xem các biến số có các giá trị thế nào.

sapply(df, function(x) length(unique(x)))
             Pregnancies                  Glucose 
                      17                      136 
           BloodPressure            SkinThickness 
                      47                       51 
                 Insulin                      BMI 
                     186                      248 
DiabetesPedigreeFunction                      Age 
                     517                       52 
                 Outcome 
                       2 

Giờ nhìn qua tóm tắt data mình đang có

sapply(df, function(x) summary(x))
        Pregnancies  Glucose BloodPressure SkinThickness   Insulin
Min.       0.000000   0.0000       0.00000       0.00000   0.00000
1st Qu.    1.000000  99.0000      62.00000       0.00000   0.00000
Median     3.000000 117.0000      72.00000      23.00000  30.50000
Mean       3.845052 120.8945      69.10547      20.53646  79.79948
3rd Qu.    6.000000 140.2500      80.00000      32.00000 127.25000
Max.      17.000000 199.0000     122.00000      99.00000 846.00000
             BMI DiabetesPedigreeFunction      Age   Outcome
Min.     0.00000                0.0780000 21.00000 0.0000000
1st Qu. 27.30000                0.2437500 24.00000 0.0000000
Median  32.00000                0.3725000 29.00000 0.0000000
Mean    31.99258                0.4718763 33.24089 0.3489583
3rd Qu. 36.60000                0.6262500 41.00000 1.0000000
Max.    67.10000                2.4200000 81.00000 1.0000000

Đường cong ROC


pred <- prediction(df$Glucose, df$Outcome )
perf <- performance(pred,"tpr","fpr")
plot(perf,col="black")
abline(a=0, b=1, col="#8AB63F")

TPR = Sensitivity và FPR = 1 - Sensitivity và TNR = Specificity.

Đường thẳng màu xanh để xác định trực quan diện tích dưới đường thằng này bằng 0.5. Đường cong nằm trên đường thẳng này có AUC lớn hơn 0.5. Nếu gần 0.5 thì khả năng của xét nghiệm (predictor) không cao, gần 1 thì khả năng (accurarcy) của xét nghiệm cao. AUC không chỉ điểm cutoff cho phép chúng ta chẩn đoán bệnh khi giá trị biến số xét nghiệm nằm trên hoặc dưới điểm cutoff này. Trong trường hợp này là Glucose.

UAC trong data này là ~ 0.79. Có nghĩa là sử dụng Glucose để chẩn đoán bệnh diabetes có độ chính xác là 79%.

auc<- performance( pred,  c("auc"))
unlist(slot(auc , "y.values"))
[1] 0.7881306

Điểm cutoff của biến số Glucose


testy <- performance(pred,"tpr","fpr")

Sử dụng str() function, testy:

plot(testy@alpha.values[[1]], testy@x.values[[1]], type='n',  
     xlab='Cutoff points of the Glucose', 
     ylab='sensitivity or specificity')
lines(testy@alpha.values[[1]], testy@y.values[[1]], 
      type='s', col="#1A425C", lwd=2)
lines(testy@alpha.values[[1]], 1-testy@x.values[[1]], 
      type='s', col="#8AB63F", lwd=2)
legend(1.11,.85, c('sensitivity', 'specificity'), 
       lty=c(1,1), col=c("#1A425C", "#8AB63F"), cex=.9, bty='n')

Khi sensitivity giảm thì specificity tăng và cutpoint là điểm khả dĩ chấp nhận về cả 2 chỉ số này. Không thể nói điểm cutoff nào là tốt nhất. Điều đó tùy thuộc vào ý nghĩa sensitivity hay specificity quan trọng như thế nào trong việc chẩn đoán một bệnh cụ thể nào đó. Từ đó người ta xác định điểm cutoff trong sự ưu tiên cho tiêu chí nào.

OptimalCutpoints package

# library(OptimalCutpoints)

optimal.cutpoint.Youden <- optimal.cutpoints(X = "Glucose", status = "Outcome", tag.healthy = 0, 
methods = "Youden", data = df, pop.prev = NULL,  
control = control.cutpoints(), ci.fit = FALSE, conf.level = 0.95, trace = FALSE)

summary(optimal.cutpoint.Youden)

Call:
optimal.cutpoints.default(X = "Glucose", status = "Outcome", 
    tag.healthy = 0, methods = "Youden", data = df, pop.prev = NULL, 
    control = control.cutpoints(), ci.fit = FALSE, conf.level = 0.95, 
    trace = FALSE)

Area under the ROC curve (AUC):  0.788 (0.755, 0.822) 

CRITERION: Youden
Number of optimal cutoffs: 1

                     Estimate
cutoff            124.0000000
Se                  0.7014925
Sp                  0.7320000
PPV                 0.5838509
NPV                 0.8206278
DLR.Positive        2.6175095
DLR.Negative        0.4077971
FP                134.0000000
FN                 80.0000000
Optimal criterion   0.4334925
#plot(optimal.cutpoint.Youden)

Điểm cutoff chẩn đoán diabetes trong data này là 124. Nếu Glucose lớn hơn 124 mg/dL là chẩn đoán tiểu đường, < 124 là không tiểu đường.

Tương ứng với điểm cutoff này là Sensitivity = 70%, specificity = 73% và các giá trị khác như PPV, NPV, FP, FN.

plot(optimal.cutpoint.Youden)

Tính AUC với Mann-Whitney U test

Wilcoxon hoặc Kruskall-Wallis test cho AUC tương đương.

wt <-wilcox.test(data=df, df$Glucose ~ df$Outcome)
1 - wt$statistic/(sum(df$Outcome==1)*sum(df$Outcome==0))
        W 
0.7881306 

p-value của the Mann-Whitney U test cho chúng ta biết sự khác biệt của AUC so với 0.5 có ý nghĩa thống kê hay không.

wt <-wilcox.test(data=df, df$Glucose ~ df$Outcome)
wt$p.value
[1] 1.200727e-39

p < 0.001, cho thấy Glucose có thể là chỉ số tin cậy để chẩn đoán bệnh tiểu đường, với điểm cutoff là 124 mg/dL.


LS0tDQp0aXRsZTogIkdp4bqjaSBUaMOtY2ggxJDGsOG7nW5nIENvbmcgUk9DIFRyb25nIFjDqXQgTmdoaeG7h20gQ2jhuqluIMSQb8OhbiINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNCiMjIENodeG6qW4gYuG7iyBkYXRhDQoNCkRhdGEgdOG6o2kgdOG7qyB0cmFuZyB3ZWIgS2FnZ2xlIHbhu4EgbcOheSB0w61uaC4gQuG7mSBk4buvIGxp4buHdSBkaWFiZXRlcyBjw7Mgbmhp4buBdSBiaeG6v24gc+G7kS4gVHJvbmcgxJHDsywgdMO0aSBz4butIGThu6VuZyBiaeG6v24gc+G7kSBHbHVjb3NlIChtZy9kTCkgdsOgIE91dGNvbWUgKGNo4bqpbiDEkW/DoW4gZGlhYmV0ZXMpIMSR4buDIHjDonkgZOG7sW5nIFJPQyB2w6AgdMOtbmggY8OhYyB0aMO0bmcgc+G7kSBsacOqbiBxdWFuIG5oxrAgQVVDLCBzZW5zaXRpdml0eSwgc3BlY2lmaWNpdHksIGN1dG9mZiBwb2ludHMgLi4uIA0KDQpgYGB7cn0NCmxpYnJhcnkodGlkeXZlcnNlKQ0KbGlicmFyeShyZWFkcikNCmxpYnJhcnkocmVhZHhsKQ0KbGlicmFyeShST0NSKQ0KbGlicmFyeShPcHRpbWFsQ3V0cG9pbnRzKQ0KZGYgPC0gcmVhZC5jc3YoIkU6L1N0YXRpc3RpY3MvUi9EYXRhL2RpYWJldGVzLmNzdiIpDQpoZWFkKGRmKQ0KYGBgDQoNCktp4buDbSB0cmEgY8OhYyBiaeG6v24gc+G7kSBj4bunYSBkYXRhc2V0IHhlbSBjw7Mgbmjhu69uZyBOQSAobWlzc2luZykgaGF5IGtow7RuZw0KDQpgYGB7cn0NCnNhcHBseShkZixmdW5jdGlvbih4KSBzdW0oaXMubmEoeCkpKQ0KDQoNCmBgYA0KDQpOaMawIHbhuq15IGzDoCBkYXRhc2V0IMSRw6MgaG/DoG4gaOG6o28gxJHhu4MgcGjDom4gdMOtY2guIFRp4bq/cCB04bulYyB4ZW0gY8OhYyBiaeG6v24gc+G7kSBjw7MgY8OhYyBnacOhIHRy4buLIHRo4bq/IG7DoG8uDQoNCmBgYHtyfQ0Kc2FwcGx5KGRmLCBmdW5jdGlvbih4KSBsZW5ndGgodW5pcXVlKHgpKSkNCmBgYA0KDQpHaeG7nSBuaMOsbiBxdWEgdMOzbSB04bqvdCBkYXRhIG3DrG5oIMSRYW5nIGPDsw0KDQpgYGB7cn0NCnNhcHBseShkZiwgZnVuY3Rpb24oeCkgc3VtbWFyeSh4KSkNCmBgYA0KDQojIyAqKsSQxrDhu51uZyBjb25nIFJPQyoqIA0KDQpgYGB7cn0NCiNsaWJyYXJ5KFJPQ1IpDQpwcmVkIDwtIHByZWRpY3Rpb24oZGYkR2x1Y29zZSwgZGYkT3V0Y29tZSApDQpwZXJmIDwtIHBlcmZvcm1hbmNlKHByZWQsInRwciIsImZwciIpDQpwbG90KHBlcmYsY29sPSJibGFjayIpDQphYmxpbmUoYT0wLCBiPTEsIGNvbD0iIzhBQjYzRiIpDQpgYGANCg0KVFBSID0gU2Vuc2l0aXZpdHkgdsOgIEZQUiA9IDEgLSBTZW5zaXRpdml0eSB2w6AgVE5SID0gU3BlY2lmaWNpdHkuDQoNCsSQxrDhu51uZyB0aOG6s25nIG3DoHUgeGFuaCDEkeG7gyB4w6FjIMSR4buLbmggdHLhu7FjIHF1YW4gZGnhu4duIHTDrWNoIGTGsOG7m2kgxJHGsOG7nW5nIHRo4bqxbmcgbsOgeSBi4bqxbmcgMC41LiDEkMaw4budbmcgY29uZyBu4bqxbSB0csOqbiDEkcaw4budbmcgdGjhurNuZyBuw6B5IGPDsyBBVUMgbOG7m24gaMahbiAwLjUuIE7hur91IGfhuqduIDAuNSB0aMOsIGto4bqjIG7Eg25nIGPhu6dhIHjDqXQgbmdoaeG7h20gKHByZWRpY3Rvcikga2jDtG5nIGNhbywgZ+G6p24gMSB0aMOsIGto4bqjIG7Eg25nIChhY2N1cmFyY3kpIGPhu6dhIHjDqXQgbmdoaeG7h20gY2FvLiBBVUMga2jDtG5nIGNo4buJIMSRaeG7g20gY3V0b2ZmIGNobyBwaMOpcCBjaMO6bmcgdGEgY2jhuqluIMSRb8OhbiBi4buHbmgga2hpIGdpw6EgdHLhu4sgYmnhur9uIHPhu5EgeMOpdCBuZ2hp4buHbSBu4bqxbSB0csOqbiBob+G6t2MgZMaw4bubaSDEkWnhu4NtIGN1dG9mZiBuw6B5LiBUcm9uZyB0csaw4budbmcgaOG7o3AgbsOgeSBsw6AgR2x1Y29zZS4NCg0KVUFDIHRyb25nIGRhdGEgbsOgeSBsw6AgXH4gMC43OS4gQ8OzIG5naMSpYSBsw6Agc+G7rSBk4bulbmcgR2x1Y29zZSDEkeG7gyBjaOG6qW4gxJFvw6FuIGLhu4duaCBkaWFiZXRlcyBjw7MgxJHhu5kgY2jDrW5oIHjDoWMgbMOgIDc5JS4NCg0KYGBge3J9DQphdWM8LSBwZXJmb3JtYW5jZSggcHJlZCwgIGMoImF1YyIpKQ0KdW5saXN0KHNsb3QoYXVjICwgInkudmFsdWVzIikpDQpgYGANCg0KIyMgKirEkGnhu4NtIGN1dG9mZiBj4bunYSBiaeG6v24gc+G7kSBHbHVjb3NlKioNCg0KYGBge3J9DQojIGxpYnJhcnkoUk9DUikNCnRlc3R5IDwtIHBlcmZvcm1hbmNlKHByZWQsInRwciIsImZwciIpDQpgYGANCg0KU+G7rSBk4bulbmcgKnN0cigpKiBmdW5jdGlvbiwgdGVzdHk6DQoNCi0gICBhbHBoYS52YWx1ZXM6IEN1dG9mZg0KDQotICAgeC52YWx1ZXM6IFNwZWNpZmljaXR5IG9yIFRydWUgTmVnYXRpdmUgUmF0ZQ0KDQotICAgeS52YWx1ZXM6IFNlbnNpdGl2aXR5IG9yIFRydWUgUG9zaXRpdmUgUmF0ZQ0KDQpgYGB7cn0NCnBsb3QodGVzdHlAYWxwaGEudmFsdWVzW1sxXV0sIHRlc3R5QHgudmFsdWVzW1sxXV0sIHR5cGU9J24nLCAgDQogICAgIHhsYWI9J0N1dG9mZiBwb2ludHMgb2YgdGhlIEdsdWNvc2UnLCANCiAgICAgeWxhYj0nc2Vuc2l0aXZpdHkgb3Igc3BlY2lmaWNpdHknKQ0KbGluZXModGVzdHlAYWxwaGEudmFsdWVzW1sxXV0sIHRlc3R5QHkudmFsdWVzW1sxXV0sIA0KICAgICAgdHlwZT0ncycsIGNvbD0iIzFBNDI1QyIsIGx3ZD0yKQ0KbGluZXModGVzdHlAYWxwaGEudmFsdWVzW1sxXV0sIDEtdGVzdHlAeC52YWx1ZXNbWzFdXSwgDQogICAgICB0eXBlPSdzJywgY29sPSIjOEFCNjNGIiwgbHdkPTIpDQpsZWdlbmQoMS4xMSwuODUsIGMoJ3NlbnNpdGl2aXR5JywgJ3NwZWNpZmljaXR5JyksIA0KICAgICAgIGx0eT1jKDEsMSksIGNvbD1jKCIjMUE0MjVDIiwgIiM4QUI2M0YiKSwgY2V4PS45LCBidHk9J24nKQ0KYGBgDQoNCktoaSBzZW5zaXRpdml0eSBnaeG6o20gdGjDrCBzcGVjaWZpY2l0eSB0xINuZyB2w6AgY3V0cG9pbnQgbMOgIMSRaeG7g20ga2jhuqMgZMSpIGNo4bqlcCBuaOG6rW4gduG7gSBj4bqjIDIgY2jhu4kgc+G7kSBuw6B5LiBLaMO0bmcgdGjhu4MgbsOzaSDEkWnhu4NtIGN1dG9mZiBuw6BvIGzDoCB04buRdCBuaOG6pXQuIMSQaeG7gXUgxJHDsyB0w7l5IHRodeG7mWMgdsOgbyDDvSBuZ2jEqWEgc2Vuc2l0aXZpdHkgaGF5IHNwZWNpZmljaXR5IHF1YW4gdHLhu41uZyBuaMawIHRo4bq/IG7DoG8gdHJvbmcgdmnhu4djIGNo4bqpbiDEkW/DoW4gbeG7mXQgYuG7h25oIGPhu6UgdGjhu4MgbsOgbyDEkcOzLiBU4burIMSRw7MgbmfGsOG7nWkgdGEgeMOhYyDEkeG7i25oIMSRaeG7g20gY3V0b2ZmIHRyb25nIHPhu7EgxrB1IHRpw6puIGNobyB0acOqdSBjaMOtIG7DoG8uDQoNCiMjIyBPcHRpbWFsQ3V0cG9pbnRzIHBhY2thZ2UNCg0KYGBge3J9DQojIGxpYnJhcnkoT3B0aW1hbEN1dHBvaW50cykNCg0Kb3B0aW1hbC5jdXRwb2ludC5Zb3VkZW4gPC0gb3B0aW1hbC5jdXRwb2ludHMoWCA9ICJHbHVjb3NlIiwgc3RhdHVzID0gIk91dGNvbWUiLCB0YWcuaGVhbHRoeSA9IDAsIA0KbWV0aG9kcyA9ICJZb3VkZW4iLCBkYXRhID0gZGYsIHBvcC5wcmV2ID0gTlVMTCwgIA0KY29udHJvbCA9IGNvbnRyb2wuY3V0cG9pbnRzKCksIGNpLmZpdCA9IEZBTFNFLCBjb25mLmxldmVsID0gMC45NSwgdHJhY2UgPSBGQUxTRSkNCg0Kc3VtbWFyeShvcHRpbWFsLmN1dHBvaW50LllvdWRlbikNCg0KYGBgDQrEkGnhu4NtIGN1dG9mZiBjaOG6qW4gxJFvw6FuIGRpYWJldGVzIHRyb25nIGRhdGEgbsOgeSBsw6AgMTI0LiBO4bq/dSBHbHVjb3NlIGzhu5tuIGjGoW4gMTI0IG1nL2RMIGzDoCBjaOG6qW4gxJFvw6FuIHRp4buDdSDEkcaw4budbmcsIDwgMTI0IGzDoCBraMO0bmcgdGnhu4N1IMSRxrDhu51uZy4gDQoNClTGsMahbmcg4bupbmcgduG7m2kgxJFp4buDbSBjdXRvZmYgbsOgeSBsw6AgU2Vuc2l0aXZpdHkgPSA3MCUsIHNwZWNpZmljaXR5ID0gNzMlIHbDoCBjw6FjIGdpw6EgdHLhu4sga2jDoWMgbmjGsCBQUFYsIE5QViwgRlAsIEZOLg0KDQoNCmBgYHtyfQ0KcGxvdChvcHRpbWFsLmN1dHBvaW50LllvdWRlbikNCmBgYA0KDQoNCiMjIFTDrW5oIEFVQyB24bubaSBNYW5uLVdoaXRuZXkgVSB0ZXN0DQoNCldpbGNveG9uIGhv4bq3YyBLcnVza2FsbC1XYWxsaXMgdGVzdCBjaG8gQVVDIHTGsMahbmcgxJHGsMahbmcuDQoNCmBgYHtyfQ0Kd3QgPC13aWxjb3gudGVzdChkYXRhPWRmLCBkZiRHbHVjb3NlIH4gZGYkT3V0Y29tZSkNCjEgLSB3dCRzdGF0aXN0aWMvKHN1bShkZiRPdXRjb21lPT0xKSpzdW0oZGYkT3V0Y29tZT09MCkpDQpgYGANCg0KcC12YWx1ZSBj4bunYSB0aGUgTWFubi1XaGl0bmV5IFUgdGVzdCBjaG8gY2jDum5nIHRhIGJp4bq/dCBz4buxIGtow6FjIGJp4buHdCBj4bunYSBBVUMgc28gduG7m2kgMC41IGPDsyDDvSBuZ2jEqWEgdGjhu5FuZyBrw6ogaGF5IGtow7RuZy4NCg0KYGBge3J9DQp3dCA8LXdpbGNveC50ZXN0KGRhdGE9ZGYsIGRmJEdsdWNvc2UgfiBkZiRPdXRjb21lKQ0Kd3QkcC52YWx1ZQ0KDQpgYGANCg0KcCA8IDAuMDAxLCBjaG8gdGjhuqV5IEdsdWNvc2UgY8OzIHRo4buDIGzDoCBjaOG7iSBz4buRIHRpbiBj4bqteSDEkeG7gyBjaOG6qW4gxJFvw6FuIGLhu4duaCB0aeG7g3UgxJHGsOG7nW5nLCB24bubaSDEkWnhu4NtIGN1dG9mZiBsw6AgMTI0IG1nL2RMLg0KDQoNCi0tLS0NCg==