THUCHANH 19/07
Võ Thị Kim Oanh
2023-07-19
Đường link: https://docs.google.com/document/d/1S29fHZiC2Bc-z-JQSQT4s7XTEB7_CZD8/edit
1 Trích xuất dữ liệu từ một dataframe có sẵn
Trong một số trường hợp chúng ta có thể không sử dụng hết các biến số hoặc tất cả các quan sát trong 1 file dữ liệu, chúng ta có thể trích dữ liệu từ 1 dataframe có sẵn bằng những câu lệnh đơn giản, trong nội dung này chúng ta sử dụng bộ số liệu CPS1988.csv để làm ví dụ minh họa. Trước hết chúng ta sơ bộ file dữ liệu có tên CPS1988.csv bằng câu 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.
CPS1988: câu lệnh cho biết chi tiết các biến của đối tượng CPS1988.
library(readxl)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# *Găn dữ liệu Crab cho biến CPS1988*
CPS1988 <- read_excel("D:/RStudio/Crab.xlsx", sheet = 1)
str(CPS1988)## tibble [173 x 6] (S3: tbl_df/tbl/data.frame)
## $ color : num [1:173] 3 4 2 4 4 3 2 4 3 4 ...
## $ spine : num [1:173] 3 3 1 3 3 3 1 2 1 3 ...
## $ width : num [1:173] 28.3 22.5 26 24.8 26 23.8 26.5 24.7 23.7 25.6 ...
## $ satell: num [1:173] 8 0 9 0 4 0 0 0 0 0 ...
## $ weight: num [1:173] 3050 1550 2300 2100 2600 2100 2350 1900 1950 2150 ...
## $ y : num [1:173] 1 0 1 0 1 0 0 0 0 0 ...1.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:
#Trích suất quan sát ở dòng 1,3,100
CPS1<-CPS1988[c(1, 3, 100),]
head(CPS1)## # A tibble: 3 x 6
## color spine width satell weight y
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3 3 28.3 8 3050 1
## 2 2 1 26 9 2300 1
## 3 4 3 28.5 1 3000 1#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),]
head(CPS2)## # A tibble: 5 x 6
## color spine width satell weight y
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3 3 28.3 8 3050 1
## 2 4 3 22.5 0 1550 0
## 3 2 1 26 9 2300 1
## 4 4 3 25.6 0 2150 0
## 5 4 3 28.5 1 3000 11.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.
#Trích xuất từ cột 2 đến cột 4 và gắn cho biến CPS3
CPS3<-CPS1988[c(2:4)]
head(CPS3)## # A tibble: 6 x 3
## spine width satell
## <dbl> <dbl> <dbl>
## 1 3 28.3 8
## 2 3 22.5 0
## 3 1 26 9
## 4 3 24.8 0
## 5 3 26 4
## 6 3 23.8 01.3 Trích xuất dữ liệu theo hàng và theo 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ục D:/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.
# Trích xuất dữ liệu từ cho các biến từ cột 2 đến 4 và lấy dòng 1,3 5 và 7 của đối tượng CPS1988, gắn cho biến CPS4.
CPS4 <-CPS1988[c(1, 3, 5, 7), c(2:4)]
head(CPS4)## # A tibble: 4 x 3
## spine width satell
## <dbl> <dbl> <dbl>
## 1 3 28.3 8
## 2 1 26 9
## 3 3 26 4
## 4 1 26.5 01.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 y = 1 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ên
CPS1988.csv trong thư mục D:/PhanmemR và đặt tên CPS1988.
CPS5 <- subset(CPS1988, CPS1988$y == 1): trích xuất các quan sát của biến y == 1 trong đối tượng CPS1988 và đặt tên CPS5.
y <- CPS1988$y
CPS5 <- subset(CPS1988, CPS1988$y == 1)
dim(CPS5) #Số các quan sát của biến y = 1## [1] 111 6head(CPS5)## # A tibble: 6 x 6
## color spine width satell weight y
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3 3 28.3 8 3050 1
## 2 2 1 26 9 2300 1
## 3 4 3 26 4 2600 1
## 4 3 3 28.2 11 3050 1
## 5 3 1 26 14 2300 1
## 6 2 1 27.1 8 2950 1Kế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 y = 1, và (2) có education ≥ 3 với đặt tên CPS6, ta thực hiện câu lệnh sau:
CPS6 <- subset(CPS1988, y == 1 & spine == 2): trích xuất các quan sát của biến y = 1 và các quan sát của biến spine = 2 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
# Trích xuất các quan sát của biến y = 1 và các quan sát của biến spine =2 trong đối tượng bike gắn vào với tên bike6
CPS6 <- subset(CPS1988, CPS1988$y == 1 & CPS1988$spine == 2)
dim(CPS6) #Số các quan sát của biến y = 1 và spine =2## [1] 7 6head (CPS6)## # A tibble: 6 x 6
## color spine width satell weight y
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3 2 23.2 4 1950 1
## 2 2 2 25 3 2300 1
## 3 4 2 26 3 2200 1
## 4 3 2 24.7 4 2550 1
## 5 2 2 24.5 6 1950 1
## 6 3 2 25 6 2250 1Hoặ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, y == 1 | width > 23): câu lệnh trích xuất các quan sát của biến y = 1 hoặc các quan sát của biến width > 23 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, CPS1988$y == 1 | CPS1988$width > 23)
dim(CPS7) #Số các quan sát của biến y = 1 và width >23## [1] 167 6head (CPS7)## # A tibble: 6 x 6
## color spine width satell weight y
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3 3 28.3 8 3050 1
## 2 2 1 26 9 2300 1
## 3 4 3 24.8 0 2100 0
## 4 4 3 26 4 2600 1
## 5 3 3 23.8 0 2100 0
## 6 2 1 26.5 0 2350 01.5 Mã hoá 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.
hoặc af <- read_excel(“D:/RStudio/afterlife.xlsx”, sheet = 1)
af: câu lệnh cho biết chi tiết đối tượng af.
# Đọc file afterlife.xlsx và đặt tên af.
af <- read_excel("D:/RStudio/afterlife.xlsx", sheet = 1)
head(af)## # A tibble: 6 x 2
## Sex Believe
## <dbl> <dbl>
## 1 2 2
## 2 2 4
## 3 1 4
## 4 2 4
## 5 2 6
## 6 1 6Dữ 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ế giớ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.
Dữ liệu bị khuyết.
Không thể chọn.
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 thay 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 thay 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) #Số quan sát của afmahoa## [1] 1260 2head (afmahoa)## # A tibble: 6 x 2
## Sex Believe
## <dbl> <dbl>
## 1 2 2
## 2 2 4
## 3 1 4
## 4 2 4
## 5 2 2
## 6 2 2#Thay những giá trị quan sát ≤ 2 trong cột Believe của afmahoa bằng 1
afmahoa$Believe[afmahoa$Believe <=2] <-1
#Thay những giá trị quan sát > 2 trong cột Believe của afmahoa bằng 0.
afmahoa$Believe[afmahoa$Believe >2] <- 0
head(afmahoa)## # A tibble: 6 x 2
## Sex Believe
## <dbl> <dbl>
## 1 2 1
## 2 2 0
## 3 1 0
## 4 2 0
## 5 2 1
## 6 2 1Sau 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.
2 Phụ lục 2. Thực hành bảng ngẫu nhiên và suy diễn thống kê trên R
2.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: 𝐵1, 𝐵2, … , 𝐵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 đó 𝑛ij là số lần quan sát được cặp thuộc tính 𝐴i , 𝐵j, còn gọi là tần số của 𝐴i , 𝐵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ụ.
2.2 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)## # A tibble: 6 x 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
table(afmahoa$Sex)##
## 1 2
## 519 741# Lập bảng tần suất cho biến Sex
table(afmahoa$Sex)/sum(table(afmahoa$Sex))##
## 1 2
## 0.4119048 0.5880952# Lập bảng tần số cho biến Sex
table(afmahoa$Believe)##
## 0 1
## 428 832# Lập bảng tần suất cho biến Sex
table(afmahoa$Believe)/sum(table(afmahoa$Believe))##
## 0 1
## 0.3396825 0.66031752.3 Bảng 2 chiều
table(afmahoa$Sex, afmahoa$Believe)##
## 0 1
## 1 188 331
## 2 240 5012.4 Tính OddsRatio
#Gọi packages DescTools
library(DescTools)
#Tạo vecto v
v <-c(189, 104, 10845, 10933)
#Tạo ma trận cấp 2 đặt tên x
x <-matrix(v, 2,2)
# Vecto y là những người nhồi máu cơ tim
y<-c(189, 104)
#Tính Oddsratio
OddsRatio(x, conf.level = 0.95, y = NULL, method = "wald")## odds ratio lwr.ci upr.ci
## 1.832054 1.440042 2.3307802.5 Một số bài toán ước lượng
2.5.1 Ước lượng tỷ lệ một tổng thể
Công thức ước lượng
Để 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:
# Gọi pakages DescTools để sử dụng hàm BinomCI
library(DescTools)
# Đọc file afterlife.xlsx và đặt tên af.
af <- read_excel("D:/RStudio/afterlife.xlsx", sheet = 1)
#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ố 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.68646772.6 Ước lượng sự khác biệt tỷ lệ hai tổng thể
# Gọi pakages DescTools để sử dụng hàm BinomCI
library(DescTools)
# Đọc file afterlife.xlsx và đặt tên af.
af <- read_excel("D:/RStudio/afterlife.xlsx", sheet = 1)
# Trích xuất những quan sát mà giá trị nhỏ hơn 5 của biến Believe
af3 <- subset(af, Believe <5)
table(af3)## Believe
## Sex 1 2 3 4
## 1 79 252 118 70
## 2 79 422 160 80# Tổng tần số giới tính nam
n1 <- table(af3$Sex)[1]
# Tổng tần số giới tính nữ
n2 <- table(af3$Sex)[2]
# Tổng tần số người nam "Tin có thế giới khác sau khi chết"
k1 <- table(af3$Sex, af3$Believe)[1,1] + table(af3$Sex, af3$Believe)[1,2]
# Tổng tần số người nữ "Tin có thế giới khác sau khi chết"
k2 <- table(af3$Sex, af3$Believe)[2,1] + table(af3$Sex, af3$Believe)[2,2]
# Câu lệnh ước lượng và 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.091688692.7 Kiểm định tính độc lập
2.7.1 Kiểm định tính độc lập cho hai biến định danh
# Kiểm định sự độc lập của 2 biến định danh
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
data1 <- as.data.frame(rbind(ungho, phandoi)) # Tạo bảng ngẫu nhiên 2x3
# Đặt tên cho bảng ngẫu nhiên
names(data1) <- c( "Thap", "Trung bình", "Cao")
data1 # Bảng ngẫu nhiên## Thap Trung bình Cao
## ungho 182 213 203
## phandoi 154 138 110chisq.test(data1)##
## Pearson's Chi-squared test
##
## data: data1
## X-squared = 7.8782, df = 2, p-value = 0.019472.7.2 Kiểm định tính độc lập cho hai biến định tính có thang đo thứ bậc (CMH test)
#KIEM DINH TINH DOC LAP CUA 2 BIEN DINH TINH COS THANG DO THU BAC
cadem <- c(81,15,9) # Tạo vecto có tên cadem
cangay <- c(118,28,10) # Tạo vecto có tên cangay
# Tạo bảng ngẫu nhiên 2 x 3
data2 <- as.matrix(rbind(cadem, cangay))
data2 # Bảng ngẫu nhiên## [,1] [,2] [,3]
## cadem 81 15 9
## cangay 118 28 103 Phụ lục 3. Hồi quy bằng phần mềm R và Eviews
3.1 Hồi quy theo phương pháp OLS
3.1.1 Hồi quy tuyến tính đơn giản
# *Găn dữ liệu Crab cho biến CPS1988*
CPS1988 <- read_excel("D:/RStudio/Crab.xlsx", sheet = 1)
#Gọi biến weight trong CPS1988 và đặt tên weight
weight <- CPS1988$weight
#Gọi biến color trong CPS1988 và đặt tên color
color <- CPS1988$color
plot(weight~color, pch = 16)
plot## function (x, y, ...)
## UseMethod("plot")
## <bytecode: 0x000001cb460bf0c8>
## <environment: namespace:base>cor (weight, color)## [1] -0.2507772cor.test(weight, color)##
## Pearson's product-moment correlation
##
## data: weight and color
## t = -3.3876, df = 171, p-value = 0.0008751
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3855516 -0.1055253
## sample estimates:
## cor
## -0.2507772#Hồi quy tuyến tính biến weight theo color và đặt tên reg
reg <- lm(weight ~ color, data = CPS1988)
reg #Kết quả hồi quy##
## Call:
## lm(formula = weight ~ color, data = CPS1988)
##
## Coefficients:
## (Intercept) color
## 3057.8 -180.4# Hiển thị chi tiết hồi quy
summary(reg)##
## Call:
## lm(formula = weight ~ color, data = CPS1988)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1216.46 -416.46 -86.02 413.98 2683.54
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3057.80 188.09 16.257 < 2e-16 ***
## color -180.44 53.27 -3.388 0.000875 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 560.2 on 171 degrees of freedom
## Multiple R-squared: 0.06289, Adjusted R-squared: 0.05741
## F-statistic: 11.48 on 1 and 171 DF, p-value: 0.0008751fitted(reg)## 1 2 3 4 5 6 7 8
## 2516.461 2336.016 2696.906 2336.016 2336.016 2516.461 2696.906 2336.016
## 9 10 11 12 13 14 15 16
## 2516.461 2336.016 2336.016 2516.461 2516.461 2155.572 2516.461 2696.906
## 17 18 19 20 21 22 23 24
## 2516.461 2516.461 2155.572 2516.461 2516.461 2696.906 2516.461 2336.016
## 25 26 27 28 29 30 31 32
## 2155.572 2155.572 2516.461 2516.461 2155.572 2516.461 2696.906 2696.906
## 33 34 35 36 37 38 39 40
## 2516.461 2516.461 2516.461 2155.572 2516.461 2516.461 2155.572 2516.461
## 41 42 43 44 45 46 47 48
## 2336.016 2696.906 2696.906 2516.461 2336.016 2336.016 2516.461 2516.461
## 49 50 51 52 53 54 55 56
## 2516.461 2516.461 2155.572 2516.461 2696.906 2516.461 2516.461 2516.461
## 57 58 59 60 61 62 63 64
## 2516.461 2336.016 2516.461 2155.572 2336.016 2336.016 2336.016 2516.461
## 65 66 67 68 69 70 71 72
## 2516.461 2516.461 2516.461 2516.461 2516.461 2516.461 2336.016 2336.016
## 73 74 75 76 77 78 79 80
## 2336.016 2516.461 2336.016 2516.461 2336.016 2516.461 2336.016 2516.461
## 81 82 83 84 85 86 87 88
## 2336.016 2155.572 2336.016 2336.016 2516.461 2155.572 2516.461 2155.572
## 89 90 91 92 93 94 95 96
## 2155.572 2516.461 2516.461 2516.461 2155.572 2516.461 2336.016 2516.461
## 97 98 99 100 101 102 103 104
## 2336.016 2336.016 2336.016 2336.016 2516.461 2516.461 2516.461 2516.461
## 105 106 107 108 109 110 111 112
## 2516.461 2336.016 2516.461 2516.461 2516.461 2516.461 2516.461 2516.461
## 113 114 115 116 117 118 119 120
## 2336.016 2516.461 2516.461 2155.572 2155.572 2336.016 2516.461 2336.016
## 121 122 123 124 125 126 127 128
## 2336.016 2696.906 2336.016 2516.461 2516.461 2516.461 2336.016 2155.572
## 129 130 131 132 133 134 135 136
## 2516.461 2516.461 2696.906 2516.461 2516.461 2155.572 2516.461 2516.461
## 137 138 139 140 141 142 143 144
## 2516.461 2336.016 2516.461 2516.461 2516.461 2516.461 2336.016 2516.461
## 145 146 147 148 149 150 151 152
## 2516.461 2155.572 2516.461 2516.461 2336.016 2516.461 2516.461 2516.461
## 153 154 155 156 157 158 159 160
## 2516.461 2516.461 2155.572 2516.461 2516.461 2516.461 2336.016 2336.016
## 161 162 163 164 165 166 167 168
## 2516.461 2516.461 2516.461 2155.572 2516.461 2516.461 2336.016 2336.016
## 169 170 171 172 173
## 2336.016 2336.016 2696.906 2155.572 2516.461#Hồi quy tuyến tính biến weight theo color và satell và đặt tên reg
reg1 <- lm(weight ~ color + satell, data = CPS1988)
summary(reg1) #Hiển thị chi tiết hồi quy##
## Call:
## lm(formula = weight ~ color + satell, data = CPS1988)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1017.90 -367.90 -33.24 332.10 2454.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2721.87 191.45 14.217 < 2e-16 ***
## color -134.66 51.22 -2.629 0.00934 **
## satell 61.13 13.05 4.686 5.68e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 528.7 on 170 degrees of freedom
## Multiple R-squared: 0.1701, Adjusted R-squared: 0.1603
## F-statistic: 17.42 on 2 and 170 DF, p-value: 1.311e-073.2 Hồi quy theo phương pháp ML
3.2.1 Mô hình Logitics
3.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:
𝑙𝑜𝑔 �1 −𝜋(𝜋𝑥()𝑥)Æ = 𝛼$ + 𝛽x + 𝛽%𝑥% + ⋯ + 𝛽!𝑥!
Trong đó E(H) #E(H) được gọi là khả năng thành công và hàm 𝑙𝑜𝑔 ç #EE(H()H)é được gọi là hàm logit của 𝜋(𝑥) .
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.
# *Găn dữ liệu Crab cho biến CPS1988*
crab <- read_excel("D:/RStudio/Crab.xlsx", sheet = 1)
head(crab) # Xem 6 quan sát đầu## # A tibble: 6 x 6
## color spine width satell weight y
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3 3 28.3 8 3050 1
## 2 4 3 22.5 0 1550 0
## 3 2 1 26 9 2300 1
## 4 4 3 24.8 0 2100 0
## 5 4 3 26 4 2600 1
## 6 3 3 23.8 0 2100 0Trong 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
# Câu lệnh hồi quy mô hình logit của biến y theo width,weight, color VÀ spine
width <- crab$width
spine <- crab$spine
logit <-glm(crab$y~ width + weight + color + spine, family = binomial (link = "logit"), data = crab)
summary(logit)##
## Call:
## glm(formula = crab$y ~ width + weight + color + spine, family = binomial(link = "logit"),
## data = crab)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.0078436 3.8034845 -1.842 0.0654 .
## width 0.2732525 0.1893304 1.443 0.1489
## weight 0.0007949 0.0006917 1.149 0.2505
## 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: 4Từ bảng kết quả ta có mô hình hồi quy:
logit[π°(x)]=−7.5994+0,2733.width+0,7949.weight−0,5915.color+0,2717.spine
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): → 𝑄𝑢𝑖𝑐𝑘 → 𝐸𝑠𝑡𝑖𝑚𝑎𝑡𝑒 𝑒𝑞𝑢𝑎𝑡𝑖𝑜𝑛 → (𝑔õ: 𝑦 𝑐 𝑥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 e0,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.
# câu lệnh tính e mũ cho các hệ số hồi quy mẫu trong mô hình
exp(logit$coefficients)## (Intercept) width weight color spine
## 0.0009047575 1.3142320426 1.0007952608 0.5534806048 1.3121719135# Câu lệnh thực hiện dự báo
dubao <- 1/(1+(exp(-logit$coef[1]-logit$coef[2]*25-logit$coef[3]*5-
logit$coef[4]*3)))
dubao #Kết quả dự báo## (Intercept)
## 0.12486533.2.1.2 Kiểm định cho mô hình logistic
3.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: Z (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:
# Cài đặt các gói pakages và gọi lệnh:
library(mdscore)## Warning: package 'mdscore' was built under R version 4.3.1## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectwald.test(logit, terms = 2)## $W
## [1] 2.082992
##
## $pvalue
## [1] 0.148948
##
## attr(,"class")
## [1] "wald.test"3.2.1.2.2 Kiểm định sự phù hợp của mô hình
crab <- read_excel("D:/RStudio/Crab.xlsx", sheet = 1) # Đọc tập tin
logit <-glm(crab$y~ width + weight + color + spine, family = binomial (link = "logit"), data = crab) # Hồi quy
anova(logit) # Kiểm định thống kê LR và kết quả## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: crab$y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 172 225.76
## width 1 31.3059 171 194.45
## weight 1 1.5608 170 192.89
## color 1 4.9669 169 187.93
## spine 1 1.2644 168 186.66