library(readxl)
## Warning: package 'readxl' was built under R version 4.2.3
library(epitools)
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:/PTDLDT1")
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 ="yellow", fill = "purple") + 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 Ướ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)

3 Hồi quy

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

  • Mã hóa biến Price
summary(data$price)
##      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+09
price <- data$price
price <- cut(data$price, breaks = c(min(price),mean(price),max(price)), labels=c("Thấp","Cao"))
table(price)
## price
##   Thấp    Cao 
## 111009  42418
  • Mã hóa biến Baths
summary(data$baths)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.000   3.000   2.956   5.000 403.000
baths <- data$baths
baths <- cut(data$baths, breaks = c(min(baths),mean(baths),max(baths)),
labels=c("It","Nhieu "))
table(baths)
## baths
##     It Nhieu  
##  28058  88546
  • Mã hóa biến Bedrooms
summary(data$bedrooms)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.000   3.000   3.259   5.000  68.000
bedrooms <- data$bedrooms
bedrooms <- cut(data$bedrooms, breaks = c(min(bedrooms),mean(bedrooms),max(bedrooms)), labels=c("It","Nhieu"))
table(bedrooms)
## bedrooms
##    It Nhieu 
## 72988 61563
  • Mã hóa biến Area_in_Marla
summary(data$Area_in_Marla)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##     0.00     4.90     7.60    12.33    12.00 16000.00
Area_in_Marla<- data$Area_in_Marla
Area_in_Marla <- cut(data$Area_in_Marla, breaks = c(min(Area_in_Marla),mean(Area_in_Marla),max(Area_in_Marla)), labels=c("Nhỏ","To"))
table(Area_in_Marla)
## Area_in_Marla
##    Nhỏ     To 
## 115969  37449

3.2 Mô hình hồi quy cho biến phụ thuộc Purpose

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

3.2.1 Mô hình 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  
## -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 ***
## 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 ***
## priceCao        20.76076   51.69003   0.402    0.688    
## bathsNhieu      -0.02784    0.02184  -1.275    0.202    
## bedroomsNhieu   -0.29317    0.02607 -11.246  < 2e-16 ***
## Area_in_MarlaTo -3.50672    0.03972 -88.289  < 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

Đố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\]

3.2.2 Mô hình 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"))
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
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  
## -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 ***
## 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 ***
## priceCao         7.913991  10.420065   0.759    0.448    
## bathsNhieu      -0.005444   0.012589  -0.432    0.665    
## bedroomsNhieu   -0.168852   0.014413 -11.715  < 2e-16 ***
## Area_in_MarlaTo -2.035376   0.021167 -96.159  < 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

Đố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\]

3.2.3 Mô hình 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  
## -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 ***
## 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 ***
## priceCao         8.46391    5.21796   1.622    0.105    
## bathsNhieu      -0.01597    0.01203  -1.328    0.184    
## bedroomsNhieu   -0.14346    0.01266 -11.336  < 2e-16 ***
## Area_in_MarlaTo -2.60711    0.03476 -74.995  < 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

Đố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\]

3.3 Lựa chọn mô hình hồi quy phù hợp cho biến phụ thuộc Purpose

3.3.1 AIC - Akaike Information Criterion

AIC(mhlogit)
## [1] 75201.83
AIC(mhprobit)
## [1] 75268.78
AIC(mhcloglog)
## [1] 75095.45

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.

3.3.2 Deviance

deviance(mhlogit)
## [1] 75179.83
deviance(mhprobit)
## [1] 75246.78
deviance(mhcloglog)
## [1] 75073.45

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.

3.3.3 Brier Score

BrierScore(mhlogit)
## [1] 0.1027722
BrierScore(mhprobit)
## [1] 0.1028198
BrierScore(mhcloglog)
## [1] 0.1026701

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