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

Ứ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 đồ 
bw$smoke = as.factor(bw$smoke) 
bs$low = as.factor(bw$low) 
library(DescTools) 
Desc(bw$smoke ~ bw$low) 
Summary: 
 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 
# 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 
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 
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) 
 bw$low 
 0 1 Sum 
factor(bw$race) 
1 freq 73 23 96 
 p.row 76.0% 24.0% . 
 p.col 56.2% 39.0% . 
2 freq 15 11 26 
 p.row 57.7% 42.3 % . 
 p.col 11.5% 18.6% . 
3 freq 42 25 67 
 p.row 62.7% 37.3% . 
 p.col 32.3% 42.4% . 
> Desc(factor(bw$race) ~ bw$smoke) 
 bw$smoke 
 0 1 Sum 
factor(bw$race) 
1 freq 44 52 96 
 p.row 45.8% 54.2% . 
 p.col 38.3% 70.3% . 
2 freq 16 10 26 
 p.row 61.5% 38.5% . 
 p.col 13.9% 13.5% . 
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: 
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! 
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 
Thuật toán tìm mô hình tối ưu 
Á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 
Cột 1 	các biến số tiên lượng 
p!=0 	xác suất biến số có liên quan với y 
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 𝛃 
nVar	 số biến tiên lượng của một mô hình 
BIC	 Bayesian information criterion (càng thấp càng tốt) 
post prob	 Xác suất hậu định (posterior probability). 
Model 1	 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	 BMA 'đề nghị' biến lwt, ptl, và ht 
Model 3 	được nhận dạng với các biến race và smoke 
 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 
Kết quả BMA qua biểu đồ 
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? 
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) ” . 
Chọn mô hình" tối ưu"