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

ob = read.csv ("C:/Users/ADMIN/OneDrive/Statistical courses/Benh vien Cho Ray-Aug2019/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

Chọn các biến liên quan cần phân tích, đặt tên data mới tên là dat

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.2

Khai thác dữ liệu

library(table1)
## Warning: package 'table1' was built under R version 3.6.1
## 
## Attaching package: 'table1'
## The following objects are masked from 'package:base':
## 
##     units, units<-
#Mô tả dữ liệu theo gender
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)
## Warning: package 'GGally' was built under R version 3.6.1
## Loading required package: ggplot2
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
ggpairs(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`.

Mô hình hồi quy tuyến tính đơn giả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

Hoặc một cách khác để có kết quả nhanh hơn.

summary (lm(pcfat ~ gender, data=ob))
## 
## 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)
## Warning: package 'ggfortify' was built under R version 3.6.1
autoplot(m)

#Chú ý biểu đồ khi xuất bản lên RMarkdown không giống như biểu đồ gốc

Phân bố pcfat theo gender

library(visreg)
## Warning: package 'visreg' was built under R version 3.6.1
visreg(m)

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 (liên quan) giữa biến kết cục và các biến tiên lượng

visreg(m0)

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 hai
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.003711
anova(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)

Task 4: Tìm mô hình tối ưu dùng BMA

library(BMA)
## Warning: package 'BMA' was built under R version 3.6.1
## Loading required package: survival
## Loading required package: leaps
## Warning: package 'leaps' was built under R version 3.6.1
## Loading required package: robustbase
## Warning: package 'robustbase' was built under R version 3.6.1
## 
## Attaching package: 'robustbase'
## The following object is masked from 'package:survival':
## 
##     heart
## Loading required package: inline
## Warning: package 'inline' was built under R version 3.6.1
## Loading required package: rrcov
## Warning: package 'rrcov' was built under R version 3.6.1
## Scalable Robust Estimators with High Breakdown Point (version 1.4-7)
# Định nghĩa bế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)

Task 5: Chia dữ liệu thành 2 nhóm nhỏ: development và validation

library(caret)
## Warning: package 'caret' was built under R version 3.6.1
## Loading required package: lattice
## 
## Attaching package: 'caret'
## The following object is masked from 'package:survival':
## 
##     cluster
sample = createDataPartition(ob$pcfat, p=0.6, list=F)
dev = ob[sample, ]
val = ob[-sample, ]

Kiểm tra bộ dữ liệu dev và val

dim (dev) ; dim (val)
## [1] 731  11
## [1] 486  11
#dev có 731 hàng 11 biến số, val có 486 hàng và 11 biến số

Huấ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 
## -17.6820  -2.4532   0.1328   2.6358  15.4526 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -21.078385   5.324512  -3.959 8.28e-05 ***
## genderM     -10.765872   0.335250 -32.113  < 2e-16 ***
## age           0.040269   0.009021   4.464 9.33e-06 ***
## bmi           3.630391   0.458246   7.922 8.76e-15 ***
## `I(bmi^2)`   -0.053656   0.009685  -5.540 4.23e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.015 on 726 degrees of freedom
## Multiple R-squared:  0.6903, Adjusted R-squared:  0.6885 
## F-statistic: 404.5 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.05122595

Kiể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
## 3  34.0 36.72967
## 8  28.0 25.62818
## 9  21.1 21.45873
## 14 24.1 26.41889
## 15 33.0 37.78287
## 17 40.7 35.04582

Tính chỉ số Rsquared

defaultSummary(validation)
##      RMSE  Rsquared       MAE 
## 3.7239987 0.7300235 2.9855669