https://drive.google.com/file/d/14PaYgcdovuM76lnAjnB0PaaloArnznMl/view?usp=share_link

Đọc file CPS1988 và đặt tên CPS1988

# Đọc file CPS1988 và đặt tên CPS1988
CPS1 <- CPS1988[c(1, 3, 100),] # Trích xuất dữ liệu quan sát ở dòng 1, 3, 100.
CPS1    # Xem 6 giá trị đầ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 giá trị đầ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

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ừ cột 2 đến cột 4 và đặt tên CPS3
head(CPS3) # Xem 6 giá trị đầ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

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 giá trị đầ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ừ dataframe gốc ban đầu.

Trích xuất dữ liệu từ ộ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ộttiê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ố quan sát của biến parttime = yes
## [1] 2524    7
head(CPS5)  # Xem 6 giá trị đầu tiên của CPS4
##      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 giá trị đầ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 giá trị đầ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

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.

library(readxl)
# Đọc file afterlife.xlsx và đặt tên af
af <- read_excel("~/Documents/D/PhanmemR/afterlife.xlsx",  sheet = "Sheet1")
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
## 
## 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)
## Loading required package: vcd
## Loading required package: grid
## Loading required package: gnm
## 
## 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

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

7.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: 0x7ff6d811ec18>
## <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.