Phân tích dữ liệu và ứng dụng - Đánh giá mô hình hồi qui Logistic

Tuan V. Nguyen
Senior Principal Research Fellow, Garvan Institute of Medical Research 
Professor, UNSW School of Public Health and Community Medicine
Professor of Predictive Medicine, University of Technology Sydney
Adj. Professor of Epidemiology and Biostatistics,
School of Medicine Sydney, University of Notre Dame Australia 
Phân tích dữ liệu và ứng dụng | Đại học Dược Hà Nội | 12/6 to 17/6/2019 © Tuan V. Nguyen
Ba tiêu chí cho một mô hình tiên lượng
• Discrimination – phân định 
• Calibration – chính xác 
• Reclassification – tái phân nhóm 
Discrimination – phân định
Hai thước đo độ tin cậy của một mô hình
• Sensitivity – độ nhạy 
• Specificity – độ đặc hiệu 
Độ nhạy 
• Trong số những người mắc bệnh, bao nhiêu % có tiên lượng 
dương tính? 
• Gold standard – mắc bệnh trong thực tế 
• Test result – mô hình tiên lượng 
Độ đặc hiệu 
• Trong số những người không mắc bệnh, bao nhiêu % có tiên 
lượng âm tính 
• Gold standard – mắc bệnh trong thực tế 
• Test result – mô hình tiên lượng 
Ví dụ 
Tiên lượng Bệnh (n=177) Không bệnh (n=81)
+ve (trên 20%) 120 20
-ve (dưới 20%) 57 61
Tổng số 177 81
Ví dụ 
Sensitivity = 120 / 177 = 0.68
Specificity = 61 / 81 = 0.75
Tiên lượng Bệnh (n=177) Không bệnh (n=81)
+ve (trên 20%) 120 20
-ve (dưới 20%) 57 61
Tổng số 177 81
Tỉ lệ dương tính giả (false +ve)
Sensitivity (dương tính thật) = 120 / 177 = 0.68
Specificity = 61 / 81 = 0.75
Dương tính giả = 1 – 0.75 = 0.25
Tiên lượng Bệnh (n=177) Không bệnh (n=81)
+ve (trên 20%) 120 20
-ve (dưới 20%) 57 61
Tổng số 177 81
ROC curve 
• Receiver operating characteristic (ROC) curve 
• Đo lường khả năng phân định (power of discrimination) của một xét 
nghiệm hay mô hình tiên lượng
• Thường dùng cho các kết quả biến liên tục 
• Y-axis (trục tung): true positive (sensitivity)
• X-axis (trục hoành): false positive (1-specificity)
Một ví dụ về một xét nghiệm hoàn hảo
microalbuminuria) the median 24 h microalbuminuria was
40mg/day and the range went from 31 to 60mg/day. Because
of the way renal dysfunction was defined (that is, a 24 h
microalbuminuria 430mg/day), there was no overlapping
between individuals with and without renal dysfunction.
In other words, in this study a 24 h microalbuminuria
430mg/day perfectly discriminates sick and healthy people
(true-negative rate¼ 100%; true-positive rate: 100%; accu-
racy: 100%) (Figure 1).
Figure 2 shows that with UACR there is overlapping
between healthy individuals and patients with renal dysfunc-
tion indicating that this indicator does not perfectly discri-
minate the two groups. Indeed, by the 30mg/g cutoff there is
a large proportion of healthy individuals and patients with
renal dysfunction correctly classified as disease-negative
(true-negative rate or specificity) and disease-positive (true-
positive rate or sensitivity), but also a proportion of healthy
individuals and patients with renal dysfunction incorrectly
classified as disease-positive and negative. Moving the cutoff
to the right (for example, to a UACR of 40mg/g,) decreases
the false-positive rate (higher specificity) but also increases
the false-negative rate (lower sensitivity). Vice versa, moving
the cutoff point to the left (UACR¼ 20mg/g) decreases the
false-negative rate (higher sensitivity) but also increases the
false-positive rate (lower specificity). Thus, the choice of the
cutoff critically affects the diagnostic performance of the test.
ROC curve analysis
Clinical practice commonly demands ‘yes or no’ decisions.
For this reason, we frequently need to convert a continuous
diagnostic test into a dichotomous test. In this vein,
we consider an individual as affected/unaffected by renal
dysfunction depending on whether he/she had a 24 h
microalbuminuria o or 430mg/day. ROC curves analysis
assesses the discrimination performance of quantitative tests
throughout the whole range of their possible values and
helps to identify the optimal cutoff value.
An ROC curve is a graph plotting the combination of
sensitivity (true-positive rate) and the complement to
specificity (that is, 1-specificity, false-positive rate) across a
series of cutoff values covering the whole range of values of a
given disease marker. Because sensitivity and specificity are
both unaffected by the disease prevalence,1 also ROC curve
analysis and accuracy are also independent of the proportion
of the diseased. To construct an ROC curve for UACR, we
consider all possible UACR cutoffs throughout the whole
range of this measurement in patients and controls. By
plotting the pairs of sensitivity and 1-specificity correspond-
ing to each UACR cutoff, we obtain an ROC curve (Figure 3,
dotted line). A test with high discrimination has an ROC
curve approaching the upper left corner of the graph.
Therefore, the closer the ROC plot is to the upper left corner,
the higher the accuracy of the test. The closer the curve to the
diagonal line (also called reference line or chance line) of the
graph, the lower the accuracy of the test.
The overall discrimination performance of a given test is
measured by calculating the area under the ROC curve
(AUC). The AUC may be considered as a global estimate of
diagnostic accuracy. The AUC may take values ranging from
0.5 (no discrimination) to 1 (perfect discrimination). A
test has (at least some) discriminatory power if the 95%
confidence interval of the AUC does not include 0.50. In our
instance, the area under the ROC curve provided by UACR is
0.754 (95% CI: 0.681–0.826) (Figure 3, top panel), a figure
higher than that of diagnostic indifference (AUC: 0.50). An
area of 0.754 implies that in a hypothetical experiment in
which we randomly select pairs of healthy individuals and
24h microalbuminuria (mg/24 h)
6050403020100
True negatives True positives
24 h microalbuminuria cutoff
Patients with renal dysfunction
Healthy subjects
Figure 1 | Scattergram of 24 h microalbuminuria in patients
with renal dysfunction (gray circles) and in healthy subjects
(green circles). Because of the way we defined renal dysfunction,
there is no overlapping between individuals with and without the
disease, indicating that a 24 h microalbuminuria 430mg/day
perfectly discriminates between sick and healthy people.
Urine albumin:creatinine ratio
(mg/g of creatinine)
False negatives
True negatives False positives
True positives
Patients with renal dysfunction
Healthy subjects
6050403020100
Urine albumin:creatinine ratio cutoff
Figure 2 | Scattergram of urine albumin:creatinine ratio
(UACR) in patients with renal dysfunction and in healthy
subjects. Because there is no perfect agreement between urine
albumin:creatinine ratio and 24 h microalbuminuria, a value of
UACR of 30mg/g does not perfectly distinguish between sick and
healthy people. The standard UACR cutoff (30mg/g) and the
additional UACR cutoffs (20 and 40mg/g) are indicated by the
broken lines (see also text).
Kidney International (2009) 76 , 252–256 253
G Tripepi et al.: Diagnostic methods abc o f ep idemio logy
Nhưng trong thực tế 
microalbuminuria) the median 24 h microalbuminuria was
40mg/day and the range went from 31 to 60mg/day. Because
of the way renal dysfunction was defined (that is, a 24 h
microalbuminuria 430mg/day), there was no overlapping
between individuals with and without renal dysfunction.
In other words, in this study a 24 h microalbuminuria
430mg/day perfectly discriminates sick and healthy people
(true-negative rate¼ 100%; true-positive rate: 100%; accu-
racy: 100%) (Figure 1).
Figure 2 shows that with UACR there is overlapping
between healthy individuals and patients with renal dysfunc-
tion indicating that this indicator does not perfectly discri-
minate the two groups. Indeed, by the 30mg/g cutoff there is
a large proportion of healthy individuals and patients with
renal dysfunction correctly classified as disease-negative
(true-negative rate or specificity) and disease-positive (true-
positive rate or sensitivity), but also a proportion of healthy
individuals and patients with renal dysfunction incorrectly
classified as disease-positive and negative. Moving the cutoff
to the right (for example, to a UACR of 40mg/g,) decreases
the false-positive rate (higher specificity) but also increases
the false-negative rate (lower sensitivity). Vice versa, moving
the cutoff point to the left (UACR¼ 20mg/g) decreases the
false-negative rate (higher sensitivity) but also increases the
false-positive rate (lower specificity). Thus, the choice of the
cutoff critically affects the diagnostic performance of the test.
ROC curve analysis
Clinical practice commonly demands ‘yes or no’ decisions.
For this reason, we frequently need to convert a continuous
diagnostic test into a dichotomous test. In this vein,
we consider an individual as affected/unaffected by renal
dysfunction depending on whether he/she had a 24 h
microalbuminuria o or 430mg/day. ROC curves analysis
assesses the discrimination performance of quantitative tests
throughout the whole range of their possible values and
helps to identify the optimal cutoff value.
An ROC curve is a graph plotting the combination of
sensitivity (true-positive rate) and the complement to
specificity (that is, 1-specificity, false-positive rate) across a
series of cutoff values covering the whole range of values of a
given disease marker. Because sensitivity and specificity are
both unaffected by the disease prevalence,1 also ROC curve
analysis and accuracy are also independent of the proportion
of the diseased. To construct an ROC curve for UACR, we
consider all possible UACR cutoffs throughout the whole
range of this measurement in patients and controls. By
plotting the pairs of sensitivity and 1-specificity correspond-
ing to each UACR cutoff, we obtain an ROC curve (Figure 3,
dotted line). A test with high discrimination has an ROC
curve approaching the upper left corner of the graph.
Therefore, the closer the ROC plot is to the upper left corner,
the higher the accuracy of the test. The closer the curve to the
diagonal line (also called reference line or chance line) of the
graph, the lower the accuracy of the test.
The overall discrimination performance of a given test is
measured by calculating the area under the ROC curve
(AUC). The AUC may be considered as a global estimate of
diagnostic accuracy. The AUC may take values ranging from
0.5 (no discrimination) to 1 (perfect discrimination). A
test has (at least some) discriminatory power if the 95%
confidence interval of the AUC does not include 0.50. In our
instance, the area under the ROC curve provided by UACR is
0.754 (95% CI: 0.681–0.826) (Figure 3, top panel), a figure
higher than that of diagnostic indifference (AUC: 0.50). An
area of 0.754 implies that in a hypothetical experiment in
which we randomly select pairs of healthy individuals and
24h microalbuminuria (mg/24 h)
6050403020100
True negatives True positives
24 h microalbuminuria cutoff
Patients with renal dysfunction
Healthy subjects
Figure 1 | Scattergram of 24 h microalbuminuria in patients
with renal dysfunction (gray circles) and in healthy subjects
(green circles). Because of the way we defined renal dysfunction,
there is no overlapping between individuals with and without the
disease, indicating that a 24 h microalbuminuria 430mg/day
perfectly discriminates between sick and healthy people.
Urine albumin:creatinine ratio
(mg/g of creatinine)
False negatives
True negatives False positives
True positives
Patients with renal dysfunction
Healthy subjects
6050403020100
Urine albumin:creatinine ratio cutoff
Figure 2 | Scattergram of urine albumin:creatinine ratio
(UACR) in patients with renal dysfunction and in healthy
subjects. Because there is no perfect agreement between urine
albumin:creatinine ratio and 24 h microalbuminuria, a value of
UACR of 30mg/g does not perfectly distinguish between sick and
healthy people. The standard UACR cutoff (30mg/g) and the
additional UACR cutoffs (20 and 40mg/g) are indicated by the
broken lines (see also text).
Kidney International (2009) 76 , 252–256 253
G Tripepi et al.: Diagnostic methods abc o f ep idemio logy
Ví dụ về tiên lượng ung thư tiền liệt tuyến
PSA level (ng/mL) Sensitivity Specificity
≥ 1 1.00 0.21
≥ 2 1.00 0.48
≥ 3 1.00 0.60
≥ 4 0.99 0.73
≥ 5 0.96 0.76
≥ 6 0.94 0.79
≥ 7 0.90 0.83
≥ 8 0.90 0.88
≥ 9 0.68 0.90
≥ 10 0.54 0.93
≥ 11 0.47 0.94
≥ 12 0.30 0.95
≥ 13 0.23 0.96
≥ 14 0.17 0.97
≥ 15 0.11 0.97
Morgan TO, et al. Age-specific reference ranges for serum prostate specific antigen in black men. N Engl J Med 1996;335:304-310 
ROC curve
!
Diện tích dưới đường biểu diễn (area under the 
curve - AUC)
!
Diễn giải AUC
AUC Meaning
>0.90 Excellent test
0.80 to 0.90 Good
0.70 to 0.80 Fair
0.60 to 0.70 Poor
0.50 to 0.60 Fail
Ý nghĩa thật của AUC 
• Định nghĩa: Probability that a randomly selected pair of healthy 
individual and patient, the test result will be higher in the 
patient than in the healthy individual (xác suất mà một cặp 
bệnh nhân và người bình thường được chọn, và bệnh nhân có 
giá trị tiên lượng cao hơn người bình thường)
• Khó hiểu! 
Một ví dụ đơn giản 
• Chúng ta có 7 người được theo dõi 5 năm 
• 4 người mắc bệnh ung thư tiền liệt tuyến, và giá trị PSA là: 
8, 2, 6, 3
• 3 người không mắc bệnh, giá trị PSA: 
3, 2, 6
Tổ hợp 
• 4 giá trị PSA của 4 bệnh nhân 
• 3 giá trị PSA của 3 người không bệnh
• Tổng số cặp có thể: 4 x 3 = 12
• AUC = số cặp mà PSA bệnh nhân cao hơn PSA của người không 
mắc bệnh chia cho 12. 
Bệnh Không bệnh Chú ý 
8 3 Concordant (8 > 3) 
2 2 Concordant (8 > 2)
6 6 Concordant (8 > 6) 
3
Nếu chúng bắt cặp bệnh nhân 1 với 3 người trong nhóm không bệnh:
Bệnh Không bệnh Chú ý 
8 3 Concordant (8 > 3) 
2 2 Concordant (8 > 2)
6 6 Concordant (8 > 6) 
3
Nếu chúng bắt cặp bệnh nhân 1 với 3 người trong nhóm không bệnh:
Đến bệnh nhân thứ 2 
Bệnh Không bệnh Chú ý 
8 3 Discordant (2 < 3) 
2 2 Tie (2 = 2)
6 6 Discordant (2 < 6) 
3
Bệnh Không bệnh Chú ý 
8 3 Concordant (6 > 3) 
2 2 Concordant (6 > 2)
6 6 Tie (6 = 6) 
3
Đến bệnh nhân thứ 3 
và bệnh nhân thứ 4 
Bệnh Không bệnh Chú ý 
8 3 Tie (3 = 3) 
2 2 Concordant (3 > 2)
6 6 Discordant (3 < 6) 
3
AUC (hay c-index)
• AUC = (concordant + 0.5 ties) / tổng số cặp 
• Chúng ta có:
– 6 cặp concordance
– 3 cặp discordance
– 3 cặp ties (trùng) 
– AUC = (6 + 1.5) / 12 = 0.625
Ví dụ: tính AUC cho mô hình hồi qui logistic
• Mục tiêu: tiên lượng thu nhập trên $50K
• Yếu tố: Age, Education.Years, Occupation, Race, Sex
• Câu hỏi: 4 yếu tố này tiên lượng thu nhập chính xác như thế nào? 
• Trả lời: AUC 
Package "pROC" 
• pROC package 
(https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-12-
77) 
• Có thể dùng để tính AUC và vẽ biểu đồ ROC 
• Các bước cần thiết 
– Xây dựng mô hình m
– Tính giá trị tiên lượng dựa trên tham số của m
– Tính AUC 
# Đọc dữ liệu income 
inc = read.csv("~/Dropbox/_Conferences and Workshops/TDTU 
2018/Datasets/adult.data.csv")
# Fit mô hình hồi qui logistic
m = glm(Income ~ Age + Education.Years + Occupation + Race + Sex, 
family=binomial, data=inc)
# Dùng pROC 
library(pROC)
# Tính xác suất tiên lượng 
pred = predict(m, type="response")
# Tạo ra object roc, tính AUC 
roccurve = roc(inc$Income ~ pred)
auc(roccurve)
ci(roccurve)
plot(roccurve)
ci.thresholds(roccurve)
Specificity
S
en
si
tiv
ity
1.0 0.8 0.6 0.4 0.2 0.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
> auc(roccurve)
Area under the curve: 0.8241
> ci(roccurve)
95% CI: 0.819-0.8291 (DeLong)
> plot.roc(roccurve)
Dùng ModelGood để tìm giá trị 
"tối ưu" 
library(ModelGood)
r = Roc(Inc50 ~ Age + 
Education.Years + Occupation 
+ Race + Sex, data=inc)
plot(r)
click.Roc(r)
Calibration – chính xác 
Đánh giá độ chính xác - calibration
• AUC không phản ảnh mô hình tiên lượng có chính xác hay 
không 
• Chính xác: giá trị tiên lượng gần với giá trị thực tế
• Calibration = phản ảnh độ chính xác 
Residual – phần dư 
• Residual là hiệu số của 
– tình trạng bệnh của một cá nhân, và 
– xác suất tiên lượng (predicted probability) do mô hình tiên lượng 
Residual (cá nhân i) = Yi - Pi
Residual của mô hình hồi qui logistic
Bệnh nhân i Yi Pi Residual 
1 0 2.31 -2.31
2 0 1.91 -1.91
3 1 98.11 1.89
4 1 79.58 20.42
5 0 4.21 -4.21
6 1 98.81 1.19
7 1 64.72 35.28
. . . .
. . . .
Chỉ số Brier 
• Brier score = residual bình phương
B = 1N Yi − Pi( )
2∑
• Brier score càng thấp càng tốt
Residual của mô hình hồi qui logistic
Cá nhân i Yi Pi Residual 
(Yi – Pi)
(Yi – Pi)2
1 0 2.31 -2.31 5.34
2 0 1.91 -1.91 3.65
3 1 98.11 1.89 3.57
4 1 79.58 20.42 416.98
5 0 4.21 -4.21 17.72
6 1 98.81 1.19 1.42
7 1 64.72 35.28 1244.68
. . . . .
. . . . .
Total S
Brier Score = S / N
Một cách thể hiện calibration 
Nhóm nguy cơ
(xác suất) 
Số cá nhân trong 
thực tế 
Số cá nhân tiên 
lượng 
1 (0.00 – 0.20) 2 5
2 (0.21 – 0.40) 5 6
3 (0.41 – 0.60) 10 8
4 (0.61 – 0.80) 16 17
5 (0.81 – 1.00) 52 51
Calibration Curve
Harrell 2001
# Đọc dữ liệu income 
inc = read.csv("~/Dropbox/_Conferences and Workshops/TDTU 
2018/Datasets/adult.data.csv")
# Fit mô hình hồi qui logistic
f = lrm(Inc50 ~ Age + Education.Years + Occupation + Race + Sex, 
data=inc)
f
# Tính xác suất tiên lượng 
pred.logit = predict(f)
phat = 1/(1+exp(-pred.logit))
# Calibration
val.prob(phat, inc$Inc50, m=20, cex=.5)
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
Predicted Probability
A
ct
ua
l P
ro
ba
bi
lit
y
Ideal
Logistic calibration
Nonparametric
Grouped observations
Dxy 
C (ROC) 
R2 
D 
U 
Q 
Brier 
Intercept
Slope 
Emax 
E90 
Eavg 
S:z 
S:p 
0.648
0.824
0.346
0.263
0.000
0.263
0.136
0.000
1.000
0.022
0.010
0.004
0.326
0.744
> f
Logistic Regression Model
lrm(formula = Income ~ Age + Education.Years + Occupation + Race + 
Sex, data = inc)
Model Likelihood Discrimination Rank Discrim. 
Ratio Test Indexes Indexes 
Obs 32561 LR chi2 8569.01 R2 0.346 C 0.824 
<=50K 24720 d.f. 21 g 1.754 Dxy 0.648 
>50K 7841 Pr(> chi2) <0.0001 gr 5.775 gamma 0.649 
max |deriv| 3e-05 gp 0.236 tau-a 0.237 
Brier 0.136 
Dùng package "givitiR"
f = lrm(Inc50 ~ Age + 
Education.Years + Occupation + 
Race + Sex, data=inc)
pred.logit = predict(f)
phat = 1/(1+exp(-pred.logit))
library(givitiR)
cb = givitiCalibrationBelt(o = 
inc$Inc50, e=phat, devel = 
"external") 
plot(cb) 
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
GiViTI Calibration Belt
e
o
Polynomial degree: 4
p-value: <0.001
n: 32561
95% NEVER0.82 - 0.96
80% NEVER0.81 - 0.96
Confidence
level
Under
the bisector
Over
the bisector
# Đọc dữ liệu income 
inc = read.csv("~/Dropbox/_Conferences and Workshops/TDTU 
2018/Datasets/adult.data.csv")
# Fit mô hình hồi qui logistic
f = lrm(Inc50 ~ Age + Education.Years + Occupation + Race + Sex, 
data=inc, subset=1:10000)
# Tính xác suất tiên lượng 
pred.logit = predict(f, inc[1:10000,])
phat = 1/(1+exp(-pred.logit))
# Calibration
val.prob(phat, inc$Inc50[10001:20000], m=20, cex=.5)
> val.prob(phat, inc$Inc50[10001:20000], m=20, cex=.5) 
Dxy C (ROC) R2 D D:Chi-sq 
0.1919587 0.5959793 -2.1099648 -0.9051756 -9050.7558517 
D:p U U:Chi-sq U:p Q 
NA 0.9289382 9291.3817905 0.0000000 -1.8341138 
Brier Intercept Slope Emax E90 
0.2549178 -0.1039013 0.2164229 0.2640170 0.2611439 
Eavg S:z S:p 
0.2323711 133.6106877 0.0000000 
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
Predicted Probability
A
ct
ua
l P
ro
ba
bi
lit
y
Ideal
Logistic calibration
Nonparametric
Grouped observations
Dxy 
C (ROC) 
R2 
D 
U 
Q 
Brier 
Intercept
Slope 
Emax 
E90 
Eavg 
S:z 
S:p 
 0.192
 0.596
 -2.110
 -0.905
 0.929
 -1.834
 0.255
 -0.104
 0.216
 0.264
 0.261
 0.232
133.611
 0.000
Dùng package "ModelGood"
library(ModelGood)
m = glm(Inc50 ~ Age + Education.Years + Occupation + Race + Sex, 
data=inc, family=binomial)
calPlot2(m)
Đánh giá mô hình tiên lượng: Tóm lược 
• Hai tiêu chí để đánh giá một mô hình tiên lượng hồi qui logistic: 
discrimination và calibration
• Discrimination: dùng AUC và ROC 
• Calibration (goodness-of-fit): Brier score và giá trị "predicted-
observed"