TH-PTDLĐT
Xuân Huyền
2023-07-19
1 Trích xuất dữ liệu
1.1 Trích xuất dữ liệu từ một dataframe có sẵn
setwd("C:/RRR")
library(readxl)
C <- read_excel("Crab.xlsx")
C## # A tibble: 173 × 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 0
## 7 2 1 26.5 0 2350 0
## 8 4 2 24.7 0 1900 0
## 9 3 1 23.7 0 1950 0
## 10 4 3 25.6 0 2150 0
## # ℹ 163 more rows1.2 Trích xuất dữ liệu theo dòng (hàng)
C <- read_excel("Crab.xlsx")
C1<-C[c(1, 3, 100),] # Trích xuất quan sát ở dòng 1, 3, 100
C1 # Xem 6 giá trị đầu tiên của C1## # A tibble: 3 × 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 dòng 3, dòng 1, dòng 100
C2 <- C[c(1:3, 10, 100),]
C2 # Xem 6 giá trị đầu tiên của C2## # A tibble: 5 × 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.3 Trích xuất dữ liệu theo cột
C3 <- C[c(2:4)] # Trích xuất từ cột 2 đến cột 4 và đặt tên MP3
head(C3) # Xem 6 giá trị đầu tiên của d3## # A tibble: 6 × 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.4 Trích xuất dữ liệu theo hàng và cột
C4 <- C[c(1,3,5,7), c(2:4)]
head(C4) # Xem 6 giá trị đầu tiên của C4## # A tibble: 4 × 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.5 Trích xuất dữ liệu từ một dataframe có sẵn bằng lệnh subset
C5 <- subset(C, y == "1")
dim(C5) # Số các quan sát của biến doanhthu = cao## [1] 111 6head(C5) # Xem 6 giá trị đầu tiên của d5## # A tibble: 6 × 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 1C6 <- subset(C, y == "1" & color == "5")
head(C6) # Xem 6 giá trị đầu tiên của C6## # A tibble: 6 × 6
## color spine width satell weight y
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 5 3 25.8 3 2000 1
## 2 5 3 24.7 5 2100 1
## 3 5 3 27 3 2250 1
## 4 5 3 29.3 12 3225 1
## 5 5 3 22.5 4 1475 1
## 6 5 3 25.8 10 2250 1C7 <- subset(C, y == "1" | weight > 2600)
head(C7) # Xem 6 giá trị đầu tiên của C7## # A tibble: 6 × 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 25.8 0 2650 0
## 5 3 3 28.2 11 3050 1
## 6 3 1 26 14 2300 11.6 Mã hóa dữ liệu
setwd("C:/RRR")
library(xlsx)## Warning: package 'xlsx' was built under R version 4.2.3# Đọc file afterlife.xlsx và đặt tên af
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# Chỉ giữ lại những quan sát mà giá trị biến Belive < 5 và đặt tên afmahoa
afmahoa <- subset(af, af$Believe < 5)
dim(afmahoa) # Số quan sát của afmahoa## [1] 1260 2afmahoa$Believe[afmahoa$Believe <= 2] <- 1
# Trong cột biến Belive của afmahoa thay những giá trị <= 2 bằng 1
afmahoa$Believe[afmahoa$Believe > 2] <- 0
# Trong cột biến Belive của afmahoa thay những giá trị > 2 bằng 0
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 12 Thực hành bảng ngẫu nhiên và thống kê suy diễn trên R
2.1 Bảng ngẫu nhiên
2.1.1 Bảng tần số và tần suất
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# 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# Lập bảng tần số cho biến Belive và kết quả
table(afmahoa$Believe)##
## 0 1
## 428 832# Lập bảng tần suất cho biến Belive và kết quả
table(afmahoa$Believe)/sum(table(afmahoa$Believe))##
## 0 1
## 0.3396825 0.66031752.1.2 Bảng 2 chiều
table(afmahoa$Sex, afmahoa$Believe)##
## 0 1
## 1 188 331
## 2 240 5012.2 Tính OddsRatio
library(DescTools) # Gọi packages DescTools## Warning: package 'DescTools' was built under R version 4.2.3v <- c(189, 104, 10845, 10933) # Tạo vector v
x <- matrix (v, 2, 2) # Tạo ma trận cấp 2 đặt tên x
y <- c(189, 104) # Vector 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.3307803 Một số bài toán ước lượng
3.1 Ước lượng tỷ lệ một tổng thể
library(DescTools) # Gọi packages DescTools để sử dụng hàm BinomCI
# Trích xuất biến Belive 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 Belive 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 Belive
n <- length(af2) # Tổng tần số của biến Belive (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.68646773.2 Ước lượng sự khác biệt tỷ lệ hai tổng thể
# Gọi packages DescTools để sử dụng hàm BinomDiffCI
library(DescTools)
# Đọc file afterlife.xlsx và đặt tên af
af <- read_excel("afterlife.xlsx")
# Trích xuất những quan sát mà giá trị nhỏ hơn 5 của biến Belive
af1<-subset(af, Believe<5)
table(af1)## Believe
## Sex 1 2 3 4
## 1 79 252 118 70
## 2 79 422 160 80n1<-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 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.091688694 Kiểm định tính độc lập
4.1 Kiểm định tính độc lập cho hai biến định danh
ungho<-c(182, 212, 203) # Tạo vector có tên ủng hộ
phandoi<- c(154, 138, 110) # Tạo vector có tên phản đối
data1 <- as.data.frame(rbind(ungho, phandoi)) # Tạo bảng ngẫu nhiên 2x3
names(data1) <- c('Thap','Trung binh','cao') # Đặt tên cho bảng ngẫu nhiên
data1 # Bảng ngẫu nhiên## Thap Trung binh cao
## ungho 182 212 203
## phandoi 154 138 110chisq.test(data1) # Câu lệnh thực hiện kiểm định##
## Pearson's Chi-squared test
##
## data: data1
## X-squared = 7.8476, df = 2, p-value = 0.019774.2 Kiểm định tính độc lập cho 2 biến định tính có thang đo thứ bậc (CMH test)
install.packages (“vcdExtra”)
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.3cadem<-c(81,15,9) # Tạo vector có tên ca đêm
cangay<-c(118, 28, 10) # Tạo vector có tên ca ngày
# Tạo bảng ngẫu nhiên 2x3
data2 <- as.matrix(rbind(cadem, cangay)) # Bảng ngẫu nhiên
# Câu lệnh kiểm định X=(0,1), Y(97.5, 24.5, 5.5)
CMHtest(data2, rscores = c(0,1), cscores = 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.628525 Hồi quy theo phương pháp OLS
5.1 Hồi quy tuyến tính đơn giản
# Đọc file Crab và đặt tên C
C <- read_excel("Crab.xlsx")
# Gọi biến weight trong C và đặt tên weight
weight<-C$weight
# Gọi biến width trong C và đặt tên width
width <- C$width
plot(weight~width, pch = 16)
plot## function (x, y, ...)
## UseMethod("plot")
## <bytecode: 0x000001a4604252c0>
## <environment: namespace:base># Câu lệnh tính hệ số tương quan giữa weight và width
cor(weight, width)## [1] 0.8868715# Kiểm định giả thuyết tương quan
cor.test(weight, width)##
## Pearson's product-moment correlation
##
## data: weight and width
## t = 25.102, df = 171, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8501665 0.9149979
## sample estimates:
## cor
## 0.8868715# Hồi quy tuyến tính y theo x và đặt tên regreg <- lm(weight~width, data = C)
reg # Kết quả hàm hồi quy##
## Call:
## lm(formula = weight ~ width, data = C)
##
## Coefficients:
## (Intercept) width
## -3944.0 242.6summary(reg) # Hiển thị chi tiết hồi quy##
## Call:
## lm(formula = weight ~ width, data = C)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1186.00 -116.15 -5.13 150.76 1015.50
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3944.019 255.027 -15.46 <2e-16 ***
## width 242.642 9.666 25.10 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.4 on 171 degrees of freedom
## Multiple R-squared: 0.7865, Adjusted R-squared: 0.7853
## F-statistic: 630.1 on 1 and 171 DF, p-value: < 2.2e-16fitted(reg) ## 1 2 3 4 5 6 7 8
## 2922.756 1515.431 2364.679 2073.508 2364.679 1830.866 2486.000 2049.244
## 9 10 11 12 13 14 15 16
## 1806.602 2267.622 1952.187 2316.150 2898.491 1151.468 2364.679 2631.585
## 17 18 19 20 21 22 23 24
## 2170.565 3092.605 2049.244 2704.378 1685.280 2122.036 1515.431 2534.528
## 25 26 27 28 29 30 31 32
## 2316.150 2413.207 3019.813 2558.792 2728.642 2097.772 3165.398 2316.150
## 33 34 35 36 37 38 39 40
## 2291.886 2291.886 2534.528 1806.602 2558.792 2728.642 1733.809 2825.699
## 41 42 43 44 45 46 47 48
## 2728.642 2388.943 2777.170 3335.247 2971.284 3068.341 2898.491 2122.036
## 49 50 51 52 53 54 55 56
## 2971.284 3408.040 2049.244 2777.170 2704.378 1612.488 2291.886 2922.756
## 57 58 59 60 61 62 63 64
## 2655.849 2413.207 2801.435 2243.358 2631.585 2000.715 2607.321 2364.679
## 65 66 67 68 69 70 71 72
## 2849.963 3335.247 3092.605 2413.207 2486.000 2413.207 2267.622 1636.752
## 73 74 75 76 77 78 79 80
## 1636.752 2219.093 1927.923 1612.488 2364.679 2219.093 2291.886 2146.301
## 81 82 83 84 85 86 87 88
## 2000.715 2728.642 1661.016 2340.414 2316.150 2607.321 2971.284 2243.358
## 89 90 91 92 93 94 95 96
## 1758.073 1879.394 3262.455 2558.792 2534.528 3019.813 1661.016 3092.605
## 97 98 99 100 101 102 103 104
## 2243.358 2486.000 2000.715 2971.284 2898.491 2000.715 2728.642 2049.244
## 105 106 107 108 109 110 111 112
## 2170.565 2680.113 2437.471 3092.605 2194.829 2486.000 2801.435 2607.321
## 113 114 115 116 117 118 119 120
## 2291.886 2122.036 3796.268 1806.602 3165.398 1394.110 2122.036 2607.321
## 121 122 123 124 125 126 127 128
## 1830.866 3383.776 2413.207 1927.923 2704.378 2219.093 2947.020 1515.431
## 129 130 131 132 133 134 135 136
## 2413.207 2097.772 2000.715 2146.301 2849.963 2316.150 2825.699 2097.772
## 137 138 139 140 141 142 143 144
## 2947.020 2655.849 2122.036 2728.642 4184.495 3456.569 3092.605 1952.187
## 145 146 147 148 149 150 151 152
## 2316.150 2122.036 3747.739 3213.926 1879.394 3335.247 2752.906 2413.207
## 153 154 155 156 157 158 159 160
## 1661.016 1612.488 2000.715 2049.244 2922.756 1855.130 1830.866 3286.719
## 161 162 163 164 165 166 167 168
## 2486.000 2364.679 2898.491 2291.886 2486.000 2316.150 1903.658 2413.207
## 169 170 171 172 173
## 2388.943 3092.605 2849.963 2607.321 2000.7155.2 Hồi quy tuyến tính đa biến
reg1 <- lm(weight~width + color, data = C)
reg1 # Kết quả hàm hồi quy##
## Call:
## lm(formula = weight ~ width + color, data = C)
##
## Coefficients:
## (Intercept) width color
## -3867.31 241.37 -12.61summary(reg1) # Hiển thị chi tiết hồi quy##
## Call:
## lm(formula = weight ~ width + color, data = C)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1191.28 -118.87 -8.52 148.38 1019.09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3867.31 301.93 -12.809 <2e-16 ***
## width 241.37 10.05 24.027 <2e-16 ***
## color -12.61 26.42 -0.477 0.634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268 on 170 degrees of freedom
## Multiple R-squared: 0.7868, Adjusted R-squared: 0.7843
## F-statistic: 313.7 on 2 and 170 DF, p-value: < 2.2e-166 Hồi quy theo phương pháp ML
6.1 Mô hình Logistic
6.1.1 Ước lượng mô hình và dự báo
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=binomail(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.
library(readxl)
# Đọc file crab.excel và đặt tên crab
C <- read_excel("Crab.xlsx")
head(C) # Xem 6 quan sát đầu## # A tibble: 6 × 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 0