https://docs.google.com/document/d/1ewxaj6_REpFTahzebKXMGtHsOiUIB7nhrLi6f7PKHhI/edit?usp=sharing

1 Phụ lục 2. Thực hành bảng ngẫu nhiên và suy diễn thống kê trên R

1.1 Bảng ngẫu nhiên

Cấu trúc xác suất cho các bảng ngẫu nhiên

Dữ liệu định tính là tần số suất hiện các biểu hiện của các biến. Cho X và Y là hai biến định tính, X có k biểu hiện: A1, 𝐴A2, … ,Am và Y có m biểu hiện: B1, B2, … , 𝐵m . Chúng ta có thể sử dụng một bảng gồm k hàng và m cột để thể hiện kết quả có thể xảy ra từ việc khảo sát:

Bảng này được gọi là bảng ngẫu nhiên hai chiều 𝑘 × 𝑚, trong đó nij là số lần quan sát được cặp thuộc tính (Ai, Bj) còn gọi là tần số của (Ai, Bj). Một bảng ngẫu nhiên hai biến được gọi là bảng hai chiều; một bảng ngẫu nhiên 3 biến gọi là bảng ngẫu nhiên 3 chiều. Trong mục này, chúng ta sử dụng bộ dữ liệu afmahoa(bộ dữ liệu đã xử lý được mã hóa từ file dữ liệu afterlife.xlsx) để làm ví dụ.

1.1.1 Bảng tần số và tần suất

Lập bảng tần số và tần suất cho biến Sex

Để lập bảng tần số và tần suất cho biến Sex, ta thực hiện câu lệnh sau:

af<-read.xlsx(“afterlife.xlsx”, sheetIndex = 1, header = TRUE): đọc file afterlife.xlsx và đặt tên af.

afmahoa<- subset(af, Believe <5): câu lệnh trích xuất những quan sát mà giá trị của biến Believe <5 và đặt tên afmahoa.

table(afmahoa$Sex): câu lệnh lập bảng tần số cho biến Sex.

table(afmahoa\(Sex/(sum(afmahoa\)Sex)): câu lệnh lập bảng tần suất cho biến Sex

library(readxl)
## Warning: package 'readxl' was built under R version 4.2.3
af <- read_excel("D:/afterlife.xlsx")
af
## # A tibble: 1,517 × 2
##      Sex Believe
##    <dbl>   <dbl>
##  1     2       2
##  2     2       4
##  3     1       4
##  4     2       4
##  5     2       6
##  6     1       6
##  7     1       5
##  8     2       2
##  9     2       2
## 10     2       5
## # ℹ 1,507 more rows
afmahoa <- subset(af, Believe <5)
afmahoa
## # A tibble: 1,260 × 2
##      Sex Believe
##    <dbl>   <dbl>
##  1     2       2
##  2     2       4
##  3     1       4
##  4     2       4
##  5     2       2
##  6     2       2
##  7     2       2
##  8     1       2
##  9     2       3
## 10     2       2
## # ℹ 1,250 more rows
head(afmahoa) # Xem 6 quan sát đầu afmahoa
## # A tibble: 6 × 2
##     Sex Believe
##   <dbl>   <dbl>
## 1     2       2
## 2     2       4
## 3     1       4
## 4     2       4
## 5     2       2
## 6     2       2
# Lập bảng tần số cho biến sex và kết quả
table(afmahoa$Sex)
## 
##   1   2 
## 519 741

Kết quả trên được trình bày lại như sau:

Sex 1(Nam) 2(Nữ)

Tần số (Số người) 519 741

# Lập bảng tần suất cho biến sex và kết quả
table(afmahoa$Sex)/sum(table(afmahoa$Sex))
## 
##         1         2 
## 0.4119048 0.5880952

Chúng ta có thể trình bày kết quả trên như sau:

Sex 1(Nam) 2(Nữ) Tần số (Số người) 579 741 Tần suất (Tỷ lệ) 0,4119048 0,5880952

Bằng các câu lệnh tương tự chúng ta có thể lập bảng tần số và tần suất cho biến Believe và ta được kết quả.

head(afmahoa) # Xem 6 quan sát đầu afmahoa
## # A tibble: 6 × 2
##     Sex Believe
##   <dbl>   <dbl>
## 1     2       2
## 2     2       4
## 3     1       4
## 4     2       4
## 5     2       2
## 6     2       2
# Lập bảng tần số cho biến Believe và kết quả
table(afmahoa$Believe)
## 
##   1   2   3   4 
## 158 674 278 150
# Lập bảng tần suất cho biến Believe và kết quả
table(afmahoa$Believe)/sum(table(afmahoa$Believe))
## 
##         1         2         3         4 
## 0.1253968 0.5349206 0.2206349 0.1190476

1.1.2 Bảng hai chiều

Để tạo bảng hai chiều chúng ta sử dụng câu lệnh sau: table(afmahoa\(Sex, afmahoa\)Believe): câu lệnh tạo bảng hai chiều cho 2 biến Sex và Believe.

table(afmahoa$Sex, afmahoa$Believe) # Bảng hai chiều
##    
##       1   2   3   4
##   1  79 252 118  70
##   2  79 422 160  80

1.2 Tính OddsRatio

Để tính OddsRatio của hàng 1 và hàng 2 trong bảng hai chiều, ta thực hiện câu lệnh sau:

B1. install.packages(“DescTools”): cài đặt gói DescTools.

B2. library(DecsTools): gọi gói DescTools.

B3. OddsRatio(data, conf.level = NULL, y = NULL, method = “wald”): câu lệnh tính giá trị OddsRatio; trong đó: data là bảng hai chiều; conf.level: độ tin cậy mặc nhiên 95%; y là véc tơ các giá trị cần tính giá trị OddsRatio; method = “wald”: phương pháp wald.

Ta lấy ví dụ sau để minh họa cho việc tính OddsRatio như sau:

    Group
                    Myocardial Infarction (nhồi máu cơ tim)
                        Yes    No     Total
    Placebo (giả dược) 189    10845   11034
      Aspirin           104   10933   11037

Ta lấy ví dụ sau để minh họa cho việc tính OddsRatio như sau:

Ta thực hiện các lệnh sau:

B1. library(DecsTools): gọi packages DescTools.

B2. v <-c(189, 104, 10845, 10933): câu lệnh tạo véc tơ v.

B3.data <-matrix(v, 2): câu lệnh tạo ma trận cấp 2 (bảng 2 chiều) và đặt tên data. data: câu lệnh xem data.

B4. y<-c(189, 104): câu lệnh tạo véc tơ y.

B5. OddsRatio(data, conf.level = NULL, y = NULL, method = “wald”): câu lệnh tính OddsRatio

library(DescTools) # Gọi package DesTools
## Warning: package 'DescTools' was built under R version 4.2.3
v <- c(189,104,10845,10933) # Tạo vecto v
x <- matrix(v,2,2) # Tạo ma trận cấp 2 đặt tên x
y <- c(189,104) # Vecto y là những người nhồi máu cơ tim
# Tính OddsRatio
OddsRatio(x,conf.level=0.95, y = NULL, method = "wald")
## odds ratio     lwr.ci     upr.ci 
##   1.832054   1.440042   2.330780

Kết quả tỷ lệ chênh (OddsRatio) là 1,832054. Giá trị này cho chúng ta biết tỷ lệ người bị bệnh nhồi máu cơ tim đối với nhóm người dùng Placebo (giả dược) cao hơn 83,2% so với nhóm người sử dụng thuốc Aspirin.

1.3 Một số bài toán ước lượng

1.3.1 Ước lượng tỷ lệ một tổng thể

Để làm bài toán ước lượng tỷ lệ bằng R, ta thực hiện câu lệnh sau:

B1. install.packages(“DescTools”): cài đặt gói DescTools.

B2. library(DecsTools): gọi gói DescTools.

B3. BinomCI(k, n, conf.level = NULL, method = “wald”): câu lệnh ước lượng tỷ lệ một tổng thể; trong đó, k: số phần tử có tính chất A, n: số phần tử mẫu, conf.level = NULL: độ tin cậy mặc định là 95%, method = wald: phương pháp wald.

Chúng ta lấy dữ liệu từ file có tên afterlife.xlsx dể minh họa cho các câu lệnh trên.

Để ước lượng tỷ lệ những người tin rằng “Tồn tại một thế giới khác sau khi chết” ta làm như sau:

library(DescTools) #Gọi package DecsTools để sử dụng hàm BinomCI
# Đọc file afterlife.xlsx và đặt tên là af
af <- read_excel("D:/afterlife.xlsx")
# Trích xuất biến Believe trong af và đặt tên af1
af1 <- af$Believe
# Trích xuất những giá trị nhỏ hơn 5 của biến Believe trong af1
af2 <- af1[af1<5]
k <- length(af2[af2<3]) # Tổng tần số những người tin của biến Believe
n <- length(af2)        # Tổng tần số của biến Believe(cỡ mẫu)
# Câu lệnh ước lượng và kết quả
BinomCI(k,n, conf.level= 0.95, method = "wald")
##            est    lwr.ci    upr.ci
## [1,] 0.6603175 0.6341672 0.6864677

Ta được khoảng ước lượng tỷ lệ những người tin “Tồn tại thế gới khác sau khi chết” với độ tin cậy 95% là: 0.6341672 ≤ 𝑝 ≤ 0.6864677

1.3.2 Ước lượng sự khác biệt tỷ lệ hai tổng thể

Để thực hiện ước lượng sự khác biệt tỷ lệ hai tổng thể bằng R, ta thực hiện các câu lẹnh sau:

B1. install.packages(“DescTools”): câu lệnh cài đặt gói DescTools.

B2. library(DecsTools): câu lệnh gọi gói DescTools. B3. BinomDiffCI(k2, n2, k1, n1, conf.level = NULL, method = “wald”): câu lệnh ước lượng sự khác biệt hai tỷ lệ tổng thể, ki: số phần tử có tính chất A của mẫu i, ni: số phần tử mẫu thứ i, conf.level = NULL: độ tin cậy mặc định là 95%, method = wald: phương pháp wald.

Chúng ta lấy dữ liệu từ file có tên afterlife.xlsx dể minh họa cho các câu lệnh trên.

Bài toán: ước lượng sự khác biệt tỷ lệ Nam và Nữ tin “Tồn tại thế giới khác sau khi chết”

library(DescTools) #Gọi package DecsTools để sử dụng hàm BinomCI
# Đọc file afterlife.xlsx và đặt tên là af
af <- read_excel("D:/afterlife.xlsx")
# Trích xuất những quan sát mà giá trị nhỏ hơn 5 của biến Believe
af1 <- subset(af, Believe<5)
table(af1)
##    Believe
## Sex   1   2   3   4
##   1  79 252 118  70
##   2  79 422 160  80
n1 <- table(af1$Sex)[1] #Tỏng tần số giới tính Nam
n2 <- table(af1$Sex)[2] #Tỏng tần số giới tính Nữ
# Tổng tần số người Nam "Tin có thê giới khác sau chết
k1 <- table(af1$Sex,af1$Believe)[1,1] + table(af1$Sex,af1$Believe)[1,2]
# Tổng tần số người Nữ "Tin có thê giới khác sau chết
k2 <- table(af1$Sex,af1$Believe)[2,1] + table(af1$Sex,af1$Believe)[2,2]
#Câu lệnh ước lượng và kết quả
BinomDiffCI(k1,n2,k1,n1,conf.level=0.95, method = "wald")
##           est.2     lwr.ci    upr.ci
## [1,] -0.1910713 -0.2457635 -0.136379

Kết quả khoảng ước lượng sự khác biệt tỷ lệ Nam và Nữ tin “Tồn tại thế gới khác sau khi chết” với độ tin cậy 95% là: −0,01499184 ≤ 𝑝1 − 𝑝2 ≤0,09168869

1.4 Kiểm định tính độc lập

1.4.1 Kiểm định tính độc lập cho hai biến định danh

Bài toán: Giả sử X và Y là hai biến quan sát định danh. Qua điều tra, biến X có k thuộc tính (hay k dấu hiệu): A1, A2,…, Ak; biến Y có m thuộc tính: B1, B2,…, Bm, nhận được bảng sau, với nij là số lần xuất hiện cặp (Ai,Bj) và ∑ij 𝑛ij = 𝑛.

Với mức ý nghĩa 𝛼, hãy xác minh xem X và Y có độc lập hay không. Việc xác minh này được gọi là kiểm định về tính độc lập.

Giả thuyết H0: X và Y độc lập nhau,

Tiêu chuẩn bác giả thuyết H0:

Hoặc p - value = Prob(chi-square) < a

Để thực hiện kiểm định tính độc lập cho hai biến định danh bằng R ta làm như sau.

B1. Nhập dữ liệu (bảng ngẫu nhiên) cho đề bài.

B2. Thực hiện kiểm định tính độc lập 2 biến định danh bởi câu lệnh: chisq.test( ).

Chúng ta minh họa kiểm định trên thông qua ví dụ sau:

Để xác minh xem sự ủng hộ của người dân trong nước về một sắc thuế mới có phụ thuộc vào mức thu nhập của họ hay không, tiến hành điều tra ngẫu nhiên 1000 công dân, người ta có số liệu sau: Thu nhập Thái độ Thấp Trung bình Cao Ủng hộ 182 213 203 Phản đối 154 138 110 Với mức ý nghĩa 5%, dựa vào điều tra, hãy cho kết luận về việc này.

Gọi Y là thái độ của công dân đối với sắc thuế mới, với 2 thuộc tính: “ủng hộ”, “phản đối”; X là mức thu nhập của công dân, với 3 thuộc tính: “thấp”, “trung bình”, “cao”. Bài toán ở đây là kiểm định về tính độc lập của X và Y:

Giả thuyeesy H0: X,Y độc lập; đối thuyết H1:X,Y không độc lập.

Ta tiến hành kiểm định tính độc lập X và Y bằng R theo các bước trên ta được kết quả như sau:

#Kiểm định sự độc lập của hai biến
ungho <- c(182,213,203) # Tạo  vecto có tên ủng hộ
phandoi <- c(154,138,110) # Tạo vecto có tên phản đối
data <- as.data.frame(rbind(ungho, phandoi)) # Tạo bảng ngẫu nhiên 2x3
names(data)<- c('Thap', 'Trung binh', 'Cao') # Đặt tên cho bảng ngẫu nhiên
data # Bảng ngẫu nhiê
##         Thap Trung binh Cao
## ungho    182        213 203
## phandoi  154        138 110

Từ bảng kết quả trên, ta có 𝑝 − 𝑣𝑎𝑙𝑢𝑒 = 0,01947 < 0,05: Bác bỏ H0

Kết luận, sự ủng hộ của người dân trong nước về một sắc thuế mới có phụ thuộc vào mức thu nhập của họ, với mức ý nghĩa 5%.

Chú ý: Đối với bài toán này, trên Eviews thao tác như sau:

  • Nhập dữ liệu của X và của Y

  • Từ cửa số Group, chọn: View→ 𝑁 − 𝑤𝑎𝑦 𝑇𝑎𝑏𝑢𝑙𝑎𝑡𝑖𝑜𝑛 → 𝑂𝐾

  • Căn cứ vào cột Prob. (p-value) Của một trong hai thống kê: Pearson hoặc Liklihood ratio để kết luận.

Nếu p-value < 𝛼 thì bác H0 (X và Y độc lập)

Với ví dụ trên, Eviews cho kết quả:

Test Statistics df Value Prob

Pearson X2 2 7.878212 0.0195

Likelihood Ratio G2 2 7.874335 0.0195

1.4.2 Kiểm định tính độc lập cho hai biến định tính có thang đo thứ bậc (CMH test)

Khi biến dòng X và biến cột Y là thứ tự, một liên kết “xu hướng” là khá phổ biến. Đó là khi mức X tăng lên, các phản ứng đối với Y có xu hướng tăng lên các cấp độ cao hơn, hoặc các phản ứng đối với Y có khuynh hướng giảm xuống các cấp độ thấp hơn. Người ta có thể sử dụng một tham số để mô tả mối liên kết xu hướng thứ bậc như vậy. Phân tích phổ biến nhất là cho điểm và định tính mức độ của xu hướng tuyến tính hoặc tương quan.

Phần tiếp theo trình bày một thống kê kiểm định nhạy cảm với các xu hướng tuyến tính tích cực hoặc tiêu cực trong mối quan hệ giữa X và Y. Nó sử dụng các thông tin quan trong dữ liệu.

  • Tiến hành gán các điểm số 𝑢$ ≤ 𝑢% ≤ ⋯ ≤ 𝑢! cho các hàng, và 𝑣$ ≤ 𝑣% ≤ ⋯ ≤ 𝑣& cho các cột. Điểm có cùng thứ tự như các cấp độ định tính và được cho là đơn điệu. Nguyên tắc gán điểm là: điểm số phản ánh khoảng cách giữa các loại, với khoảng cách lớn hơn giữa các loại được sắp xếp xa nhau hơn. Khi đó bảng thống kê 𝑘 × 𝑚 có dạng:

Giả thuyết 𝐻, : X và Y độc lập (thống kê) Tiêu chuẩn bác giả thuyết H0

𝑀2 ≥ 𝜒12(𝛼) ; ( 𝑀2 = (𝑛 − 1)𝑟2)

Hoặc 𝑝 − 𝑣𝑎𝑙𝑢𝑒 = 𝑃𝑟𝑜𝑏(𝑐ℎ𝑖 − 𝑠𝑞𝑢𝑎𝑟𝑒) < 𝛼 Để thực hiện kiểm định tính độc lập cho biến định tính có thang đo thứ bậc ta thực hiện các bước sau:

B1. Cài đặt và gọi packages (“vcdExtra”)

B2. Thực hiện kiểm định bằng câu lệnh: CMHtest(data, rscores = c(), cscores = c()), trong đó data: bảng ngẫu nhiên, rscores = c() là điểm cho tương ứng theo biến X (hàng), cscores = c() là điểm cho tương ứng theo biến Y (cột).

Chúng ta minh họa kiểm định trên bằng ví dụ sau: Theo dõi ngẫu nhiên một số sản phẩm về chất lượng: Loại 1, Loại 2, Loại 3 được sản xuất ở các ca: ngày, đêm, có kết quả sau:

Chất lượng sản phẩm Số sản phẩm Loại 1 Loại 2 Loại 3 Ca ngày 118 28 10 Ca đêm 81 15 9

Với mức ý nghĩa 5%, có thể cho rằng chất lượng sản phẩm phụ thuộc vào ca sản xuất hay không?

Ta cần kiểm định về tính độc lập của X (ca sản xuất) và Y (chất lượng sản phẩm), giả thuyết H0 : X và Y độc lập nhau.

Đối với X, ta gán 0 cho ca đêm, 1 cho ca ngày. Đối với Y ta sử dụng điểm midrank : Loại 3 được gán số 5,5 ; loại 2 được gán số 24,5; loại 1 được gán số 97,5.

Các bước kiểm định tính độc lập cho 𝑋 và 𝑌 được thực hiện như sau:

#KIEM DINH SU DOC LAP CHO HAI BIEN DINH TINH CO THANG DO THU BAC
cadem <- c(81,15,9) # Tạo vecto có tên ca đêm
cangay <- c(118,28,10) # Tạo vecto có tên ca ngày
# Tạo bảng ngẫu nhiên 2x3
data <- as.matrix(rbind(cadem,cangay))
data #Bảng ngẫu nhiên
##        [,1] [,2] [,3]
## cadem    81   15    9
## cangay  118   28   10
library(vcdExtra)
## Warning: package 'vcdExtra' was built under R version 4.2.3
## Loading required package: vcd
## Warning: package 'vcd' was built under R version 4.2.3
## Loading required package: grid
## Loading required package: gnm
## Warning: package 'gnm' was built under R version 4.2.3
# Câu lệnh kiểm định X=(0,1), Y=(97.5,24.5,5.5)
CMHtest(data, rscores = c(0,1), cscore= c(97.5,24.5,5.5))
## Cochran-Mantel-Haenszel Statistics 
## 
##                  AltHypothesis    Chisq Df    Prob
## cor        Nonzero correlation 0.025707  1 0.87262
## rmeans  Row mean scores differ 0.025707  1 0.87262
## cmeans  Col mean scores differ 0.928770  2 0.62852
## general    General association 0.928770  2 0.62852

Từ bảng kết quả trên ta có 𝑝 − 𝑣𝑎𝑙𝑢𝑒 = 0,87260 > 0,05: Chấp nhận H0,

Kết luận, chất lượng sản phẩm không phụ thuộc vào ca sản xuất, với mức ý nghĩa 5%.