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
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
testy <- performance(pred,"tpr","fpr")
Sử dụng str() function, testy:
alpha.values: Cutoff
x.values: Specificity or True Negative Rate
y.values: Sensitivity or True Positive Rate
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.
# 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)
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.