GAM_Return in education

Tài liệu buổi thực hành ngày 13/4/2022: Xử lý số liệu thô với R: Bộ số liệu Điều tra mức sống dân cư 2014

Speaker: PGS.TS. Ngô Hoàng Long, Trường ĐH Sư phạm Hà Nội TS. Trịnh Thị Hường, Trường ĐH Thương Mại

Chi tiết tại: https://sites.google.com/view/tkud/home?authuser=1

Tài liệu thực hành có thể download tại đây. (Chọn chuột phải tại chữ “tại đây”, chọn open new tab)

Hoặc copy link này: https://drive.google.com/drive/folders/1QoUDwzN49S0ohWOU3ztD__3uFXOfOqrA?usp=sharing

1. Tìm hiểu chủ đề nghiên cứu

Chúng tôi mô phỏng lại quá trình tiến hành xử lý số liệu thô của bài nghiên cứu “TÁC ĐỘNG CỦA GIÁO DỤC ĐẾN THU NHẬP TẠI VIỆT NAM GIAI ĐOẠN 2014-2020: KẾT QUẢ TỪ MÔ HÌNH HỒI QUY CỘNG TÍNH TỔNG QUÁT GAM” đã được chấp nhận đăng trên tạp Tạp chí kinh tế phát triển năm 2022.

Tóm tắt: Nghiên cứu này phân tích mối quan hệ giữa giáo dục và thu nhập cá nhân của người lao động tại Việt Nam trong các năm 2014, 2014, 2014 và 2020. Chúng tôi sử dụng dữ liệu cá nhân bao gồm thu nhập bình quân theo giờ, bằng cấp giáo dục cao nhất, số năm đào tạo và các thông tin nhân khẩu học từ các bộ số liệu Điều tra mức sống Dân cư của Tổng cục Thống kê Việt Nam. Kết quả từ mô hình hồi quy cộng tính tổng quát (A generalized additive model, GAM) thể hiện mối quan hệ phi tuyến và tích cực giữa số năm đi học và thu nhập theo giờ. Trong đó, lợi tức từ giáo dục của người lao động tăng 1 năm đào tạo ở trình độ cao là lớn hơn so với lợi tức từ tăng 1 năm đào tạo của các cá nhân ở trình độ thấp. Chúng tôi sử dụng biểu đồ xác suất q-q và tiêu chuẩn xác định chéo để kiểm chứng sự phù hợp của mô hình GAM so với mô hình hàm thu nhập Mincer. Kết quả nghiên cứu cho thấy vai trò quan trọng của việc đầu tư cho giáo dục, đặc biệt là đầu tư cho giáo dục ở trình độ cao.

Trong phần thực hành này, chúng tôi minh họa quá trình phân tích cho dữ liệu năm 2014.

2. Một số vấn đề về gói mgcv và hàm gam() trong R

Tìm hiểu hàm mgcv

require(mgcv)
#?gam
n <- 5000
eg <- gamSim(2,n=n,scale=.5)

## Bivariate smoothing example

attach(eg)

2.1. GAM in 2D

ind<-sample(1:n,200,replace=FALSE)

b5<-gam(y~s(x),data=data,
        knots=list(x=data$x[ind],z=data$z[ind]))
## various visualizations
plot(b5)

2.2 GAM in 3D

ind<-sample(1:n,200,replace=FALSE)
b5<-gam(y~s(x,z,k=40),data=data,
        knots=list(x=data$x[ind],z=data$z[ind]))
ind<-sample(1:n,200,replace=FALSE)
b5<-gam(y~s(x,z,k=40),data=data,
        knots=list(x=data$x[ind],z=data$z[ind]))
## various visualizations
vis.gam(b5,theta=30,phi=30)

plot(b5)

plot(b5,scheme=1,theta=50,phi=20)

plot(b5,scheme=2)

summary(b5)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## y ~ s(x, z, k = 40)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.300046   0.007093    42.3   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F p-value    
## s(x,z) 21.18  27.61 23.05  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.114   Deviance explained = 11.8%
## GCV = 0.25266  Scale est. = 0.25154   n = 5000

str(b5)

## List of 46
##  $ coefficients     : Named num [1:40] 0.3 0.0146 -0.0623 0.0664 -0.0424 ...
##   ..- attr(*, "names")= chr [1:40] "(Intercept)" "s(x,z).1" "s(x,z).2" "s(x,z).3" ...
##  $ residuals        : num [1:5000] -0.129 0.778 0.63 -0.298 -0.511 ...
##  $ fitted.values    : num [1:5000] 0.4381 0.439 0.2862 0.0372 0.2088 ...
##  $ family           :List of 11
##   ..$ family    : chr "gaussian"
##   ..$ link      : chr "identity"
##   ..$ linkfun   :function (mu)  
##   ..$ linkinv   :function (eta)  
##   ..$ variance  :function (mu)  
##   ..$ dev.resids:function (y, mu, wt)  
##   ..$ aic       :function (y, n, mu, wt, dev)  
##   ..$ mu.eta    :function (eta)  
##   ..$ initialize:  expression({  n <- rep.int(1, nobs)  if (family$link == "inverse")  mustart <- y + (y == 0) * sd(y) * 0.01  else | __truncated__
##   ..$ validmu   :function (mu)  
##   ..$ valideta  :function (eta)  
##   ..- attr(*, "class")= chr "family"
##  $ linear.predictors: num [1:5000] 0.4381 0.439 0.2862 0.0372 0.2088 ...
##  $ deviance         : num 1252
##  $ null.deviance    : num 1420
##  $ iter             : int 1
##  $ weights          : num [1:5000] 1 1 1 1 1 1 1 1 1 1 ...
##  $ prior.weights    : num [1:5000] 1 1 1 1 1 1 1 1 1 1 ...
##  $ df.null          : int 4999
##  $ y                : num [1:5000] 0.309 1.217 0.916 -0.261 -0.302 ...
##  $ converged        : logi TRUE
##  $ sig2             : num 0.252
##  $ edf              : Named num [1:40] 1 0.995 0.926 0.924 0.955 ...
##   ..- attr(*, "names")= chr [1:40] "(Intercept)" "s(x,z).1" "s(x,z).2" "s(x,z).3" ...
##  $ edf1             : Named num [1:40] 1 1.017 0.998 0.991 0.994 ...
##   ..- attr(*, "names")= chr [1:40] "(Intercept)" "s(x,z).1" "s(x,z).2" "s(x,z).3" ...
##  $ hat              : num [1:5000] 0.00395 0.00351 0.00373 0.00559 0.00482 ...
##  $ R                : num [1:40, 1:40] -70.7 0 0 0 0 ...
##  $ boundary         : logi FALSE
##  $ sp               : Named num 4.68
##   ..- attr(*, "names")= chr "s(x,z)"
##  $ nsdf             : int 1
##  $ Ve               : num [1:40, 1:40] 5.03e-05 -3.86e-20 -8.04e-20 -8.86e-20 3.81e-20 ...
##  $ Vp               : num [1:40, 1:40] 5.03e-05 -6.58e-20 -1.39e-19 -1.52e-19 7.23e-20 ...
##  $ rV               : num [1:40, 1:40] 4.88e-19 -3.11e-04 -4.52e-04 -9.10e-04 -5.54e-04 ...
##  $ mgcv.conv        :List of 7
##   ..$ full.rank      : int 40
##   ..$ rank           : int 40
##   ..$ fully.converged: logi FALSE
##   ..$ hess.pos.def   : logi TRUE
##   ..$ iter           : int 4
##   ..$ score.calls    : int 18
##   ..$ rms.grad       : num 5.17e-17
##  $ gcv.ubre         : Named num 0.253
##   ..- attr(*, "names")= chr "GCV.Cp"
##  $ aic              : num 7313
##  $ rank             : int 40
##  $ gcv.ubre.dev     : num 0.253
##  $ scale.estimated  : logi TRUE
##  $ method           : chr "GCV"
##  $ smooth           :List of 1
##   ..$ :List of 27
##   .. ..$ term          : chr [1:2] "x" "z"
##   .. ..$ bs.dim        : num 40
##   .. ..$ fixed         : logi FALSE
##   .. ..$ dim           : int 2
##   .. ..$ p.order       : num 0
##   .. ..$ by            : chr "NA"
##   .. ..$ label         : chr "s(x,z)"
##   .. ..$ xt            : NULL
##   .. ..$ id            : NULL
##   .. ..$ sp            : Named num -1
##   .. .. ..- attr(*, "names")= chr "s(x,z)"
##   .. ..$ drop.null     : num 0
##   .. ..$ S             :List of 1
##   .. .. ..$ : num [1:39, 1:39] 24.69 2.19 -2.65 2.32 1.03 ...
##   .. ..$ UZ            : num [1:203, 1:40] 0.89 -3.71 -6.3 3.97 3.49 ...
##   .. ..$ Xu            : num [1:200, 1:2] -0.499 -0.489 -0.487 -0.479 -0.475 ...
##   .. ..$ df            : num 39
##   .. ..$ shift         : num [1:2(1d)] 0.5 0.49
##   .. ..$ rank          : num 37
##   .. ..$ null.space.dim: num 2
##   .. ..$ plot.me       : logi TRUE
##   .. ..$ side.constrain: logi TRUE
##   .. ..$ repara        : logi TRUE
##   .. ..$ S.scale       : num 41.1
##   .. ..$ vn            : chr [1:2] "x" "z"
##   .. ..$ first.para    : num 2
##   .. ..$ last.para     : num 40
##   .. ..$ first.sp      : num 1
##   .. ..$ last.sp       : num 1
##   .. ..- attr(*, "class")= chr [1:2] "tprs.smooth" "mgcv.smooth"
##   .. ..- attr(*, "qrc")=List of 4
##   .. .. ..$ qr   : num [1:40, 1] 3.4428 -0.0858 -0.0289 -0.0109 -0.0107 ...
##   .. .. ..$ rank : int 1
##   .. .. ..$ qraux: num 1.01
##   .. .. ..$ pivot: int 1
##   .. .. ..- attr(*, "class")= chr "qr"
##   .. ..- attr(*, "nCons")= int 1
##  $ formula          :Class 'formula'  language y ~ s(x, z, k = 40)
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##  $ var.summary      :List of 2
##   ..$ x: num [1:3] 3.65e-05 5.02e-01 1.00
##   ..$ z: num [1:3] 4.85e-05 4.90e-01 9.99e-01
##  $ cmX              : Named num [1:40] 1.00 -3.52e-19 6.99e-18 5.02e-18 1.79e-18 ...
##   ..- attr(*, "names")= chr [1:40] "(Intercept)" "" "" "" ...
##  $ model            :'data.frame':   5000 obs. of  3 variables:
##   ..$ y: num [1:5000] 0.309 1.217 0.916 -0.261 -0.302 ...
##   ..$ x: num [1:5000] 0.7432 0.6189 0.5054 0.9429 0.0788 ...
##   ..$ z: num [1:5000] 0.794 0.686 0.235 0.399 0.722 ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language y ~ 1 + x + z
##   .. .. ..- attr(*, "variables")= language list(y, x, z)
##   .. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:3] "y" "x" "z"
##   .. .. .. .. ..$ : chr [1:2] "x" "z"
##   .. .. ..- attr(*, "term.labels")= chr [1:2] "x" "z"
##   .. .. ..- attr(*, "order")= int [1:2] 1 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(y, x, z)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "numeric" "numeric"
##   .. .. .. ..- attr(*, "names")= chr [1:3] "y" "x" "z"
##  $ control          :List of 19
##   ..$ nthreads    : num 1
##   ..$ irls.reg    : num 0
##   ..$ epsilon     : num 1e-07
##   ..$ maxit       : num 200
##   ..$ trace       : logi FALSE
##   ..$ mgcv.tol    : num 1e-07
##   ..$ mgcv.half   : num 15
##   ..$ rank.tol    : num 1.49e-08
##   ..$ nlm         :List of 6
##   .. ..$ ndigit           : num 7
##   .. ..$ gradtol          : num 1e-06
##   .. ..$ stepmax          : num 2
##   .. ..$ steptol          : num 1e-04
##   .. ..$ iterlim          : num 200
##   .. ..$ check.analyticals: logi FALSE
##   ..$ optim       :List of 1
##   .. ..$ factr: num 1e+07
##   ..$ newton      :List of 5
##   .. ..$ conv.tol: num 1e-06
##   .. ..$ maxNstep: num 5
##   .. ..$ maxSstep: num 2
##   .. ..$ maxHalf : num 30
##   .. ..$ use.svd : logi FALSE
##   ..$ outerPIsteps: num 0
##   ..$ idLinksBases: logi TRUE
##   ..$ scalePenalty: logi TRUE
##   ..$ efs.lspmax  : num 15
##   ..$ efs.tol     : num 0.1
##   ..$ keepData    : logi FALSE
##   ..$ scale.est   : chr "fletcher"
##   ..$ edge.correct: logi FALSE
##  $ terms            :Classes 'terms', 'formula'  language y ~ 1 + x + z
##   .. ..- attr(*, "variables")= language list(y, x, z)
##   .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:3] "y" "x" "z"
##   .. .. .. ..$ : chr [1:2] "x" "z"
##   .. ..- attr(*, "term.labels")= chr [1:2] "x" "z"
##   .. ..- attr(*, "order")= int [1:2] 1 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(y, x, z)
##   .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "numeric" "numeric"
##   .. .. ..- attr(*, "names")= chr [1:3] "y" "x" "z"
##  $ pred.formula     :Class 'formula'  language ~x + z
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "full")=Class 'formula'  language ~y + x + z
##   .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##  $ pterms           :Classes 'terms', 'formula'  language y ~ 1
##   .. ..- attr(*, "variables")= language list(y)
##   .. ..- attr(*, "factors")= int(0) 
##   .. ..- attr(*, "term.labels")= chr(0) 
##   .. ..- attr(*, "order")= int(0) 
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(y)
##   .. ..- attr(*, "dataClasses")= Named chr "numeric"
##   .. .. ..- attr(*, "names")= chr "y"
##  $ assign           : int 0
##  $ offset           : num [1:5000] 0 0 0 0 0 0 0 0 0 0 ...
##  $ df.residual      : num 4978
##  $ min.edf          : num 3
##  $ optimizer        : chr "magic"
##  $ call             : language gam(formula = y ~ s(x, z, k = 40), data = data, knots = list(x = data$x[ind],      z = data$z[ind]))
##  - attr(*, "class")= chr [1:3] "gam" "glm" "lm"

Một số minh họa về cơ sở spline tại:[tại đây] (https://rpubs.com/Lucie/806781)

3. Quá trình phân tích số liệu

3.1. Đọc số liệu sạch

rm(list = ls()) #remove all environment

library(tidyverse)
library(nortest)
require(goft)
library(fitdistrplus)
library(tableone)
library(flexsurv) # on CRAN
require(here)
# Load data
#setwd("D:/TaphuanVIASM2022/SeminarTKUD/HuongTrinh/13April")
load(here("EDUC14V3.RData"))
EDUC14V3 <- na.omit(EDUC14V3)

dim(EDUC14V3)

## [1] 7155   24

head(EDUC14V3)

##              ID           IDind  hsalary Ownership.type PROVINCE MEMBER.CODE
## 1   01_1_4_8_15   01_1_4_8_15_2 38.35227          state       01           2
## 2   01_1_7_6_15   01_1_7_6_15_3 44.43182          state       01           3
## 3   01_1_7_6_15   01_1_7_6_15_4 20.26515        private       01           4
## 4 01_1_16_20_14 01_1_16_20_14_2 22.06439        private       01           2
## 5 01_1_22_19_13 01_1_22_19_13_3 38.35227        foreign       01           3
## 6 01_1_22_19_13 01_1_22_19_13_2 33.61742        foreign       01           2
##   Gender Marital Age hsize hsize15_59 hsizeNo DeRatio URBAN Ethnicity
## 1 Female Married  33     5          2       3   150.0     1      Kinh
## 2   Male Married  30     5          2       3   150.0     1      Kinh
## 3 Female Married  30     5          2       3   150.0     1      Kinh
## 4   Male Married  43     5          3       2    66.7     1      Kinh
## 5 Female Married  32     6          3       3   100.0     1      Kinh
## 6   Male Married  39     6          3       3   100.0     1      Kinh
##        Highest.EDU               Highest.VACATION YEDU            AREA
## 1 No qualification                             No   16 Red River Delta
## 2 No qualification                             No   16 Red River Delta
## 3   Primary school Professional vocational school   14 Red River Delta
## 4   Primary school                             No   12 Red River Delta
## 5 No qualification                             No   16 Red River Delta
## 6   Primary school Middle-level vocational school   14 Red River Delta
##   Income_avg Incomes tentinh       sp GDPP
## 1       9876  592600  HA NOI 35.80185 4113
## 2       5769  346160  HA NOI 35.80185 4113
## 3       5769  346160  HA NOI 35.80185 4113
## 4       2676  160600  HA NOI 35.80185 4113
## 5       4666  336000  HA NOI 35.80185 4113
## 6       4666  336000  HA NOI 35.80185 4113

3.2. Thống kê mô tả

hsalary: Trung bình thu nhập theo giờ của người lao động (nghìn đồng/giờ), được thu thập từ tiền công và phúc lợi của người lao động trong vòng 12 tháng trước thời điểm khảo sát.

Biểu đồ tần suất của hsalary và logarit hsalary

hist(log(EDUC14V3$hsalary))

hist(x = log(EDUC14V3$hsalary), freq = FALSE, xlim = c(0, 6))
lines(x = density(x = log(EDUC14V3$hsalary)), col = "red")

Kiểm tra phân phối của hsalary tuân theo phân phối chuẩn

ks.test(EDUC14V3$hsalary, "pnorm")

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  EDUC14V3$hsalary
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided

# Kolmogorov-Smirnov test: The null hypothesis states that there is no difference between the two distributions

#H0 : Tổng thể tuân theo phân phối lý thuyết;
#H1 : Tổng thể KHÔNG tuân theo phân phối lý thuyết.
# p-value < 0.05 => bác bỏ H0, chấp nhận H1.

qqnorm(log(EDUC14V3$hsalary))
qqline(log(EDUC14V3$hsalary))

#-----Fiting gamma distribution without log
gamma_test(EDUC14V3$hsalary)

## 
##  Test of fit for the Gamma distribution
## 
## data:  EDUC14V3$hsalary
## V = -4.1228, p-value = 0.003554

#p-value = 0.004 > 0.01
# Chấp nhận Ho => EDUC14V3$hsalary tuân theo phân phối gamma
gammafit14 <- fitdistrplus::fitdist(EDUC14V3$hsalary, "gamma")

qqcomp(gammafit14, main = "fitting hsalary to gamma dist", legendtext = c("2014"), fitcol = "blue")

#============== THONG KE MO TA ==============

catVars <- c("Gender","URBAN","Marital","AREA","Ethnicity", "Ownership.type")
myVars  <- c("Age", "hsize", "YEDU", "DeRatio", "hsalary", "Gender", "URBAN", 
             "Marital", "AREA", "Ethnicity", "Ownership.type")

tab14 <- CreateTableOne(vars = myVars, data = EDUC14V3 ,
                        factorVars = catVars)
print(tab14)

##                                     
##                                      Overall      
##   n                                   7155        
##   Age (mean (SD))                    35.42 (10.95)
##   hsize (mean (SD))                   4.36 (1.51) 
##   YEDU (mean (SD))                    9.54 (5.06) 
##   DeRatio (mean (SD))                56.24 (58.59)
##   hsalary (mean (SD))                21.50 (11.11)
##   Gender = Female (%)                 2936 (41.0) 
##   URBAN = 0 (%)                       4390 (61.4) 
##   Marital = Married (%)               5208 (72.8) 
##   AREA (%)                                        
##      Central Highlands                 326 ( 4.6) 
##      Mekong River Delta               1416 (19.8) 
##      Northern-Coastal Central         1626 (22.7) 
##      Northern Midlands and Mountains   855 (11.9) 
##      Red River Delta                  1713 (23.9) 
##      Southeastern Area                1219 (17.0) 
##   Ethnicity = Minor (%)                706 ( 9.9) 
##   Ownership.type (%)                              
##      private                          4567 (63.8) 
##      state                            1934 (27.0) 
##      foreign                           654 ( 9.1)

#====== Bieu do hslary nam 2014 ======

ggplot(EDUC14V3, aes(x = Highest.EDU, y = hsalary)) +
  geom_boxplot(outlier.shape = NA)+
  ggtitle("hsalary by education levels") +
  ylim(c(5,70))

ggplot(EDUC14V3[EDUC14V3$Highest.VACATION != "No", ], aes(x = Highest.VACATION, y = hsalary)) +
  geom_boxplot(outlier.shape = NA)+
  ggtitle("hsalary by Vocational education levels") +
  ylim(c(5,70))

3.3 Mô hình hồi quy

Mô hình hồi quy bội với log(hsalary)

#-----------Regression--------------

reg1.14 <- lm(log(hsalary) ~ YEDU + Marital + DeRatio + Ethnicity + Gender + Age + 
                Ownership.type + URBAN +  AREA, data = EDUC14V3)

summary(reg1.14)

## 
## Call:
## lm(formula = log(hsalary) ~ YEDU + Marital + DeRatio + Ethnicity + 
##     Gender + Age + Ownership.type + URBAN + AREA, data = EDUC14V3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.78616 -0.27784  0.04803  0.32215  1.35734 
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                          2.365e+00  3.740e-02  63.236  < 2e-16 ***
## YEDU                                 3.788e-02  1.401e-03  27.035  < 2e-16 ***
## MaritalMarried                       1.535e-01  1.420e-02  10.807  < 2e-16 ***
## DeRatio                              6.195e-04  9.605e-05   6.449 1.20e-10 ***
## EthnicityMinor                      -1.028e-01  2.023e-02  -5.082 3.83e-07 ***
## GenderFemale                        -1.101e-01  1.153e-02  -9.548  < 2e-16 ***
## Age                                  3.064e-03  5.858e-04   5.230 1.75e-07 ***
## Ownership.typestate                  1.343e-01  1.533e-02   8.761  < 2e-16 ***
## Ownership.typeforeign                2.055e-01  2.045e-02  10.052  < 2e-16 ***
## URBAN0                              -1.826e-01  1.196e-02 -15.265  < 2e-16 ***
## AREAMekong River Delta               1.372e-02  2.886e-02   0.475   0.6346    
## AREANorthern-Coastal Central        -3.603e-02  2.848e-02  -1.265   0.2058    
## AREANorthern Midlands and Mountains -1.798e-04  3.059e-02  -0.006   0.9953    
## AREARed River Delta                  5.764e-02  2.867e-02   2.010   0.0444 *  
## AREASoutheastern Area                2.791e-01  2.945e-02   9.476  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4641 on 7140 degrees of freedom
## Multiple R-squared:  0.3067, Adjusted R-squared:  0.3053 
## F-statistic: 225.6 on 14 and 7140 DF,  p-value: < 2.2e-16

Ước lượng mô hình GAM, phân phối của hsalry là gamma và hàm liên kết là log

#YEDU: số năm đào tạo/số năm đi học

reg2.14 <- gam(hsalary ~ s(YEDU, k = 10,bs = "cr") + 
                 #k: the dimension of the basis used to represent the smooth term
                 # bs =  "cr" for cubic regression spline
                 Marital + DeRatio + Ethnicity + Gender + Age +
                 Ownership.type + URBAN +  AREA, data = EDUC14V3 ,
               family = Gamma(link = "log"),
               method = "REML")
summary(reg2.14)

## 
## Family: Gamma 
## Link function: log 
## 
## Formula:
## hsalary ~ s(YEDU, k = 10, bs = "cr") + Marital + DeRatio + Ethnicity + 
##     Gender + Age + Ownership.type + URBAN + AREA
## 
## Parametric coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                          2.822e+00  3.140e-02  89.862  < 2e-16 ***
## MaritalMarried                       1.197e-01  1.320e-02   9.073  < 2e-16 ***
## DeRatio                              5.009e-04  8.949e-05   5.597 2.26e-08 ***
## EthnicityMinor                      -9.093e-02  1.885e-02  -4.824 1.43e-06 ***
## GenderFemale                        -1.133e-01  1.076e-02 -10.531  < 2e-16 ***
## Age                                  3.853e-03  5.467e-04   7.047 1.99e-12 ***
## Ownership.typestate                  1.014e-01  1.488e-02   6.815 1.02e-11 ***
## Ownership.typeforeign                2.083e-01  1.903e-02  10.950  < 2e-16 ***
## URBAN0                              -1.538e-01  1.117e-02 -13.777  < 2e-16 ***
## AREAMekong River Delta               2.641e-03  2.684e-02   0.098   0.9216    
## AREANorthern-Coastal Central        -4.157e-02  2.646e-02  -1.571   0.1162    
## AREANorthern Midlands and Mountains  1.121e-02  2.845e-02   0.394   0.6935    
## AREARed River Delta                  6.603e-02  2.667e-02   2.476   0.0133 *  
## AREASoutheastern Area                2.578e-01  2.738e-02   9.416  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##           edf Ref.df     F p-value    
## s(YEDU) 8.004  8.644 106.1  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.348   Deviance explained = 31.7%
## -REML =  25433  Scale est. = 0.1858    n = 7155

plot(reg2.14)

3.4 So sánh độ phù hợp của mô hình qua tiêu chuẩn xác định chéo (Cross validation)

#================ Ve hinh fitted===================

#================Chay Cross validation=========================
f.CV.mean.EDUC <- function(year)
{
  # require(mgcv)
  don <- EDUC14V3
  validation <- don[sample(1:nrow(don), 1000,replace = FALSE),] 
  train <- don[don$ID%in%validation$ID == FALSE,] #Set the training set
  
  reg.lm.log <- lm(log(hsalary) ~ YEDU + Marital + DeRatio + Ethnicity + Gender +
                     Age + Ownership.type + URBAN + AREA, data = train) 
  
  newpred.lm.log <- predict(reg.lm.log, type = "response", newdata = validation) 
  
  RSS.lm.log <- mean((validation$hsalary - (exp(newpred.lm.log + (summary(reg.lm.log)$sigma)^2/2))) ^2)
  ##
  
  reg.Gamma <- mgcv::gam(hsalary ~ s(YEDU) + Marital + DeRatio + Ethnicity + Gender + Age +
                           Ownership.type + URBAN +  AREA,
                         family = Gamma(link = "log"), data = train) 
  newpred.Gamma <- predict(reg.Gamma,type = "response", newdata = validation) 
  RSS.Gamma.log <- mean((validation$hsalary - newpred.Gamma)^2)
  
  return(c(RSS.lm.log, RSS.Gamma.log))
}

#--------------------------------------------------------------

# #-----YEAR 2014---
# 
CV14 <- data.frame(LM = NA, GAM = NA)

for (i in 1:10) # Nghiên cứu thực tế là 10.000 lượt!!
{
  CV14[i, ] <- f.CV.mean.EDUC(2014)
}
# write.csv(CV14, file = "CV14.CSV")

boxplot(CV14[, c("LM", "GAM")],
        main = "CV14")

So sánh hai trung bình

t.test(CV14$LM, CV14$GAM)

## 
##  Welch Two Sample t-test
## 
## data:  CV14$LM and CV14$GAM
## t = 0.93029, df = 18, p-value = 0.3645
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.923692  4.981180
## sample estimates:
## mean of x mean of y 
##  83.36914  81.84039

3.5 Vẽ đồ thị minh họa so sánh kết quả ước lượng LM và GAM

Vẽ hai đường hồi quy và có dịch chuyển về tỉ lệ thực của hsalary.

#-------------------------------------------------------------------------
#------------------ UOC LUONG GAM ---------------------------------------
f.smooth.GAMLog.link <- function( model.gam)
{
  DON.new <- NULL
  
  n <- length(EDUC14V3[, "ID"])
  DON <- EDUC14V3
  DON.new <- data.frame(YEDU = seq(min(DON$YEDU), max(DON$YEDU),length.out = n),
                        Marital = factor(rep("Married", n)),
                        DeRatio = rep(60, n),
                        Ethnicity = factor(rep("Kinh", n)),
                        Gender = factor(rep("Male", n)),
                        Age = rep(36, n),
                        Ownership.type = factor(rep("private", n)),
                        URBAN = factor(rep("1", n)),
                        AREA = factor(rep("Red River Delta", n)))
  
  pre.GAMLog <- predict(model.gam, type = "link", newdata = DON.new , se.fit = TRUE)
  seq.YEDU <- DON.new$YEDU
  e.fit.GAMLog <- exp(pre.GAMLog$fit)
  fit.up95.GAMLog <- exp(pre.GAMLog$fit - 1.96*pre.GAMLog$se.fit)
  fit.low95.GAMLog <- exp(pre.GAMLog$fit + 1.96*pre.GAMLog$se.fit)
  return(c(seq.YEDU,e.fit.GAMLog, fit.up95.GAMLog, fit.low95.GAMLog))
}


income.2014.GAMGAMLog <- matrix(unlist(f.smooth.GAMLog.link( reg2.14)), 
                                ncol = 4,byrow = FALSE)

# plot "UOC LUONG GAM"
plot(income.2014.GAMGAMLog[,1], income.2014.GAMGAMLog[,2], type = "n", lwd = 3,
     main = "UOC LUONG GAM",xlab="YEDU", 
     ylab = "Hsalary", ylim = c(10,80), xlim = c(0,25))

# 2014
polygon(c(income.2014.GAMGAMLog[,1], rev(income.2014.GAMGAMLog[,1])), 
        c(income.2014.GAMGAMLog[,3], rev(income.2014.GAMGAMLog[,4])), col = "grey70",
        border = NA)

lines(income.2014.GAMGAMLog[,1],income.2014.GAMGAMLog[,2], lwd = 2,col = "green")

legend("topleft", legend = c("2014"), col = c("green"), lty = 1, cex = 1, lwd = 2)

# ==================== UOC LUONG LM =====================================
# =======================================================================
##----- Function -----
f.smooth.LM.link <- function(model.LM)
{ 
  DON.new <- NULL
  n <- length(EDUC14V3[,"ID"])
  DON <- EDUC14V3
  DON.new <- data.frame(YEDU = seq(min(DON$YEDU), max(DON$YEDU), length.out=n),
                        Marital = factor(rep("Married", n)),
                        DeRatio = rep(60, n),
                        Ethnicity = factor(rep("Kinh", n)),
                        Gender = factor(rep("Male", n)),
                        Age = rep(36, n),
                        Ownership.type = factor(rep("private", n)),
                        URBAN = factor(rep("1", n)),
                        AREA = factor(rep("Red River Delta", n)))
  
  pre.LM <- predict(model.LM, type = "response", newdata = DON.new, se.fit = TRUE)
  seq.YEDU <- DON.new$YEDU
  e.fit.LM <- exp(pre.LM$fit + (summary(model.LM)$sigma)^2/2)
  fit.up95.LM <- exp(pre.LM$fit + (summary(model.LM)$sigma)^2/2 - 1.96*pre.LM$se.fit)
  fit.low95.LM <- exp(pre.LM$fit + (summary(model.LM)$sigma)^2/2 + 1.96*pre.LM$se.fit)
  return(c(seq.YEDU, e.fit.LM, fit.up95.LM, fit.low95.LM))
}


income.2014.LM <- matrix(unlist(f.smooth.LM.link( reg1.14)), 
                         ncol = 4,byrow = FALSE)


#
plot(income.2014.LM[,1], income.2014.LM[,2], type = "n", lwd = 3,
     main = "UOC LUONG LM",xlab = "YEDU", 
     ylab = "Hsalary", ylim = c(10,80), xlim = c(0,25))
#2014
polygon(c(income.2014.LM[,1], rev(income.2014.LM[,1])), 
        c(income.2014.LM[,3],rev(income.2014.LM[,4])), col = "grey70",
        border = NA)
lines(income.2014.LM[,1], income.2014.LM[,2], lwd = 2, col = "green")
legend("topleft", legend = c("2014"), col = c("green"),lty = 1,cex = 1, lwd = 2)

Xuất kết quả ước lượng mô hình hồi quy GAM cho bài nghiên cứu

#=======================================================================
#=================Bang he so hoi quy===============================
Coef.OLS = data.frame(name = c(names(coef(reg1.14 )), "Rsquared"),
                      coef = NA,sd = NA,  row.names = "name")

gen.OLS <- function(model.OLS)
{
  p.v.OLS = 2*pnorm(abs(summary(model.OLS)$coefficients[, 1]/summary(model.OLS)$coefficients[, 2]), lower.tail=FALSE)
  star.OLS = ifelse(p.v.OLS > 0.1, " ",
                    ifelse(p.v.OLS <= 0.1 & p.v.OLS > 0.05, ".",
                           ifelse(p.v.OLS <= 0.05 & p.v.OLS > 0.01, "*",
                                  ifelse(p.v.OLS <= 0.01 & p.v.OLS > 0.001, "**", "***"))))
  Coef.OLS.tempt = Coef.OLS
  
  Coef.OLS.tempt$coef <- c(paste(round(coefficients(model.OLS), 3), star.OLS),
                           round(summary(reg1.14) $ adj.r.squared, 3))
  
  Coef.OLS.tempt$sd <- c(paste0("(", c(round(summary(model.OLS)$coefficients[, 2], 3)), ")"), "")
  
  return(Coef.OLS.tempt)
}

Coef.OLS.All <- data.frame(Coef14 = gen.OLS(reg1.14))
Coef.OLS.All

##                                     Coef14.coef Coef14.sd
## (Intercept)                           2.365 ***   (0.037)
## YEDU                                  0.038 ***   (0.001)
## MaritalMarried                        0.153 ***   (0.014)
## DeRatio                               0.001 ***       (0)
## EthnicityMinor                       -0.103 ***    (0.02)
## GenderFemale                          -0.11 ***   (0.012)
## Age                                   0.003 ***   (0.001)
## Ownership.typestate                   0.134 ***   (0.015)
## Ownership.typeforeign                 0.206 ***    (0.02)
## URBAN0                               -0.183 ***   (0.012)
## AREAMekong River Delta                  0.014     (0.029)
## AREANorthern-Coastal Central           -0.036     (0.028)
## AREANorthern Midlands and Mountains         0     (0.031)
## AREARed River Delta                     0.058 *   (0.029)
## AREASoutheastern Area                 0.279 ***   (0.029)
## Rsquared                                  0.305

Coef.GAM= data.frame(name = c(names(coef(reg2.14))[1:14], "Rsquared"),
                     coef = NA, sd = NA, row.names="name")

gen.GAM <- function(model.GAM)
{
  p.v.OLS = 2*pnorm(abs(summary(model.GAM)$coefficients[, 1]/summary(model.GAM)$coefficients[, 2]),lower.tail = FALSE)
  star.OLS = ifelse(p.v.OLS > 0.1, " ",
                    ifelse(p.v.OLS <= 0.1 & p.v.OLS > 0.05, ".",
                           ifelse(p.v.OLS <= 0.05 & p.v.OLS > 0.01, "*",
                                  ifelse(p.v.OLS <= 0.01 & p.v.OLS > 0.001, "**", "***"))))
  Coef.GAM.tempt = Coef.GAM
  Coef.GAM.tempt$coef <- c(paste0(round(summary(model.GAM)$p.coeff[1:14], 3), star.OLS),
                            round(summary(reg2.14)$ r.sq, 3))
  Coef.GAM.tempt$sd <- c(paste0("(", round(summary(model.GAM)$se[1:14], 3),")"), " ")
  
  return( Coef.GAM.tempt)
}

Coef.GAM.All <- data.frame(Coef14 = gen.GAM(reg2.14))
Coef.GAM.All

##                                     Coef14.coef Coef14.sd
## (Intercept)                               2.822   (0.031)
## MaritalMarried                             0.12   (0.013)
## DeRatio                                   0.001       (0)
## EthnicityMinor                           -0.091   (0.019)
## GenderFemale                             -0.113   (0.011)
## Age                                       0.004   (0.001)
## Ownership.typestate                       0.101   (0.015)
## Ownership.typeforeign                     0.208   (0.019)
## URBAN0                                   -0.154   (0.011)
## AREAMekong River Delta                    0.003   (0.027)
## AREANorthern-Coastal Central             -0.042   (0.026)
## AREANorthern Midlands and Mountains       0.011   (0.028)
## AREARed River Delta                       0.066   (0.027)
## AREASoutheastern Area                     0.258   (0.027)
## Rsquared                                  0.348

#print(Coef.GAM.All)

4. Tìm hiểu thêm về cơ sở đã sử dụng để ước lượng hàm s(YEDU)

coef(reg2.14)

##                         (Intercept)                      MaritalMarried 
##                        2.8219118847                        0.1197270995 
##                             DeRatio                      EthnicityMinor 
##                        0.0005008955                       -0.0909271974 
##                        GenderFemale                                 Age 
##                       -0.1133153210                        0.0038530024 
##                 Ownership.typestate               Ownership.typeforeign 
##                        0.1013906495                        0.2083414397 
##                              URBAN0              AREAMekong River Delta 
##                       -0.1538349876                        0.0026412792 
##        AREANorthern-Coastal Central AREANorthern Midlands and Mountains 
##                       -0.0415683655                        0.0112113573 
##                 AREARed River Delta               AREASoutheastern Area 
##                        0.0660276766                        0.2577829853 
##                           s(YEDU).1                           s(YEDU).2 
##                       -0.3882141449                       -0.0090997686 
##                           s(YEDU).3                           s(YEDU).4 
##                       -0.0770564170                        0.1866618578 
##                           s(YEDU).5                           s(YEDU).6 
##                        0.0734711911                        0.3049757697 
##                           s(YEDU).7                           s(YEDU).8 
##                        0.4381662997                        0.4308176101 
##                           s(YEDU).9 
##                        0.5462551373

TRÂN TRỌNG MỜI ĐẠI BIỂU THAM DỰ VÀ TRÂN TRỌNG CẢM ƠN!