Saturday, November 7, 2015

Vẽ đường biểu diễn ROC

ROC là viết tắt từ chữ "Receiver Operating Characteristic", thường được dùng để đánh giá một phương pháp hay một mô hình tiên lượng. Đường biểu diễn ROC gồm có trục tung là độ nhậy, và trục hoành là xác suất dương tính giả (tức lấy 1 trừ cho độ đặc hiệu). Do đó, có thể xem đường biểu diễn ROC là một thước đo dung hoà giữa độ nhậy và tỉ lệ dương tính giả của một mô hình tiên lượng.


Cách vẽ đường biểu diễn ROC gồm ba bước như sau:

  •       Bước 1 là tạo ra một mô hình tiên lượng (cũng có thể xem là "rule"), rồi tính xác suất tiên lượng bằng cách dùng mô hình đó; 
  •       Bước 2 là tính độ nhậy (sensitivity) và độ đặc hiệu (specificity) cho mỗi giá trị tiên lượng;
  •       Bước 3 là vẽ đường biểu diễn.


Tất cả 3 bước đó có thể triển khai trong package có tên là pROC.


dat = read.table("~/Google Drive/_QA Book (Vietnamese)/Data/osteo data.txt", header=T, na.strings=".")

dat = na.omit(dat)

> head(dat)
  id gender age wt   bmi fnbmd lsbmd   lsz   lst lstbs fx pfx fall
1 10 Female  62 72 24.06  0.84  1.14  1.51 -0.36  1.27  0   0    0
2 27 Female  64 85 30.48  0.86  1.09  1.13 -0.76  1.04  0   0    0
3 28 Female  76 48 20.50  0.65  0.89 -0.45 -2.45  1.28  0   0    1
4 33 Female  75 52 21.37  0.83  0.84 -0.86 -2.81  1.17  0   0    0
5 37 Female  60 60 23.15  0.79  0.82 -1.19 -3.03  1.10  1   0    0
6 51 Female  68 65 24.46  0.79  1.34  3.20  1.32  1.11  0   0    0

# Ước tính tham số của mô hình hồi qui logistic

model1 = glm(fx ~ age + fnbmd, family=binomial, data=dat)

# Tính xác suất tiên lượng và đưa vào dataframe dat

dat$pred = predict(model1, type="response")


# Tính AUC
library(pROC)
auc = roc(dat$fx, dat$pred)
auc

> auc

Call:
roc.default(response = dat$fx, predictor = dat$pred)

Data: dat$pred in 1465 controls (dat$fx 0) < 530 cases (dat$fx 1).
Area under the curve: 0.6786

# Tính khoảng tin cậy 95%

> ci(auc)
95% CI: 0.6518-0.7053 (DeLong)

# Vẽ đường biểu diễn ROC

plot.roc(smooth(auc), col="blue")


# Có thể tô màu dưới dường biểu diễn

plot(smooth(auc), print.auc=TRUE, auc.polygon=TRUE)




plot(smooth(auc), print.auc=TRUE, auc.polygon=TRUE, grid=c(0.1, 0.2), grid.col=c("green", "red"), max.auc.polygon=TRUE, auc.polygon.col="yellow", print.thres=TRUE)



# Mô hình 2, thêm fall  

model2 = glm(fx ~ age + fnbmd+lstbs+fall, family=binomial, data=dat)

dat$pred2 = predict(model2, type="response")

auc2 = roc(dat$fx, dat$pred2)

> auc2

Call:
roc.default(response = dat$fx, predictor = dat$pred2)

Data: dat$pred2 in 1465 controls (dat$fx 0) < 530 cases (dat$fx 1).
Area under the curve: 0.6942


plot.roc(smooth(auc2), add=TRUE, col="red")

 


# Tính điểm tối ưu
ci.coords(dat$fx, dat$pred2, x="best", input = "threshold", ret=c("specificity", "ppv", "tp"))

> ci.coords(dat$fx, dat$pred2, x="best", input = "threshold", ret=c("specificity", "ppv", "tp"))
  |========================================================================| 100%
95% CI (2000 stratified bootstrap replicates):
                                2.5%      50%    97.5%
threshold best: specificity   0.5113   0.7201   0.7734
threshold best: ppv           0.3644   0.4279   0.4688
threshold best: tp          277.0000 315.0000 418.0000




No comments:

Post a Comment