File code:

https://drive.google.com/file/d/11_uXVovoKhgdxnbltdrQQh_WciB2e-Jy/view?usp=sharing

7.1. Trích xuất dữ liệu theo dòng (hàng)

Có thể trích xuất ra một dataframe con bằng cách chỉ định cụ thể quan sát nào muốn lấy theo vị trí của hàng. Câu lệnh như sau:

CPS1988<-read.csv(“crab.csv”, header =TRUE): đọc file có tên CPS1988.csv trong thư mục D:/PhanmemR và đặt tên CPS1988.

CPS1<-CPS1988[c(1, 3, 100),]:trích xuất các quan sát dòng 1, 3, 100 của đối tượng CPS1988 và đặt tên CPS1.

CPS2<-CPS1988[c(1: 3, 10, 100),]:trích xuất các quan sát dòng từ 1 đến dòng 3, dòng 10, dòng 100 của đối tượng CPS1988 và đặt tên CPS2.

Kết quả của hai câu lệnh trên:

# Đọc file CPS1988 và đặt tên CPS1988.
CPS1988 <- read.csv("CPS1988.csv", header = TRUE)
CPS1<-CPS1988[c(1, 3, 100),]# Trích xuất quan sát dòng 1, 3, 100
CPS1 # Xem 6 quan sát đầu tiên của CPS1
##       wage education experience ethnicity smsa    region parttime
## 1   354.94         7         45      cauc  yes northeast       no
## 3   370.37         9          9      cauc  yes northeast       no
## 100 688.51        18         16      cauc   no northeast       no
## Trích xuất quan sát từ dòng 1 đến 3, dòng 10, dòng 100 
CPS2<-CPS1988[c(1: 3, 10, 100),]
CPS2 # Xem 6 quan sát đầu tiên của CPS2
##        wage education experience ethnicity smsa    region parttime
## 1    354.94         7         45      cauc  yes northeast       no
## 2    123.46        12          1      cauc  yes northeast      yes
## 3    370.37         9          9      cauc  yes northeast       no
## 10  1643.83        14         18      cauc  yes northeast       no
## 100  688.51        18         16      cauc   no northeast       no

7.2. Trích xuất dữliệu theo cột

Trong một số trường hợp chúng ta chỉ quan tâm đến một số biến chứ không phải toàn bộ các biến. Để trích xuất dữ liệu theo cột ta thực hiện lệnh sau:

CPS1988<-read.csv(“crab.csv”, header =TRUE): đọc file có tên CPS1988.csv trong thư mục D:/PhanmemR và đặt tên CPS1988.

CPS3<-CPS1988[c(2:4)]:lấy cột từ 2 đến cột 4 của đối tượng CPS1988 và đặt tên CPS3

head(CPS3):xem 6 quan sát đầu tiên của đối tượng CPS3.

CPS3<-CPS1988[c(2:4)] #Trích xuất từ 2 đến cột 4 của đối tượng CPS1988 và đặt tên CPS3
head(CPS3)# Xem 6 quan sát đầu tiên của CPS3
##   education experience ethnicity
## 1         7         45      cauc
## 2        12          1      cauc
## 3         9          9      cauc
## 4        11         46      cauc
## 5        12         36      cauc
## 6        16         22      cauc

7.3. Trích xuất dữ liệu theo hàng và cột

Để trích xuất dữ liệu theo hàng và cột ta thực hiện câu lệnh:

CPS1988<-read.csv(“crab.csv”, header =TRUE): câu lệnh đọc file có tên CPS1988.csv trong thư mụcD:/PhanmemR và đặt tên CPS1988.

CPS4 <-CPS1988[c(1, 3, 5, 7), c(2:4)]:câu lệnh trích xuất dữ liệu từ cột 2 đến cột 4, dòng 2 đến dòng 4 của đối tượng CPS1988 và đặt tên CPS4.

head(CPS4):câu lệnh xem 6 quan sát đầu tiên của CPS4.

CPS4 <-CPS1988[c(1, 3, 5, 7), c(2:4)]
head(CPS4)# Xem 6 quan sát đầu tiên của CPS4
##   education experience ethnicity
## 1         7         45      cauc
## 3         9          9      cauc
## 5        12         36      cauc
## 7         8         51      cauc

Kết quả câu lệnh trên tạo ra một dataframe gồm các dòng 1, 3, 5, 7 và chỉ lấy các biến ở cột 2, 3, 4 từ data frame gốc ban đầu.

7.4. Trích xuất dữ liệu từ một dataframe có sẵn bằng lệnh subset

Trong nhiều tình huống, chúng ta cần trích xuất dữ liệu cho các quan sát theo một tiêu chí nào đó. Chẳng hạn, từ tập dữ liệu CPS1988.csv chúng ta muốn trích xuất tất cảcác quan sát ứng với parttime=yes mà thôi với tên gọi CPS5 chẳng hạn, ta thực hiện câu lệnh sau:

CPS1988<-read.csv(“crab.csv”, header =TRUE): đọc file có tênCPS1988.csv trong thư mục D:/PhanmemR và đặt tên CPS1988.

CPS5 <-subset(CPS1988, parttime == “yes”):trích xuất các quan sát của biến parttime = yes trong đối tượng CPS1988 và đặt tên CPS5.

CPS5 <-subset(CPS1988, parttime == "yes")
dim(CPS5)# Số các quan sát của biến pảttime = yes
## [1] 2524    7
head(CPS5)# Xem 6 quan sát đầu tiên của CPS5
##      wage education experience ethnicity smsa    region parttime
## 2  123.46        12          1      cauc  yes northeast      yes
## 38 166.67        12         51      cauc   no northeast      yes
## 41 222.93        11         41      cauc   no northeast      yes
## 43 277.78        11         40      cauc  yes northeast      yes
## 55 284.90        16          2      cauc  yes northeast      yes
## 64 109.04        12          1      cauc  yes northeast      yes

Kết quả câu lệnh cho thấy, có 2524 quan sát có làm bán thời gian được trích xuất.

Tương tự, chúng ta cũng có thể trích xuất dữ liệu theo một số tiêu chí nào đó.

Ví dụ, cũng từ dữ liệu CPS1988 ta trích xuất dữ liệu cho các quan sát các tiêu chí mà: (1) ứng với parttime =yes, và (2) có education ≥ 3 với đặt tên CPS6, ta thực hiện câu lệnh sau:

CPS6 <-subset(CPS1988, parttime == “yes” & education == 5):trích xuất các quan sát của biến parttime = yes và các quan sát của biến education = 5 trong đối tượng CPS1988 với tên CPS6.

head(CPS6):câu lệnh xem 6 quan sát đầu tiên của CPS6.

CPS6 <-subset(CPS1988, parttime == "yes" & education == 5)
head(CPS6)# Xem 6 quan sát đầu tiên của CPS6
##         wage education experience ethnicity smsa  region parttime
## 8024  118.71         5         56      afam  yes midwest      yes
## 15108  98.24         5         46      afam  yes   south      yes
## 16841 130.58         5         56      afam   no   south      yes
## 17806 178.33         5         13      cauc  yes   south      yes
## 20135  74.07         5         56      cauc  yes   south      yes
## 21488 115.74         5         53      cauc  yes   south      yes

Hoặc chúng ta cũng có thể trích xuất dữ liệu cho các quan sát trong CPS1988 theo các tiêu chí (1) ứng với education <2 hoặc (2) experience >50 và đặt tên CPS7, ta thực hiện câu lệnh:

CPS7 <-subset(CPS1988, parttime == “yes” | experience > 50):câu lệnh trích xuất các quan sát của biến parttime = yes hoặc các quan sát của biến experience > 50 trong đối tượng CPS1988 và đặt tên CPS7.

head(CPS7):câu lệnh xem 6 quan sát đầu tiên của CPS7.

CPS7 <-subset(CPS1988, parttime == "yes" | experience > 50)
head(CPS7)# Xem 6 quan sát đầu tiên của CPS7
##      wage education experience ethnicity smsa    region parttime
## 2  123.46        12          1      cauc  yes northeast      yes
## 7  284.90         8         51      cauc  yes northeast       no
## 38 166.67        12         51      cauc   no northeast      yes
## 41 222.93        11         41      cauc   no northeast      yes
## 43 277.78        11         40      cauc  yes northeast      yes
## 55 284.90        16          2      cauc  yes northeast      yes

8. Mã hóa dữ liệu

Trong mục này chúng ta sửdụng dữliệu của file afterlife.xlsx trong thư mục D:/PhanmemRđểlàm ví dụminh họa.

Trước tiên chúng ta mô tả sơ bộ về file dữ liệu afterlife.xlsx bằng câu lệnh sau:

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

af: câu lệnh cho biết chi tiết đối tượng af.

setwd("D:/PhanmemR")
# Đọc file afterlife.xlsx và đặt tên là af
library(readxl)
af <- read_excel("afterlife.xlsx")
head(af) # Xem 6 quan sát đầu
## # A tibble: 6 × 2
##     Sex Believe
##   <dbl>   <dbl>
## 1     2       2
## 2     2       4
## 3     1       4
## 4     2       4
## 5     2       6
## 6     1       6

Dữ liệu file afterlife.xlsx có hai biến định tính.

Biến thứ nhất: Sex (giới tính) được mã hóa như sau: 1 là nam và 2 là nữ.

Biến thứ 2: Believe (niềm tin) của người được khảo sát “Sau khi chết có tồn tại một thếg iới khác” được mã hóa như sau:

  1: Tin là có, khẳng định chắc chắn.

  2: Tin là có, nhưng, nhưng không chắc chắn.

  3: Tin là không có, nhưng không chắc chắn.

  4: Tin là không có, khẳng định chắc chắn.

  5: Dữ liệu bị khuyết.

  6: Không thể chọn.

  7: Không trả lời.

Trong dữ liệu này chúng ta chỉ giữ lại những quan sát mà biến Believe < 5, sau đó thay những quan sát có giá trị 1 và 2 của biến Believe thành giá trị 1 (những quan sát Tin là có), đồng thời những quan sát có giá trị 3 và 4 của biến Believe thành giá trị 0 (những quan sát Tin là không có).

Để thực hiện công việc trên ta thực hiện các lệnh sau:

af<-read.xlsx(“afterlife.xlsx”, sheetIndex =1, header =TRUE) :đọc file tên 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.

afmahoa\(Believe[afmahoa\)Believe <=2] <-1: câu lệnh thấy những giá trị quan sát ≤2 trong cột Believe của afmahoa bằng 1.

afmahoa\(Believe[afmahoa\)Believe >2] <-0:câu lệnh thấy những giá trị quan sát >2 trong cột Believe của afmahoa bằng 0.

# Chỉ giữ lại những quan sát mà giá trị của biến Believe <5 và đặt tên afmahoa.
afmahoa<-subset(af, Believe <5)
dim(afmahoa)  # Xem 6 quan sát đầu afmahoa
## [1] 1260    2
afmahoa$Believe[afmahoa$Believe <=2] <-1
# Trong cột biến Believe của afmahoa thấy những giá trị ≤2 bằng 1.
afmahoa$Believe[afmahoa$Believe >2] <-0
# Trong cột biến Believe của afmahoa thấy những giá trị >2 bằng 0
head(afmahoa)
## # A tibble: 6 × 2
##     Sex Believe
##   <dbl>   <dbl>
## 1     2       1
## 2     2       0
## 3     1       0
## 4     2       0
## 5     2       1
## 6     2       1

Sau khi thực hiện các câu lệnh trên, kết quả cho thấy đối tượng afmahoa chỉ còn 1260 quan sát, biến Believe chỉ nhận một trong hai giá trị: 0, 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. 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: \(𝐴_{1}\),\(𝐴_{2}\),…,\(𝐴_{k}\) và Y có m biểu hiện: \(B_{1}\),\(B_{2}\),…,\(B_{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:

X \(B_{1}\) \(B_{2}\) … \(B_{m}\)
\(𝐴_{1}\) \(𝑛_{11}\) \(𝑛_{12}\) … \(𝑛_{1m}\)
\(𝐴_{2}\) \(𝑛_{21}\) \(𝑛_{22}\) … \(𝑛_{2m}\)
: : : … :
\(𝐴_{k}\) \(𝑛_{k1}\) \(𝑛_{k2}\) … \(𝑛_{km}\)

Bảng này được gọi là bảng ngẫu nhiên hai chiều 𝑘×𝑚, trong đó \(𝑛_{ij}\)là số lần quan sát được cặp thuộc tính (\(𝐴_{i}\),\(B_{j}\)), còn gọi là tần số của (\(𝐴_{i}\),\(B_{j}\)). 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. 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.

head(afmahoa) # Xem 6 quan sát đầu afmahoa
## # A tibble: 6 × 2
##     Sex Believe
##   <dbl>   <dbl>
## 1     2       1
## 2     2       0
## 3     1       0
## 4     2       0
## 5     2       1
## 6     2       1
# 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) 519 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       1
## 2     2       0
## 3     1       0
## 4     2       0
## 5     2       1
## 6     2       1
# Lập bảng tần số cho biến Believe và kết quả
table(afmahoa$Believe)
## 
##   0   1 
## 428 832
# Lập bảng tần suất cho biến Believe và kết quả
table(afmahoa$Believe)/sum(table(afmahoa$Believe))
## 
##         0         1 
## 0.3396825 0.6603175

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

Sex 0(Không tin) 1(Tin)
Tần số (Số Người) 428 832
Tần suất (Tỷlệ) 0,3396825 0,6603175

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)
##    
##       0   1
##   1 188 331
##   2 240 501

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

Sex Believe
0(Không tin) 1(Tin)
1(Nam) 188 331
2(Nữ) 240 501

2.Tính OddsRatioD

Để 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 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 packages DescTools
## Warning: package 'DescTools' was built under R version 4.2.3
## 
## Attaching package: 'DescTools'
## The following object is masked from 'package:car':
## 
##     Recode
v <- c(189, 104, 10845, 10933) # Tạo véc tơ v
x <- matrix(v, 2, 2) # Tạo ma trận cấp 2 đặt tên là x
y<-c(189, 104) # Verto Y là những người nhồi máu cơ tim 
# Tính OddsRatio

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.

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

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

Công thức ước lượng

\(f- u_\frac{\alpha}{2}\) ≤𝑝≤ \(f + u_{\frac{\alpha}{2}}\)

Để làm bài toán ước lươnhj tỷ lrrj 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(DescTools): 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 để minh họa

Để ướ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 DescTools để sử dụng hàm BinomCI
#Đọc file afterlife và đặt tên af
af <- read_excel("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]
# Tổng tần số những người tin của biến Believe
k <- length(af2[af2 < 3])
# Tổng tần số của biến Believe (cỡ mẫu)
n <- length(af2)
# 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ế giới khác sau khi chết” với độ tin cậy 95% là: 3.6341672 <= p <= 0.6864677

5.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ài đặt gói DescTools.

B2. library(DescTools): 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ể; trong đó, 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 để minh họa

Để ướ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:

# Gọi backage DescTools để sử dụng hàm BinomCI
library(DescTools)
# Đọc file afterlife.xlsx vầ đặt tên af
af <- read_excel("afterlife.xlsx")
# Trích cuấ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 khi 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 khi chết"
k2 <- table(af1$Sex, af1$Believe)[2,1] + table(af1$Sex, af1$Believe) [2,2]
# Câu lệnh ước lượng kết quả
BinomDiffCI(k2, n2, k1, n1, conf.level = 0.95, method = "wald")
##           est.2      lwr.ci     upr.ci
## [1,] 0.03834843 -0.01499184 0.09168869

kết quả khoảng ước lượng sự khác biệt tỷ lệ Nam và Nữ tin “Tồn tại thế giới khác nhau sau khi chết” vưới độ tin cậy 95% là -0,01499184 <= p1−p2<=0,09168869

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

Để thực hiện kiểm định tính độc lập cho 2 biến định danh

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: chisqtest()

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àng điều tra ngẫu nhiên 1000 công dân, ta có số liệu sau:

Với mức ý nghĩa 5% dựa vào điều tra, 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ứuc 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ả thuyết 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 được kết quả sau

# Kiểm định sự độc lập của 2 biến định danh
ungho <- c(182, 213, 203) # Tạo vetor có tên ủng hộ
phandoi <- c(154, 138, 110) # Tạo vetor có tên phản đối
data <- as.data.frame(rbind(ungho, phandoi)) # Tạo bảng tần số ngẫu nhiên 2x3
names(data) <- c("Thap", "Trung bình", "Cao") # Đặt tên cho bảng ngẫu nhiên
# Bảng ngẫu nhiên
data
##         Thap Trung bình Cao
## ungho    182        213 203
## phandoi  154        138 110
chisq.test(data) # câu lệnh thực hiện kiểm định
## 
##  Pearson's Chi-squared test
## 
## data:  data
## X-squared = 7.8782, df = 2, p-value = 0.01947

Từ bảng kết quả ta có p-value =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%

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

Để 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ụ: Theo dõi ngẫu nhiên một số sản phẩm về chất lượng: loại 1, loại2, loại 3 được sản xuất ở các ca: ngày, đêm, có kết quả sau:

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 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 độc lập cho X và Y được thực hiện như sau:

#Kiểm định sự độc lập cho 2 biến định tính có thang đo thứ bậc
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
## 
## Attaching package: 'vcdExtra'
## The following object is masked from 'package:carData':
## 
##     Burt
cadem <- c(81, 15, 9) # Tạo vector có tên ca đêm
cangay <- c(118, 28, 10) # Tạo vetor có ten ca ngay
# Tạo bảng ngẫu nhiên 2x3
data <- as.matrix(rbind(cadem, cangay))
# Bảng ngẫu nhiên
data
##        [,1] [,2] [,3]
## cadem    81   15    9
## cangay  118   28   10
# 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 trên có p-value=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%.

Phụ lục 3. Hồi quy bằng phần mềm R và Eviews

1. Hồi quy theo phương pháp OLS

**1.1 Hồi quy tuyến tính đơn giản

Chúng ta xét mô hình hồi quy 2 biến số: y=β1+β2x+u

Trong phần này, chúng ta lấy bộ dữ liệu CPS1988.csv làm ví dụ minh họa cho nội dung trên.

Trước hết chúng ta đánh giá mối quan hệ giữa biến education và wage

Từ số liệu ta nhận thấy khi education tăng lên thì wage cũng tăng, ta có đồ thị phân tán của wage theo education bằng câu lệnh sau:

plot(wage~education, pch = 16): đồ thị phân tác của wage theo education

# Gọi biến weight trong CPS1988 và đặt tên weight
wage <- CPS1988$wage
# Gọi biến color trong CPS1988 và đặt tên color
education <- CPS1988$education
plot(wage~education, pch = 16)

plot
## function (x, y, ...) 
## UseMethod("plot")
## <bytecode: 0x000001f71847d5e8>
## <environment: namespace:base>

Dựa vào đồ thị phân tán ta cũng có kết quả như nhận xét ban đầu

Tuy nhiên để đánh giá chính xác mối quan hệ tuyến tính giữa wage và education ta dùng hệ số tương quan để đánh giá.

Hệ số tương quan giữa wage và education được tính bằng hàm câu lệnh sau:

cor(wage,education): câu lệnh tính hệ số tương quan giữa wage và education

cor(wage, education)
## [1] 0.301644

Kết quả câu lệnh cho thấy hệ số tương quan giữa wage và education là 0,301644.

Chúng ta có thể kiểm định giả thuyết tương quan giữa wage và education (tức là hệ số tương quan giữa wage và education bằng 0) với hàm cor.test() có sẵn trong R

cor.test(wage, education)
## 
##  Pearson's product-moment correlation
## 
## data:  wage and education
## t = 53.085, df = 28153, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2909884 0.3122247
## sample estimates:
##      cor 
## 0.301644

Ta có p-value < 2.2e-16 < 0,05, ta bác bỏ giả thuyết hệ số tương quan giữa wage và education bằng 0.

Để ước lượng mô hình hồi quy đơn bằng phần mềm R, chúng ta thực hiện câu lệnh sau:

reg<-lm(y~x,data=CPS1988): câu lệnh hồi quy tuyến tính y theo x và đặt tên reg; data=CPS1988: dữ liệu chính là đối tượng CPS1988

Hàm lm(viết tắt từu linear model) trong R có thể tính toán các giá trị của β1^ và β2^, cũng như σ̂ 2 một cách nhanh gọn.

# Hồi quy tuyến tính y theo x và đặt tên reg
reg<- lm(wage~education, data=CPS1988)
reg #kết quả hồi quy
## 
## Call:
## lm(formula = wage ~ education, data = CPS1988)
## 
## Coefficients:
## (Intercept)    education  
##      -12.83        47.18

Trong lệnh lm(wage~education) cho ta thấy wage là một hàm của education. Kết quả tính toán của lm cho thấy β1 = -12,83 và β2 = 47,18. Nói cách khác với hai tham số này chúng ta có thể tính wage cho bất kỳ education nào trong mẫu dựa vào phương trình tuyến tính

wage=−12,83+47,18.education

Thật ra hàm lm còn cung cấp cho chúng ta nhiều thông tin khác và các thông tin đó được lưu trong đối tượng reg. Để xem các thông tin đó chúng ta thực hiện

summary(reg) #Hiển thị chi tiết hồi quy
## 
## Call:
## lm(formula = wage ~ education, data = CPS1988)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
##  -786.0  -264.9   -54.8   174.5 18035.1 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -12.8281    11.8970  -1.078    0.281    
## education    47.1810     0.8888  53.085   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 432.4 on 28153 degrees of freedom
## Multiple R-squared:  0.09099,    Adjusted R-squared:  0.09096 
## F-statistic:  2818 on 1 and 28153 DF,  p-value: < 2.2e-16

Lệnh summary(reg) cầu R liệt kê các tính toán trong gồm các phần sau:

  1. Phần một mô tả chi tiết phần dư (Residuals) của hàm hồi quy

  2. Phần hai trình bày ước số của b1 và B2cùng với sai số chuẩn và giá trị của kiểm định t.

Giá trị kiểm định t cho B2là 53,085 với trị số B2, cho thấy B2 không phải bằng 0.

  1. Cột Coefficients: tên biến

  2. Cột Estimate: giá trị của các hệ số hồi quy

  3. Cột Std. Error: độ lệch chuẩn của các hệ số hồi quy

  4. Cột t value giá trị thống kê t, tiêu chuẩn kiểm định của bài toán kiểm định giả thuyết: B1= 0

  5. Cột Pr(>|t|): p_value của bài toán kiểm định giả thuyết:

  6. Cột Pr(>|t|): p_value của bài toán kiểm định giả thuyết: \(𝛽_i\) = 0

  7. Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

Dấu hiệu cho biết các biến độc lập có nghĩa hay không, với các mức ý nghĩa cho trước: ’*’: 1 ‰; ’’: 1%; “0.05 “: 5%; “0.1”: 1%.

  1. Phần ba của kết quả cho chúng ta thông tin về phương sai của phần dư (residual mean square). Ở đây, \(𝜎^2\) = 432.4. Trong kết quả này còn có kiểm định F, cũng chỉ là một kiểm định xem có quả thật \(𝛽_2\) = 0, tức có ý nghĩa tương tự như kiểm định t trong phần trên. Nói chung, trong trường hợp phân tích hồi qui tuyến tính đơn giản (với một yếu tố) chúng ta không cần phải quan tâm đến kiểm định F.

Ngoài ra, phần ba còn cho chúng ta một thông tin quan trọng, đó là trị số \(R^2\) hay hệ số xác định bội (coefficient of determination). Trị số R2 trong ví dụ này là 0.09099. Một hệ số cũng cần đề cập ở đây là hệ số điều chỉnh xác định bội (mà trong kết quả trên R gọi là “Adjusted R-squared”).

Với lệnh fitted() chúng ta có thể tính toán \(𝑦_i\) cho từng cá nhân như sau:

fitted(reg)

** 1.2.Hồi quy tuyến tính đa biến**

Chúng ta xét mô hình hồi quy k biến số:

\(𝑦 = 𝛽_1 + 𝛽_2𝑥_2 + ⋯ + 𝛽_i𝑥_i + ⋯𝛽_k𝑥_k + u\)

Để thực hiện hồi quy mô hình này ta thực hiện câu lệnh sau:

hoiquy <- lm(y~x2 + x3 + xi + … + xk): câu lệnh hồi quy tuyến tính y theo các biến x2, x3, …, xk.

Chúng ta tiếp tục sử dụng bộ số liệu CPS1988.csv để minh họa cho hồi quy đa biến. Ta thực hiện hồi quy wage theo education và experience ta được kết quả như sau:

reg <- lm(wage~education + experience, data = CPS1988) summary(reg)

Kết quả hàm hồi quy mẫu của wage theo education và experience:

\(wage\) = −385,0834 + 60,8964. education + 10,6057.experience

2. Hồi quy theo phương pháp ML

2.1 Mô hình Logitics

2.1.1 Ước lượng mô hình và dự báo

Giả sử chúng ta cần ước lượng mô hình sau:

Việc ước lượng các hệ số của mô hình logit được thực hiện bằng phương pháp ước lượng hợp lí cực đại (maximum likelihood). Với R chúng ta có thể thực hiện ước lượng này lệnh glm (viết tắt của generalized linear models – các mô hình tuyến tính tổng quát) như sau:

logit <-glm(y~x1 + x2+ …+ xk, family = binomial(link = “logit”), data = object): câu lệnh hồi quy logit của hàm y theo x1, x2, …, xk và đặt tên logit.

Chúng ta lấy dữ liệu từ tập tin crab.csv minh họa cho hồi quy logit.

Trước tiên chúng ta mô tả các biến trong file crab.

# Đọc file crab.csv và đặt tên crab
crab <- read_excel("D:/PhanmemR/Crab.xlsx")

Trong bộ dữ liệu gồm các biến sau:

sat: số vệ tinh xung quanh cua móng ngựa, nhận các giá trị từ 0 đến 15

𝑦: biến nhị phân nhận giá trị 1 nếu 𝑠𝑎𝑡 >= 1, nhận giá trị 0 nếu 𝑠𝑎𝑡 = 0

weight: trọng lượng cua móng ngựa

width: độ rộng mai cua móng ngựa

color: màu sắc cua móng ngựa

spine: số gai của cua móng ngựa

Để hồi quy mô hình Logit của biến 𝑦 theo các biến weight, width, color, spine bằng R ta thực hiện các câu lệnh sau:

crab <- read.csv(“crab.csv” , header = TRUE): câu lệnh đọc file crab.csv và đặt tên crab.

logit <-glm(y~x1 + x2+ …+ xk, family = binomial(link = “logit”), data = object): câu lệnh hồi quy mô hình Logit của biến y theo weight, width, color, spine và đặt tên logit; data = crab: dữ liệu chính là đối tượng crab.

summary(logit): câu lệnh hiển thị chi tiết kết quả hồi quy mô hình logit.

logit <-glm(y~weight+width+color+spine, family = binomial(link = "logit"), data = crab)
summary(logit)
## 
## Call:
## glm(formula = y ~ weight + width + color + spine, family = binomial(link = "logit"), 
##     data = crab)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1928  -0.9465   0.4919   0.8682   1.9892  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -7.0078436  3.8034845  -1.842   0.0654 .
## weight       0.0007949  0.0006917   1.149   0.2505  
## width        0.2732525  0.1893304   1.443   0.1489  
## color       -0.5915286  0.2417241  -2.447   0.0144 *
## spine        0.2716837  0.2410407   1.127   0.2597  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 225.76  on 172  degrees of freedom
## Residual deviance: 186.66  on 168  degrees of freedom
## AIC: 196.66
## 
## Number of Fisher Scoring iterations: 4

Từ bảng kết quả ta có mô hình hồi quy:

\(𝑙𝑜𝑔𝑖𝑡[𝜋°(𝑥)] = −7.5994 + 0,2733 . 𝑤𝑖𝑑𝑡ℎ + 0,7949. 𝑤𝑒𝑖𝑔ℎ𝑡 − 0,5915. 𝑐𝑜𝑙𝑜𝑟 + 0,2717. 𝑠𝑝𝑖𝑛𝑒\)

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

B1. Nhập số liệu cho các biến Y, X1,…, Xm, trong đó biến đáp ứng Y là nhị phân (dữ liệu được gán 0 hoặc 1).

B2. H.quy ước lượng cho (1):

Trong ví dụ trên, Eviews cho kết quả]

Chúng ta có thể diễn giải kết quả trên như sau. Ta có hệ số của biến 𝑤𝑖𝑑𝑡ℎ là 0,2733, có nghĩa là độ rộng tăng 1 đơn vị thì tỷ lệ cược (Odds) được nhân lên \(𝑒^{0.2733}\) = 1,31423. Tương tự cho các hệ số của các biến còn lại, Tương tự, chúng ta có thể tính Odds cho các biến còn lại bằng câu lệnh sau:

exp(logit$coefficients): câu lệnh tính e mũ cho các hệ số hồi quy mẫu trong mô hình.

exp(logit$coef)
##  (Intercept)       weight        width        color        spine 
## 0.0009047575 1.0007952608 1.3142320426 0.5534806048 1.3121719135

Kết quả cho chúng ta các giá trị của \(𝑒Zf\) và \(𝑒x2\); chẳng hạn \(𝑒x% = 2,2143187226\).

Ngoài ra, dựa trên mô hình hồi quy logit ở trên chúng ta có thể dự báo xác suất có vệ tinh của con cua cái với 𝑤𝑖𝑑𝑡ℎ = 25, 𝑤𝑒𝑖𝑔ℎ𝑡 = 3, 𝑐𝑜𝑙𝑜𝑟 = 5 𝑣à 𝑠𝑝𝑖𝑛𝑒 = 3 như sau bằng câu lệnh sau:

dubao <- 1/(1+(exp(-logit\(coef[1]-logit\)coef[2]25-logit\(coef[3]*3- logit\)coef[4]5-logit$coef[5]*3))): câu lệnh thực hiện dự báo và đặt tên dubao

dubao: câu lệnh xem kết quả dự báo.

Kết quả sau khi thự hiện câu lệnh

1/(1+(exp(-logit$coef[1]-logit$coef[2]*25-logit$coef[3]*3- logit$coef[4]*5-logit$coef[5]*3)))
##  (Intercept) 
## 0.0002457858

2.1.2. Kiểm định cho mô hình logistic

2.1.2.1.Kiểm định sự ảnh hưởng của biến dự báo (giải thích) lên biến đáp ứng bằng kiểm định Wald

Thống kê Wald: \(𝑍^2 = (\frac{\beta}{se(\beta)})\) (có phân phối xấp xỉ Chi – bình phương khi mẫu lớn). Kiểm định này gọi là kiểm định Wald.

Để thực hiện kiểm định Wald bằng phần mềm R ta làm như sau: B1. Cài đặt và gọi các packages sau:

library(mdscore)
## Warning: package 'mdscore' was built under R version 4.2.3
## Loading required package: MASS
## Warning: package 'MASS' was built under R version 4.2.3

B2. Hồi quy mô hình logistic.

B3. Thực hiện kiểm định Wald bằng câu lệnh: wald.test(model , 𝑡𝑒𝑟𝑚); trong đó model (tên mô hình hồi quy), terms (số tham số trong mô hình)

Chúng ta sử dụng tập dữ liệu crab.csv tiến hành kiểm định Wald bằng R, ta được kết quả như sau:

# Đọc tập tin crab.csv và đặt ten crab
Crab <- read_excel("D:/PhanmemR/Crab.xlsx")
# Hồi quy y theo width
logit <- glm(y~width, family = binomial(link="logit"), data = crab)
wald.test(logit,terms = 2)
## $W
## [1] 23.88748
## 
## $pvalue
## [1] 1.02134e-06
## 
## attr(,"class")
## [1] "wald.test"

Kết quả kiểm định Wald, giá trị thống kê \(𝑍^2\)= 23,88748, 𝑝 − 𝑣𝑎𝑙𝑢𝑒 = 1,02134 × \(10^{-6}\) < 0,05. Kết luận: với mức ý nghĩa 5%, sự xuất hiện của vệ tinh thực sự phụ thuộc vào chiều rộng của mai của con cua cái.

2.1.2.2. Kiểm định sự phù hợp của mô hình

Đối với mô hình hồi quy logistic, để kiểm định giả thuyết \(𝐻_0\), người ta thường dùng thống kê tỷ số hợp lý (LR: \(𝑃𝑠𝑒𝑢𝑑𝑜 − 𝑅^2)\)

\[ 𝐿𝑅 = −2(𝐿_0 − 𝐿) \]

Dựa trên bộ dữ liệu crab.csv ta hồi quy logitics y theo with, bài toán kiểm định giả thuyết : “Mô hình không phù hợp với dữ liệu điều tra” được thực hiện trên R như sau

Crab <- read_excel("D:/PhanmemR/Crab.xlsx")
ogit <- glm(y~width, family = binomial(link="logit"), data = crab)
Anova(logit)
## Analysis of Deviance Table (Type II tests)
## 
## Response: y
##       LR Chisq Df Pr(>Chisq)    
## width   31.306  1  2.204e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ta có kết quả: Thống kê LR= 31,306 và $ p_{value} = 2,204 × 10^{-8} < 0,05$ : Bác bỏ �, Kết luận, mô hình phù hợp với mức ý nghĩa 5%.