Ứng dụng 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 ứng dụng hình hồi qui logistic
• Đánh giá mối liên quan, ảnh hưởng (association)
• Hiệu chỉnh (adjustment) cho yếu tố nhiễu
• Tiên lượng (prediction)
Ứng dụng 1: Đánh giá mối liên quan, ảnh hưởng
Cân nặng của trẻ sơ sinh và mẹ hút thuốc lá
• Câu hỏi nghiên cứu: có mối liên quan giữa mẹ hút thuốc lá và trọng lượng
sơ sinh?
• Nghiên cứu cắt ngang trên 189 bà mẹ và trẻ sơ sinh
•
Biến outcome: low, biến tiên lượng: smoke
bw = read.csv("~/Dropbox/_Conferences and Workshops/Dai hoc Duoc 6- 2019/Datasets/birthwt.csv")
head(bw)
id low age lwt race smoke ptl ht ui ftv bwt
1 85 0 19 182 2 0 0 0 1 0 2523
2 86 0 33 155 3 0 0 0 0 3 2551
3 87 0 20 105 1 1 0 0 0 1 2557
Hiển thị mối liên quan bằng biểu đồ
Summary:
bw$smoke = as.factor(bw$smoke)
bs$low = as.factor(bw$low)
library(DescTools)
Desc(bw$smoke ~ bw$low)
estimate lwr.ci upr.ci' odds ratio 2.022 1.081 3.783 rel. risk (col1) 1.258 1.013 1.561 rel. risk (col2) 0.622 0.409 0.945 bw$low 0 1 Sum bw$smoke
0 freq 86 29
115
p.row 74.8% 25.2%
.
p.col 66.2% 49.2%
.
1 freq 44 30
74
p.col 33.8% 50.8%
.
Sum freq 130 59
189
p.row . .
.
p.col . .
.
Mô hình hồi qui logistic
Mô hình liên quan giữa smoke và low như sau: gọi P là xác suất trọng lượng thấp
(cid:230) (cid:246)
(cid:247) = a + b smoke
log
(cid:231)
P 1- P
Ł ł
# Tri n khai b ng R ể ằ
m = glm(low ~ smoke, family=binomial, data=bw) summary(m)
# tính odds ratio library(epiDisplay) logistic.display(m)
Kết quả phân tích
Mô hình liên quan giữa smoke và low như sau: gọi P là xác suất trọng lượng thấp
(cid:230) (cid:246)
(cid:247) = a + b smoke
log
(cid:231)
P 1- P
Ł ł
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.0871 0.2147 -5.062 4.14e-07 ***
smoke1 0.7041 0.3196 2.203 0.0276 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(cid:230) (cid:246)
log
(cid:247) = - 1.087 + 0.704 · smoke
(cid:231)
P 1- P
Ł ł
Kết quả phân tích
> logistic.display(m)
Logistic regression predicting low : 1 vs 0 OR(95%CI) P(Wald's test) P(LR-test) smoke: 1 vs 0 2.02 (1.08,3.78) 0.028 0.027 Log-likelihood = -114.9023 No. of observations = 189 AIC value = 233.8046
Diễn giải: có mối liên quan giữa bà mẹ hút thuốc lá và sinh con thiếu cân. Tính trung bình, những bà mẹ hút thuốc lá có odds sinh con thiếu cân cao 2 lần (KTC 95%: 1.1 đến 3.8) so với những bà mẹ không hút thuốc lá.
Ứng dụng 2: Hiệu chỉnh cho yếu tố nhiễu
Yếu tố nhiễu (confounding factors)
• Câu hỏi nghiên cứu: có mối liên quan giữa bà mẹ hút thuốc lá và
sinh con thiếu cân ?
• Yếu tố nhiễu: có liên quan đến yếu tố nguy cơ (risk factor) và
outcome
• Mối liên quan giữa hút thuốc lá và thiếu cân có thể do yếu tố
chủng tộc và độ tuổi?
• Câu hỏi mới: mối liên quan giữa hút thuốc lá và thiếu cân có độc
lập với yếu tố chủng tộc và độ tuổi?
Yếu tố nhiễu: chủng tộc
> Desc(factor(bw$race) ~ bw$low)
> Desc(factor(bw$race) ~ bw$smoke)
bw$low 0 1 Sum factor(bw$race)
bw$smoke 0 1 Sum factor(bw$race)
1 freq 73 23 96 p.row 76.0% 24.0% . p.col 56.2% 39.0% .
1 freq 44 52 96 p.row 45.8% 54.2% . p.col 38.3% 70.3% .
2 freq 15 11 26 p.row 57.7% 42.3% . p.col 11.5% 18.6% .
2 freq 16 10 26 p.row 61.5% 38.5% . p.col 13.9% 13.5% .
3 freq 42 25 67 p.row 62.7% 37.3% . p.col 32.3% 42.4% .
3 freq 55 12 67 p.row 82.1% 17.9% . p.col 47.8% 16.2% .
Mô hình hiệu chỉnh cho yếu tố nhiễu
• Câu hỏi mới: mối liên quan giữa hút thuốc lá và thiếu cân có độc
lập với yếu tố chủng tộc và độ tuổi?
• Giải pháp: phân tích mô hình hồi qui logistic
• Mô hình mới:
b
(cid:230) (cid:246)
(cid:247) = a
log
1smoke + b
2race + b
3age
· (cid:231)
P 1- P
Ł ł
m = glm(low ~ smoke + factor(race) + age, family = binomial, data=bw)
logistic.display(m)
Kết quả phân tích
> m = glm(low ~ smoke + factor(race) + age, family = binomial, data=bw)
> logistic.display(m)
OR lower95ci upper95ci Pr(>|Z|)
smoke1 3.0058203 1.4499820 6.231081 0.003087348
factor(race)2 2.7494834 1.0453178 7.231924 0.040385185
factor(race)3 2.8769483 1.2983141 6.375061 0.009239817
age 0.9657186 0.9045206 1.031057 0.296337134
Diễn giải: Sau khi hiệu chỉnh cho yếu tố chủng tộc và tuổi, những bà mẹ hút thuốc lá có odds sinh con thiếu cân tăng gấp 3 lần (KTC 95%: 1.45 đến 6.23) so với những bà mẹ không hút thuốc lá. Những bà mẹ da đen (OR 2.7; KTC95: 1.05 – 7.23) và Hispanics (OR 2.88; 95% CI: 1.30 – 6.37) có odds sinh con thiếu cân cao hơn những bà mẹ da trắng. Độ tuổi có mẹ không có liên quan có ý nghĩa thống kê với sinh con thiếu cân.
Yếu tố nào có liên quan đến sinh con thiếu cân?
• Có 8 yếu tố nguy cơ trong nghiên cứu: age, race, lwt,
smoke, ptl, ht, ui, ftv
• Câu hỏi mới: yếu tố nào có liên quan đến sinh con thiếu cân?
• Có ít nhất 2^8 = 256 mô hình khả dĩ !
• Giải pháp: tìm mô hình tối ưu (model selection, feature selection)
• Phương pháp: Bayesian model averaging (BMA)
Tiêu chí cho một mô hình đa biến tối ưu
• Đơn giản (simplicity): ít biến số
• Đầy đủ (adequacy): mô tả dữ liệu một cách thỏa đáng
• Có ý nghĩa thực tế (practicality): phải được yểm trợ bằng lí
thuyết hay có ý nghĩa sinh học
Thách thức: Tìm mô hình với ít biến số nhất nhưng giải
thích nhiều dữ liệu nhất (nguyên tắc parsimony)
Tiêu chuẩn thống kê: AIC
• AIC = Akaike Information Criterion.
AIC = Deviance + 2x(thông số trong mô hình)
• Mô hình tối ưu nhất có AIC thấp nhất. AIC “phạt” những mô hình có
nhiều thông số.
• Nhiều thuật toán thống kê đã được phát triển đi tìm một mô hình với
k thông số sao cho AIC thấp nhất.
• Tìm một mô hình tối ưu là một thách thức rất lớn!
Thuật toán tìm mô hình tối ưu
• Stepwise (forward) algorithm
– Bắt đầu với mô hình 1 thông số
– Từng bước thêm vào 1 thông số cho đến khi deviance không còn giảm
nữa.
• Backward algorithm
– Bắt đầu với mô hình gồm k thông số.
– Loại bỏ từng bước những thông số không có ý nghĩa thống kê
cho đến khi deviance thấp nhất.
• AIC based algorithm: Tìm tất cả các mô hình với AIC thấp nhất
• Bayesian model average (BMA). Tìm tổ các thông số sao cho BIC
thấp nhất
Áp dụng BMA vào nghiên cứu thiếu cân
# Gọi package "BMA". Nếu chưa có, cần phải cài đặt trước library(BMA)
# Chuẩn bị dữ liệu, biến y và biến x, loại bỏ biến id (cột số 1) và y (cột 2) và bwt (cột xvars = bw[, -c(1, 2, 11)]
yvar = bw[, 2]
# Tìm mô hình tối ưu bằng BMA
m = bic.glm(xvars, yvar, strict=F, OR=20, glm.family="binomial")
summary(m)
imageplot.bma(m)
Kết quả phân tích BMA
> summary(m)
84 models were selected
Best 5 models (cumulative posterior probability = 0.2531 ):
p!=0 EV SD model 1 model 2 model 3 model 4 model 5
Intercept 100 -0.390128 1.575728 1.45068 1.09291 -2.32488 -0.35754 1.06795
age 10.4 -0.004815 0.018070 . . . . .
lwt 54.8 -0.008473 0.009253 -0.01865 -0.01707 . -0.01535 -0.01692
race 44.3 0.212462 0.280368 . . 0.55898 0.48955 .
smoke 52.1 0.484523 0.552668 . . 1.11668 1.08002 .
ptl 41.2 0.291512 0.410590 . 0.72560 . . .
ht 59.7 1.011382 0.999519 1.85551 1.85604 . 1.74427 1.96157
ui 30.0 0.263111 0.470489 . . . . 0.93000
ftv 2.0 -0.001015 0.024588 . . . . .
nVar 2 3 2 4 3
BIC -753.82285 -753.75940 -753.62525 -753.44086 -753.11035
post prob 0.058 0.056 0.052 0.048 0.040
BMA phân tích tất cả 84 mô hình, nhưng chỉ báo cáo 5 mô hình tốt nhất
Diễn giải kết quả BMA
p!=0 EV SD
Intercept 100 -0.390128 1.575728
age 10.4 -0.004815 0.018070
lwt 54.8 -0.008473 0.009253
race 44.3 0.212462 0.280368
smoke 52.1 0.484523 0.552668
ptl 41.2 0.291512 0.410590
ht 59.7 1.011382 0.999519
ui 30.0 0.263111 0.470489
ftv 2.0 -0.001015 0.024588
xác suất biến số có liên quan với y
Cột 1 các biến số tiên lượng p!=0 EV expected value: giá trị trung bình của hệ số hồi qui logistic (𝛃) SD standard deviation: độ lệch chuẩn của 𝛃
p!=0 EV SD model 1 model 2 model 3 model 4 model 5 Intercept 100 -0.390128 1.575728 1.45068 1.09291 -2.32488 -0.35754 1.06795 age 10.4 -0.004815 0.018070 . . . . . lwt 54.8 -0.008473 0.009253 -0.01865 -0.01707 . -0.01535 -0.01692 race 44.3 0.212462 0.280368 . . 0.55898 0.48955 . smoke 52.1 0.484523 0.552668 . . 1.11668 1.08002 . ptl 41.2 0.291512 0.410590 . 0.72560 . . . ht 59.7 1.011382 0.999519 1.85551 1.85604 . 1.74427 1.96157 ui 30.0 0.263111 0.470489 . . . . 0.93000 ftv 2.0 -0.001015 0.024588 . . . . . nVar 2 3 2 4 3 BIC -753.82285 -753.75940 -753.62525 -753.44086 -753.11035 post prob 0.058 0.056 0.052 0.048 0.040
số biến tiên lượng của một mô hình
Bayesian information criterion (càng thấp càng tốt)
nVar BIC post prob Model 1 Xác suất hậu định (posterior probability). bao gồm biến lwt (cân nặng của mẹ) và ht (cao huyết áp). Mô hình này có xác
suất hậu định ~6% và được xem là tốt nhất trong các mô hình
Model 2 Model 3 BMA 'đề nghị' biến lwt, ptl, và ht được nhận dạng với các biến race và smoke
Kết quả BMA qua biểu đồ
Models selected by BMA
age
lwt
race
smoke
ptl
ht
ui
ftv
1
2
3
4
5
6
7 8
10
12 14
17
20 23
27 31
36
42 48 56 65 76
Model #
Phương pháp BMA
• Không đưa ra một mô hình duy nhất, mà đề nghị hơn 1 mô hình tốt
nhất
• Mô hình nào tốt nhất còn tuỳ thuộc vào kiến thức chuyên ngành
• BMA được 'chứng minh' là phương pháp tốt nhất khi chọn mô hình
tối ưu
• Dựa vào kết quả, chúng ta chọn mô hình nào để đánh giá mối liên quan giữa các yếu tố nguy cơ và nguy cơ sinh con thiếu cân?
Chọn mô hình" tối ưu"
• Xây dựng mô hình đa biến vừa là một nghệ thuật, vừa là một khoa
học
• Một mô hình đẹp là mô hình mô tả sát với thực tế.
• Một mô hình phản ánh 100% thực tế thì đó không còn là “mô hình”
nữa, hay quá phức tạp không thể ứng dụng được.
• Ngược lại một mô hình chỉ mô tả thực tế khoảng 1% thì cũng không
thể sử dụng được.
•
“Mô hình nào cũng sai so với thực tế, nhưng trong số các mô hình sai đó, có một vài mô hình có ích (George Box)”.