BTVN PTDLĐT

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

Tập dữ liệu “SmokeBan” trong gói AER của R là một tập dữ liệu nghiên cứu về tác động của việc áp đặt lệnh cấm hút thuốc tại nơi làm việc đối với việc hút thuốc của những người làm việc. Dưới đây là mô tả và định dạng của tập dữ liệu này:

smoker: cá nhân có phải là người hút thuốc hiện tại hay không.

ban: có áp dụng lệnh cấm hút thuốc tại nơi làm việc hay không.

age: Tuổi của cá nhân

education: cho biết trình độ giáo dục cao nhất đã đạt được: bỏ học trung học (hs), tốt nghiệp trung học, một số đại học, tốt nghiệp đại học, thạc sĩ (hoặc cao hơn).

afam: cho biết liệu cá nhân có phải là người Mỹ gốc Phi-Ameirca hay không.

hispanic: cá nhân có phải là người Mỹ gốc Tây Ban Nha hay không.

gender: giới tính của cá nhân.

library(readxl)

## Warning: package 'readxl' was built under R version 4.2.3

data1 <- read_excel("C:/Users/dell/Downloads/SmokeBan.xlsx")
View(data1)

library(AER)

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

## Loading required package: car

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

## Loading required package: carData

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

## Loading required package: lmtest

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

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.2.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.2.3

## Loading required package: survival

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.2.3

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

## Warning: package 'tibble' was built under R version 4.2.3

## Warning: package 'tidyr' was built under R version 4.2.3

## Warning: package 'readr' was built under R version 4.2.3

## Warning: package 'purrr' was built under R version 4.2.3

## Warning: package 'dplyr' was built under R version 4.2.3

## Warning: package 'stringr' was built under R version 4.2.3

## Warning: package 'forcats' was built under R version 4.2.3

## Warning: package 'lubridate' was built under R version 4.2.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::recode() masks car::recode()
## ✖ purrr::some()   masks car::some()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(DT)

## Warning: package 'DT' was built under R version 4.2.3

library(DescTools)

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

## 
## Attaching package: 'DescTools'
## 
## The following object is masked from 'package:car':
## 
##     Recode

library(ggplot2)

Lập bảng tần số

# Thay đổi giá trị của biến 
data1$smoker<- recode(data1$smoker, "yes" = "có hút", "no" = "không hút")
data1$ban <- recode(data1$ban, "yes" = "cấm", "no" = "không cấm")

Bảng tần số biến ban

table(data1$ban)

## 
##       cấm không cấm 
##      6098      3902

Bảng tần suất

table(data1$ban)/sum(table(data1$ban))

## 
##       cấm không cấm 
##    0.6098    0.3902

Bảng tần số số người hút thuốc

table(data1$smoker)

## 
##    có hút không hút 
##      2423      7577

Bảng tần suất

table(data1$smoker)/sum(table(data1$smoker))

## 
##    có hút không hút 
##    0.2423    0.7577

Vẽ đồ thị

library(ggplot2)

Biến ban

data1 |> ggplot(aes(x = ban, y = after_stat(count))) +
geom_bar(fill = 'red') +geom_text(aes(label = scales::percent(after_stat(count/sum(count)))), stat ='count', color = 'blue', vjust = - .5) +labs(x = 'Ban', y = 'Số nơi cấm')

Nhận xét: Trong tổng số 10000 quan sát:

Số nơi làm việc cấm hút thuốc là 61%

Số nơi làm việc cho hút thuốc là 39%

Biến Smoker

data1 |> ggplot(aes(x = smoker, y = after_stat(count))) +
geom_bar(fill = 'blue') +geom_text(aes(label = scales::percent(after_stat(count/sum(count)))), stat ='count', color = 'purple', vjust = - .5) +labs(x = 'Số người hút thuốc', y = 'số người')

Nhận xét:

Số người đang hút thuốc là 24%

Số người không hút thuốc là 76%

Ước lượng tỷ lệ

Ước lượng tỷ lệ cá nhân hút thuốc ở nơi làm việc cấm hút

smoke<- table(data1[data1$smoker == "có hút",]$ban)
smoke

## 
##       cấm không cấm 
##      1293      1130

prop.test(smoke["cấm"], sum(smoke), p = 0.4)

## 
##  1-sample proportions test with continuity correction
## 
## data:  smoke["cấm"] out of sum(smoke), null probability 0.4
## X-squared = 179.74, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.4
## 95 percent confidence interval:
##  0.5135283 0.5536361
## sample estimates:
##        p 
## 0.533636

Ta ước lượng tỷ lệ với giả thuyết (HO): Tỷ lệ người hút thuốc ở nơi cấm hút là 40%

Kết quả cho p - value là 2.2e-16. Ta bác bỏ giả thuyết tỷ lệ người hút thuốc ở nơi cấm hút là 40%

Ước lượng tỷ lệ cá nhân không hút thuốc ở nơi làm việc cấm hút thuốc

nosmoke<- table(data1[data1$smoker == "không hút",]$ban)
nosmoke

## 
##       cấm không cấm 
##      4805      2772

prop.test(nosmoke["cấm"], sum(nosmoke), p = 0.6)

## 
##  1-sample proportions test with continuity correction
## 
## data:  nosmoke["cấm"] out of sum(nosmoke), null probability 0.6
## X-squared = 36.689, df = 1, p-value = 1.385e-09
## alternative hypothesis: true p is not equal to 0.6
## 95 percent confidence interval:
##  0.6231788 0.6449964
## sample estimates:
##        p 
## 0.634156

Ta ước lượng tỷ lệ với giả thuyết (HO): Tỷ lệ người không hút thuốc ở nơi cấm hút là 40%

Kết quả cho p - value là 1.385e-09-16. Ta bác bỏ giả thuyết tỷ lệ người không hút thuốc ở nơi làm việ c cấm hút là 60%

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

da_hut <- data1[data1$smoker == 'có hút',]
da_khonghut <- data1[data1$smoker == 'không hút',]

hut_cam <- da_hut[da_hut$ban == 'cấm',]
khonghut_cam <- da_khonghut[da_khonghut$ban == 'cấm',]

a <- c(nrow(da_hut), nrow(da_khonghut))
a

## [1] 2423 7577

b <- c(nrow(hut_cam), nrow(khonghut_cam))
b

## [1] 1293 4805

prop.test(b,a)

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  b out of a
## X-squared = 77.54, df = 1, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.12342380 -0.07761622
## sample estimates:
##   prop 1   prop 2 
## 0.533636 0.634156

Ta kiểm định giả thuyết H0: Không có sự khác biệt đáng kể giữa tỉ lệ người hút thuốc và người không ở hút thuốc ở nơi làm việc cấm hút thuốc

Vì p-value = 2.2e-16 < 0.05. Vì vậy bác bỏ giả thuyết H0. Nghĩa là có sự khác biệt đáng kể giữa tỉ lệ người hút thuốc và người không ở hút thuốc ở nơi làm việc cấm hút thuốc

Với độ tin cậy 95%, khoảng chênh lệch giữa tỉ ̣lệ người hút thuốc và người không ở hút thuốc ở nơi làm việc cấm hút thuộc trong khoảng từ 0.533636 đến 0.63415

Ước lượng Relative Risk và Odd Ratio

library(DescTools)
library(epitools)

## 
## Attaching package: 'epitools'

## The following object is masked from 'package:survival':
## 
##     ratetable

da <- table(data1$smoker, data1$ban)
addmargins(da)

##            
##               cấm không cấm   Sum
##   có hút     1293      1130  2423
##   không hút  4805      2772  7577
##   Sum        6098      3902 10000

RelRisk(da)

## [1] 0.8414901

da <- table(data1$smoker, data1$ban)
riskratio(da)

## $data
##            
##              cấm không cấm Total
##   có hút    1293      1130  2423
##   không hút 4805      2772  7577
##   Total     6098      3902 10000
## 
## $measure
##            risk ratio with 95% C.I.
##              estimate     lower     upper
##   có hút    1.0000000        NA        NA
##   không hút 0.7844602 0.7447901 0.8262433
## 
## $p.value
##            two-sided
##             midp.exact fisher.exact   chi.square
##   có hút            NA           NA           NA
##   không hút          0 1.795861e-18 1.050602e-18
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

Nhận xét: Kết quả trên cho ta thấy tỉ lệ rủi ro tương đối (risk ratio) của những người được hút thuốc ở những nơi làm việc cho phép hút là 0.8413901 tức là tỷ lệ những người được hút thuốc ở những nơi làm việc cho phép gấp 0.8414901 lần tỷ lệ những người không được hút thuốc ở những nơi làm việc không cho phép

epitab(da, method = 'oddsratio')

## $tab
##            
##              cấm        p0 không cấm        p1 oddsratio     lower     upper
##   có hút    1293 0.2120367      1130 0.2895951 1.0000000        NA        NA
##   không hút 4805 0.7879633      2772 0.7104049 0.6601155 0.6017953 0.7240874
##            
##                  p.value
##   có hút              NA
##   không hút 1.795861e-18
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

Kết quả tỷ lệ chênh lệch (Oddsratio) là 0.6601155 tức tỉ lệ những người hút thuốc ở nơi làm việc cho phép hút thuốc bé hơn khoảng 33.98% so với tỉ lệ những người không hút thuốc ở nơi không cho phép hút thuốc