Nhiệm vụ tuần 7

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## Warning: package 'ggplot2' was built under R version 4.3.3

## Warning: package 'psych' was built under R version 4.3.3

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

## Warning: package 'DescTools' was built under R version 4.3.3

## 
## Attaching package: 'DescTools'

## The following objects are masked from 'package:psych':
## 
##     AUC, ICC, SD

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.0
## ✔ readr     2.1.5     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ psych::%+%()    masks ggplot2::%+%()
## ✖ psych::alpha()  masks ggplot2::alpha()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## 
## Attaching package: 'scales'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## 
## The following object is masked from 'package:readr':
## 
##     col_factor
## 
## 
## The following objects are masked from 'package:psych':
## 
##     alpha, rescale

## Warning: package 'AER' was built under R version 4.3.3

## Loading required package: car

## Warning: package 'car' was built under R version 4.3.3

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.3.3

## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:purrr':
## 
##     some
## 
## The following object is masked from 'package:DescTools':
## 
##     Recode
## 
## The following object is masked from 'package:psych':
## 
##     logit
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## Loading required package: lmtest

## Warning: package 'lmtest' was built under R version 4.3.3

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.3.3

## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Loading required package: sandwich

## Warning: package 'sandwich' was built under R version 4.3.3

## Loading required package: survival
## 
## Attaching package: 'survival'
## 
## The following object is masked from 'package:epitools':
## 
##     ratetable

## Warning: package 'datarium' was built under R version 4.3.3

## Warning: package 'Ecdat' was built under R version 4.3.3

## Loading required package: Ecfun

## Warning: package 'Ecfun' was built under R version 4.3.3

## 
## Attaching package: 'Ecfun'
## 
## The following object is masked from 'package:DescTools':
## 
##     BoxCox
## 
## The following object is masked from 'package:base':
## 
##     sign
## 
## 
## Attaching package: 'Ecdat'
## 
## The following object is masked from 'package:carData':
## 
##     Mroz
## 
## The following object is masked from 'package:datasets':
## 
##     Orange

## Warning: package 'ISLR' was built under R version 4.3.3

d <- Wages
datatable(d)

1 Giới thiệu bộ dữ liệu Wages

Bộ dữ liệu Wages được lấy từ package Ecdat, chứa thông tin về mức lương, trình độ học vấn và các đặc điểm nhân khẩu học của người lao động tại Hoa Kỳ. Bộ dữ liệu Wages chứa 10 biến với tổng cộng 4165 quan sát.

Dữ liệu này thường được sử dụng trong phân tích kinh tế lao động và các mô hình kinh tế lượng để kiểm định tác động của giới tính, chủng tộc, học vấn và tình trạng hôn nhân đến thu nhập.

• Ý nghĩa các biến trong bộ dữ liệu

Dưới đây là giải thích chi tiết về các biến trong bộ dữ liệu:

Tên biến	Kiểu dữ liệu	Ý nghĩa
education	Numeric	Số năm học chính quy mà cá nhân đã hoàn thành.
experience	Numeric	Số năm kinh nghiệm làm việc.
sex	Factor	Giới tính của cá nhân: male hoặc female.
ethnic	Factor	Chủng tộc: white, black, hispanic, other.
married	Factor	Tình trạng hôn nhân: married hoặc single.
occupation	Factor	Nhóm nghề nghiệp: manager, sales, clerical, service, other.
industry	Factor	Ngành nghề: manufacturing, construction, trade, services, other.
union	Factor	Thành viên công đoàn: yes hoặc no.
wage	Numeric	Mức lương theo giờ (USD).

Loại biến và Số lượng

Loại biến	Số lượng	Danh sách biến
Định tính	6	sex, ethnic, married, occupation, industry, union
Định lượng	3	education, experience, wage

• Kết luận

Bộ dữ liệu Wages là nguồn dữ liệu tiêu biểu và dễ tiếp cận cho các bài tập kinh tế lượng cơ bản về phân tích mức lương và tác động của các yếu tố nhân khẩu học đến thu nhập. Dữ liệu bao gồm nhiều biến định tính có thể dùng cho phân tích thống kê mô tả, kiểm định độc lập, hoặc xây dựng mô hình hồi quy lương theo từng nhóm đặc điểm.

---

2 Phân tích mối quan hệ giữa việc tham gia công đoàn và các đặc điểm nhân khẩu học của người lao động tại Hoa Kỳ

2.1 Phân tích mối quan hệ giữa biến union và các biến định tính khác

2.1.1 Biến union và biến sex

Bảng tuần suất chéo

a <- table(d$union,d$sex)
addmargins(a)

##      
##       female male  Sum
##   no     370 2279 2649
##   yes     99 1417 1516
##   Sum    469 3696 4165

Bảng tuần suất chéo theo tỷ lệ hàng

a1 <-prop.table(a, margin=1)

Trực quan hóa

a2 <- as.data.frame(a1)
colnames(a2) <- c("union", "sex", "Proportion")

ggplot(a2, aes(x = sex, y = Proportion, fill = union)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(
    title = "Tỷ lệ tham gia công đoàn theo tình trạng giới tính",
    x = "Tình trạng giới tính",
    y = "Tỷ lệ",
    fill = "Tình trạng tham gia công đoàn"
  ) +
  theme_minimal() +
  scale_y_continuous(labels = scales::percent) +
  scale_x_discrete(labels = c("female" = "Nữ", "male" = "Nam")) +
  scale_fill_manual(values = c("no" = "pink", "yes" = "blue")) +
  geom_text(aes(label = scales::percent(Proportion, accuracy = 0.1)),
            position = position_dodge(width = 0.9),
            vjust = -0.25)

Mối quan hệ giữa hai biến

chisq.test(a)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  a
## X-squared = 52.63, df = 1, p-value = 4.027e-13

• Relative risk

TH1: Không tham gia công đoàn làm biến tham chiếu

riskratio(a)

## $data
##        
##         female male Total
##   no       370 2279  2649
##   yes       99 1417  1516
##   Total    469 3696  4165
## 
## $measure
##      risk ratio with 95% C.I.
##       estimate    lower    upper
##   no  1.000000       NA       NA
##   yes 1.086446 1.064604 1.108737
## 
## $p.value
##      two-sided
##         midp.exact fisher.exact   chi.square
##   no            NA           NA           NA
##   yes 3.819167e-14 4.377756e-14 2.760744e-13
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

TH2: Tham gia công đoàn làm biến tham chiếu

riskratio(a, rev="b")

## $data
##        
##         male female Total
##   yes   1417     99  1516
##   no    2279    370  2649
##   Total 3696    469  4165
## 
## $measure
##      risk ratio with 95% C.I.
##       estimate    lower    upper
##   yes 1.000000       NA       NA
##   no  2.138867 1.729222 2.645556
## 
## $p.value
##      two-sided
##         midp.exact fisher.exact   chi.square
##   yes           NA           NA           NA
##   no  3.819167e-14 4.377756e-14 2.760744e-13
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

• Oddratio

oddsratio(a)

## $data
##        
##         female male Total
##   no       370 2279  2649
##   yes       99 1417  1516
##   Total    469 3696  4165
## 
## $measure
##      odds ratio with 95% C.I.
##       estimate    lower    upper
##   no  1.000000       NA       NA
##   yes 2.320732 1.848487 2.939374
## 
## $p.value
##      two-sided
##         midp.exact fisher.exact   chi.square
##   no            NA           NA           NA
##   yes 3.819167e-14 4.377756e-14 2.760744e-13
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"

2.1.2 Biến union và biến industry

Bảng tuần suất chéo

b <- table(d$union,d$married)
addmargins(b)

##      
##         no  yes  Sum
##   no   581 2068 2649
##   yes  192 1324 1516
##   Sum  773 3392 4165

Bảng tuần suất chéo theo tỷ lệ hàng

b1 <-prop.table(b, margin=1)

Trực quan hóa

b2 <- as.data.frame(b1)
colnames(b2) <- c("union", "married", "Proportion")

ggplot(b2, aes(x = married, y = Proportion, fill = union)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(
    title = "Tỷ lệ tham gia công đoàn theo tình trạng hôn nhân",
    x = "Tình trạng hôn nhân",
    y = "Tỷ lệ",
    fill = "Tình trạng tham gia công đoàn"
  ) +
  theme_minimal() +
  scale_y_continuous(labels = scales::percent) +
  scale_x_discrete(labels = c("no" = "Chưa kết hôn", "yes" = "Đã kết hôn")) +
  scale_fill_manual(values = c("no" = "orange", "yes" = "seagreen")) +
  geom_text(aes(label = scales::percent(Proportion, accuracy = 0.1)),
            position = position_dodge(width = 0.9),
            vjust = -0.25)

Mối quan hệ giữa hai biến

chisq.test(b)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  b
## X-squared = 54.181, df = 1, p-value = 1.828e-13

• Relative risk

TH1: Không tham gia công đoàn làm biến tham chiếu

riskratio(b)

## $data
##        
##          no  yes Total
##   no    581 2068  2649
##   yes   192 1324  1516
##   Total 773 3392  4165
## 
## $measure
##      risk ratio with 95% C.I.
##       estimate    lower    upper
##   no  1.000000       NA       NA
##   yes 1.118717 1.088005 1.150296
## 
## $p.value
##      two-sided
##         midp.exact fisher.exact   chi.square
##   no            NA           NA           NA
##   yes 3.863576e-14 5.029762e-14 1.339327e-13
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

TH2: Tham gia công đoàn làm biến tham chiếu

riskratio(b, rev="b")

## $data
##        
##          yes  no Total
##   yes   1324 192  1516
##   no    2068 581  2649
##   Total 3392 773  4165
## 
## $measure
##      risk ratio with 95% C.I.
##       estimate    lower    upper
##   yes 1.000000       NA       NA
##   no  1.731778 1.489884 2.012945
## 
## $p.value
##      two-sided
##         midp.exact fisher.exact   chi.square
##   yes           NA           NA           NA
##   no  3.863576e-14 5.029762e-14 1.339327e-13
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

• Oddratio

oddsratio(b)

## $data
##        
##          no  yes Total
##   no    581 2068  2649
##   yes   192 1324  1516
##   Total 773 3392  4165
## 
## $measure
##      odds ratio with 95% C.I.
##       estimate    lower    upper
##   no  1.000000       NA       NA
##   yes 1.936146 1.624601 2.316163
## 
## $p.value
##      two-sided
##         midp.exact fisher.exact   chi.square
##   no            NA           NA           NA
##   yes 3.863576e-14 5.029762e-14 1.339327e-13
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"

2.2 Mô hình hồi quy logistic

2.2.1 Biến union và biến sex

Câu hỏi nghiên cứu: Xác suất người lao động tại Hoa kỳ là người tham gia công đoàn phụ thuộc vào giới tính là bao nhiêu ?

📌 Thành phần mô hình:

• Random component:
  \(Y \sim \text{Bernoulli}(p)\)

• Systematic component:
  \(\eta = \beta_0 + \beta_1 \times \text{sex}\)

• Link function:
  \(g(p) = \log\left(\frac{p}{1-p}\right)\)

📌 Phương trình mô hình:

\[ \log \left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_1 \times \text{sex}_i \]

Trong đó:

• \(p_i\) là xác suất người lao động thứ \(i\) tham gia công đoàn
  

• \(\text{sex}_i\) là biến giới tính: male hoặc female
  

• \(\beta_0\) là hệ số chặn (intercept)
  

• \(\beta_1\) là hệ số ảnh hưởng của giới tính đến xác suất tham gia công đoàn

• Ước lượng mô hình:

d$union_bin <- ifelse(d$union == "yes", 1, 0)   

logit_model <- glm(union_bin ~ sex, data = d, family = binomial(link = "logit"))
summary(logit_model)

## 
## Call:
## glm(formula = union_bin ~ sex, family = binomial(link = "logit"), 
##     data = d)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.3184     0.1132 -11.651  < 2e-16 ***
## sexmale       0.8432     0.1181   7.139 9.37e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 5461.8  on 4164  degrees of freedom
## Residual deviance: 5404.3  on 4163  degrees of freedom
## AIC: 5408.3
## 
## Number of Fisher Scoring iterations: 4

Phương trình hồi quy:
\(\log \left(\frac{p_i}{1-p_i}\right) = -1.3184 + 0.8432 \times \text{sex}_i\)

• Dự báo: Sử dụng quy ước ngưỡng 0.5 cho phần dự báo.

Dự báo với type = response

Đây là xác suất dự báo một người lao động là nam sẽ tham gia công đoàn theo mô hình logistic.

new_worker <- data.frame(sex = "male")
predict(logit_model, newdata = new_worker, type = "response")

##         1 
## 0.3833874

• Giá trị dự báo: 38.34%. Tức là nếu người lao động đó là nam thì mô hình logistic dự báo rằng xác suất người đó tham gia công đoàn là 38.84%.
  
• Với ngưỡng quy ước là 0.5 thì giá trị dự báo 38.84% cho thấy người lao động nam này không tham gia công đoàn.

Dự báo với type = link
👉 Ý nghĩa:

• Đây là log odds (logit) dự báo cho người nam.
  
• Nó tính theo công thức:

\[ \text{logit}(p) = \beta_0 + \beta_1 \times \text{sex} \]

new_worker <- data.frame(sex = "male")
predict(logit_model, newdata = new_worker, type = "link")

##          1 
## -0.4751948

• Giải thích giá trị:

• -0.4751948 là log odds của việc tham gia công đoàn nếu là nam.

• Tính xác suất từ log odds này:

\[ p = \frac{e^{-0.4751948}}{1 + e^{-0.4751948}} \approx 0.3834 \]

\(\rightarrow\) đúng bằng giá trị ở trên với type = "response"

Dự báo với type = terms

new_worker <- data.frame(sex = "male")
predict(logit_model, newdata = new_worker, type = "terms")

##          sex
## 1 0.09494726
## attr(,"constant")
## [1] -0.5701421

👉 Ý nghĩa:

Đây là phần đóng góp (effect) của từng biến vào giá trị logit dự báo.

Gồm:

• Effect của biến sex = 0.09494726
  
• Intercept (constant) = -0.5701421

Giải thích cụ thể:

\[ \text{Logit}(p) = \text{Intercept} + \text{Effect của sex} \]

Cộng lại:

\[ -0.5701421 + 0.09494726 = -0.4751948 \]

\(\rightarrow\) Chính là giá trị type = "link" phía trên.

2.3 Mô hình hồi quy probit

2.3.1 Biến union và biến sex

Câu hỏi nghiên cứu: Xác suất người lao động tại Hoa kỳ là người tham gia công đoàn phụ thuộc vào giới tính là bao nhiêu ?

📌 Thành phần mô hình:

• Random component:
  \(Y \sim \text{Bernoulli}(p)\)

• Systematic component:
  \(\eta = \beta_0 + \beta_1 \times \text{sex}\)

• Link function:
  \(g(p) = \Phi^{-1}(p)\)

📌 Phương trình mô hình:

\[ \Phi^{-1}(p_i) = \beta_0 + \beta_1 \times \text{sex}_i \]

Trong đó:

• \(\Phi^{-1}\) là hàm phân phối tích lũy chuẩn nghịch đảo (inverse CDF)
  

• Các ký hiệu khác như mô hình Logistic (áp dụng cho trường hợp này là mô hình Probit, với \(p_i\) là xác suất người lao động thứ \(i\) tham gia công đoàn, \(\text{sex}_i\) là biến giới tính, \(\beta_0\) là hệ số chặn và \(\beta_1\) là hệ số ảnh hưởng của giới tính).

• Ước lượng mô hình:

probit_model <- glm(union_bin ~ sex, data = d, family = binomial(link = "probit"))
summary(probit_model)

## 
## Call:
## glm(formula = union_bin ~ sex, family = binomial(link = "probit"), 
##     data = d)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.80265    0.06518 -12.313  < 2e-16 ***
## sexmale      0.50606    0.06847   7.391 1.46e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 5461.8  on 4164  degrees of freedom
## Residual deviance: 5404.3  on 4163  degrees of freedom
## AIC: 5408.3
## 
## Number of Fisher Scoring iterations: 4

Phương trình hồi quy probit:
\[ \Phi^{-1}(p_i) = -0.80265 + 0.50606 \times \text{sex}_i \]

• Dự báo:

Dự báo với typr = response

new_worker <- data.frame(sex = "male")
predict(probit_model, newdata = new_worker, type = "response")

##         1 
## 0.3833874

Ý nghĩa:

• Xác suất tham gia công đoàn nếu là nam là 38.34%.

Nếu quy ước ngưỡng 0.5:

\[ 0.3834 < 0.5 \]

\(\rightarrow\) Dự đoán không tham gia công đoàn.

Dự báo với type = link

new_worker <- data.frame(sex = "male")
predict(probit_model, newdata = new_worker, type = "link")

##          1 
## -0.2965961

Ý nghĩa:

• Đây là log-probit (inverse CDF chuẩn) dự báo cho người nam.

Nó tương đương:

\[ \Phi^{-1}(p_{\text{nam}}) = -0.80265 + 0.50606 \times 1 = -0.2965961 \]

\(\rightarrow\) Nếu bạn lấy \(\Phi(-0.2966) \approx 0.3834\)

\(\Phi(-0.2966) \approx 0.3834\) đúng với type="response"

Dự báo với type=terms

new_worker <- data.frame(sex = "male")
predict(probit_model, newdata = new_worker, type = "terms")

##          sex
## 1 0.05698465
## attr(,"constant")
## [1] -0.3535807

Ý nghĩa:

• Effect của biến sexmale là +0.057
  
• Intercept (constant) là -0.3536

Tổng log-probit dự báo:

\[ -0.3536 + 0.057 = -0.2966 \]

2.4 Mô hình hồi quy cloglog

2.4.1 Biến union và biến sex

Câu hỏi nghiên cứu: Xác suất người lao động tại Hoa kỳ là người tham gia công đoàn phụ thuộc vào giới tính là bao nhiêu ?

📌 Thành phần mô hình:

• Random component:
  \(Y \sim \text{Bernoulli}(p)\)

• Systematic component:
  \(\eta = \beta_0 + \beta_1 \times \text{sex}\)

• Link function:
  \(g(p) = \log(-\log(1-p))\)

📌 Phương trình mô hình:

\[ \log \left( -\log(1-p_i) \right) = \beta_0 + \beta_1 \times \text{sex}_i \]

• Ước lượng mô hình:

cloglog_model <- glm(union_bin ~ sex, data = d, family = binomial(link = "cloglog"))
summary(cloglog_model)

## 
## Call:
## glm(formula = union_bin ~ sex, family = binomial(link = "cloglog"), 
##     data = d)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.4393     0.1007 -14.287  < 2e-16 ***
## sexmale       0.7126     0.1042   6.836 8.17e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 5461.8  on 4164  degrees of freedom
## Residual deviance: 5404.3  on 4163  degrees of freedom
## AIC: 5408.3
## 
## Number of Fisher Scoring iterations: 5

• Dự báo:

new_worker <- data.frame(sex = "male")
predict(cloglog_model, newdata = new_worker, type = "response")

##         1 
## 0.3833874

new_worker <- data.frame(sex = "male")
predict(cloglog_model, newdata = new_worker, type = "link")

##          1 
## -0.7266742

new_worker <- data.frame(sex = "male")
predict(cloglog_model, newdata = new_worker, type = "terms")

##          sex
## 1 0.08024237
## attr(,"constant")
## [1] -0.8069165