Phân tích dữ liệu và ứng dụng - Mô hình hồi qui logistic (logistic regression)

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
Mô hình hồi qui logistic (logistic regression)
• Ví dụ dẫn nhập 
• Khái niệm odds, logit, và mô hình hồi qui logistic
• Ước tính và R 
The Challenger shuttle disaster
Flight Temp Damage
STS-1 66 0
STS-2 70 1
STS-3 69 0
STS-4 80
STS-5 68 0
STS-6 67 0
STS-7 72 0
STS-8 73 0
STS-9 70 0
STS 41B 57 1
STS 41C 63 1
STS 41D 70 1
STS 41G 78 0
STS 51A 67 0
STS 51C 53 1
STS 51D 67 0
Flight Temp Damage
STS 51B 75 0
STS 51G 70 0
STS 51F 81 0
STS 51I 76 0
STS 51J 79 0
STS 61A 75 1
STS 61B 76 0
STS 61C 58 1
Temp = c(66, 70, 69, 80, 68, 67, 72, 73, 70, 57, 63, 
70, 78, 67, 53, 67, 75, 70, 81, 76, 79, 75, 76, 58)
Damage = c(0, 1, 0, ., 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 
0, 0, 0, 0, 0, 1, 0, 1) 
Nghiên cứu bệnh tiểu đường (type 2 diabetes)
• Nghiên cứu cắt ngang trên 3155 cá nhân
• Outcome: chẩn đoán tiểu đường
• Yếu tố nguy cơ: tuổi, giới tính, huyết áp, vòng eo, vòng mông, cân nặng, 
tỉ trọng cơ thể (BMI)
• Câu hỏi 1: BMI có liên quan đến bệnh tiểu đường? 
• Câu hỏi 2: Yếu tố nào có liên quan đến tiểu đường? 
• Câu hỏi 2: Có thể xây dựng mô hình tiên lượng nguy cơ mắc bệnh? 
> db = read.csv("~/Dropbox/_Conferences and Workshops/Dai hoc Duoc 6-
2019/Datasets/Diabetes data.csv", header=T)
> head(db)
id age gender height weight waist hip sysbp diabp active hypertension
1 1 76 Female 163 53 90 93 160 90 0 1
2 1 40 Female 149 51 74 94 100 60 0 0
3 1 51 Female 151 55 91 100 120 80 0 0
4 1 43 Female 158 62 78 96 120 80 1 0
5 2 72 Female 148 47 91 95 130 60 1 0
6 2 44 Male 155 48 69 86 120 80 0 0
bmi whr diabetes
1 19.95 0.97 IFG
2 22.97 0.79 Normal
3 24.12 0.91 Normal
4 24.84 0.81 Normal
5 21.46 0.96 IFG
6 19.98 0.80 Normal
Gian lận thẻ tín dụng (credit card)
• Nghiên cứu cắt ngang trên 284807 transactions
• Outcome: gian lận (yes / no)
• Yếu tố nguy cơ: Time, Amount, V1-V28
• Câu hỏi: Có thể xây dựng mô hình tiên lượng gian lận 
> cc = read.csv("~/Dropbox/_Conferences and Workshops/Dai hoc Duoc 6-
2019/Datasets/Credit card data.csv", header=T)
> head(cc, 3)
Time V1 V2 V3 V4 V5 V6
1 0 -1.359807 -0.07278117 2.5363467 1.3781552 -0.33832077 0.46238778
2 0 1.191857 0.26615071 0.1664801 0.4481541 0.06001765 -0.08236081
3 1 -1.358354 -1.34016307 1.7732093 0.3797796 -0.50319813 1.80049938
V7 V8 V9 V10 V11 V12
1 0.23959855 0.09869790 0.3637870 0.09079417 -0.5515995 -0.61780086
2 -0.07880298 0.08510165 -0.2554251 -0.16697441 1.6127267 1.06523531
3 0.79146096 0.24767579 -1.5146543 0.20764287 0.6245015 0.06608369
V13 V14 V15 V16 V17 V18
1 -0.9913898 -0.3111694 1.4681770 -0.4704005 0.2079712 0.02579058
2 0.4890950 -0.1437723 0.6355581 0.4639170 -0.1148047 -0.18336127
3 0.7172927 -0.1659459 2.3458649 -2.8900832 1.1099694 -0.12135931
V19 V20 V21 V22 V23 V24
1 0.403993 0.25141210 -0.01830678 0.2778376 -0.1104739 0.06692808
2 -0.145783 -0.06908314 -0.22577525 -0.6386720 0.1012880 -0.33984648
3 -2.261857 0.52497973 0.24799815 0.7716794 0.9094123 -0.68928096
V25 V26 V27 V28 Amount Class
1 0.1285394 -0.1891148 0.133558377 -0.02105305 149.62 0
2 0.1671704 0.1258945 -0.008983099 0.01472417 2.69 0
3 -0.3276418 -0.1390966 -0.055352794 -0.05975184 378.66 0
Đặc tính của các nghiên cứu 
• Outcome (dependent) variable: biến nhị phân (binary variable), chỉ có 
2 giá trị
• Predictor (independent) variables: đa dạng (nhị phân, biến liên tục) 
Không thể dùng mô hình hồi qui tuyến tính!
Ứng dụng của mô hình hồi qui logistic
• Mô tả mối liên quan giữa biến outcome và biến tiên lượng 
• Kiểm soát các biến nhiễu (Controlling for confounders) 
• Phát triển mô hình tiên lượng (Developing prognostic models)
Ông "tổ" của mô hình hồi qui logistic
Professor David R. Cox
Imperial College, London
1970
Khi nào cần sử dụng mô hình hồi qui logistic
• Logistic regression: 
– outcome là biến phân loại (thường có 2 giá trị yes/no) 
– biến tiên lượng có thể là biến phân loại hay liên tục
• Mô hình hồi qui tuyến tính (Linear regression) 
– biến outcome là biến liên tục 
– biến tiên lượng có thể là biến phân loại hay liên tục
Vài khái niệm cơ bản 
Risk, probability và odds
• Risk: probability (P) of an event [during a period] – xác suất của 
một biến cố trong một thời gian 
• Odds: xác suất biến cố xảy ra chia cho xác suất biến cố không xảy 
ra: 
• n =5 bệnh nhân, 1 bệnh nhân bị đột quị: 
P = 1/ 5 = 0.20
Odds = 0.2 / 0.8 = 0.25 
Odds = P
1− P
Probability và odds
• P = 1/5 = 0.2 or 20%
• Odds = (P) / (1-P)
• Odds = 0.2 / 0.8 hay 1:4
Probability, odds, và logit
• Probability: từ 0 đến 1
• Odds: biến liên tục 
– Khi P = 0.5, odds = 1
• Logit = log odds 
logit p( ) = log p1− p
"
#
$
%
&
'
Mô hình hồi qui logistic dựa trên logit
• Gọi X là biến tiên lượng 
• Gọi P là xác suất của một biến cố (outcome) 
• Mô hình hồi qui logistic phát biểu rằng: 
logit p( ) =α +βX
log p1− p
"
#
$
%
&
'=α +βXhay
Mô hình hồi qui logistic
Điều này cũng có nghĩa là: 
log p1− p
"
#
$
%
&
'=α +βX
p = e
α+βX
1+ eα+βX
Mối liên quan giữa X, p và logit(p)Logistic Regression Model 
0
1
x
P
(x
)
P (x) = exp[0+1x]1+exp[0+1x]log[
P (x)
1P (x) ] = 0 + 1x
!
6
!
4
!
2
0
2
4
6
8
x
lo
g 
[ P
(x
) 
/ (
 1
 !
 P
(x
) 
) 
]
linear form nonlinear form
31
log p1− p
"
#
$
%
&
' α +βX p e
α+βX
1+ eα+βX
Ý nghĩa của tham số mô hình logistic
• α là log odds của biến outcome khi X = 0 
• β là log odds ratio (tỉ số) liên quan với một đơn vị tăng của X 
• Odds ratio = exp(β )
log p1− p
"
#
$
%
&
'=α +βX
Nghiêm chỉnh hơn!
• Mô hình hồi qui logistic:
π i = Pr Yi =1| Xi = xi( ) =
exp β0 +βi xi( )
1+ exp β0 +βi xi( )
• Hay viết dưới dạng logit
logit π i( ) = log
π i
1−π i
!
"
#
$
%
&= β0 +β1xi1 +β2xi2 +...
Giả định mô hình hồi qui logistic 
• Mô hình cung cấp một sự "đại diện" tiêu biểu giữa outcome và X 
• Outcomes độc lập với nhau
• Biến tiên lượng không có sai số ngẫu nhiên
Advantages of logistic regression model
• Xác suất của outcome có thể thay đổi với giá trị của biến tiên 
lượng 
• Hệ số có thể diễn giải như là log odds ratio 
• Có thể áp dụng cho nhiều mô hình nghiên cứu
• Nhiều software có thể dùng để ước tính tham số 
Ước tính tham số 
• The maximum likelihood estimator (MLE) for (β0, β1) is obtained by 
finding ( ) that maximizes 
L β0,β1( ) = π iyi
i=1
N
∏ 1−π i( )
ni−yi
=
exp yi β0 +β1xi( )( )
1+ exp β0 +β1xi( )i=1
N
∏
• This is implemented in R program called “glm” and “lrm” 
Hàm glm trong base R 
• Công thức chung: 
m = glm(outcome ~ riskfactor, family = binomial)
outcome có giá trị (0, 1)
riskfactor – bất cứ biến nào 
• Có có khoảng tin cậy 95% OR:
library(epiDisplay)
logistic.display(m)
Hàm lrm trong rms package 
• rms package – Frank Harrell (regression modeling strategies) 
• lrm – logistic regression model 
lrm(outcome ~ factor1 + factor2 + ...)
Phân tích với R
Tai nạn phi thuyền Challenger
temp = c(66, 70, 69, 80, 68, 67, 72, 73, 70, 57, 63, 70, 78, 
67, 53, 67, 75, 70, 81, 76, 79, 75, 76, 58)
damage = c(0, 1, 0, NA, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 
0, 0, 0, 0, 1, 0, 1) 
dat = data.frame(temp, damage)
# hiển thị mối liên quan 
library(ggplot2)
p = ggplot(data=dat, aes(x = temp, y = damage))
p + geom_point(alpha = 0.15) + geom_smooth(method = "glm", 
method.args = list(family = "binomial"))
0.00
0.25
0.50
0.75
1.00
60 70 80
temp
da
m
ag
e
library(ggplot2)
p = ggplot(data=dat, aes(x = temp, y = damage))
p + geom_point(alpha = 0.15) + geom_smooth(method = "glm", 
method.args = list(family = "binomial"))
Tai nạn phi thuyền Challenger
temp = c(66, 70, 69, 80, 68, 67, 72, 73, 70, 57, 63, 70, 78, 
67, 53, 67, 75, 70, 81, 76, 79, 75, 76, 58)
damage = c(0, 1, 0, NA, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 
0, 0, 0, 0, 1, 0, 1) 
dat = data.frame(temp, damage)
# mô hình hồi qui logistic
m = glm(damage ~ temp, family=binomial, data=dat)
summary(m)
Kết quả 
Coefficients:
Estimate Std. Error z value Pr(>|z|) 
(Intercept) 15.0429 7.3786 2.039 0.0415 *
temp -0.2322 0.1082 -2.145 0.0320 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 28.267 on 22 degrees of freedom
Residual deviance: 20.315 on 21 degrees of freedom
(1 observation deleted due to missingness)
AIC: 24.315
Kết quả: ý nghĩa 
Coefficients:
Estimate Std. Error z value Pr(>|z|) 
(Intercept) 15.0429 7.3786 2.039 0.0415 *
temp -0.2322 0.1082 -2.145 0.0320 *
Mô hình là: 
Nói cách khác: log(odds ratio) = -0.23
Odds ratio = exp(-0.23) = 0.79
log p
1− p
⎛
⎝
⎜
⎞
⎠
⎟=15.04−0.23× temp
Tính toán OR – odds ratio 
> library(epiDisplay)
> logistic.display(m)
Logistic regression predicting damage 
OR(95%CI) P(Wald's test) P(LR-test)
temp (cont. var.) 0.79 (0.64,0.98) 0.032 0.005
Log-likelihood = -10.1576
No. of observations = 23
AIC value = 24.3152
Interpretation: Each degree F increase in temperature was associated with a 21% decrease in the 
odds of being damage, and this association was statistically significant (P = 0.032)
Hiển thị kết quả mô hình hồi qui logistic
Có thể dùng package "GGally" để vẽ biểu đồ odds ratio cho yếu tố nguy cơ
library(GGally); library(ggplot2)
ggcoef(m, exponentiate=T, exclude_intercept=T, vline_color = 
"red", errorbar_color = "blue", errorbar_height = 0.10)
ggcoef(m, exponentiate = TRUE, exclude_intercept=T, mapping = 
aes(x = estimate, y = term, size = p.value)) + 
scale_size_continuous(trans = "reverse")
temp
0.6 0.7 1.0
estimate
te
rm
p.value
0.0319561
ggcoef(m, exponentiate = TRUE, exclude_intercept=T, mapping = 
aes(x = estimate, y = term, size = p.value)) + 
scale_size_continuous(trans = "reverse")
Mô hình hồi qui logistic: tóm tắt
• Tất cả mô hình nghiên cứu (cắt ngang, đoàn hệ, bệnh chứng)
• Mô hình hồi qui logistic rất có ích cho đánh giá mối liên quan
– mô tả mối liên quan giữa biến outcome và yếu tố nguy cơ 
– phát triển mô hình tiên lượng 
• Biến outcome: nhị phân (yes/no)
• Biến tiên lượng (predictor variables): tất cả các dạng liên tục, phân nhóm
• Mô hình logit(p) = α + βx 
Mô hình hồi qui logistic: tóm tắt
• Mô hình logit(p) = α + βx 
• Hàm R
m = glm(y ~ x, family="binomial", data = xxx)
• Hiển thị kết quả (odds ratio)
library(GGally)
ggcoef(m, exponentiate=T, exclude_intercept=T, vline_color = 
"red", errorbar_color = "blue", errorbar_height = 0.10)
Thực hành ...
WHR và tiểu đường 
Câu hỏi: mối liên quan giữa WHR (tỉ số eo-mông) và tiểu đường
> db = read.csv("~/Dropbox/_Conferences and Workshops/Dai hoc Duoc 6-
2019/Datasets/Diabetes data.csv", header=T)
> head(db)
id age gender height weight waist hip sysbp diabp active hypertension
1 1 76 Female 163 53 90 93 160 90 0 1
2 1 40 Female 149 51 74 94 100 60 0 0
3 1 51 Female 151 55 91 100 120 80 0 0
4 1 43 Female 158 62 78 96 120 80 1 0
bmi whr diabetes
1 19.95 0.97 IFG
2 22.97 0.79 Normal
3 24.12 0.91 Normal
4 24.84 0.81 Normal
WHR và tiểu đường 
# Vì biến diabetes có 3 giá trị (Normal, IFG, Yes), chúng ta cần tạo ra biến mới chỉ có 
2 giá trị là 0=không tiểu đường, 1=tiểu đường 
db$diab = ifelse(db$diabetes=="Yes", 1, 0)
# Tìm hiểu xem có bao nhiêu nam và nữ mắc tiểu đường 
library(DescTools)
Desc(db$diab ~ db$gender)
# Vẽ biểu đồ hộp WHR và tiểu đường theo giới tính 
library(ggplot2)
p = ggplot(data=db, aes(x=factor(diab), y=whr, col=factor(diab)))
p + geom_boxplot() + geom_jitter(alpha=0.05) + facet_grid(~gender)
# Vẽ biểu đồ logistic giữa whr và diab 
p = ggplot(data=db, aes(x = whr, y = diab))
p + geom_point(alpha = 0.15) + geom_smooth(method = "glm", method.args = 
list(family = "binomial"))
# Triển khai câu hỏi nghiên cứu, ước tính tham số mô hình hồi qui logistic 
m1 = glm(diab ~ whr, family = binomial, data = db)
summary(m1) 
library(epiDisplay)
logistic.display(m1)
m2 = glm(diab ~ I(whr/0.08), family = binomial, data = db)
summary(m2)
# Khác biệt giữa nam và nữ? 
m2 = glm(diab ~ gender, family = binomial, data = db)
logistic.display(m2)
# Nhưng mối liên quan có thể bị confounded bởi giới tính. Cần hiệu chỉnh 
m3 = glm(diab ~ gender + I(whr/0.08), family = binomial, data = db)
summary(m3)
logistic.display(m3)
# Vẽ biểu đồ odds ratio 
library(GGally)
ggcoef(m3, exponentiate=T, exclude_intercept=T, vline_color = "red", 
errorbar_color = "blue", errorbar_height = 0.10)