1 Chuẩn bị môi trường

1.1 Load các package

library(carData)
library(car)
library(caret)
library(rms)
library(GGally)
library(olsrr)
library(BMA)
library(glmnet)

2 Mô hình hồi quy tuyến tính và hiệu chỉnh yếu tố nhiễu

2.1 Đọc và khảo sát dữ liệu `Salaries`

data("Salaries")

dim(Salaries)

## [1] 397   6

head(Salaries, 6)

str(Salaries)

## 'data.frame':    397 obs. of  6 variables:
##  $ rank         : Factor w/ 3 levels "AsstProf","AssocProf",..: 3 3 1 3 3 2 3 3 3 3 ...
##  $ discipline   : Factor w/ 2 levels "A","B": 2 2 2 2 2 2 2 2 2 2 ...
##  $ yrs.since.phd: int  19 20 4 45 40 6 30 45 21 18 ...
##  $ yrs.service  : int  18 16 3 39 41 6 23 45 20 18 ...
##  $ sex          : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 1 ...
##  $ salary       : int  139750 173200 79750 115000 141500 97000 175000 147765 119250 129000 ...

Dataset Salaries gồm các thông tin liên quan đến hạng giáo sư, chuyên ngành, số năm sau tiến sĩ, số năm công tác, giới tính và mức lương.

2.2 Mô hình 1: Mối liên quan giữa giới tính và lương

Mô hình:

\[ Salary = \beta_0 + \beta_1 Sex + \epsilon \]

model1 <- lm(
  salary ~ sex,
  data = Salaries
)

summary(model1)

## 
## Call:
## lm(formula = salary ~ sex, data = Salaries)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -57290 -23502  -6828  19710 116455 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   101002       4809  21.001  < 2e-16 ***
## sexMale        14088       5065   2.782  0.00567 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30030 on 395 degrees of freedom
## Multiple R-squared:  0.01921,    Adjusted R-squared:  0.01673 
## F-statistic: 7.738 on 1 and 395 DF,  p-value: 0.005667

2.2.1 Hệ số hồi quy

coef(model1)

## (Intercept)     sexMale 
##   101002.41    14088.01

2.2.2 Khoảng tin cậy 95%

confint(model1)

##                 2.5 %    97.5 %
## (Intercept) 91547.216 110457.60
## sexMale      4131.107  24044.91

2.2.3 Giá trị dự đoán

head(predict(model1), 10)

##        1        2        3        4        5        6        7        8        9       10 
## 115090.4 115090.4 115090.4 115090.4 115090.4 115090.4 115090.4 115090.4 115090.4 101002.4

2.2.4 Kiểm tra giả định của mô hình

par(
  mfrow = c(2, 2),
  mar = c(4, 4, 2, 1)
)

plot(model1)

par(mfrow = c(1, 1))

2.3 Mô hình 2: Mối liên quan giữa thời gian sau tiến sĩ và lương

Mô hình:

\[ Salary = \beta_0 + \beta_1 YearsSincePhD + \epsilon \]

model2 <- lm(
  salary ~ yrs.since.phd,
  data = Salaries
)

summary(model2)

## 
## Call:
## lm(formula = salary ~ yrs.since.phd, data = Salaries)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -84171 -19432  -2858  16086 102383 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    91718.7     2765.8  33.162   <2e-16 ***
## yrs.since.phd    985.3      107.4   9.177   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 27530 on 395 degrees of freedom
## Multiple R-squared:  0.1758, Adjusted R-squared:  0.1737 
## F-statistic: 84.23 on 1 and 395 DF,  p-value: < 2.2e-16

2.3.1 Hệ số hồi quy

coef(model2)

##   (Intercept) yrs.since.phd 
##    91718.6854      985.3421

2.3.2 Khoảng tin cậy 95%

confint(model2)

##                    2.5 %    97.5 %
## (Intercept)   86281.1714 97156.199
## yrs.since.phd   774.2636  1196.421

2.3.3 Giá trị dự đoán

head(predict(model2), 10)

##         1         2         3         4         5         6         7         8         9        10 
## 110440.19 111425.53  95660.05 136059.08 131132.37  97630.74 121278.95 136059.08 112410.87 109454.84

2.3.4 Biểu đồ mối liên quan

plot(
  Salaries$yrs.since.phd,
  Salaries$salary,
  main = "Mối liên quan giữa thời gian sau tiến sĩ và lương",
  xlab = "Số năm sau tiến sĩ",
  ylab = "Salary",
  pch = 19,
  col = "steelblue"
)

abline(
  model2,
  col = "red",
  lwd = 2
)

2.3.5 Kiểm tra giả định của mô hình

par(
  mfrow = c(2, 2),
  mar = c(4, 4, 2, 1)
)

plot(model2)

par(mfrow = c(1, 1))

2.4 Mô hình 3: Giới tính và lương sau khi hiệu chỉnh thời gian sau tiến sĩ

Mô hình:

\[ Salary = \beta_0 + \beta_1 Sex + \beta_2 YearsSincePhD + \epsilon \]

model3 <- lm(
  salary ~ sex + yrs.since.phd,
  data = Salaries
)

summary(model3)

## 
## Call:
## lm(formula = salary ~ sex + yrs.since.phd, data = Salaries)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -84167 -19735  -2551  15427 102033 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    85181.8     4748.3  17.939   <2e-16 ***
## sexMale         7923.6     4684.1   1.692   0.0915 .  
## yrs.since.phd    958.1      108.3   8.845   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 27470 on 394 degrees of freedom
## Multiple R-squared:  0.1817, Adjusted R-squared:  0.1775 
## F-statistic: 43.74 on 2 and 394 DF,  p-value: < 2.2e-16

2.4.1 Hệ số hồi quy

coef(model3)

##   (Intercept)       sexMale yrs.since.phd 
##    85181.8181     7923.6239      958.0793

2.4.2 Khoảng tin cậy 95%

confint(model3)

##                    2.5 %    97.5 %
## (Intercept)   75846.6147 94517.021
## sexMale       -1285.2949 17132.543
## yrs.since.phd   745.1244  1171.034

2.4.3 Giá trị dự đoán

head(predict(model3), 10)

##         1         2         3         4         5         6         7         8         9        10 
## 111308.95 112267.03  96937.76 136219.01 131428.62  98853.92 121847.82 136219.01 113225.11 102427.25

2.4.4 Kiểm tra giả định

par(
  mfrow = c(2, 2),
  mar = c(4, 4, 2, 1)
)

plot(model3)

par(mfrow = c(1, 1))

2.5 Mô hình 4: Rank, discipline và sex

Mô hình:

\[ Salary = \beta_0 + \beta_1 Sex + \beta_2 Discipline + \beta_3 Rank + \epsilon \]

model4 <- lm(
  salary ~ sex + discipline + rank,
  data = Salaries
)

summary(model4)

## 
## Call:
## lm(formula = salary ~ sex + discipline + rank, data = Salaries)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -66268 -14127  -1566  10813  97718 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      68224       4477  15.238  < 2e-16 ***
## sexMale           4492       3860   1.164 0.245291    
## disciplineB      13709       2295   5.972 5.25e-09 ***
## rankAssocProf    13723       3959   3.466 0.000586 ***
## rankProf         47403       3133  15.130  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22640 on 392 degrees of freedom
## Multiple R-squared:  0.4469, Adjusted R-squared:  0.4412 
## F-statistic: 79.18 on 4 and 392 DF,  p-value: < 2.2e-16

2.5.1 Hệ số hồi quy

coef(model4)

##   (Intercept)       sexMale   disciplineB rankAssocProf      rankProf 
##     68223.531      4491.801     13708.688     13723.416     47403.317

2.5.2 Khoảng tin cậy 95%

confint(model4)

##                   2.5 %   97.5 %
## (Intercept)   59421.202 77025.86
## sexMale       -3097.562 12081.16
## disciplineB    9195.780 18221.60
## rankAssocProf  5939.859 21506.97
## rankProf      41243.519 53563.11

2.5.3 Đánh giá mức độ quan trọng của từng biến

Anova(
  model4,
  type = 2
)

2.5.4 Biểu đồ hệ số hồi quy

coef_model <- coef(model4)[-1]

names(coef_model) <- c(
  "Sex: Male",
  "Discipline B",
  "Rank: AssocProf",
  "Rank: Prof"
)

barplot(
  coef_model,
  horiz = TRUE,
  main = "Ảnh hưởng của các yếu tố đến lương",
  xlab = "Hệ số hồi quy (Estimate)",
  col = c(
    "deepskyblue",
    "orange",
    "lightgreen",
    "violet"
  ),
  border = "gray30",
  las = 1,
  cex.names = 1,
  xlim = c(0, max(coef_model) * 1.2)
)

abline(
  v = 0,
  col = "red",
  lty = 2,
  lwd = 2
)

3 Đánh giá mô hình tiên lượng

3.1 Đọc dữ liệu mô phỏng béo phì

ob <- read.csv(
  "D:/RData/Sim obesity data.csv"
)

dim(ob)

## [1] 1217   12

head(ob)

str(ob)

## 'data.frame':    1217 obs. of  12 variables:
##  $ X        : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ id       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ gioi     : chr  "F" "F" "F" "M" ...
##  $ chieucao : num  150 154 154 161 156 ...
##  $ cannang  : num  49.7 61.8 46.4 61.3 59 ...
##  $ bmi      : num  22.1 26 19.5 23.5 24.1 ...
##  $ tuoi     : num  81.4 13.5 67.6 80.5 61.1 ...
##  $ bmc      : num  1117 1805 1450 2106 1452 ...
##  $ bmd      : num  0.727 1.076 0.908 1.129 0.886 ...
##  $ mo       : num  19951 21805 15245 17816 26822 ...
##  $ co       : num  32106 33838 29262 42930 30095 ...
##  $ titrongmo: num  36.6 37.4 34.8 27.8 45 ...

3.2 Xây dựng mô hình tiên lượng

Biến kết quả là titrongmo.

Các biến dự đoán gồm:

gioi
tuoi
bmi

modelob <- lm(
  titrongmo ~ gioi + tuoi + bmi,
  data = ob
)

summary(modelob)

## 
## Call:
## lm(formula = titrongmo ~ gioi + tuoi + bmi, data = ob)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.7956  -2.6381  -0.0685   2.7692  14.3410 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.137347   0.852356   9.547  < 2e-16 ***
## gioiM       -10.806252   0.254068 -42.533  < 2e-16 ***
## tuoi          0.047147   0.006852   6.881 9.53e-12 ***
## bmi           1.089355   0.038508  28.289  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.974 on 1213 degrees of freedom
## Multiple R-squared:  0.6947, Adjusted R-squared:  0.6939 
## F-statistic: 919.9 on 3 and 1213 DF,  p-value: < 2.2e-16

3.2.1 Hệ số hồi quy

coef(modelob)

##  (Intercept)        gioiM         tuoi          bmi 
##   8.13734710 -10.80625248   0.04714735   1.08935502

3.2.2 Khoảng tin cậy 95%

confint(modelob)

##                    2.5 %       97.5 %
## (Intercept)   6.46509220   9.80960201
## gioiM       -11.30471340 -10.30779156
## tuoi          0.03370369   0.06059102
## bmi           1.01380613   1.16490391

3.2.3 Giá trị tiên lượng

head(
  predict(modelob),
  10
)

##        1        2        3        4        5        6        7        8        9       10 
## 36.02263 37.08244 32.61431 26.73402 37.29437 38.16709 37.24439 36.12422 24.90434 25.34778

3.2.4 Kiểm tra giả định mô hình

par(
  mfrow = c(2, 2),
  mar = c(4, 4, 2, 1)
)

plot(modelob)

par(mfrow = c(1, 1))

4 Chia dữ liệu Development và Test bằng `caret`

Dữ liệu được chia theo tỷ lệ:

Development: 60%
Test: 40%

set.seed(123)

index <- createDataPartition(
  ob$titrongmo,
  p = 0.6,
  list = FALSE
)

dev <- ob[index, ]

test <- ob[-index, ]

4.0.1 Kích thước dữ liệu

cat("Kích thước dữ liệu gốc:\n")

## Kích thước dữ liệu gốc:

dim(ob)

## [1] 1217   12

cat("\nKích thước Development:\n")

## 
## Kích thước Development:

dim(dev)

## [1] 732  12

cat("\nKích thước Test:\n")

## 
## Kích thước Test:

dim(test)

## [1] 485  12

4.1 Xây dựng mô hình trên tập Development

modeldev <- lm(
  titrongmo ~ gioi + tuoi + bmi,
  data = dev
)

summary(modeldev)

## 
## Call:
## lm(formula = titrongmo ~ gioi + tuoi + bmi, data = dev)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.0662  -2.8517   0.0328   2.8488  14.5981 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.765301   1.100168   7.967 6.26e-15 ***
## gioiM       -11.204788   0.327981 -34.163  < 2e-16 ***
## tuoi          0.039615   0.008958   4.422 1.12e-05 ***
## bmi           1.085899   0.049618  21.885  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.961 on 728 degrees of freedom
## Multiple R-squared:  0.7072, Adjusted R-squared:  0.706 
## F-statistic: 586.2 on 3 and 728 DF,  p-value: < 2.2e-16

4.2 Dự đoán trên tập Test

pred_test <- predict(
  modeldev,
  newdata = test
)

result_test <- data.frame(
  Thuc_te = test$titrongmo,
  Du_doan = pred_test
)

head(result_test, 10)

4.2.1 So sánh giá trị thực tế và dự đoán

knitr::kable(
  head(result_test, 20),
  digits = 2,
  caption = "Giá trị thực tế và giá trị dự đoán trên tập Test"
)

Giá trị thực tế và giá trị dự đoán trên tập Test
	Thuc_te	Du_doan
1	36.62	35.96
2	37.36	37.52
4	27.81	26.28
6	38.74	38.46
7	35.85	37.38
10	34.81	25.23
11	24.90	29.46
12	32.71	31.33
14	22.55	22.55
21	26.85	24.74
25	38.74	36.93
28	32.41	35.44
30	29.83	37.16
33	25.16	26.13
36	34.69	34.31
41	33.54	34.42
43	37.83	40.21
44	32.85	37.54
47	31.38	33.82
48	23.70	26.96

4.3 Đánh giá sai số dự đoán

4.3.1 RMSE

RMSE <- sqrt(
  mean(
    (result_test$Thuc_te - result_test$Du_doan)^2
  )
)

RMSE

## [1] 4.009811

4.3.2 MAE

MAE <- mean(
  abs(
    result_test$Thuc_te - result_test$Du_doan
  )
)

MAE

## [1] 3.144478

4.3.3 R-squared trên tập Test

R2_test <- cor(
  result_test$Thuc_te,
  result_test$Du_doan
)^2

R2_test

## [1] 0.6734275

4.3.4 Tổng hợp kết quả kiểm định

performance_test <- data.frame(
  Chi_so = c(
    "RMSE",
    "MAE",
    "R-squared"
  ),
  Gia_tri = c(
    RMSE,
    MAE,
    R2_test
  )
)

knitr::kable(
  performance_test,
  digits = 4,
  caption = "Đánh giá khả năng tiên lượng trên tập Test"
)

Đánh giá khả năng tiên lượng trên tập Test
Chi_so	Gia_tri
RMSE	4.0098
MAE	3.1445
R-squared	0.6734

5 Bootstrap validation bằng package `rms`

5.1 Thiết lập `datadist`

dd <- datadist(ob)

options(
  datadist = "dd"
)

5.2 Xây dựng mô hình OLS

modelboot <- ols(
  titrongmo ~ gioi + tuoi + bmi,
  data = ob,
  x = TRUE,
  y = TRUE
)

modelboot

## Linear Regression Model
## 
## ols(formula = titrongmo ~ gioi + tuoi + bmi, data = ob, x = TRUE, 
##     y = TRUE)
## 
##                  Model Likelihood    Discrimination    
##                        Ratio Test           Indexes    
## Obs    1217    LR chi2    1443.82    R2       0.695    
## sigma3.9739    d.f.             3    R2 adj   0.694    
## d.f.   1213    Pr(> chi2)  0.0000    g        6.785    
## 
## Residuals
## 
##       Min        1Q    Median        3Q       Max 
## -12.79562  -2.63806  -0.06854   2.76916  14.34105 
## 
## 
##           Coef     S.E.   t      Pr(>|t|)
## Intercept   8.1373 0.8524   9.55 <0.0001 
## gioi=M    -10.8063 0.2541 -42.53 <0.0001 
## tuoi        0.0471 0.0069   6.88 <0.0001 
## bmi         1.0894 0.0385  28.29 <0.0001

5.3 Bootstrap 200 lần

set.seed(123)

bootstrap_result <- validate(
  modelboot,
  method = "boot",
  B = 200
)

bootstrap_result

##           index.orig training    test optimism index.corrected   Lower   Upper   n
## R-square      0.6947   0.6965  0.6936   0.0030          0.6917  0.6664  0.7202 200
## MSE          15.7399  15.5980 15.7962  -0.1982         15.9381 14.6795 17.1570 200
## g             6.7849   6.7824  6.7812   0.0012          6.7837  6.4318  7.1196 200
## Intercept     0.0000   0.0000  0.0250  -0.0250          0.0250 -1.6113  1.5853 200
## Slope         1.0000   1.0000  0.9993   0.0007          0.9993  0.9535  1.0490 200

6 Đánh giá đa cộng tuyến

6.1 Biểu đồ đa biến bằng `GGally`

ggpairs(
  Salaries[, c(
    "salary",
    "sex",
    "rank",
    "discipline",
    "yrs.since.phd",
    "yrs.service"
  )],
  title = "Mối liên quan giữa các biến trong dữ liệu Salaries"
)

6.2 Mô hình hồi quy đa biến

Đánh giá mối liên quan giữa giới tính và lương sau khi hiệu chỉnh:

rank
discipline
yrs.since.phd
yrs.service

modelTest <- lm(
  salary ~
    sex +
    rank +
    discipline +
    yrs.since.phd +
    yrs.service,
  data = Salaries
)

summary(modelTest)

## 
## Call:
## lm(formula = salary ~ sex + rank + discipline + yrs.since.phd + 
##     yrs.service, data = Salaries)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -65248 -13211  -1775  10384  99592 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    65955.2     4588.6  14.374  < 2e-16 ***
## sexMale         4783.5     3858.7   1.240  0.21584    
## rankAssocProf  12907.6     4145.3   3.114  0.00198 ** 
## rankProf       45066.0     4237.5  10.635  < 2e-16 ***
## disciplineB    14417.6     2342.9   6.154 1.88e-09 ***
## yrs.since.phd    535.1      241.0   2.220  0.02698 *  
## yrs.service     -489.5      211.9  -2.310  0.02143 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22540 on 390 degrees of freedom
## Multiple R-squared:  0.4547, Adjusted R-squared:  0.4463 
## F-statistic:  54.2 on 6 and 390 DF,  p-value: < 2.2e-16

6.2.1 Hệ số hồi quy

coef(modelTest)

##   (Intercept)       sexMale rankAssocProf      rankProf   disciplineB yrs.since.phd   yrs.service 
##    65955.2324     4783.4928    12907.5879    45065.9987    14417.6256      535.0583     -489.5157

6.2.2 Khoảng tin cậy 95%

confint(modelTest)

##                     2.5 %      97.5 %
## (Intercept)   56933.74315 74976.72156
## sexMale       -2802.90123 12369.88691
## rankAssocProf  4757.69991 21057.47589
## rankProf      36734.75104 53397.24631
## disciplineB    9811.37975 19023.87139
## yrs.since.phd    61.24805  1008.86852
## yrs.service    -906.19882   -72.83261

6.3 Kiểm tra giả định mô hình đa biến

par(
  mfrow = c(2, 2),
  mar = c(4, 4, 2, 1)
)

plot(
  modelTest,
  which = 1
)

plot(
  modelTest,
  which = 2
)

plot(
  modelTest,
  which = 3
)

plot(
  modelTest,
  which = 5
)

par(mfrow = c(1, 1))

6.4 Đánh giá đa cộng tuyến bằng VIF

vif_result <- vif(modelTest)

vif_result

##       sexMale rankAssocProf      rankProf   disciplineB yrs.since.phd   yrs.service 
##      1.030805      1.815873      3.102637      1.064105      7.518936      5.923038

6.4.1 Trình bày VIF dạng bảng

vif_table <- data.frame(
  Bien = names(vif_result),
  VIF = as.numeric(vif_result)
)

knitr::kable(
  vif_table,
  digits = 3,
  caption = "Variance Inflation Factor (VIF)"
)

Variance Inflation Factor (VIF)
Bien	VIF
sexMale	1.031
rankAssocProf	1.816
rankProf	3.103
disciplineB	1.064
yrs.since.phd	7.519
yrs.service	5.923

6.4.2 Nguyên tắc diễn giải

VIF < 5: chưa có đa cộng tuyến đáng kể.
VIF từ 5 đến 10: có dấu hiệu đa cộng tuyến.
VIF > 10: đa cộng tuyến nghiêm trọng.

7 Lựa chọn mô hình tối ưu

7.1 Đọc dữ liệu BMA

df <- read.csv(
  "D:/RData/Sim data for BMA.csv"
)

dim(df)

## [1] 2000   53

nrow(df)

## [1] 2000

ncol(df)

## [1] 53

head(df)

names(df)

##  [1] "ID"  "Sex" "y"   "X01" "X02" "X03" "X04" "X05" "X06" "X07" "X08" "X09" "X10" "X11" "X12" "X13" "X14" "X15" "X16"
## [20] "X17" "X18" "X19" "X20" "X21" "X22" "X23" "X24" "X25" "X26" "X27" "X28" "X29" "X30" "X31" "X32" "X33" "X34" "X35"
## [39] "C01" "C02" "C03" "C04" "C05" "C06" "C07" "C08" "C09" "C10" "C11" "C12" "C13" "C14" "C15"

str(df)

## 'data.frame':    2000 obs. of  53 variables:
##  $ ID : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Sex: chr  "Male" "Male" "Male" "Male" ...
##  $ y  : num  17.2 19 18.4 22.3 30.9 ...
##  $ X01: num  0.998 -0.893 -2.398 0.878 0.558 ...
##  $ X02: num  0.657 0.36 -1.097 -0.284 1.224 ...
##  $ X03: num  -2.14934 0.22864 0.03713 -0.00807 0.40362 ...
##  $ X04: num  1.15861 -0.14192 -0.74704 0.00531 0.95459 ...
##  $ X05: num  1.782 1.26 -0.561 -0.72 0.49 ...
##  $ X06: num  -0.541 2.36 0.16 -0.147 -1.049 ...
##  $ X07: num  0.272 0.998 0.397 1.063 -1.637 ...
##  $ X08: num  -0.176 -0.374 -0.193 1.303 0.179 ...
##  $ X09: num  -0.219 -0.779 -0.758 1.175 0.198 ...
##  $ X10: num  -0.973 0.547 1.544 -0.799 0.318 ...
##  $ X11: num  -0.9468 -0.2579 0.0407 0.5013 0.1396 ...
##  $ X12: num  1.334 0.483 -0.417 -0.171 -0.454 ...
##  $ X13: num  0.578 0.259 0.47 0.508 -0.769 ...
##  $ X14: num  -0.577 0.491 0.451 -0.33 -1.275 ...
##  $ X15: num  2.1303 -0.101 0.6644 -0.8729 0.0114 ...
##  $ X16: num  0.987 0.578 0.842 1.156 1.561 ...
##  $ X17: num  0.4315 1.2801 0.0873 1.2728 0.822 ...
##  $ X18: num  1.619 -0.453 -0.588 1.039 1.516 ...
##  $ X19: num  0.7806 -0.1559 -0.4908 0.375 0.0112 ...
##  $ X20: num  0.735 1.101 1.553 -1.283 1.331 ...
##  $ X21: num  0.501 0.528 0.547 0.427 0.646 ...
##  $ X22: num  0.449 0.1834 0.2882 0.3074 -0.0228 ...
##  $ X23: num  0.0287 0.6926 -0.4956 -0.8664 -1.2179 ...
##  $ X24: num  -0.973 -0.378 -0.361 -1.09 -1.272 ...
##  $ X25: num  -1.415 1.551 -0.24 0.767 -0.898 ...
##  $ X26: num  1.026 0.377 -0.394 0.363 -0.516 ...
##  $ X27: num  1.041 -2.523 1.081 0.327 -0.865 ...
##  $ X28: num  0.853 0.706 1.445 -1.314 0.877 ...
##  $ X29: num  -0.876 2.051 -0.21 0.652 0.455 ...
##  $ X30: num  -1.1261 -1.613 -0.7183 -0.9661 -0.0798 ...
##  $ X31: num  -0.38 -0.304 0.385 -0.783 0.358 ...
##  $ X32: num  0.685 -0.974 0.966 -0.563 0.72 ...
##  $ X33: num  0.0649 -1.0697 -0.7381 -0.389 -0.6131 ...
##  $ X34: num  0.581 -0.611 -1.003 -0.445 -0.639 ...
##  $ X35: num  0.9616 -0.9248 1.4258 -0.4652 -0.0227 ...
##  $ C01: chr  "No" "No" "No" "Yes" ...
##  $ C02: chr  "No" "Yes" "No" "Yes" ...
##  $ C03: chr  "No" "Yes" "No" "No" ...
##  $ C04: chr  "Yes" "Yes" "Yes" "No" ...
##  $ C05: chr  "No" "No" "Yes" "No" ...
##  $ C06: chr  "Yes" "Yes" "No" "No" ...
##  $ C07: chr  "Yes" "No" "Yes" "No" ...
##  $ C08: chr  "Yes" "No" "No" "Yes" ...
##  $ C09: chr  "No" "Yes" "Yes" "Yes" ...
##  $ C10: chr  "Yes" "No" "No" "No" ...
##  $ C11: chr  "Yes" "No" "Yes" "Yes" ...
##  $ C12: chr  "Yes" "No" "No" "No" ...
##  $ C13: chr  "Yes" "No" "Yes" "No" ...
##  $ C14: chr  "No" "No" "Yes" "No" ...
##  $ C15: chr  "No" "No" "No" "No" ...

8 Stepwise Regression

Nếu ID chỉ là mã định danh, biến này được loại khỏi mô hình.

df_stepwise <- df[
  ,
  setdiff(
    names(df),
    "ID"
  )
]

8.1 Xây dựng mô hình đầy đủ

Dấu . có nghĩa là sử dụng tất cả các biến còn lại để giải thích y.

modelStepwise <- lm(
  y ~ .,
  data = df_stepwise
)

summary(modelStepwise)

## 
## Call:
## lm(formula = y ~ ., data = df_stepwise)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.4639  -2.2447  -0.1142   2.2701  12.8027 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.517e+01  2.819e-01 124.756   <2e-16 ***
## SexMale     -1.017e+01  1.742e-01 -58.356   <2e-16 ***
## X01          3.076e+00  8.326e-02  36.946   <2e-16 ***
## X02         -2.597e+00  8.910e-02 -29.145   <2e-16 ***
## X03          2.190e+00  9.138e-02  23.968   <2e-16 ***
## X04         -1.135e-01  9.417e-02  -1.205   0.2282    
## X05         -4.327e-02  9.008e-02  -0.480   0.6310    
## X06         -5.083e-03  9.182e-02  -0.055   0.9559    
## X07         -2.054e-02  9.100e-02  -0.226   0.8214    
## X08          7.177e-04  9.033e-02   0.008   0.9937    
## X09          1.060e-01  8.988e-02   1.180   0.2383    
## X10          6.344e-04  9.023e-02   0.007   0.9944    
## X11         -1.855e-01  9.208e-02  -2.015   0.0441 *  
## X12          1.331e-01  9.268e-02   1.436   0.1510    
## X13         -1.553e-01  9.358e-02  -1.660   0.0971 .  
## X14         -6.839e-03  9.110e-02  -0.075   0.9402    
## X15         -3.314e-04  9.324e-02  -0.004   0.9972    
## X16         -7.469e-02  8.963e-02  -0.833   0.4048    
## X17          3.265e-02  8.987e-02   0.363   0.7164    
## X18          6.455e-02  9.439e-02   0.684   0.4941    
## X19         -1.240e-03  9.033e-02  -0.014   0.9890    
## X20         -9.779e-02  8.846e-02  -1.105   0.2691    
## X21         -1.498e-01  9.247e-02  -1.620   0.1054    
## X22          1.128e-01  9.140e-02   1.234   0.2174    
## X23          8.059e-02  8.979e-02   0.898   0.3695    
## X24          2.240e-02  9.119e-02   0.246   0.8060    
## X25          6.932e-02  8.971e-02   0.773   0.4398    
## X26         -2.754e-02  8.883e-02  -0.310   0.7566    
## X27         -2.427e-02  8.916e-02  -0.272   0.7855    
## X28         -8.364e-02  9.004e-02  -0.929   0.3531    
## X29          1.139e-01  9.116e-02   1.249   0.2118    
## X30         -2.541e-02  8.798e-02  -0.289   0.7728    
## X31         -7.017e-02  8.846e-02  -0.793   0.4277    
## X32         -8.466e-02  9.033e-02  -0.937   0.3487    
## X33         -5.568e-02  9.089e-02  -0.613   0.5402    
## X34         -7.954e-02  9.167e-02  -0.868   0.3857    
## X35          8.287e-02  8.510e-02   0.974   0.3302    
## C01Yes       2.683e+00  1.681e-01  15.965   <2e-16 ***
## C02Yes      -2.342e+00  1.620e-01 -14.457   <2e-16 ***
## C03Yes      -2.151e-01  1.665e-01  -1.292   0.1965    
## C04Yes      -2.887e-02  1.699e-01  -0.170   0.8651    
## C05Yes       3.710e-01  1.742e-01   2.129   0.0334 *  
## C06Yes      -8.553e-02  1.705e-01  -0.502   0.6160    
## C07Yes      -7.746e-02  1.609e-01  -0.481   0.6303    
## C08Yes       1.237e-01  1.809e-01   0.684   0.4942    
## C09Yes      -2.556e-01  1.610e-01  -1.588   0.1124    
## C10Yes       1.848e-02  1.654e-01   0.112   0.9111    
## C11Yes      -5.354e-03  1.601e-01  -0.033   0.9733    
## C12Yes       1.145e-01  1.618e-01   0.708   0.4793    
## C13Yes      -7.027e-02  1.659e-01  -0.424   0.6719    
## C14Yes       1.056e-01  1.800e-01   0.587   0.5574    
## C15Yes      -1.572e-01  1.658e-01  -0.949   0.3430    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.536 on 1948 degrees of freedom
## Multiple R-squared:  0.7552, Adjusted R-squared:  0.7488 
## F-statistic: 117.9 on 51 and 1948 DF,  p-value: < 2.2e-16

8.2 Stepwise selection bằng `olsrr`

stepwise <- ols_step_both_p(
  modelStepwise
)

stepwise

## 
## 
##                                 Stepwise Summary                                
## ------------------------------------------------------------------------------
## Step    Variable         AIC          SBC         SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------------
##  0      Base Model    13493.949    13505.151    7815.576    0.00000    0.00000 
##  1      Sex (+)       12399.331    12416.133    6720.677    0.42207    0.42178 
##  2      X01 (+)       11903.063    11925.466    6224.265    0.54952    0.54907 
##  3      X02 (+)       11621.554    11649.558    5942.687    0.60906    0.60847 
##  4      X03 (+)       11165.903    11199.509    5488.156    0.68902    0.68840 
##  5      C01 (+)       10931.706    10970.912    5254.847    0.72366    0.72297 
##  6      C02 (+)       10733.477    10778.284    5057.794    0.74999    0.74923 
##  7      C05 (+)       10730.596    10781.004    5054.962    0.75059    0.74972 
##  8      X11 (+)       10729.063    10785.072    5053.477    0.75104    0.75003 
##  9      C09 (+)       10728.149    10789.759    5052.611    0.75140    0.75027 
##  10     X32 (+)       10727.269    10794.480    5051.785    0.75176    0.75051 
##  11     X21 (+)       10726.374    10799.186    5050.949    0.75211    0.75074 
## ------------------------------------------------------------------------------
## 
## Final Model Output 
## ------------------
## 
##                           Model Summary                           
## -----------------------------------------------------------------
## R                       0.867       RMSE                   3.512 
## R-Squared               0.752       MSE                   12.333 
## Adj. R-Squared          0.751       Coef. Var             11.008 
## Pred R-Squared          0.749       AIC                10726.374 
## MAE                     2.796       SBC                10799.186 
## -----------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
##  AIC: Akaike Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
## 
##                                  ANOVA                                   
## ------------------------------------------------------------------------
##                  Sum of                                                 
##                 Squares          DF    Mean Square       F         Sig. 
## ------------------------------------------------------------------------
## Regression    74841.886          11       6803.808    548.348    0.0000 
## Residual      24666.772        1988         12.408                      
## Total         99508.657        1999                                     
## ------------------------------------------------------------------------
## 
##                                     Parameter Estimates                                     
## -------------------------------------------------------------------------------------------
##       model       Beta    Std. Error    Std. Beta       t        Sig       lower     upper 
## -------------------------------------------------------------------------------------------
## (Intercept)     35.004         0.160                 218.825    0.000     34.690    35.318 
##     SexMale    -10.162         0.172       -0.660    -59.058    0.000    -10.500    -9.825 
##         X01      3.074         0.082        0.447     37.359    0.000      2.913     3.235 
##         X02     -2.597         0.088       -0.376    -29.639    0.000     -2.769    -2.425 
##         X03      2.142         0.084        0.304     25.436    0.000      1.977     2.307 
##      C01Yes      2.713         0.165        0.184     16.412    0.000      2.389     3.037 
##      C02Yes     -2.316         0.159       -0.163    -14.557    0.000     -2.628    -2.004 
##      C05Yes      0.381         0.172        0.025      2.214    0.027      0.044     0.719 
##         X11     -0.134         0.079       -0.019     -1.693    0.091     -0.290     0.021 
##      C09Yes     -0.267         0.158       -0.019     -1.685    0.092     -0.578     0.044 
##         X32     -0.142         0.079       -0.020     -1.784    0.075     -0.297     0.014 
##         X21     -0.135         0.079       -0.019     -1.697    0.090     -0.290     0.021 
## -------------------------------------------------------------------------------------------

9 Bayesian Model Averaging - BMA

9.1 Chuẩn bị biến kết quả Y

Y <- df$y

head(Y)

## [1] 17.24257 19.04815 18.37182 22.33739 30.85075 22.97769

9.2 Chuẩn bị các biến dự đoán X

Loại:

y: biến kết quả
ID: mã định danh

X_BMA <- df[
  ,
  setdiff(
    names(df),
    c(
      "y",
      "ID"
    )
  )
]

dim(X_BMA)

## [1] 2000   51

names(X_BMA)

##  [1] "Sex" "X01" "X02" "X03" "X04" "X05" "X06" "X07" "X08" "X09" "X10" "X11" "X12" "X13" "X14" "X15" "X16" "X17" "X18"
## [20] "X19" "X20" "X21" "X22" "X23" "X24" "X25" "X26" "X27" "X28" "X29" "X30" "X31" "X32" "X33" "X34" "X35" "C01" "C02"
## [39] "C03" "C04" "C05" "C06" "C07" "C08" "C09" "C10" "C11" "C12" "C13" "C14" "C15"

9.3 Chạy Bayesian Model Averaging

modelBMA <- bicreg(
  X_BMA,
  Y,
  strict = FALSE,
  OR = 20
)

summary(modelBMA)

## 
## Call:
## bicreg(x = X_BMA, y = Y, strict = FALSE, OR = 20)
## 
## 
##   11  models were selected
##  Best  5  models (cumulative posterior probability =  0.768 ): 
## 
##            p!=0    EV         SD       model 1     model 2     model 3     model 4     model 5   
## Intercept  100.0   34.983553  0.14376     34.9874     34.8730     34.9883     34.9913     34.9941
## SexMale    100.0  -10.156922  0.17253    -10.1564    -10.1636    -10.1574    -10.1515    -10.1567
## X01        100.0    3.071367  0.08248      3.0713      3.0708      3.0687      3.0702      3.0726
## X02        100.0   -2.597639  0.08766     -2.5982     -2.5892     -2.5920     -2.5973     -2.6004
## X03        100.0    2.144985  0.08518      2.1434      2.1396      2.1388      2.1482      2.1448
## X04          3.8   -0.005217  0.03138       .           .           .           .           .    
## X09          0.0    0.000000  0.00000       .           .           .           .           .    
## X11          6.3   -0.009428  0.04141       .           .         -0.1492       .           .    
## X12          0.0    0.000000  0.00000       .           .           .           .           .    
## X13          5.0   -0.007138  0.03614       .           .           .         -0.1437       .    
## X16          0.0    0.000000  0.00000       .           .           .           .           .    
## X18          0.0    0.000000  0.00000       .           .           .           .           .    
## X20          4.4   -0.005649  0.03081       .           .           .           .           .    
## X21          4.6   -0.006225  0.03308       .           .           .           .           .    
## X22          0.0    0.000000  0.00000       .           .           .           .           .    
## X23          0.0    0.000000  0.00000       .           .           .           .           .    
## X25          0.0    0.000000  0.00000       .           .           .           .           .    
## X28          0.0    0.000000  0.00000       .           .           .           .           .    
## X29          0.0    0.000000  0.00000       .           .           .           .           .    
## X31          2.7   -0.002921  0.02171       .           .           .           .           .    
## X32          5.0   -0.006894  0.03496       .           .           .           .         -0.1387
## X34          0.0    0.000000  0.00000       .           .           .           .           .    
## X35          0.0    0.000000  0.00000       .           .           .           .           .    
## C01Yes     100.0    2.718417  0.16567      2.7203      2.7222      2.7288      2.7120      2.7122
## C02Yes     100.0   -2.309810  0.15942     -2.3094     -2.3049     -2.3184     -2.3116     -2.3124
## C03Yes       3.0   -0.006978  0.04901       .           .           .           .           .    
## C05Yes      12.0    0.045644  0.13728       .          0.3804       .           .           .    
## C08Yes       0.0    0.000000  0.00000       .           .           .           .           .    
## C09Yes       4.8   -0.013091  0.06799       .           .           .           .           .    
## C12Yes       0.0    0.000000  0.00000       .           .           .           .           .    
## C15Yes       0.0    0.000000  0.00000       .           .           .           .           .    
##                                                                                                  
## nVar                                         6           7           7           7           7   
## r2                                         0.750       0.751       0.750       0.750       0.750 
## BIC                                    -2726.9033  -2724.1080  -2722.8254  -2722.3446  -2722.3446
## post prob                                  0.485       0.120       0.063       0.050       0.050

9.4 Xác suất biến khác 0

modelBMA$probne0

##  [1] 100.0 100.0 100.0 100.0   3.8   0.0   6.3   0.0   5.0   0.0   0.0   4.4   4.6   0.0   0.0   0.0   0.0   0.0   2.7
## [20]   5.0   0.0   0.0 100.0 100.0   3.0  12.0   0.0   4.8   0.0   0.0

9.5 Xác suất hậu nghiệm của các mô hình

modelBMA$postprob

##  [1] 0.48546840 0.11999515 0.06318981 0.04968775 0.04968775 0.04773667 0.04586228 0.04406156 0.03753909 0.02952163
## [11] 0.02724989

9.6 Xác định mô hình có posterior probability cao nhất

best_BMA <- which.max(
  modelBMA$postprob
)

best_BMA

## [1] 1

9.6.1 Posterior probability của mô hình tốt nhất

modelBMA$postprob[
  best_BMA
]

## [1] 0.4854684

9.6.2 Các biến thuộc mô hình tốt nhất

modelBMA$which[
  best_BMA,
]

## SexMale     X01     X02     X03     X04     X09     X11     X12     X13     X16     X18     X20     X21     X22     X23 
##    TRUE    TRUE    TRUE    TRUE   FALSE   FALSE   FALSE   FALSE   FALSE   FALSE   FALSE   FALSE   FALSE   FALSE   FALSE 
##     X25     X28     X29     X31     X32     X34     X35  C01Yes  C02Yes  C03Yes  C05Yes  C08Yes  C09Yes  C12Yes  C15Yes 
##   FALSE   FALSE   FALSE   FALSE   FALSE   FALSE   FALSE    TRUE    TRUE   FALSE   FALSE   FALSE   FALSE   FALSE   FALSE

10 LASSO Regression

10.1 Chuẩn bị dữ liệu

Loại biến:

y
ID

X_LASSO <- df[
  ,
  setdiff(
    names(df),
    c(
      "y",
      "ID"
    )
  )
]

10.2 Chuyển X thành matrix

glmnet yêu cầu dữ liệu X ở dạng matrix.

xvars <- as.matrix(
  X_LASSO
)

yvars <- df$y

dim(xvars)

## [1] 2000   51

length(yvars)

## [1] 2000

10.3 Cross-validation LASSO

set.seed(123)

modelLASSO <- cv.glmnet(
  x = xvars,
  y = yvars,
  alpha = 1,
  standardize = TRUE
)

10.4 Lambda tối ưu

10.4.1 Lambda minimum

modelLASSO$lambda.min

## [1] 0.13093

10.4.2 Lambda 1-SE

modelLASSO$lambda.1se

## [1] 0.5285668

10.5 Hệ số tại `lambda.min`

coef_min <- coef(
  modelLASSO,
  s = "lambda.min"
)

coef_min

## 52 x 1 sparse Matrix of class "dgCMatrix"
##               lambda.min
## (Intercept) 31.907732543
## Sex          .          
## X01          2.856654239
## X02         -2.057386592
## X03          1.859464568
## X04         -0.040351421
## X05         -0.092950418
## X06         -0.004008194
## X07          .          
## X08          0.020597834
## X09          .          
## X10          0.001990350
## X11          .          
## X12          .          
## X13         -0.109713310
## X14          .          
## X15          .          
## X16          .          
## X17          0.038193293
## X18          .          
## X19         -0.073484577
## X20         -0.075867038
## X21         -0.010084846
## X22         -0.106851383
## X23          .          
## X24          .          
## X25          .          
## X26          .          
## X27          .          
## X28         -0.039921229
## X29          .          
## X30          .          
## X31          .          
## X32         -0.018972496
## X33          .          
## X34          .          
## X35          .          
## C01          .          
## C02          .          
## C03          .          
## C04          .          
## C05          .          
## C06          .          
## C07          .          
## C08          .          
## C09          .          
## C10          .          
## C11          .          
## C12          .          
## C13          .          
## C14          .          
## C15          .

10.6 Hệ số tại `lambda.1se`

coef_1se <- coef(
  modelLASSO,
  s = "lambda.1se"
)

coef_1se

## 52 x 1 sparse Matrix of class "dgCMatrix"
##             lambda.1se
## (Intercept)  31.939746
## Sex           .       
## X01           2.263692
## X02          -1.252123
## X03           1.234748
## X04           .       
## X05           .       
## X06           .       
## X07           .       
## X08           .       
## X09           .       
## X10           .       
## X11           .       
## X12           .       
## X13           .       
## X14           .       
## X15           .       
## X16           .       
## X17           .       
## X18           .       
## X19           .       
## X20           .       
## X21           .       
## X22           .       
## X23           .       
## X24           .       
## X25           .       
## X26           .       
## X27           .       
## X28           .       
## X29           .       
## X30           .       
## X31           .       
## X32           .       
## X33           .       
## X34           .       
## X35           .       
## C01           .       
## C02           .       
## C03           .       
## C04           .       
## C05           .       
## C06           .       
## C07           .       
## C08           .       
## C09           .       
## C10           .       
## C11           .       
## C12           .       
## C13           .       
## C14           .       
## C15           .

10.7 Các biến được chọn tại `lambda.min`

selected_min <- rownames(
  coef_min
)[
  as.vector(
    coef_min != 0
  )
]

selected_min <- setdiff(
  selected_min,
  "(Intercept)"
)

selected_min

##  [1] "X01" "X02" "X03" "X04" "X05" "X06" "X08" "X10" "X13" "X17" "X19" "X20" "X21" "X22" "X28" "X32"

10.8 Các biến được chọn tại `lambda.1se`

selected_1se <- rownames(
  coef_1se
)[
  as.vector(
    coef_1se != 0
  )
]

selected_1se <- setdiff(
  selected_1se,
  "(Intercept)"
)

selected_1se

## [1] "X01" "X02" "X03"

10.9 Số lượng biến được lựa chọn

cat(
  "Số biến tại lambda.min:",
  length(selected_min),
  "\n"
)

## Số biến tại lambda.min: 16

cat(
  "Số biến tại lambda.1se:",
  length(selected_1se),
  "\n"
)

## Số biến tại lambda.1se: 3

10.10 Biểu đồ Cross-validation LASSO

plot(
  modelLASSO,
  main = "Cross-validation của mô hình LASSO"
)

10.11 Đường đi hệ số LASSO

plot(
  modelLASSO$glmnet.fit,
  xvar = "lambda",
  label = FALSE,
  main = "Đường đi hệ số LASSO"
)

11 Tóm tắt kết quả LASSO

cat(
  "============================================\n"
)

## ============================================

cat(
  "KẾT QUẢ LASSO\n"
)

## KẾT QUẢ LASSO

cat(
  "============================================\n\n"
)

## ============================================

cat(
  "Lambda min:",
  modelLASSO$lambda.min,
  "\n"
)

## Lambda min: 0.13093

cat(
  "Lambda 1se:",
  modelLASSO$lambda.1se,
  "\n\n"
)

## Lambda 1se: 0.5285668

cat(
  "Biến được chọn tại lambda.min:\n"
)

## Biến được chọn tại lambda.min:

print(
  selected_min
)

##  [1] "X01" "X02" "X03" "X04" "X05" "X06" "X08" "X10" "X13" "X17" "X19" "X20" "X21" "X22" "X28" "X32"

cat(
  "\nBiến được chọn tại lambda.1se:\n"
)

## 
## Biến được chọn tại lambda.1se:

print(
  selected_1se
)

## [1] "X01" "X02" "X03"

cat(
  "\nSố biến tại lambda.min:",
  length(selected_min),
  "\n"
)

## 
## Số biến tại lambda.min: 16

cat(
  "Số biến tại lambda.1se:",
  length(selected_1se),
  "\n"
)

## Số biến tại lambda.1se: 3

12 Tổng kết

Trong bài thực hành này, các phương pháp được sử dụng gồm:

Hồi quy tuyến tính đơn biến.
Hồi quy tuyến tính đa biến.
Hiệu chỉnh yếu tố nhiễu.
Kiểm tra giả định mô hình.
Chia dữ liệu Development/Test bằng caret.
Đánh giá khả năng tiên lượng bằng RMSE, MAE và R-squared.
Bootstrap validation bằng rms.
Đánh giá đa cộng tuyến bằng VIF.
Stepwise Regression bằng olsrr.
Bayesian Model Averaging bằng BMA.
LASSO Regression bằng glmnet.

Kết thúc bài thực hành.

Phân tích mô hình hồi quy tuyến tính và lựa chọn mô hình tối ưu

Linear Regression - Stepwise - BMA - LASSO

Chính Dương

2026-08-17

1 Chuẩn bị môi trường

1.1 Load các package

2 Mô hình hồi quy tuyến tính và hiệu chỉnh yếu tố nhiễu

2.1 Đọc và khảo sát dữ liệu Salaries

2.2 Mô hình 1: Mối liên quan giữa giới tính và lương

2.2.1 Hệ số hồi quy

2.2.2 Khoảng tin cậy 95%

2.2.3 Giá trị dự đoán

2.2.4 Kiểm tra giả định của mô hình

2.3 Mô hình 2: Mối liên quan giữa thời gian sau tiến sĩ và lương

2.3.1 Hệ số hồi quy

2.3.2 Khoảng tin cậy 95%

2.3.3 Giá trị dự đoán

2.3.4 Biểu đồ mối liên quan

2.3.5 Kiểm tra giả định của mô hình

2.4 Mô hình 3: Giới tính và lương sau khi hiệu chỉnh thời gian sau tiến sĩ

2.4.1 Hệ số hồi quy

2.4.2 Khoảng tin cậy 95%

2.4.3 Giá trị dự đoán

2.4.4 Kiểm tra giả định

2.5 Mô hình 4: Rank, discipline và sex

2.5.1 Hệ số hồi quy

2.5.2 Khoảng tin cậy 95%

2.5.3 Đánh giá mức độ quan trọng của từng biến

2.5.4 Biểu đồ hệ số hồi quy

3 Đánh giá mô hình tiên lượng

3.1 Đọc dữ liệu mô phỏng béo phì

3.2 Xây dựng mô hình tiên lượng

3.2.1 Hệ số hồi quy

3.2.2 Khoảng tin cậy 95%

3.2.3 Giá trị tiên lượng

3.2.4 Kiểm tra giả định mô hình

4 Chia dữ liệu Development và Test bằng caret

4.0.1 Kích thước dữ liệu

4.1 Xây dựng mô hình trên tập Development

4.2 Dự đoán trên tập Test

4.2.1 So sánh giá trị thực tế và dự đoán

4.3 Đánh giá sai số dự đoán

4.3.1 RMSE

4.3.2 MAE

4.3.3 R-squared trên tập Test

4.3.4 Tổng hợp kết quả kiểm định

5 Bootstrap validation bằng package rms

5.1 Thiết lập datadist

5.2 Xây dựng mô hình OLS

5.3 Bootstrap 200 lần

6 Đánh giá đa cộng tuyến

6.1 Biểu đồ đa biến bằng GGally

6.2 Mô hình hồi quy đa biến

6.2.1 Hệ số hồi quy

6.2.2 Khoảng tin cậy 95%

6.3 Kiểm tra giả định mô hình đa biến

6.4 Đánh giá đa cộng tuyến bằng VIF

6.4.1 Trình bày VIF dạng bảng

6.4.2 Nguyên tắc diễn giải

7 Lựa chọn mô hình tối ưu

7.1 Đọc dữ liệu BMA

8 Stepwise Regression

8.1 Xây dựng mô hình đầy đủ

8.2 Stepwise selection bằng olsrr

9 Bayesian Model Averaging - BMA

9.1 Chuẩn bị biến kết quả Y

9.2 Chuẩn bị các biến dự đoán X

9.3 Chạy Bayesian Model Averaging

9.4 Xác suất biến khác 0

9.5 Xác suất hậu nghiệm của các mô hình

9.6 Xác định mô hình có posterior probability cao nhất

9.6.1 Posterior probability của mô hình tốt nhất

9.6.2 Các biến thuộc mô hình tốt nhất

10 LASSO Regression

10.1 Chuẩn bị dữ liệu

10.2 Chuyển X thành matrix

10.3 Cross-validation LASSO

10.4 Lambda tối ưu

10.4.1 Lambda minimum

10.4.2 Lambda 1-SE

2.1 Đọc và khảo sát dữ liệu `Salaries`

4 Chia dữ liệu Development và Test bằng `caret`

5 Bootstrap validation bằng package `rms`

5.1 Thiết lập `datadist`

6.1 Biểu đồ đa biến bằng `GGally`

8.2 Stepwise selection bằng `olsrr`

10.5 Hệ số tại `lambda.min`

10.6 Hệ số tại `lambda.1se`

10.7 Các biến được chọn tại `lambda.min`

10.8 Các biến được chọn tại `lambda.1se`