Phân tích dữ liệu và ứng dụng - Ứng dụng mô hình hồi qui tuyến tính - Ứng dụng 3: mô hình tiên lượng (prediction)

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
Nội dung
• Mô hình hồi qui tuyến tính đa biến
• Ứng dụng 1: đánh giá mối liên quan (association / assessment)
• Ứng dụng 2: hiệu chỉnh cho yếu tố nhiễu (adjustment)
• Ứng dụng 3: mô hình tiên lượng (prediction)
Xây dựng mô hình tiên lượng tỉ trọng mỡ 
• Để đo tỉ trọng mỡ (pcfat), cần phải có máy DXA (đắt tiền)
• Có thể xây dựng một mô hình tiên lượng pcfat chỉ cần dùng các yếu tố 
'thường qui' (có thể thu thập từ bệnh nhân)
• Các biến thuờng qui: giới tính (gender), chiều cao (height), cân nặng 
(weight), tỉ trọng cơ thể (bmi), tuổi (age)
• Giải pháp: mô hình hồi qui tuyến tính 
Xây dựng mô hình tiên lượng tỉ trọng mỡ 
• Bước 1: Phân tích khai thác (exploratory analysis)
• Bước 2: Tìm các biến liên quan (có giá trị thống kê)
• Bước 3: Chia dữ liệu thành 2 nhóm: development và validation
• Bước 4: Phát triển mô hình dựa vào biến bước 2 trên nhóm development
• Bước 5: Kiểm tra mô hình ở bước 5 trên nhóm validation
Bước 1: Phân tích mô tả / khai thác 
# Các biến có thể liên quan
dat = ob[, c("gender", "weight", "height", "bmi", "age", "pcfat")]
library(GGally)
ggpairs(dat)
Bước 2: Dùng BMA tìm biến liên quan 
> head(ob)
id gender height weight bmi age WBBMC wbbmd 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
library(BMA)
yvar = ob[, ("pcfat")]
xvars = ob[, c("gender", "height", "weight", "bmi", "age")] 
bma = bicreg(xvars, yvar, strict=FALSE, OR=20)
summary(bma)
Dùng BMA tìm mô hình tối ưu 
> summary(bma)
3 models were selected
Best 3 models (cumulative posterior probability = 1 ): 
p!=0 EV SD model 1 model 2 model 3 
Intercept 100.0 5.26146 4.582901 7.958e+00 -7.928e-01 8.137e+00
genderM 100.0 -11.25139 0.429659 -1.144e+01 -1.143e+01 -1.081e+01
height 31.4 0.01759 0.028494 . 5.598e-02 . 
weight 39.2 0.03102 0.042611 7.921e-02 . . 
bmi 100.0 1.01265 0.111625 8.942e-01 1.089e+00 1.089e+00
age 100.0 0.05259 0.008048 5.497e-02 5.473e-02 4.715e-02
nVar 4 4 3 
r2 0.697 0.696 0.695 
BIC -1.423e+03 -1.423e+03 -1.422e+03
post prob 0.392 0.314 0.294 
Diễn giải kết quả BMA: p!=0 
> summary(bma)
3 models were selected
Best 3 models (cumulative posterior probability = 1 ): 
p!=0 EV SD model 1 model 2 model 3 
Intercept 100.0 5.26146 4.582901 7.958e+00 -7.928e-01 8.137e+00
genderM 100.0 -11.25139 0.429659 -1.144e+01 -1.143e+01 -1.081e+01
height 31.4 0.01759 0.028494 . 5.598e-02 . 
weight 39.2 0.03102 0.042611 7.921e-02 . . 
bmi 100.0 1.01265 0.111625 8.942e-01 1.089e+00 1.089e+00
age 100.0 0.05259 0.008048 5.497e-02 5.473e-02 4.715e-02
• Dùng 5 biến cung cấp (gender, height, weight, bmi và age) BMA tìm được 
3 mô hình 'tối ưu'. 
• "p!=0" có nghĩa là xác suất biến số có 'ảnh hưởng'. Vd: xác suất mà height 
có ảnh hưởng đến pcfat là 31.4%, weight là 39.2%, và age là 100% 
Bước 3: Chia dữ liệu thành 2 nhóm 
Chúng ta chia dữ liệu thành 2 nhóm: development (60%) và validation 
(40%) 
Tổng số cỡ mẫu 
(n = 1217)
Development / Training
(n = 1217*0.6)
Validation / Testing 
(n = 1217*0.4)
Xây dựng mô hình tiên lượng Kiểm tra mô hình tiên lượng
Bước 3: Chia dữ liệu thành 2 nhóm 
Tổng số cỡ mẫu 
(n = 1217)
Development / Training
(n = 1217*0.6)
Validation / Testing 
(n = 1217*0.4)
Dùng phương pháp "thủ công"
rows = nrow(ob)
prop = 0.6 
upper = floor(prop*rows)
permutation = ob[sample(rows), ]
dev = permutation[1:upper, ]
val = permutation[(upper+1):rows, ]
Dùng "caret", hàm createDataPartition
library(caret)
sample = createDataPartition(ob$pcfat, 
p=0.6, list=F)
dev = ob[sample, ]
val = ob[-sample, ]
Bước 4: Xây dựng mô hình (training)
rows = nrow(ob)
prop = 0.6 
upper = floor(prop*rows)
permutation = ob[sample(rows), ]
dev = permutation[1:upper, ]
val = permutation[(upper+1):rows, ]
# Xây dựng mô hình dùng dữ liệu của dev 
m = lm(pcfat ~ gender + age + bmi + weight, data=dev)
summary(m)
Bước 4: Xây dựng mô hình (training)
> m = lm(pcfat ~ gender + age + bmi + weight, data=dev) 
> summary(m)
Coefficients:
Estimate Std. Error t value Pr(>|t|) 
(Intercept) 8.11642 1.13548 7.148 2.15e-12 ***
genderM -12.03375 0.44972 -26.758 < 2e-16 ***
age 0.06165 0.00955 6.456 1.97e-10 ***
bmi 0.73709 0.10090 7.305 7.32e-13 ***
weight 0.13856 0.03629 3.818 0.000146 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.997 on 725 degrees of freedom
Multiple R-squared: 0.6885, Adjusted R-squared: 0.6868 
F-statistic: 400.6 on 4 and 725 DF, p-value: < 2.2e-16
Bước 5: Kiểm tra mô hình (validation)
m = lm(pcfat ~ gender + age + bmi + weight, data = dev)
# Kiểm tra mô hình dùng dữ liệu của val 
# Dùng hàm predict 
val$pred = predict(m, newdata = val)
val$resid = val$pred-val$pcfat
# Vẽ residuals vs giá trị tiên lượng 
plot(val$resid ~ val$pred) 
# Tính RMSE – residual mean square error
sum(val$resid^2) / (nrow(val)-5)
# Tính R-square 
cor(val$pred, val$pcfat)^2 
20 30 40 50
-1
5
-1
0
-5
0
5
10
15
val$pred
va
l$
re
si
d
val$pred = predict(m, newdata = val)
val$resid = val$pred-val$pcfat
> sum(val$resid^2) / (nrow(val)-5)
[1] 17.97093
> cor(val$pred, val$pcfat)^2 
[1] 0.6565799
Bước 5: Kiểm tra mô hình (validation)
val$pred = predict(m, newdata = val)
val$resid = val$pred-val$pcfat
> sum(val$resid^2) / (nrow(val)-5)
[1] 17.97093
> cor(val$pred, val$pcfat)^2 
[1] 0.6565799
Training và testing mô hình qua "caret"
library(caret)
# Chia mẫu thành development và validation 
sample = createDataPartition(ob$pcfat, p=0.6, list=F)
dev = ob[sample, ]
val = ob[-sample, ]
# Huấn luyện mô hình: dùng hàm "train"
control = trainControl(method="cv", number=10)
training = train(pcfat ~ gender + age + bmi + weight, data=dev, 
method="lm", trControl=control, metric="Rsquared")
summary(training)
Kết quả training
> summary(training)
Coefficients:
Estimate Std. Error t value Pr(>|t|) 
(Intercept) 7.108911 1.113591 6.384 3.08e-10 ***
genderM -11.987929 0.430910 -27.820 < 2e-16 ***
age 0.054309 0.009447 5.749 1.32e-08 ***
bmi 0.830378 0.101628 8.171 1.36e-15 ***
weight 0.123183 0.036436 3.381 0.000762 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.895 on 726 degrees of freedom
Multiple R-squared: 0.7152, Adjusted R-squared: 0.7136 
F-statistic: 455.7 on 4 and 726 DF, p-value: < 2.2e-16
Kiểm tra mô hình tiên lượng (caret)
# Kiểm tra mô hình (val), tính giá trị tiên lượng 
pred = predict(training, newdata=val)
model.values = data.frame(obs=val$pcfat, pred)
plot(pred ~ ob$pcfat, pch=16)
defaultSummary(model.values)
> defaultSummary(model.values)
RMSE Rsquared 
4.1844599 0.6483567
10 20 30 40
15
20
25
30
35
40
45
val$pcfat
pr
ed
Tóm tắt
Mô hình hồi qui tuyến tính đa biến
• Biến outcome: liên tục; biến tiên lượng: liên tục và phân nhóm, định tính
• Ứng dụng 1: đánh giá các mối liên quan
– đánh giá ý nghĩa thống kê
– đánh giá tầm quan trọng 
• Ứng dụng 2: hiệu chỉnh cho các 'covariates'
– lưu ý đến ảnh hưởng tương tác 
• Ứng dụng 3: tiên lượng 
– training và testing 
Các hàm R đã học 
• Hàm căn bản
m = lm(y ~ x1 + x2 + x3, data = xx)
• Hiển thị mô hình
library(visreg)
visreg(m, xvar="age", gg=T) 
• Đánh giá tầm quan trọng
library(relaimpo)
calc.relimp(m, type="lmg", rela=T, rank=T)
• Tìm biến liên quan 
library(BMA)
bma = bicreg(xvars, yvar, strict=FALSE, OR=20)
Các hàm R đã học 
• Chia data thành 2 nhóm 
library(caret)
sample = createDataPartition(ob$pcfat, p=0.6, list=F)
dev = ob[sample, ]
val = ob[-sample, ]
• Training mô hình 
control = trainControl(method="cv", number=10)
training = train(y ~ ., data=dev, method="lm", 
trControl=control, metric="Rsquared")
• Kiểm tra mô hình 
pred = predict(training, newdata=val)
model.values = data.frame(obs=val$pcfat, pred)
defaultSummary(model.values)