Nhập dữ liệu và hiển thị, chọn biến liên quan

ob = read.csv("C:/Users/Dell/Desktop/Dataset/Obesity data.csv")
dat = ob[, c('gender', 'age', "height", 'weight', 'bmi', 'pcfat')]
head(ob,6)
##   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
library(GGally) #GGally là package chuyên dụng cho phân tích tương quan đa biến
## 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`.

#Nếu không customize để có màu, syntax mặc định ggpairs(dat)

Task 2: Thử tạo một hàm tuyến tính đơn giản

m = lm(pcfat~gender, data = ob) #Hồi quy tuyến tính pcfat theo gender
summary(m) #Hàm summary tóm tắt các thông số cho hồi quy tuyết tính, ở đây là 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
#Phân tích hồi qui tuyến tính cho thấy nam có 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

Mapping linear model using visreg package

library(visreg)
visreg(m)

Difference between male and female after ajusting by BMI

n = lm(pcfat ~ gender + bmi + bmi:gender, data = ob)
visreg(n, xvar = "bmi", by ="gender", overlay = T)

#visreg làm hàm trong visreg library
#xvar: chọn biến biểu diễn trục hoành ( do trong mô hình tuyến tính có 2 biến tiên lượng -> hàm yêu cầu cần chị rõ trục hoành là biến nào)
#by ="x", x là biến còn lại trong mô hình, ta có thể chọn để phân thành các mô hình với các giá trị khác nhau của biến -> biến này thường hữu hạn về giá trị và thường sẽ là biến phân loại ( ví dụ gender chỉ 2 giá trị là 0 - 1 hay là nam - nữ)
#overlay = T/F - True/False : các hình vễ chồng lên nhau hay tách thành các hình vẽ khác nhau

#Sau hiệu chỉnh bmi theo giới, ta thấy không có hiện tượng tương tác - interaction

m1 = lm(pcfat ~ bmi, data = ob)
m2 = lm(pcfat ~ bmi + I(bmi^2), data = ob)
m3 = lm(pcfat ~ bmi + I(bmi^2) + I(bmi^3), data = ob)
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ổ hiển thị

par(mfrow = c(1,3)) 
#Hàm par(ở đây có thể chọn (mfcol, mfrow) = c(a,b) với a,b là số
#Patrition the screen into 3 windows, 
visreg(m1)
visreg(m2)
visreg(m3)

library(ggplot2)
p = ggplot(data = ob, aes(x = bmi, y =pcfat, fill = gender,col = gender))+geom_point()+ geom_smooth()
#geom_point(): chọn đối tượng hình học biểu diễn là điểm
#geom_smooth(): vẽ đường thẳng biểu hiện giá trị trung vị và vùng cùng màu 2 bên là khoảng dao động 95%
p
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

#Như hình trên, visualizatiion cho ta thất có vẻ dữ liệu là một hàm bậc 2, từ đó ta chuyển tuyến tính theo kiểu bậc 2

Thử dùng visref2d

x1 = lm(pcfat ~ gender + age + bmi + (bmi^2),data = ob)
visreg2d(x1,'gender', 'bmi')

visreg2d(x1,'age', 'bmi')

n1 = lm(pcfat ~ bmi, data = ob)
summary(n1)
## 
## Call:
## lm(formula = pcfat ~ bmi, data = ob)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -19.612  -4.181   1.392   4.690  18.241 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.39889    1.36777   6.141 1.11e-09 ***
## bmi          1.03619    0.06051  17.123  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.45 on 1215 degrees of freedom
## Multiple R-squared:  0.1944, Adjusted R-squared:  0.1937 
## F-statistic: 293.2 on 1 and 1215 DF,  p-value: < 2.2e-16

#R-squared : 0.19

n2 = lm(pcfat ~ bmi + I(bmi^2), data = ob)
summary(n2)
## 
## Call:
## lm(formula = pcfat ~ bmi + I(bmi^2), data = ob)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -19.503  -4.290   1.512   4.633  18.002 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -16.31919    6.50871  -2.507 0.012296 *  
## bmi           3.20632    0.56205   5.705 1.46e-08 ***
## I(bmi^2)     -0.04675    0.01204  -3.883 0.000109 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.413 on 1214 degrees of freedom
## Multiple R-squared:  0.2043, Adjusted R-squared:  0.203 
## F-statistic: 155.8 on 2 and 1214 DF,  p-value: < 2.2e-16

#R-Squared: 0.2 => tốt hơn lần đầu

n3 = lm(pcfat ~ bmi + gender + I(bmi^2), data = ob)
summary(n3)
## 
## Call:
## lm(formula = pcfat ~ bmi + gender + I(bmi^2), data = ob)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.4652  -2.5166   0.0536   2.6760  14.8865 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -19.739284   4.024626  -4.905 1.06e-06 ***
## bmi           3.677757   0.347640  10.579  < 2e-16 ***
## genderM     -11.108667   0.250710 -44.309  < 2e-16 ***
## I(bmi^2)     -0.054377   0.007444  -7.305 5.02e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.964 on 1213 degrees of freedom
## Multiple R-squared:  0.6961, Adjusted R-squared:  0.6954 
## F-statistic: 926.3 on 3 and 1213 DF,  p-value: < 2.2e-16
R-Squared: 0.7 = > rất tốt, R-squared càng cao thì khả năng các số đo được dùng có thể giải thích cho khác biệt phương sai của biến càng cao
Ta có thể việt lại thành phương trình
pcfat = -19.7 + 3.68*bmi - 11.1(Nam) -0.054*bmi*bmi
#Nam, bmi = 20
-19.7 + 3.68*20 -11.1 - 0.054*20*20
## [1] 21.2
final = lm(pcfat ~ gender + age + bmi + I(bmi^2),data = ob)
summary(final) 
## 
## Call:
## lm(formula = pcfat ~ gender + age + bmi + I(bmi^2), data = ob)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.8829  -2.5326   0.0884   2.5596  15.3824 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -18.889287   3.959182  -4.771 2.06e-06 ***
## genderM     -10.863196   0.249341 -43.568  < 2e-16 ***
## age           0.044058   0.006736   6.541 8.99e-11 ***
## bmi           3.471971   0.343248  10.115  < 2e-16 ***
## I(bmi^2)     -0.051226   0.007335  -6.984 4.72e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.898 on 1212 degrees of freedom
## Multiple R-squared:  0.7065, Adjusted R-squared:  0.7055 
## F-statistic: 729.3 on 4 and 1212 DF,  p-value: < 2.2e-16
#Nam, age = 60, bmi = 23
-18.89 - 10.87 + 0.05*60 + 3.47*23 -0.051*23*23
## [1] 26.071

LÀM LẠI BIỂU ĐỒ TRONG SÁCH

Phân tích tương quan đa biến

dat1 = ob[,-c(1,2)] #loại bỏ 2 biến số 1 và 2 (id,gender) trong object ob
library(GGally) #Gọi thư viện/package 'GGally'
ggcorr(dat1, label = TRUE) 

#ggcorr() là hàm trong thư viện GGally-(tương quan đa biến) nhằm vẽ biểu đồ thể hiện mức độ tương quan, mặc định là phần biệt qua màu
#'label = T/F(True/False)'là hàm để thêm correlation coefficient dạng number vào biểu đồ
ggpairs(dat1)

Dùng thư viện psych thay cho GGally

library(psych)
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
pairs.panels(dat1) 

#ưu thế của thư viện là còn có vẽ luôn epsilon(số dư-sai số ngẫu) của hồi quy tuyến tính thực tế(<> dự báo)

#Đường thẳng màu đỏ là 1 dạng KIỂM ĐỊNH GIẢ ĐỊNH, đường đỏ trên hình ta thấy không thẳng hàng, chứng tỏ có 1 sự sai lệch của dự báo so với thực tế

head(ob) #mặc định là 6 dòng
##   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

Vẽ các biểu đồ, để đánh giá tổng quan trước về dữ liệu, trước khi bước đến phân tích, chọn mô hình kiểm định

#Vẽ biểu đồ tương quan 2 biến
plot(ob$pcfat ~ ob$age, pch = 16, col = 'Green') #plot là hàm base, pch/pchisq: độ tương phản màu

#Qua hình trên ta thấy, tuổi tăng có vẻ pcfat cũng tăng theo => khả năng có tương quan theo hồi quy tuyến tính ( y =aX )

Vẽ biểu đồ hộp ( box plot)

boxplot(ob$pcfat) #boxplot là hàm base, mô tả trung vị, bách phân vị 75,25, outlier

boxplot(ob$pcfat ~ ob$gender, col = "blue", border = "green")#phân ra 2 giới

#Ảnh hưởng của tuổi đến tỉ trọng mỡ (pcfat)
a1 = lm(pcfat ~ age, data = ob)
summary(a1)
## 
## Call:
## lm(formula = pcfat ~ age, data = ob)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.7114  -4.0069   0.8378   4.9514  18.2192 
## 
## 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
#Mô hình hồi quy tuyến tính trên có  α = 25.58 và β = 0.13 trong phương trình Y =  α + β(age) + ε
# Giữa các giá trị của biến X ( ở đây là age) có sự khác biệt là 0.12% về tỉ trọng mỡ (p<0.0001-làm tròn từ 2.2e-16) => tuổi tăng 1 đơn vị thì tỉ trọng mỡ tăng 0.12% ( p < 0.0001)
# R-Squared: 0.09 -> tức chỉ có 9% phương sai ( dao động về giá trị) của biến Y (ở đây là tỉ trọng mỡ) có thể được giải thích bằng biến X 
plot(ob$pcfat ~ ob$age, col = "blue", pch = 16 )#Biểu độ tương quan (phân tán) của pcfat theo age
a1 = lm(pcfat ~ age, data = ob)#Tạo mô hình/phương trình tuyến tính pcfat theo age lưu vào biến a1
abline(a1, col ="red")#Vẽ biểu đồ biến/mô hình/phương trình a1

#Đường đỏ chính là biểu diễn của phương trình Y(pcfat) = 25.6(α) + 0.13(β)*X(tuổi)

KIỂM TRA GIẢ ĐỊNH QUA “autoplot” ( Bước cuối cùng)

m1 = lm(pcfat ~ age, data = ob) #Mô hình hồi quy tuyến tính
library(ggfortify) #Gọi hàm ggfortify
autoplot(m1)

#Ta thấy đường dự báo ( nét đứt ) và thực tế (xanh) khác nhau

MÔ HÌNH HỒI QUI TUYẾN TÍNH ĐA BIẾN Công thức chung : lm (y ~ x1 + x2 + x3 + …)

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
## The following object is masked from 'package:psych':
## 
##     cushny
## Loading required package: inline
## Loading required package: rrcov
## Scalable Robust Estimators with High Breakdown Point (version 1.4-7)
yvar = ob[,c("pcfat")]
xvar = ob[,c('gender', 'height', 'weight', 'bmi', 'age')]
bma = bicreg(xvar, yvar, strict = FALSE, OR = 20)
summary(bma)
## 
## Call:
## bicreg(x = xvar, y = yvar, 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
## 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
#nVar: số lượng biến trong mô hình
#r2: R-squared ?
#post prob: xác suất hậu định - là đại lượng quan trọng nhất
#Trên bảng ta thấy, với tổng 31 mô hình có thể có, qua kiểm định Bayesian ta chọn được 3 mô hình tốt ưu nhất, và ta đặt tổng xác suất hậu định tích lũy của 3 mô hình này là 1. Ta xét tiếp 3 mô hình ấy, ta thấy model 1 có post prob là cao nhất với 39.2% -> tức 39,2% model 1 sẽ xuất hiện trong việc lập 3 mô hình tích lũy tối nhất

#Vẽ biểu đồ

imageplot.bma(bma)#tạo hình ảnh của những mô hình được chọn lựa bởi bayesian