Thực tập buổi 3
Trần Ngọc Nguyên
8/6/2019
Mô hình hồi quy tuyến tính
Package cần làm quen/sử dụng: BMA, visreg, GGally, table1, caret, relaimpo Dữ liệu thực hành: “obesity data.csv” Mục tiêu: đánh giá mối liên quan giữa pcfat với các biến gender, age, height, weight, bmi
Task 1
ob = read.csv("C:\\Users\\Nguyen\\Desktop\\Workshop NCKH 2019\\Dataset\\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# Chọn biến liên quan
dat = ob[, c("gender", "age", "height", "weight" , "bmi", "pcfat")]
head(dat)## gender age height weight bmi pcfat
## 1 F 53 150 49 21.8 37.3
## 2 M 65 165 52 19.1 16.8
## 3 F 64 157 57 23.1 34.0
## 4 F 56 156 53 21.8 33.8
## 5 M 54 160 51 19.9 14.8
## 6 F 52 153 47 20.1 32.2Task 2: Khai thác dữ liệu (mô tả)
# Mô tả dữ liệu theo gender
library(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-table1(~age + height + weight + bmi + pcfat | gender, data=ob)| F (n=862) |
M (n=355) |
Overall (n=1217) |
|
|---|---|---|---|
| age | |||
| Mean (SD) | 48.6 (16.4) | 43.7 (18.8) | 47.2 (17.3) |
| Median [Min, Max] | 49.0 [14.0, 85.0] | 44.0 [13.0, 88.0] | 48.0 [13.0, 88.0] |
| height | |||
| Mean (SD) | 153 (5.55) | 165 (6.73) | 157 (7.98) |
| Median [Min, Max] | 153 [136, 170] | 165 [146, 185] | 155 [136, 185] |
| weight | |||
| Mean (SD) | 52.3 (7.72) | 62.0 (9.59) | 55.1 (9.40) |
| Median [Min, Max] | 51.0 [34.0, 95.0] | 62.0 [38.0, 95.0] | 54.0 [34.0, 95.0] |
| bmi | |||
| Mean (SD) | 22.3 (3.05) | 22.7 (3.04) | 22.4 (3.06) |
| Median [Min, Max] | 22.1 [15.2, 37.1] | 22.5 [14.5, 34.7] | 22.2 [14.5, 37.1] |
| pcfat | |||
| Mean (SD) | 34.7 (5.19) | 24.2 (5.76) | 31.6 (7.18) |
| Median [Min, Max] | 34.7 [14.6, 48.4] | 24.6 [9.20, 39.0] | 32.4 [9.20, 48.4] |
library(GGally)## Loading required package: ggplot2## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2ggpairs(dat, mapping = aes(color = gender))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Phân tích hồi quy tuyến tính cho thấy nam có ít tỷ trọng mỡ thấp hơn nữ 10.5% (sai số chuẩn 0.34%) và sự khác biệt này có ý nghĩa thống kê (p<0.0001). Yếu tố giới tính giải thích khoảng 44% những khác biệt về tỷ trọng mỡ giữa các cá nhân.
Task 3: Mô hình hồi quy tuyến tính đơn giản
Khác biệt gữa nam và nữ
m = lm(pcfat ~ gender, data=ob)
summary(m)##
## Call:
## lm(formula = pcfat ~ gender, data = ob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.0724 -3.2724 0.1484 3.6276 14.8439
##
## 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 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 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# Kiểm tra giả định của mô hình
library(ggfortify)
library(ggplot2)
autoplot(m)
# Phân bố pcfat theo gender
library(visreg)
visreg(m)
Phân tích hồi quy tuyến tính cho thấy nam có ít tỷ trọng mỡ thấp hơn nữ 10.5% (sai số chuẩn 0.34%) và sự khác biệt này có ý nghĩa thống kê (p<0.0001). Yếu tố giới tính giải thích khoảng 44% những khác biệt về tỷ trọng mỡ giữa các cá nhân.
Khác biệt gữa nam và nữ, hiệu chỉnh cho BMI
m0 = lm(pcfat ~ gender + bmi, data=ob)
summary(m0)##
## Call:
## lm(formula = pcfat ~ gender + bmi, data = ob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.4709 -2.4780 0.1773 2.6903 15.1761
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.01035 0.85880 10.49 <2e-16 ***
## genderM -11.06631 0.25599 -43.23 <2e-16 ***
## bmi 1.15303 0.03809 30.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.049 on 1214 degrees of freedom
## Multiple R-squared: 0.6828, Adjusted R-squared: 0.6822
## F-statistic: 1306 on 2 and 1214 DF, p-value: < 2.2e-16# Biểu đồ tương quan giữa biến kết cục và các biến tiên lượng
visreg(m0)
Phân tích hồi quy tuyến tính cho thấy ở cùng BMI, nam có ít tỷ trọng mỡ thấp hơn nữ 11% (sai số chuẩn 0.26%) và sự khác biệt này có ý nghĩa thống kê (p<0.0001). Phân tích hồi quy tuyến tính cũng cho thấy cùng giới tính, khi BMI tăng 1kg/m2 thì tỷ trọng mỡ tăng 1.15% (sai số chuẩn 0.04%) và sự khác biệt này có ý nghĩa thống kê (p<0.0001). Yếu tố giới tính và BMI giải thích khoảng 68% những khác biệt về tỷ trọng mỡ giữa các cá nhân.
Tương quan với bmi
# Hồi quy bậc nhất
m1 = lm(pcfat ~ bmi, data=ob)
m1##
## Call:
## lm(formula = pcfat ~ bmi, data = ob)
##
## Coefficients:
## (Intercept) bmi
## 8.399 1.036# Hồi quy bậc 2
m2 = lm(pcfat ~ bmi + I(bmi^2), data=ob)
m2##
## Call:
## lm(formula = pcfat ~ bmi + I(bmi^2), data = ob)
##
## Coefficients:
## (Intercept) bmi I(bmi^2)
## -16.31919 3.20632 -0.04675# Hồi quy bậc 3
m3 = lm(pcfat ~ bmi + I(bmi^2) + I(bmi^3), data=ob)
m3##
## Call:
## lm(formula = pcfat ~ bmi + I(bmi^2) + I(bmi^3), data = ob)
##
## Coefficients:
## (Intercept) bmi I(bmi^2) I(bmi^3)
## -68.046558 9.799399 -0.320765 0.003711anova(m1, m2, m3)## Analysis of Variance Table
##
## Model 1: pcfat ~ bmi
## Model 2: pcfat ~ bmi + I(bmi^2)
## Model 3: pcfat ~ bmi + I(bmi^2) + I(bmi^3)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1215 50541
## 2 1214 49921 1 620.14 15.1175 0.0001065 ***
## 3 1213 49758 1 162.30 3.9565 0.0469163 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Chia thành 3 cửa số và hiển thị mô hình m1, m2, m3
par(mfrow=c(1,3))
visreg(m1)
visreg(m2)
visreg(m3)
Viết vài dòng nhận xét các kết quả phân tích. Bạn chọn mô hình nào?
Tương quan với tuổi, hiệu chỉnh theo BMI - biểu đồ nhiệt
m4 = lm(pcfat ~ age + bmi + I(bmi^2), data = ob)
visreg2d(m4, "age", "bmi")
Tiên lượng pcfat
m5 = lm(pcfat ~ gender + age + bmi + I(bmi^2), data = ob)
ob$predicted = predict(m5, newdata=ob)
head(ob)## id gender height weight bmi age bmc bmd fat lean pcfat predicted
## 1 1 F 150 49 21.8 53 1312 0.88 17802 28600 37.3 34.79014
## 2 2 M 165 52 19.1 65 1309 0.84 8381 40229 16.8 20.73820
## 3 3 F 157 57 23.1 64 1230 0.84 19221 36057 34.0 36.79828
## 4 4 F 156 53 21.8 56 1171 0.80 17472 33094 33.8 34.92231
## 5 5 M 160 51 19.9 54 1681 0.98 7336 40621 14.8 21.43289
## 6 6 F 153 47 20.1 52 1358 0.91 14904 30068 32.2 32.49256Task 4: Tìm mô hình tối ưu dùng BMA
library(BMA)## Loading required package: survival## Loading required package: leaps## Loading required package: robustbase##
## Attaching package: 'robustbase'## The following object is masked from 'package:survival':
##
## heart## Loading required package: inline## Loading required package: rrcov## Scalable Robust Estimators with High Breakdown Point (version 1.4-7)# Định nghĩa biến y và x.
y = ob$pcfat
x = ob[, c("gender", "age", "height", "weight", "bmi")]
# Tìm biến liên quan
bma = bicreg(x, y, strict=FALSE, OR=20)
summary(bma)##
## Call:
## bicreg(x = x, y = y, strict = FALSE, OR = 20)
##
##
## 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
## age 100.0 0.05259 0.008048 5.497e-02 5.473e-02 4.715e-02
## 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
##
## 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# Hình ảnh các mô hình có thể có khi tiên lượng pcfat
imageplot.bma(bma)
Bạn chọn mô hình nào?
Task 5:
Chia dữ liệu thành 2 nhóm nhỏ: development và validation
library(caret)## Loading required package: lattice##
## Attaching package: 'caret'## The following object is masked from 'package:survival':
##
## clustersample = createDataPartition(ob$pcfat, p=0.6, list=F)
dev = ob[sample, ]
val = ob[-sample, ]
dim(dev); dim(val)## [1] 731 12## [1] 486 12Huấn luyện (xây dựng) mô hình
control = trainControl(method="cv", number=10)
training = train(pcfat ~ gender + age + bmi + I(bmi^2), data=dev, method="lm", trControl=control, metric="Rsquared")
summary(training)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.1788 -2.3742 0.1027 2.4671 14.1412
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -22.967251 5.353875 -4.290 2.03e-05 ***
## genderM -11.074593 0.314226 -35.244 < 2e-16 ***
## age 0.035987 0.008386 4.291 2.02e-05 ***
## bmi 3.893305 0.467374 8.330 4.02e-16 ***
## `I(bmi^2)` -0.060775 0.010083 -6.028 2.65e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.793 on 726 degrees of freedom
## Multiple R-squared: 0.7149, Adjusted R-squared: 0.7133
## F-statistic: 455 on 4 and 726 DF, p-value: < 2.2e-16Đánh giá tầm quan trọng của từng biến
library(relaimpo)## Loading required package: MASS## Loading required package: boot##
## Attaching package: 'boot'## The following object is masked from 'package:lattice':
##
## melanoma## The following object is masked from 'package:robustbase':
##
## salinity## The following object is masked from 'package:survival':
##
## aml## Loading required package: survey## Loading required package: grid## Loading required package: Matrix##
## Attaching package: 'survey'## The following object is masked from 'package:graphics':
##
## dotchart## Loading required package: mitools## This is the global version of package relaimpo.## If you are a non-US user, a version with the interesting additional metric pmvd is available## from Ulrike Groempings web site at prof.beuth-hochschule.de/groemping.m = lm(pcfat ~ gender + age + bmi + I(bmi^2), data=ob)
calc.relimp(m, type="lmg", rela=T, rank=T)## Response variable: pcfat
## Total response variance: 51.5935
## Analysis based on 1217 observations
##
## 4 Regressors:
## gender age bmi I(bmi^2)
## Proportion of variance explained by model: 70.65%
## Metrics are normalized to sum to 100% (rela=TRUE).
##
## Relative importance metrics:
##
## lmg
## gender 0.64634359
## age 0.06186434
## bmi 0.15442688
## I(bmi^2) 0.13736518
##
## Average coefficients for different model sizes:
##
## 1X 2Xs 3Xs 4Xs
## gender -10.51634414 -10.717838213 -10.88911021 -10.86319562
## age 0.12768705 0.092389959 0.06189251 0.04405818
## bmi 1.03619023 1.759526891 2.52740016 3.47197075
## I(bmi^2) 0.02152954 -0.001260437 -0.02421181 -0.05122595Kiểm định (testing) mô hình
# Tính giá trị tiên lượng
pred = predict(training, newdata=val)
plot(pred ~ val$pcfat, pch=16)
# Data mô tả giá trị kết cục quan sát trên thực tế và kết cục có được từ mô hình tiên lượng
validation = data.frame(obs=val$pcfat, pred)
head(validation)## obs pred
## 7 35.3 37.93766
## 9 21.1 21.41446
## 19 35.5 29.38222
## 21 42.6 39.72023
## 31 25.2 22.73233
## 32 42.0 37.61378# Tính chỉ số Rsquared
defaultSummary(validation)## RMSE Rsquared MAE
## 4.0664825 0.6923407 3.1888860