Bước 1: Phân tích mô tả/khai thác (exploratory analysis)

Nhập dữ liệu

ob = read.csv("C:/Users/Dell/Desktop/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 các biến có thể liên quan ( dựa vào kinh nghiệm chuyên ngành )

tl = ob[, c('gender', 'age', "height", 'weight', 'bmi', 'pcfat')]
head(tl)

##   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

Đánh giá ‘các biến liên quan theo kinh nghiệm’ dựa vào tương quan đa biến để xác nhận lại

library(GGally)

## Loading required package: ggplot2

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

ggpairs(tl)# có thể thêm 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`.

Ta thấy coefficient of correlation của các biến với biến pcfat đều khá ổn, > 0.05 => pcfat và 4 biến còn lại có mối tương quan với nhau -> 4 biến age, height, weight, bmi và biến phân loại là gender đều có thể có tác dụng trong mô hình tiên lượng
Tương quan đa biến đúng ra các biến phải liên tục, vậy sao lại nghĩ dùng cả gender
Có cơ sở sinh lý bệnh và thực ra có thể test xem gender và pcfat có liên quan nhau không vì pcfat là liên tục, gender là phân loại, ta đang mong muốn xây dựng tiên lượng cho pcfat -> thử đặt mô hình tuyến tính với pcfat phụ thuộc gender
lm(pcfat ~ gender, data = ob),rồi summary() hay anova() rồi xét RMSE hay tốt nhất là R-squared -> kết quả R-squared là 0.44 => khác biệt nam nữ giải thích 44% những khác biệt về phương sai của tỉ trọng mỡ

Tất cả bước trên mới chỉ nằm ở phân tích mô tả/khai thác

Bước 2: Tìm các biến liên quan tối ưu bằng thống kê, ở đây dùng phương pháp bayesian

Gọi thư viện ‘BMA’ : thư viện thực hiện Bayesian

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

#Mục đích hiện nay là tiên lượng cho pcfat
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)
#strict = F : vẫn trả về các mô hình khác mà thỏa mãn 1/OR của best model, còn để "True" sẽ loại những cái không thỏa
#Tuy là False những các mô hình vẫn thỏa điều kiện 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

Vậy có 3 mô hình phù hợp và các biến tương ứng, biến nào không có trong mô hình nào thì giá trị tương ứng trong mô hình đó là ‘.’
Việc chọn mô hình nào phù hợp ngoài giá trị thống kế post prob còn phải phù hợp với kinh nghiệm chuyên ngành, giả dụ weight mình có đo trực tiếp không hay nghe lời khai của bệnh nhân, vậy khả năng ứng dụng trên thực tế có hiệu quả lâm sàng, không sai số khi thu nhận giá trị nhiều, không tốn tài nguyên công sức quá mức

Vẽ hình các mô đình được chọn cho dễ hình dung

imageplot.bma(bma)

Bước 3: Chia dữ liệu thành 2 nhóm: development và validation

Gọi thư viện caret

library(caret) # Đây là thư viện cực kì quan trong trong xây dựng mô hình, kể cả là logistic regression, vì cũng ứng dụng trong chia data để training và testing

## Loading required package: lattice

## 
## Attaching package: 'caret'

## The following object is masked from 'package:survival':
## 
##     cluster

Chia dữ liệu dùng hàm createDataPartition

sample = createDataPartition(ob$pcfat, p=0.6, list=F) #sample sẽ chiếm 0.6 dữ liệu
dev = ob[sample, ] #dev sẽ là dữ liệu từ ob với các thông số sample phía trên
val = ob[-sample, ] #val sẽ là dữ liệu từ ob mà không bao gồm 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ố

Bước 4: Huấn luyện (xây dựng) mô hình

Ở đây ta dùng gói caret thay cho làm thủ công

control = trainControl (method = "cv", number = 10)
#method = 'cv" : cv ở đây là ám chỉ phương pháp crossvalidaion
#number ở đây là số lần làm, có thể set tùy thích
training = train (pcfat ~ gender + age + bmi + weight, data= dev, method = "lm", trControl = control, metric = "Rsquared")

summary(training)

## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.2229  -2.5569   0.0139   2.6727  15.6470 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.749287   1.127707   7.758 2.92e-14 ***
## genderM     -11.374452   0.457747 -24.849  < 2e-16 ***
## age           0.052255   0.009901   5.278 1.73e-07 ***
## bmi           0.831568   0.113666   7.316 6.79e-13 ***
## weight        0.091365   0.040111   2.278    0.023 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.075 on 726 degrees of freedom
## Multiple R-squared:  0.6846, Adjusted R-squared:  0.6829 
## F-statistic:   394 on 4 and 726 DF,  p-value: < 2.2e-16

Đánh giá tầm quan trọng của từng biến

#Cần dùng thư viện relaimpo
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 + weight, 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 weight 
## Proportion of variance explained by model: 69.66%
## Metrics are normalized to sum to 100% (rela=TRUE). 
## 
## Relative importance metrics: 
## 
##               lmg
## gender 0.59317775
## age    0.06893066
## bmi    0.24613695
## weight 0.09175463
## 
## Average coefficients for different model sizes: 
## 
##                  1X          2Xs          3Xs          4Xs
## gender -10.51634414 -11.71834412 -11.80453842 -11.44430262
## age      0.12768705   0.10445197   0.05168496   0.05496623
## bmi      1.03619023   1.50631405   1.54278433   0.89419395
## weight   0.04319324  -0.05539405  -0.06907993   0.07920690

Bước 5: Kiểm tra mô hình tiên lượng

Ở đây ta dùng caret

pred = predict(training, newdata=val)
plot(pred ~ val$pcfat, pch=16)#Vẽ biểu đồ tương quan giữa các giá trị của pred và val$pcfat( pcfat trong tập dữ liệu kiểm tra)

Biểu đồ thể hiện mối tương quan giữa kết quả quan sát trên thực tế và kết quả có được từ mô hình tiên lượng

Ta cụ thể con số để chính xác về thống kê

validation = data.frame(obs=val$pcfat, pred)
head(validation)

##     obs     pred
## 2  16.8 21.40531
## 4  33.8 34.64607
## 6  32.2 32.47519
## 7  35.3 37.53804
## 9  21.1 21.96112
## 12 28.6 25.80076

Tính chỉ số Rsquared

defaultSummary(validation)

##      RMSE  Rsquared       MAE 
## 3.7997057 0.7163455 2.9651441

R-squared bây giờ là 0.66 thay đổi không đáng kể so với R-squared có được từ training 0.7182

Vậy mô hình thu được bây giờ khá tốt, đủ để viết luận án tiến sĩ, còn ứng dụng lâm sàng thì còn khá thấp, thường phải 0.8 trở lên

Linear regression - Ứng dụng xây dựng mô hình tiên lượng có dùng Bayesian