Bài giảng Phân tích dữ liệu và ứng dụng - Bài 8a: Mô hình hồi quy 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 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)