Bài giảng Phân tích dữ liệu và ứng dụng - Bài 7a: Mô hình hồi quy tuyến tính đơn giản

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 tuyến tính • Phân tích tương quan • Mô hình hồi qui tuyến tính đơn giản • Kiểm tra giả định mô hình hồi qui tuyến tính Khái niệm tương quan (correlation) • Khi hai biến số (x và y) có liên quan với nhau • Mối liên quan có thể cùng chiều hay nghịch đảo • Ví dụ: mối liên quan giữa tiêu thụ chocolate và giải Nobel (?) Tương quan giữa 2 biến liên tục Làm sao định lượng mối liên quan? Sir Francis Galton (16/2/1822 – 17/1/1911) Galton’s conclusions: • Nature dominates: “families of reputation were much more likely than ordinary families to produce offspring of ability” • Recommended “judicious marriages during several generations” to “produce a highly gifted race of men” • His “genetic utopia”: “Bright, healthy individuals were treated and paid well, and encouraged to have plenty of children. Social undesirables were treated with reasonable kindness so long as they worked hard and stayed celibate.” Didn’t have data on “intelligence” so instead studied HEIGHT • Although a self-proclaimed genius, who wrote that he could read @2½, write/do arithmetic @4, and was comfortable with Latin texts @8, he couldn’t figure out how to model these data(!) • He went to JD Dickson, a mathematician at Cambridge, who formalized the relationship by developing what we now know as linear regression Research interest: “Those qualifications of intellect and disposition which lead to reputation” J Singer's slide Làm thế nào để mô tả mối tương quan tuyến tính? • Gọi X và Y là 2 biến ngẫu nhiên từ n quan sát • Đo lường độ biến thiên: phương sai (variance) ( ) ( )å = - - = n i i n xxx 1 2 1 var ( ) ( )å = - - = n i i n yyy 1 2 1 var • Chúng ta cần một thước đo độ "hiệp biến" giữa X và Y • Covariance là trung bình của tích số X và Y ( ) ( )( )å = -- - = n i ii yyxx n yx 11 1,cov Ước tính hệ số tương quan • Covariance có đơn vị đo lường (X * Y). • Coefficient of correlation (r) giữa X và Y là một standardized covariance – không có đơn vị đo lường • r định nghĩa như sau: ( ) ( ) ( ) ( ) yx SDSD yx yx yxr ´ = ´ = ,cov varvar ,cov Obesity data (Vietnam) • Nghiên cứu cắt ngang >1100 nam và nữ (Việt Nam) • Mục tiêu: ước tính hệ số tương quan giữa tỉ trọng cơ thể (bmi) và tỉ trọng mỡ (pcfat) ob = read.csv("~/Dropbox/_Conferences and Workshops/Dai hoc Duoc 6- 2019/Datasets/obesity data.csv") dim(ob) [1] 1217 11 > head(ob) id gender height weight bmi age bmc bmd fat lean pcfat 1 1 F 150 49 21.8 53 1312 0.88 17802 28600 37.3 2 2 M 165 52 19.1 65 1309 0.84 8381 40229 16.8 3 3 F 157 57 23.1 64 1230 0.84 19221 36057 34.0 4 4 F 156 53 21.8 56 1171 0.80 17472 33094 33.8 5 5 M 160 51 19.9 54 1681 0.98 7336 40621 14.8 6 6 F 153 47 20.1 52 1358 0.91 14904 30068 32.2 Hàm cor.test > cor.test(ob$bmi, ob$pcfat) data: ob$bmi and ob$pcfat t = 17.123, df = 1215, p-value < 2.2e-16 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.3944998 0.4850923 sample estimates: cor 0.4409183 Phân tích tương quan đa biến (GGally) • Hàm cor.test chỉ tính toán tương quan giữa 2 biến • Với k biến, chúng ta có k(k – 1)/2 hệ số tương quan • Mục tiêu: tính toán tất cả mối tương quan • Package "GGally" library("GGally") ggcorr(dat, label=T) Phân tích tương quan đa biến Dữ liệu: tất cả biến số phải là biến liên tục # Mục tiêu: tính hệ số tương quan cho nhiều biến trong ob # loại trừ 2 cột đầu tiên của ob, chỉ quan tâm đến các biến liên tục dat = ob[, -c(1,2)] # loại bỏ 2 biến số 1 và 2 (id, gender) library("GGally") ggcorr(dat, label=T) dat = ob[, -c(1,2)] # loại bỏ 2 biến số 1 và 2 (id, gender) library("GGally") ggcorr(dat, label=T) 0.6 0 -0.4 0.7 0.4 -0.1 0.8 -0.5 0.8 0 0.6 0.3 0.6 0.8 0.1 0.2 0.2 0.1 0.8 0.5 0.4 -0.5 -0.4 0.2 -0.2 0.3 0.9 0 0.7 -0.4 -0.1 0.4 -0.3 0.1 0.8 -0.4 height weight bmi age bmc bmd fat lean pcfat -1.0 -0.5 0.0 0.5 1.0 ggpairs(dat) Hàm pairs.panels trong psych dat = ob[, -c(1,2)] # Gọi package psych và dùng hàm pairs.panels library(psych) pairs.panels(dat) r và R2 • r là hệ số tương quan • R2 là hệ số xác định (coefficient of determination) phản ảnh phần trăm phương sai của y có thể giải thích bởi biến x • r(weight, BMD) = 0.33 có nghĩa là R2 = (0.33)2 = 0.11. 11% độ khác biệt về BMD có thể giải thích bằng những khác biệt về cân nặng Ý tưởng và nhu cầu mô hình hồi qui tuyến tính 2 4 6 8 10 12 14 16 20 30 40 50 60 X Y 2 4 6 8 10 12 14 16 -4 -2 0 2 X Y Làm cách nào để định lượng hóa mối liên quan giữa Y và X ? Ý tưởng mô hình hồi qui tuyến tính • Biến phụ thuộc (Y) phải là biến liên tục (vd: pcfat) • Biến tiên lượng (X) hay predictor variables: không giới hạn (vd: giới tính, tuổi) • Hồi qui tuyến tính đơn giản (simple linear regression model) – có một biến tiên lượng Mô hình hồi qui tuyến tính Mô hình: Y = a + bX + e a : intercept b : slope / gradient ε : sai số ngẫu nhiên (random error – những dao động về Y trong mỗi giá trị X) • Mối liên quan giữa X và Y là tuyến tính (linear) về tham số • X không có sai số ngẫu nhiên • Giá trị của Y độc lập với nhau (vd, Y1 không liên quan với Y2) ; • Sai số ngẫu nhiên (e): phân bố chuẩn, trung bình 0, phương sai bất biến ε ~ N(0, s2) Giả định Tham số của mô hình HQTT • Mô hình Y = a + bX + e • Chúng ta không biết α và β • Nhưng có thể dùng dữ liệu thí nghiệm / thực tế để ước tính 2 tham số đó • Ước số (estimate) của α và β là a và b Phương pháp "Least square" Y X ii bxay +=ˆ iii yyd ˆ-= yi Tìm công thức (estimator) để tính a và b sao cho tổng d2 là nhỏ nhất à Least square method = Bình phương nhỏ nhất Ước tính tham số mô hình HQTT bằng R • Mô hình hồi qui tuyến tính: y = α + β*X + ε • Triển khai bằng R: hàm "lm" lm(y~ x) Ví dụ: nghiên cứu về tỉ trọng mỡ (pcfat) • Nghiên cứu cắt ngang, 1217 nam và nữ (Việt Nam) • Câu hỏi 1: Ảnh hưởng của tuổi đến tỉ trọng mỡ ? • Câu hỏi 2: Tỉ trọng mỡ khác nhau giữa nam và nữ ob = read.csv("~/Dropbox/_Conferences and Workshops/Dai hoc Duoc 6- 2019/Datasets/obesity data.csv") > head(ob) id gender height weight bmi age bmc bmd fat lean pcfat 1 1 F 150 49 21.8 53 1312 0.88 17802 28600 37.3 2 2 M 165 52 19.1 65 1309 0.84 8381 40229 16.8 3 3 F 157 57 23.1 64 1230 0.84 19221 36057 34.0 4 4 F 156 53 21.8 56 1171 0.80 17472 33094 33.8 5 5 M 160 51 19.9 54 1681 0.98 7336 40621 14.8 6 6 F 153 47 20.1 52 1358 0.91 14904 30068 32.2 Biểu đồ 20 40 60 80 10 20 30 40 ob$age ob $p cf at F M 10 20 30 40 plot(ob$pcfat ~ ob$age, pch=16, col="blue") boxplot(ob$pcfat ~ ob$gender, col="blue", border="green") Triển khai câu hỏi nghiên cứu bằng R • Ảnh hưởng của tuổi đến tỉ trọng mỡ ? m1 = lm(pcfat ~ age, data=ob) summary(m1) • Khác biệt về tỉ trọng mỡ giữa nam và nữ ? m2 = lm(pcfat ~ gender, data=ob) summary(m2) Câu hỏi 1: tuổi và tỉ trọng mỡ > m1 = lm(pcfat ~ age, data=ob) > summary(m1) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 25.58408 0.57003 44.88 <2e-16 *** age 0.12769 0.01135 11.25 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 6.839 on 1215 degrees of freedom Multiple R-squared: 0.09431, Adjusted R-squared: 0.09357 F-statistic: 126.5 on 1 and 1215 DF, p-value: < 2.2e-16 Tuổi và tỉ trọng mỡ: diễn giải Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 25.58408 0.57003 44.88 <2e-16 *** age 0.12769 0.01135 11.25 <2e-16 *** Mô hình (giả thuyết) pcfat = α + β(age) + ε Phương trình pcfat = 25.6 + 0.13*age Diễn giải: Mỗi năm tăng độ tuổi, tỉ trọng mỡ tăng 0.13% (SE 0.011), và mối liên quan này có ý nghĩa thống kê (P < 0.0001) a = 25.6 (ước số của α) b = 0.13 (ước số của β) m1 = lm(pcfat ~ age, data=ob) plot(ob$pcfat ~ ob$age, pch=16, col="blue") abline(m1, col="red") 20 40 60 80 10 20 30 40 ob$age ob $p cf at pcfat = 25.6 + 0.13*age Khác biệt về pcfat giữa nam và nữ: diễn giải Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 34.6724 0.1826 189.9 <2e-16 *** genderM -10.5163 0.3381 -31.1 <2e-16 *** Mô hình (giả thuyết) pcfat = α + β(gender) + ε Phương trình pcfat = 34.7 – 10.5*gender(M) Diễn giải: Nam có tỉ trọng mỡ thấp hơn nữ 10.5% (SE 0.34%), và sự khác biệt này có ý nghĩa thống kê (P < 0.0001) Khác biệt gvề pcfat giữa nam và nữ: diễn giải Residual standard error: 5.362 on 1215 degrees of freedom Multiple R-squared: 0.4432, Adjusted R-squared: 0.4428 F-statistic: 967.3 on 1 and 1215 DF, p-value: < 2.2e-16 Diễn giải: Khác biệt giữa nam và nữ giải thích 44% những khác biệt về phương sai của tỉ trọng mỡ. Báo cáo qua package "stargazer" library(stargazer) stargazer(m1, type="text") # text, htmll latex > stargazer(m1, type="text") =============================================== Dependent variable: --------------------------- pcfat ----------------------------------------------- age 0.128*** (0.011) Constant 25.584*** (0.570) ----------------------------------------------- Observations 1,217 R2 0.094 Adjusted R2 0.094 Residual Std. Error 6.839 (df = 1215) F Statistic 126.519*** (df = 1; 1215) =============================================== Note: *p<0.1; **p<0.05; ***p<0.01 • Mối liên quan giữa X và Y là tuyến tính (linear) về tham số • X không có sai số ngẫu nhiên • Giá trị của Y độc lập với nhau (vd, Y1 không liên quan với Y2) ; • Sai số ngẫu nhiên (e): phân bố chuẩn, trung bình 0, phương sai bất biến (homogeneity) ε ~ N(0, s2) Giả định mô hình hồi qui tuyến tính Kiểm tra qua "plot" m1 = lm(pcfat ~ age, data=ob) par(mfrow=c(2,2)) plot(m1) 28 30 32 34 36 -2 0 0 10 20 Fitted values R es id ua ls Residuals vs Fitted 972210 858 -3 -2 -1 0 1 2 3 -3 -1 1 2 3 Theoretical Quantiles S ta nd ar di ze d re si du al s Normal Q-Q 972210858 28 30 32 34 36 0. 0 0. 5 1. 0 1. 5 Fitted values S ta nd ar di ze d re si du al s Scale-Location 972210 858 0.000 0.001 0.002 0.003 0.004 0.005 -4 -2 0 1 2 3 Leverage S ta nd ar di ze d re si du al s Cook's distance Residuals vs Leverage 858 245 972 Tóm tắt: mô hình hồi qui tuyến tính • Một phương pháp định lượng hóa một mối tương quan • Mô hình cơ bản y = α + β*X + ε • Triển khai bằng R lm(y ~ x, data=xxx)