library(readxl)
## Warning: package 'readxl' was built under R version 4.2.3
library(DescTools)
## Warning: package 'DescTools' was built under R version 4.2.3
library(ggplot2)
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
setwd("C:/PTDLDT")
data <- read_excel("For_EDA_dataset 1.xlsx")
## Warning: Expecting numeric in G3081 / R3081C7: got a date
## New names:
## • `` -> `...1`
data
str(data)
## tibble [153,430 × 15] (S3: tbl_df/tbl/data.frame)
##  $ ...1         : num [1:153430] 0 1 2 3 4 5 6 7 8 9 ...
##  $ property_type: chr [1:153430] "Flat" "Flat" "House" "House" ...
##  $ price        : num [1:153430] 10000000 6900000 16500000 43500000 7000000 34500000 27000000 7800000 50000000 40000000 ...
##  $ location     : chr [1:153430] "G-10" "E-11" "G-15" "Bani Gala" ...
##  $ city         : chr [1:153430] "Islamabad" "Islamabad" "Islamabad" "Islamabad" ...
##  $ province_name: chr [1:153430] "Islamabad Capital" "Islamabad Capital" "Islamabad Capital" "Islamabad Capital" ...
##  $ latitude     : num [1:153430] 3.37e+06 3.37e+07 3.36e+16 3.37e+13 3.35e+07 ...
##  $ longitude    : num [1:153430] 7.30e+06 7.30e+07 7.29e+07 7.32e+12 7.33e+07 ...
##  $ baths        : num [1:153430] 2 3 6 4 3 8 8 2 7 5 ...
##  $ purpose      : chr [1:153430] "For Sale" "For Sale" "For Sale" "For Sale" ...
##  $ bedrooms     : num [1:153430] 2 3 5 4 3 8 8 2 7 5 ...
##  $ date_added   : POSIXct[1:153430], format: "2019-02-04" "2019-05-04" ...
##  $ agency       : chr [1:153430] "Self" "Self" "Self" "Self" ...
##  $ agent        : chr [1:153430] "Self" "Self" "Self" "Self" ...
##  $ Area_in_Marla: num [1:153430] 4 5.6 8 40 8 32 20 6.2 20 20 ...

1 Thống kê mô tả các biến phụ thuộc Purpose

1.1 Bảng tần số, bảng tần suất

table(data$purpose)
## 
## For Rent For Sale 
##    43183   110247
table(data$purpose)/sum(table(data$purpose))
## 
##  For Rent  For Sale 
## 0.2814508 0.7185492

Biến Purpose là biến định tính, thể hiện mục đích sử dụng căn hộ (For Rent: cho thuê; For Sale:rao bán). Dựa vào kết quả thống kê trên ta thấy trong tổng số 43183 căn hộ cho thuê (chiếm 28,14% tổng số căn hộ) và 110247 căn hộ để rao bán (chiếm 71,85% tổng số căn hộ).

1.2 Đồ thị cột

ggplot(data,aes(purpose)) + geom_bar(color ="navy", fill = "pink") + ylab("Số căn hộ") + xlab("Mục đích sử dụng")

Đồ thị cột thể hiện số căn hộ được chia theo mục đích sử dụng là cho thuê hoặc rao bán. Ta có thể thấy mục đích sử dụng giữa cho thuê và rao bán có sự chênh lệch rõ rệt khi tỷ lệ rao bán (chiếm 72%) cao hơn tỷ lệ cho thuê là (chiếm 28%).

2 Mã hóa các biến định lượng sang định tính

2.1 Mã hóa biến Price

price <- cut(data$price, breaks = c(0,8700000,2000000000), labels=c("Thấp","Cao"))
table(price)
## price
##  Thấp   Cao 
## 76817 76610

Với dữ liệu gốc, biến price là biến định lượng nhận các giá trị từ 0 đến 2000000000 nên tác giả đã đặt quy ước về việc mã hoá biến price như sau: Chi phí thấp (Thấp): từ 0 đến 8700000. Chi phí cao (Cao): chi phí từ 8700000 đến 2000000000.

2.2 Mã hóa biến Baths

baths <- cut(data$baths, breaks = c(0,3,403), labels=c("It","Nhiều"))
table(baths)
## baths
##    It Nhiều 
## 59031 57573

Với dữ liệu gốc, biến baths là biến định lượng nhận các giá trị từ 0 đến 403 nên tác giả đã đặt quy ước về việc mã hoá biến baths như sau: Số lượng ít từ 0 đến 3. Số lượng nhiều từ 3 đến 403.

2.3 Mã hóa biến Bedrooms

bedrooms <- cut(data$bedrooms, breaks = c(0,3,68), labels=c("It","Nhiều"))
table(bedrooms)
## bedrooms
##    It Nhiều 
## 72988 61563

Với dữ liệu gốc, biến bedrooms là biến định lượng nhận các giá trị từ 0 đến 68 nên tác giả đã đặt quy ước về việc mã hoá biến bedrooms như sau: Số lượng ít từ 0 đến 3. Số lượng nhiều từ 3 đến 68.

2.4 Mã hóa biến Area_in_Marla

Area_in_Marla <- cut(data$Area_in_Marla, breaks = c(0,12,16000), labels=c("Nhỏ","To"))
table(Area_in_Marla)
## Area_in_Marla
##    Nhỏ     To 
## 115824  37594

Với dữ liệu gốc, biến Area_in_Marla là biến định lượng nhận các giá trị từ 0 đến 16000 nên tác giả đã đặt quy ước về việc mã hoá biến Area_in_Marla như sau: Diện tích nhỏ từ 0 đến 12. Diện tích lớn từ 12 đến 16000.

3 Ước lượng tỷ lệ cho biến Purpose

Ước lượng tỷ lệ số căn hộ được dùng để rao bán có phải là 72% không (nghĩa là chúng ta kiểm định giả thuyết \(H_{0}\): p = 0,72)

d <- data[data$purpose == "For Sale",]
prop.test(length(d$purpose),length(data$purpose),p= 0.72)
## 
##  1-sample proportions test with continuity correction
## 
## data:  length(d$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.7185492

Dựa vào kết quả ta thấy p-value = 0,2066 > 5%, chấp nhận \(H_{0}\), tức là tại mứuc ý nghĩa 5% số căn hộ với mục đích rao bán là 72%.

Khoảng ước lượng với độ tin cậy 95% của tỷ lệ số căn hộ rao bán là (0,7162903; 0,7207971)

4 Mô hình hồi quy

purpose <- data$purpose
pt <-  data$property_type

4.1 Mô hình hồi quy logit

k <- data.frame(purpose, pt, price, baths, bedrooms, Area_in_Marla)
mhlogit <- glm(data = k, formula = factor(purpose) ~  pt + price + baths + bedrooms + Area_in_Marla, family = binomial(link = "logit"))
levels(factor(purpose))
## [1] "For Rent" "For Sale"
summary(mhlogit)
## 
## Call:
## glm(formula = factor(purpose) ~ pt + price + baths + bedrooms + 
##     Area_in_Marla, family = binomial(link = "logit"), data = k)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -5.1298  -0.0660   0.0020   0.0136   3.9227  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      3.33037    0.43371   7.679 1.61e-14 ***
## ptFlat          -2.84648    0.43393  -6.560 5.39e-11 ***
## ptHouse         -1.85083    0.43354  -4.269 1.96e-05 ***
## ptLower Portion -5.59273    0.43736 -12.787  < 2e-16 ***
## ptPenthouse     -2.60185    0.47304  -5.500 3.79e-08 ***
## ptRoom          -7.15869    0.56084 -12.764  < 2e-16 ***
## ptUpper Portion -4.79993    0.43508 -11.032  < 2e-16 ***
## priceCao        13.20083    0.70942  18.608  < 2e-16 ***
## bathsNhiều      -1.30432    0.03918 -33.287  < 2e-16 ***
## bedroomsNhiều   -0.21879    0.03988  -5.487 4.09e-08 ***
## Area_in_MarlaTo -3.86518    0.08284 -46.657  < 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:  52927  on 116240  degrees of freedom
##   (37179 observations deleted due to missingness)
## AIC: 52949
## 
## Number of Fisher Scoring iterations: 12

Đối với mô hình hồi quy logit của biến phụ thuộc purpose với các biến property_type, biến price, biến baths, biến bedrooms, biến Area_in_Marla. Ta thấy tại mức ý nghĩa 5% cả 3 biến property_type, bedrooms và Area_in_Marla có ý nghĩa thống kê, 2 biến còn lại là price và baths không có ý nghĩa thống kê. Chiều hướng tác động của các biến property_type, bedrooms và Area_in_Marla đến biến phụ thuộc purpose là ngược chiều.

Vậy nên mô hình logit được xác định như sau:

MHlogit:

\[ logit(π) = log(π/1−π) = 3,63737-2,77101ptFlat-1,69364ptHouse-5,46444ptLowerPorrtion-2,42300ptPenhouse -7,07815ptRoom-4,69798ptUpperPortion-0,29317bedroomsNhieu-3,50672Area_in_MarlaTo\]

4.2 Mô hình hồi quy probit

k <- data.frame(purpose, pt, price, baths, bedrooms, Area_in_Marla)
mhprobit <- glm(data = k, formula = factor(purpose) ~  pt + price + baths + bedrooms + Area_in_Marla, family = binomial(link = "probit"))
levels(factor(purpose))
## [1] "For Rent" "For Sale"
summary(mhprobit)
## 
## Call:
## glm(formula = factor(purpose) ~ pt + price + baths + bedrooms + 
##     Area_in_Marla, family = binomial(link = "probit"), data = k)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -6.3759  -0.0288   0.0001   0.0100   4.5586  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      1.71347    0.24377   7.029 2.08e-12 ***
## ptFlat          -1.40892    0.24392  -5.776 7.64e-09 ***
## ptHouse         -0.83230    0.24365  -3.416 0.000635 ***
## ptLower Portion -3.01278    0.24545 -12.274  < 2e-16 ***
## ptPenthouse     -1.25051    0.26895  -4.650 3.33e-06 ***
## ptRoom          -3.67631    0.27963 -13.147  < 2e-16 ***
## ptUpper Portion -2.58842    0.24458 -10.583  < 2e-16 ***
## priceCao         5.96344    0.18349  32.500  < 2e-16 ***
## bathsNhiều      -0.79384    0.02334 -34.010  < 2e-16 ***
## bedroomsNhiều   -0.11783    0.02380  -4.951 7.37e-07 ***
## Area_in_MarlaTo -2.04432    0.03633 -56.278  < 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:  52972  on 116240  degrees of freedom
##   (37179 observations deleted due to missingness)
## AIC: 52994
## 
## Number of Fisher Scoring iterations: 12

Đối với mô hình hồi quy probit của biến phụ thuộc purpose với các biến property_type, biến price, biến baths, biến bedrooms, biến Area_in_Marla. Ta thấy tại mức ý nghĩa 5% cả 3 biến property_type, bedrooms và Area_in_Marla có ý nghĩa thống kê, 2 biến còn lại là price và baths không có ý nghĩa thống kê. Chiều hướng tác động của các biến property_type, bedrooms và Area_in_Marla đến biến phụ thuộc purpose là ngược chiều.

Vậy nên mô hình probit được xác định như sau:

MHprobit:

\[ probit(π) = Φ − 1(π) = 2,102495-1,565692ptFlat-0,963865ptHouse-3,176946ptLowerPorrtion-1,342918ptPenhouse -3,929479ptRoom-2,744848ptUpperPortion-0,168852bedroomsNhieu-2,035376Area_in_MarlaTo\]

4.3 Mô hình hồi quy cloglog

k <- data.frame(purpose, pt, price, baths, bedrooms, Area_in_Marla)
mhcloglog <- glm(data = k, formula = factor(purpose) ~  pt + price + baths + bedrooms + Area_in_Marla, family = binomial(link = "cloglog"))
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(mhcloglog)
## 
## Call:
## glm(formula = factor(purpose) ~ pt + price + baths + bedrooms + 
##     Area_in_Marla, family = binomial(link = "cloglog"), data = k)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -8.4904  -0.0806   0.0000   0.0000   3.7718  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      2.48644    0.37536   6.624 3.49e-11 ***
## ptFlat          -2.53062    0.37547  -6.740 1.59e-11 ***
## ptHouse         -1.96941    0.37534  -5.247 1.55e-07 ***
## ptLower Portion -4.76690    0.37881 -12.584  < 2e-16 ***
## ptPenthouse     -2.40454    0.39359  -6.109 1.00e-09 ***
## ptRoom          -6.15084    0.48853 -12.590  < 2e-16 ***
## ptUpper Portion -4.07193    0.37654 -10.814  < 2e-16 ***
## priceCao         7.00510    0.16119  43.459  < 2e-16 ***
## bathsNhiều      -0.84572    0.02553 -33.128  < 2e-16 ***
## bedroomsNhiều   -0.05855    0.02508  -2.334   0.0196 *  
## Area_in_MarlaTo -3.44831    0.08052 -42.823  < 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:  53108  on 116240  degrees of freedom
##   (37179 observations deleted due to missingness)
## AIC: 53130
## 
## Number of Fisher Scoring iterations: 15

Đối với mô hình hồi quy cloglog của biến phụ thuộc purpose với các biến property_type, biến price, biến baths, biến bedrooms, biến Area_in_Marla. Ta thấy tại mức ý nghĩa 5% cả 3 biến property_type, bedrooms và Area_in_Marla có ý nghĩa thống kê, 2 biến còn lại là price và baths không có ý nghĩa thống kê. Chiều hướng tác động của các biến property_type, bedrooms và Area_in_Marla đến biến phụ thuộc purpose là ngược chiều.

Vậy nên mô hình probit được xác định như sau:

MHprobit:

\[ cloglog(π) = log(−log(1−π)) = 2,30393-2,11722ptFlat-1,55161ptHouse-4,23572ptLowerPorrtion-1,98745ptPenhouse -5,77482ptRoom-3,55301ptUpperPortion-0,14346bedroomsNhieu-2,607111Area_in_MarlaTo\]

4.4 Lựa chọn mô hình hồi quy phù hợp

4.4.1 AIC - Akaike Information Criterion

AIC(mhlogit)
## [1] 52949.04
AIC(mhprobit)
## [1] 52994.4
AIC(mhcloglog)
## [1] 53130.08

Từ chỉ số AIC của 3 mô hình trên ta thấy mô hình cloglog có chỉ số AIC thấp nhất (75095,45). Vì thế đối với tiêu chí đánh giá AIC thì mô hình cloglog là phù hợp để xem xét sự tác động của các yếu tố đến biến phụ thuộc Purpose hơn mô hình logit và probit.

4.4.2 Deviance

deviance(mhlogit)
## [1] 52927.04
deviance(mhprobit)
## [1] 52972.4
deviance(mhcloglog)
## [1] 53108.08

Từ chỉ số deviance của 3 mô hình trên ta thấy mô hình cloglog có chỉ số deviance thấp nhất (75073,45). Vì thế đối với tiêu chí đánh giá devience thì mô hình cloglog là phù hợp để xem xét sự tác động của các yếu tố đến biến phụ thuộc Purpose hơn mô hình logit và probit.

4.4.3 Brier Score

BrierScore(mhlogit)
## [1] 0.07636095
BrierScore(mhprobit)
## [1] 0.07639882
BrierScore(mhcloglog)
## [1] 0.0765309

Từ chỉ số BrierScore của 3 mô hình trên ta thấy mô hình cloglog có chỉ số BrierScore thấp nhất (0,1026701). Vì thế đối với tiêu chí đánh giá BrierScore thì mô hình cloglog là phù hợp để xem xét sự tác động của các yếu tố đến biến phụ thuộc Purpose hơn mô hình logit và probit.

Thông qua các tiêu chí đánh giá mô hình AIC, Deviance và BrierScore thì mô hình phù hợp nhất để xem xét tác động của các yếu tố property_type, bedrooms, Area_in_Marla tới biến phụ thuộc purpose là mô hình hồi quy cloglog