Bai tap 26/07/2023
Nguyễn Thị Yến Nhi
2023-07-26
require(tidyverse)## Loading required package: tidyverse## Warning: package 'tibble' was built under R version 4.2.3## Warning: package 'dplyr' was built under R version 4.2.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorslibrary(DescTools)## Warning: package 'DescTools' was built under R version 4.2.3library(ggplot2)
library(readxl)
setwd("C:/Users/YEN NHI/Documents/PTDLDT")1 Thống kê mô tả biến purpose
1.1 Thống kê mô tả
data <- read_excel("For_EDA_dataset 1.xlsx")## Warning: Expecting numeric in G3081 / R3081C7: got a date## New names:
## • `` -> `...1`pp = data$purpose
table(pp)## pp
## For Rent For Sale
## 43183 110247table(pp)/sum(table(pp))## pp
## For Rent For Sale
## 0.2814508 0.7185492Theo kết quả trên biến purpose (mục đích) là biến định tính, thể hiện mục đích sử dụng của căn hộ là để bán hay để cho thuê (For Rent: cho thuê; For Sale: bán). Trong tổng số 153430 căn hộ, số căn hộ dùng để cho thuê là 43183 căn, số căn hộ dùng để bán là 110247.
1.2 Đồ thị
data |> ggplot(aes(x = pp, y = after_stat(count))) +
geom_bar(fill = 'blue') +geom_text(aes(label = scales::percent(after_stat(count/sum(count)))), stat ='count', color = 'red', vjust = - .5) +labs(x = 'Số căn hộ', y = 'Mục đích sử dụng')
Dựa vào đồ thị trên, ta thấy rằng số căn hộ dùng để bán chiếm tỷ lệ cao hơn với tỷ lệ là 72%. Ngược lại, số căn hộ dùng để cho thuê chiếm tỷ lệ thấp hơn với tỷ lệ 28%. Tỷ lệ căn hộ dùng để bán gấp 2,57 lần so với tỷ lệ căn hộ dùng để cho thuê.
1.3 Mã hóa dữ liệu
1.3.1 Mã hóa dữ liệu price
pr = data$price
summary(pr)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000e+00 2.000e+05 8.700e+06 1.842e+07 2.000e+07 2.000e+09pr <- cut(pr, breaks = c(min(pr),mean(pr),max(pr)), labels = c('thap','cao'))
table(pr)## pr
## thap cao
## 111009 424181.3.2 Mã hóa dữ liệu bath
bt = data$baths
summary(bt)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.000 3.000 2.956 5.000 403.000bt <- cut(bt, breaks = c(min(bt),mean(bt),max(bt)), labels = c('it','nhieu'))
table(bt)## bt
## it nhieu
## 28058 885461.3.3 Mã hóa dữ liệu bedrooms
br = data$bedrooms
summary(br)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2.000 3.000 3.259 5.000 68.000br <- cut(br, breaks = c(min(br),mean(br),max(br)), labels = c('it','nhieu'))
table(br)## br
## it nhieu
## 72988 615631.3.4 Mã hóa dữ liệu Area
ar = data$Area_in_Marla
summary(ar)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 4.90 7.60 12.33 12.00 16000.00ar <- cut(ar, breaks = c(min(ar),mean(ar),max(ar)), labels = c('thap','cao'))
table(ar)## ar
## thap cao
## 115969 374491.3.5 Dữ liệu property_type
pt = data$property_type1.4 Rủi ro tương đối
library(epitools)
dt1 <- table(pp, bt)
addmargins(dt1)## bt
## pp it nhieu Sum
## For Rent 10021 20994 31015
## For Sale 18037 67552 85589
## Sum 28058 88546 116604RelRisk(dt1)## [1] 1.5331792 Ước lượng tỷ lệ
Ước lượng tỷ lệ căn hộ với mục đích cho bán có phải là 72% hay không.
data1 <- data[pp == "For Sale",]
prop.test(length(data1$purpose), length(data$purpose), p = 0.72)##
## 1-sample proportions test with continuity correction
##
## data: length(data1$purpose) out of length(data$purpose), null probability 0.72
## X-squared = 1.5948, df = 1, p-value = 0.2066
## alternative hypothesis: true p is not equal to 0.72
## 95 percent confidence interval:
## 0.7162903 0.7207971
## sample estimates:
## p
## 0.7185492Kết quả cho thấy, p_value = 0,2066 > 5%, tức là chấp nhận \(H_0\). Tại mức ý nghĩa 5% số căn hộ với mục đích bán là 72%. Khoảng ước lượng với độ tin cậy 95% của tỷ lệ số căn hộ bán là (0,7162903;0,7207971).
3 Mô hình hồi quy
3.1 Mô hình logit
mh1 <- glm(data = data, formula = factor(pp) ~ pr + bt + br + ar + pt, family = binomial(link = "logit"))
levels(factor(pp))## [1] "For Rent" "For Sale"summary(mh1)##
## Call:
## glm(formula = factor(pp) ~ pr + bt + br + ar + pt, family = binomial(link = "logit"),
## data = data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0384 -0.1225 0.0001 0.6002 3.7279
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.63737 0.26715 13.615 < 2e-16 ***
## prcao 20.76076 51.69003 0.402 0.688
## btnhieu -0.02784 0.02184 -1.275 0.202
## brnhieu -0.29317 0.02607 -11.246 < 2e-16 ***
## arcao -3.50672 0.03972 -88.289 < 2e-16 ***
## ptFlat -2.77101 0.26726 -10.368 < 2e-16 ***
## ptHouse -1.69364 0.26679 -6.348 2.18e-10 ***
## ptLower Portion -5.46444 0.27089 -20.173 < 2e-16 ***
## ptPenthouse -2.42300 0.31946 -7.585 3.33e-14 ***
## ptRoom -7.07815 0.39661 -17.846 < 2e-16 ***
## ptUpper Portion -4.69798 0.26854 -17.495 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 134646 on 116250 degrees of freedom
## Residual deviance: 75180 on 116240 degrees of freedom
## (37179 observations deleted due to missingness)
## AIC: 75202
##
## Number of Fisher Scoring iterations: 18\[ logit(\pi) = log(\frac{\pi}{1-\pi}) = 3.63737 - 0,29317brnhieu - 3,50672arcao - 2,77101ptFlat - 1,69364ptHouse - 5,46444ptLower Portion -2,423ptPentHouse - 7,07815ptRoom - 4,69798ptUpper Portion \]
3.2 Mô hình probit
mh2 <- glm(data = data, formula = factor(pp) ~ pr + bt + br + ar + pt, family = binomial(link = "probit"))## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredlevels(factor(pp))## [1] "For Rent" "For Sale"summary(mh2)##
## Call:
## glm(formula = factor(pp) ~ pr + bt + br + ar + pt, family = binomial(link = "probit"),
## data = data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0299 -0.0659 0.0000 0.6054 4.4244
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.102495 0.165922 12.672 < 2e-16 ***
## prcao 7.913991 10.420065 0.759 0.448
## btnhieu -0.005444 0.012589 -0.432 0.665
## brnhieu -0.168852 0.014413 -11.715 < 2e-16 ***
## arcao -2.035376 0.021167 -96.159 < 2e-16 ***
## ptFlat -1.565692 0.165995 -9.432 < 2e-16 ***
## ptHouse -0.963865 0.165751 -5.815 6.06e-09 ***
## ptLower Portion -3.176946 0.167654 -18.949 < 2e-16 ***
## ptPenthouse -1.342918 0.194357 -6.910 4.86e-12 ***
## ptRoom -3.929479 0.206266 -19.051 < 2e-16 ***
## ptUpper Portion -2.744848 0.166708 -16.465 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 134646 on 116250 degrees of freedom
## Residual deviance: 75247 on 116240 degrees of freedom
## (37179 observations deleted due to missingness)
## AIC: 75269
##
## Number of Fisher Scoring iterations: 18\[ probit(\pi) = \Phi^{-1}(\pi) = 2,102495 - 0,168852brnhieu - 2,035376arcao - 1,565692ptFlat - 0,963865ptHouse - 3,176946ptLower Portion -1,342918ptPentHouse - 3,929479ptRoom - 2,744848ptUpper Portion \]
3.3 Mô hình cloglog
mh3 <- glm(data = data, formula = factor(pp) ~ pr + bt + br + ar + pt, family = binomial(link = "cloglog"))## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredlevels(factor(pp))## [1] "For Rent" "For Sale"summary(mh3)##
## Call:
## glm(formula = factor(pp) ~ pr + bt + br + ar + pt, family = binomial(link = "cloglog"),
## data = data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0601 -0.1899 0.0000 0.5981 3.4869
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.30393 0.19488 11.822 < 2e-16 ***
## prcao 8.46391 5.21796 1.622 0.105
## btnhieu -0.01597 0.01203 -1.328 0.184
## brnhieu -0.14346 0.01266 -11.336 < 2e-16 ***
## arcao -2.60711 0.03476 -74.995 < 2e-16 ***
## ptFlat -2.11722 0.19495 -10.860 < 2e-16 ***
## ptHouse -1.55161 0.19479 -7.965 1.65e-15 ***
## ptLower Portion -4.23572 0.19941 -21.241 < 2e-16 ***
## ptPenthouse -1.98745 0.21746 -9.139 < 2e-16 ***
## ptRoom -5.77482 0.34846 -16.572 < 2e-16 ***
## ptUpper Portion -3.55301 0.19638 -18.093 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 134646 on 116250 degrees of freedom
## Residual deviance: 75073 on 116240 degrees of freedom
## (37179 observations deleted due to missingness)
## AIC: 75095
##
## Number of Fisher Scoring iterations: 20\[ cloglog(\pi) = log(-log(1-\pi)) = 2,3093 - 0,14346brnhieu - 2,60711arcao - 2,11722ptFlat - 1,55161ptHouse - 4,23572ptLower Portion -1,98745ptPentHouse - 5,77482ptRoom - 3,55301ptUpper Portion \]
3.4 Lựa chọn mô hình
Để đánh giá các mô hình hồi quy trên, ta sử dụng các tiêu chí sau:
# Tiêu chí AIC - Akaike Information Criterion
aic1 <- AIC(mh1)
aic2 <- AIC(mh2)
aic3 <- AIC(mh3)
AIC <-cbind(aic1,aic2,aic3)
AIC## aic1 aic2 aic3
## [1,] 75201.83 75268.78 75095.45# Tiêu chí Deviance
de1 <- deviance(mh1)
de2 <- deviance(mh2)
de3 <- deviance(mh3)
deviance <- cbind(de1,de2,de3)
deviance## de1 de2 de3
## [1,] 75179.83 75246.78 75073.45# Tiêu chí Brier Score
bs1 <- BrierScore(mh1)
bs2 <- BrierScore(mh2)
bs3 <- BrierScore(mh3)
BrierScore <- cbind(bs1,bs2,bs3)
BrierScore## bs1 bs2 bs3
## [1,] 0.1027722 0.1028198 0.1026701Dựa vào các tiêu chí trên, ta thấy các giá trị AIC, deviance và Brier Score của mô hình cloglog là nhỏ nhất, tức là mô hình cloglog là mô hình tốt nhất trong 3 mô hình.