Bài tập 1

Lấy data

library(tidyverse)
## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr
## Conflicts with tidy packages ----------------------------------------------
## filter(): dplyr, stats
## lag():    dplyr, stats
library(hexView)
setwd("D:/Duy Anh/DSR/data_giaotrinh")
duyanh <- readEViews("ch2bt3.WF1", as.data.frame = T)
## Warning in readEViews("ch2bt3.WF1", as.data.frame = T): Skipping boilerplate variable

## Warning in readEViews("ch2bt3.WF1", as.data.frame = T): Skipping boilerplate variable
names(duyanh)
## [1] "EXP01" "GRADE" "UNION" "WAGE"

Bài 2.5

Câu 1

ols1 <- lm(WAGE ~ GRADE + EXP01, data = duyanh)
summary(ols1)
## 
## Call:
## lm(formula = WAGE ~ GRADE + EXP01, data = duyanh)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.291  -4.568  -1.152   2.137  45.765 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.85554    4.14526  -1.654    0.101    
## GRADE        1.32896    0.29352   4.528 1.70e-05 ***
## EXP01        0.30566    0.06378   4.792 5.93e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.79 on 97 degrees of freedom
## Multiple R-squared:  0.2565, Adjusted R-squared:  0.2411 
## F-statistic: 16.73 on 2 and 97 DF,  p-value: 5.732e-07

Nếu gia tăng một năm học mà số năm kinh nghiệm giữa nguyên thì mức lương trung bình tăng thêm 1.32 đơn vị ### Câu 2

mean(duyanh$EXP01)
## [1] 19.48
dubao <- data.frame(GRADE = 12, EXP01 = mean(duyanh$EXP01))
predict(ols1, dubao)
##        1 
## 15.04619

Mức lương trung bình của một người có số năm học là 12, kinh nghiệm là 19.48 năm là 15.04619 đơn vị ### Câu 3 Nếu chỉ quan tâm đến lương, vì hệ số ước lượng của GRADE lớn hơn hệ số ước lượng của EXP, nên người đó nên tăng thời gian đi học và do đó giảm số năm kinh nghiệm

Bài 2.6

Câu 1

ols2 <- duyanh %>%lm(log(WAGE) ~ GRADE + EXP01, data = .) 
summary(ols2)
## 
## Call:
## lm(formula = log(WAGE) ~ GRADE + EXP01, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.82144 -0.24795 -0.04623  0.21754  1.19179 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.262617   0.193793   6.515 3.24e-09 ***
## GRADE       0.081364   0.013722   5.929 4.67e-08 ***
## EXP01       0.018342   0.002982   6.152 1.72e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3642 on 97 degrees of freedom
## Multiple R-squared:  0.3667, Adjusted R-squared:  0.3536 
## F-statistic: 28.08 on 2 and 97 DF,  p-value: 2.388e-10

Ý nghĩa hệ số ước lượng: Khi tăng thêm 1 năm học, các yếu tố khác không đổi, mức lương trung bình tăng thêm 8% Khi tăng thêm 1 năm kinh nghiệm, các yếu tố khác không đổi, mức lương trung bình tăng thêm 2%

Câu 2: Ko thể so sánh trực tiếp hệ số xác định của hai mô hình với nhau được vì không có cùng biến phụ thuộc

So sánh 2 mô hình bằng phương pháp tái chọn mẫu

Phương pháp K-fold Cross Validation

Chọn tập training set: gồm 90 phần tử testing set: gồm 10 phần tử Tính MSE cho 2 mô hình

library(Metrics)
MSEtesting1 <- c()
MSEtesting2 <- c()
n <- 1000
for (i in 1:n) {
  id <- sample(nrow(duyanh), 10)
  training <- duyanh[-id,]
  testing <- duyanh[id,]
  ols01 <- lm(WAGE ~ GRADE + EXP01, data = training)
  ols02 <- lm(log(WAGE) ~ GRADE + EXP01, data = training)
  MSEtesting1 <- append(MSEtesting1, mse(testing$WAGE, predict(ols01,testing)))
  MSEtesting2 <- append(MSEtesting2, mse(testing$WAGE, exp(predict(ols02, testing))))
  }
kq <- data.frame(MSEtesting1, MSEtesting2)
summary(kq)
##   MSEtesting1       MSEtesting2     
##  Min.   :  4.943   Min.   :  2.322  
##  1st Qu.: 21.631   1st Qu.: 16.544  
##  Median : 32.291   Median : 30.605  
##  Mean   : 63.536   Mean   : 62.444  
##  3rd Qu.: 58.679   3rd Qu.: 57.684  
##  Max.   :447.456   Max.   :447.930
dim(kq)
## [1] 1000    2
View(kq)

MSE của mô hình 2 thấp hơn nên mô hình 2 dự báo tốt hơn

Minh họa bằng hình ảnh

library(ggplot2)
error <- c(MSEtesting1,MSEtesting2)
cat <- c(rep("Model_1",500), rep("Model_2",500))
mydf <- data.frame(error,cat)
ggplot(mydf, aes(cat, error, fill= cat)) +
  geom_boxplot(show.legend = F)

ggplot(mydf, aes(error,  fill = cat)) + geom_density(alpha = 0.3)

Phương pháp LOOCV

mse_ols1 <- c()
mse_ols2 <- c()
for (i in 1:nrow(duyanh)){
  testing_LOOCV <- duyanh[i,]
  training_LOOCV <- duyanh[-i,]
  ols_LOOCV_1 <- lm(WAGE ~ GRADE + EXP01, data = training_LOOCV)
  ols_LOOCV_2 <- lm(log(WAGE) ~ GRADE + EXP01, data =training_LOOCV)
  mse_LOOCV_1 <- mse(testing_LOOCV$WAGE, predict(ols_LOOCV_1, testing_LOOCV))
  mse_LOOCV_2 <- mse(testing_LOOCV$WAGE, exp(predict(ols_LOOCV_2, testing_LOOCV)))
  mse_ols1 <- c(mse_ols1,mse_LOOCV_1)
  mse_ols2 <- c(mse_ols2,mse_LOOCV_2)
  kq2 <- data.frame(mse_ols1,mse_ols2)
}
summary(kq2)
##     mse_ols1            mse_ols2        
##  Min.   :   0.0003   Min.   :   0.0006  
##  1st Qu.:   5.0924   1st Qu.:   2.0226  
##  Median :  13.5150   Median :  11.7928  
##  Mean   :  63.6383   Mean   :  62.5821  
##  3rd Qu.:  34.5729   3rd Qu.:  26.3419  
##  Max.   :2306.4988   Max.   :2355.2451

Giá trị trung bình của MSE của mô hình 2 thấp hơn mô hình 1. Do đó mô hình 2 dự báo tốt hơn

sum(kq2$mse_ols1==kq2$mse_ols2)
## [1] 0

Ko có trường hợp hai mse của mô hình bằng nhau

sum(kq2$mse_ols1>kq2$mse_ols2)
## [1] 61

Có 61 trường hợp mse của mô hình 1 lớn hơn