BÀI TẬP

Giới thiệu các biến định tính sẽ sử dụng

Bộ dữ liệu bao gồm 4165, 14 biến, trong đó sử dụng 2 biến định tính là:

• Biến married: Biến phân loại. Cá nhân có kết hôn không?

• Biến gender: phân loại giới tính của cá nhân đó là nam hay nữ.

Khai báo các package cần thiết

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ 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

library(epitools)
library(DescTools)
library(DT)
library(energy)
library(ggplot2)
library(AER)

## Loading required package: car
## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:DescTools':
## 
##     Recode
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
## 
## Loading required package: lmtest
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Loading required package: sandwich
## Loading required package: survival
## 
## Attaching package: 'survival'
## 
## The following object is masked from 'package:epitools':
## 
##     ratetable

Bảng dữ liệu

options(digits = 4)
data("PSID7682")
th <- PSID7682
datatable(th)

1. Bảng tần số

Lập bảng tần số cho biến married

table(th$married)

## 
##   no  yes 
##  773 3392

• Tỷ lệ kết hôn chiếm số lượng lớn, có 3392 người đã kết hôn (chiếm 81.4%), 773 người chưa kết hôn (18.6%), cho thấy là có sự chênh lệch rất lớn giữa người đã kết hôn và người chưa kết hôn (62.8%)

Lập bảng tần biến cho biến gender

table(th$gender)

## 
##   male female 
##   3696    469

• Nhận xét, nam là 3696 người (chiếm khoảng 88.7%), nữ là 469 người (chiếm 11.3%), giới tính tương đối đồng đều, với nam và nữ không chênh lệch quá nhiều (chênh lệch khoảng 77.4%).

tp <- th
tp <- table(th$married, th$gender)

addmargins(tp)

##      
##       male female  Sum
##   no   316    457  773
##   yes 3380     12 3392
##   Sum 3696    469 4165

tp <- prop.table(tp)
tp

##      
##           male   female
##   no  0.075870 0.109724
##   yes 0.811525 0.002881

• Nhận xét: Trong số nam, có khoảng 81.15% đã kết hôn, trong khi đó, chỉ có khoảng 7.58% không kết hôn. Trong số nữ, tỷ lệ kết hôn khoảng 0.29% thấp hơn rất nhiều so với nam giới, và tỷ lệ không kết hôn 10.97%.

2. Đồ thị

2.1. Đồ thị cột cho biến married

th |> ggplot(aes(married))+geom_bar(color = "blue", fill = "red")

th |> ggplot(aes(x = th$married, y = after_stat(count))) + geom_bar(fill = 'red') + geom_text(aes(label = scales::percent(after_stat(count/sum(count)))), stat = 'count', color = 'green', vjust = - .5) + theme_classic() + labs(x = 'Giới tính', y = 'Số người')

## Warning: Use of `th$married` is discouraged.
## ℹ Use `married` instead.
## Use of `th$married` is discouraged.
## ℹ Use `married` instead.

th |> ggplot(aes(married)) +
  geom_bar(aes(y = (..count..)/sum(..count..)),color = "blue", fill = "red") +
  ylab('Tỷ lệ %') + xlab('Tình trạng hôn nhân')

## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

2.2. Dồ thị cột cho biến gender

th |> ggplot(aes(gender))+geom_bar(color = "green", fill = "blue")

th |> ggplot(aes(x = th$gender, y = after_stat(count))) + geom_bar(fill = 'blue') + geom_text(aes(label = scales::percent(after_stat(count/sum(count)))), stat = 'count', color = 'green', vjust = - .5) + theme_classic() + labs(x = 'Giới tính', y = 'Số người')

## Warning: Use of `th$gender` is discouraged.
## ℹ Use `gender` instead.
## Use of `th$gender` is discouraged.
## ℹ Use `gender` instead.

th |> ggplot(aes(gender)) +
  geom_bar(aes(y = (..count..)/sum(..count..)),color = "green", fill = "blue") +
  ylab('Tỷ lệ %') + xlab('Giới tính')

2.3. Đồ thị bánh cho biến married

tp <- th
tp <- table(th$married)
tp <- th |> group_by(married) |> summarise(freq = n()) |> mutate(tp, per = freq/sum(freq))
tp |> ggplot(aes(x = '', y = per, fill = married)) +
  geom_bar(stat = 'identity') + 
  coord_polar('y')

2.4. Đồ thị bánh cho biến gender

tp <- th
tp <- table(th$gender)
tp <- th |> group_by(gender) |> summarise(freq = n()) |> mutate(tp, per = freq/sum(freq))
tp |> ggplot(aes(x = '', y = per, fill = )) +
  geom_bar(stat = 'identity') + 
  coord_polar('y')

h <- data.frame(th$gender,th$married)
ggplot(th, aes(x = th$gender, fill = th$married)) +
  geom_bar(position = "dodge") +
  labs(title = "Gender vs Marital Status",
       x = "Gender",
       y = "Count",
       fill = "Marital Status") +
  theme_minimal()

## Warning: Use of `th$gender` is discouraged.
## ℹ Use `gender` instead.

## Warning: Use of `th$married` is discouraged.
## ℹ Use `married` instead.

3. Ước lượng tỉ lệ

3.1. Ước lượng tỉ lệ cho cá nhân đã kết hôn

mar_mal <- table(th[th$married == "yes", ]$gender)
mar_mal

## 
##   male female 
##   3380     12

prop.test(mar_mal["male"], sum(mar_mal), p = 0.4)

## 
##  1-sample proportions test with continuity correction
## 
## data:  mar_mal["male"] out of sum(mar_mal), null probability 0.4
## X-squared = 5026, df = 1, p-value <2e-16
## alternative hypothesis: true p is not equal to 0.4
## 95 percent confidence interval:
##  0.9936 0.9981
## sample estimates:
##      p 
## 0.9965

prop.test(mar_mal["female"], sum(mar_mal), p = 0.6)

## 
##  1-sample proportions test with continuity correction
## 
## data:  mar_mal["female"] out of sum(mar_mal), null probability 0.6
## X-squared = 5026, df = 1, p-value <2e-16
## alternative hypothesis: true p is not equal to 0.6
## 95 percent confidence interval:
##  0.001918 0.006360
## sample estimates:
##        p 
## 0.003538

3.2. Ước lượng tỉ lệ cho cá nhân chưa kết hôn

mar_fe <- table(th[th$married == "no", ]$gender)
mar_fe

## 
##   male female 
##    316    457

prop.test(mar_fe["male"], sum(mar_fe), p = 0.4)

## 
##  1-sample proportions test with continuity correction
## 
## data:  mar_fe["male"] out of sum(mar_fe), null probability 0.4
## X-squared = 0.21, df = 1, p-value = 0.6
## alternative hypothesis: true p is not equal to 0.4
## 95 percent confidence interval:
##  0.3740 0.4445
## sample estimates:
##      p 
## 0.4088

prop.test(mar_fe["female"], sum(mar_fe), p = 0.6)

## 
##  1-sample proportions test with continuity correction
## 
## data:  mar_fe["female"] out of sum(mar_fe), null probability 0.6
## X-squared = 0.21, df = 1, p-value = 0.6
## alternative hypothesis: true p is not equal to 0.6
## 95 percent confidence interval:
##  0.5555 0.6260
## sample estimates:
##      p 
## 0.5912

4. Ước lượng chênh lệch tỷ lệ

th_yes <- th[th$married == "đakethon",]
th_no <- th[th$married == "chuakethon",]

dakethon_nam <- th_yes[th_yes$gender == 'nam',]
chuakethon_nu <- th_no[th_no$gender == 'nu' ,]

a <- c(nrow(th_yes), nrow(th_no))
a

## [1] 0 0

b <- c(nrow(dakethon_nam), nrow(chuakethon_nu))
b

## [1] 0 0

5. Ước lượng Relative risk

d <- table(th$gender,th$married)
addmargins(d)

##         
##            no  yes  Sum
##   male    316 3380 3696
##   female  457   12  469
##   Sum     773 3392 4165

RelRisk(d)

## [1] 0.08774

Nhận xét : RelRisk được tính toán cho việc kết hôn giữa nam và nữ. Nếu RelRisk = 0.08774 như trong trường hợp này, điều này có thể được hiểu là nam có một mức độ rủi ro tương đối cao hơn một chút so với nữ khi liên quan đến việc kết hôn, nhưng sự khác biệt không lớn.

Khoảng ước lượng cho Relative risk

m <- matrix(c(316 ,457 ,3380 ,12),nrow = 2)
RelRisk(m, conf.level = .95)

## rel. risk    lwr.ci    upr.ci 
##   0.08774   0.07885   0.09755

Ước lượng cho Tỷ lệ Rủi ro Tương đối giữa nhóm đã kết hôn và nhóm chưa kết hôn là khoảng 0.09.

Khoảng tin cậy 95% cho ước lượng này là từ khoảng 0.079 đến 0.1. Điều này ngụ ý rằng có 95% khả năng rằng giá trị thực sự của Tỷ lệ Rủi ro Tương đối nằm trong khoảng từ 0.079 đến 0.1.

6. Odd ratio

OddsRatio(d, conf.level = .95)

## odds ratio     lwr.ci     upr.ci 
##   0.002455   0.001368   0.004405

Odds Ratio (Tỷ lệ cơ hội): ước lượng cho tỷ lệ giữa xác suất chênh lệch giữa người nam và nữ (giữa đã kết hôn và chưa kết hôn) là 0.002455

Lwr.ci và Upr.ci: Đây là khoảng tin cậy 95% cho ước lượng của tỷ lệ cơ hội. Nó cho biết phạm vi mà ước lượng tỷ lệ cơ hội có thể nằm trong đó với mức độ tin cậy 95%. Trong trường hợp này, khoảng tin cậy này từ 0.001368 đến 0.004405.

oddsratio(d)

## $data
##         
##           no  yes Total
##   male   316 3380  3696
##   female 457   12   469
##   Total  773 3392  4165
## 
## $measure
##         odds ratio with 95% C.I.
##          estimate    lower    upper
##   male   1.000000       NA       NA
##   female 0.002496 0.001315 0.004279
## 
## $p.value
##         two-sided
##          midp.exact fisher.exact chi.square
##   male           NA           NA         NA
##   female          0            0          0
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"

epitab(d,method = "oddsratio")

## $tab
##         
##           no     p0  yes       p1 oddsratio    lower    upper p.value
##   male   316 0.4088 3380 0.996462  1.000000       NA       NA      NA
##   female 457 0.5912   12 0.003538  0.002455 0.001368 0.004405       0
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

từ bảng ước lượng trên, cho thấy: - giá trị p

7. Hồi quy

Trước tiên, tôi tiến hành kiểm định tính độc lập cho hai biến như sau:

\(H_0\) = Giữa hai biến là độc lập

\(H_1\) = Giữa hai biến có liên quan với nhau

chisq.test(table(PSID7682$gender, PSID7682$married))

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(PSID7682$gender, PSID7682$married)
## X-squared = 2170, df = 1, p-value <2e-16

Kết quả:

• Giá trị p-value = 0 < 5% cho thấy giới tính có ảnh hưởng đến kết hôn

Sau đó, để làm rõ sự ảnh hưởng của giới tính đến khả năng kết hôn, tôi tiến hành phân tích nhị phân Logit như sau:

p <- glm(factor(married) ~ gender, family = binomial(link = "logit"), data = th)
p

## 
## Call:  glm(formula = factor(married) ~ gender, family = binomial(link = "logit"), 
##     data = th)
## 
## Coefficients:
##  (Intercept)  genderfemale  
##         2.37         -6.01  
## 
## Degrees of Freedom: 4164 Total (i.e. Null);  4163 Residual
## Null Deviance:       4000 
## Residual Deviance: 2270  AIC: 2270

Ta có hàm hồi quy sau: y = 2.37 - 6.01*x (Với x là Gender (giới tính))

Nhận xét:

• Sự thay đổi giới tính nam và nữ có khác biệt lớn về khả năng kết hôn với hệ số hồi quy là - 6,01

=> Kết quả cho thấy, nếu là nam thì khả năng sống sót là cao hơn trung bình khoảng 6,01 lần so với giới tính nữ.