1 Giới thiệu
Kaggle là một trong những tổ chức lớn nhất trong lĩnh vực Data Science. Website này tổ chức rất nhiều các cuộc thi với data thực tế từ các tổ chức với sự tham gia của các chuyên gia hàng đầu trên toàn thế giới với các giải thưởng lớn.
Kaggle bike sharing demand - là một trong những cuộc thi của Kaggle cho giới Data Scientist được diễn ra vào năm 2015 với câu hỏi được đặt ra là: Dựa vào data của một hãng cung cấp cho thuê xe đạp ở Washington D.C, người tham gia cần phải dự báo được số lượng xe đạp sẽ được thuê.
Bài viết này sẽ trình bày một phương pháp dự báo lọt vào top 5% phương pháp xuất sắc nhất trong cuộc thi. Tuy không phải là bài thi suất sắc nhất nhưng sẽ giúp chúng ta nâng cao được khả năng phân tích và sử dụng các tips trong phân tích
2 Phân tích
2.1 Khám phá dữ liệu
Trước khi đi vào các phương pháp predictive modeling, chúng ta cần phải tìm hiểu những đặc tính bằng cách phân tích khám phá dữ liệu và đưa ra các câu hỏi, giả thuyết kết hợp giữa vấn đề kinh doanh và phân tích thống kê
test <- read.csv(file = "./test.csv")
train <- read.csv(file = "./train.csv")
test$registered <- 0
test$casual <- 0
test$count <- 0
data <- rbind(test, train)
save(data, file = "./data.Rda")
rm(list = ls())2.2 Data structure
- Data bao gồm 12 biến: Mùa trong năm, có phải ngày nghỉ lễ không, thời tiết, nhiệt độ, độ ẩm …
library(dplyr)
library(ggplot2)
library(ggthemes)
library(rpart)
library(rattle)
library(rpart.plot)
library(RColorBrewer)
library(lubridate)
library(randomForest)
load(file = "./data.Rda")
test <- read.csv(file = "./test.csv")
train <- read.csv(file = "./train.csv")
data %>% str## 'data.frame': 17379 obs. of 12 variables:
## $ datetime : Factor w/ 17379 levels "2011-01-20 00:00:00",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ season : int 1 1 1 1 1 1 1 1 1 1 ...
## $ holiday : int 0 0 0 0 0 0 0 0 0 0 ...
## $ workingday: int 1 1 1 1 1 1 1 1 1 1 ...
## $ weather : int 1 1 1 1 1 1 1 1 1 2 ...
## $ temp : num 10.7 10.7 10.7 10.7 10.7 ...
## $ atemp : num 11.4 13.6 13.6 12.9 12.9 ...
## $ humidity : int 56 56 56 56 56 60 60 55 55 52 ...
## $ windspeed : num 26 0 0 11 11 ...
## $ registered: num 0 0 0 0 0 0 0 0 0 0 ...
## $ casual : num 0 0 0 0 0 0 0 0 0 0 ...
## $ count : num 0 0 0 0 0 0 0 0 0 0 ...2.3 Phân phối của các biến
Một vài nhận xét được rút ra như sau
Biến Season sẽ gồm 4 giá trị tương ứng với 4 mùa
Biến weather có phân phối chủ yếu ở giá trị 1
2 biến workingday và holiday không có thông tin gì để rút ra
Các biến humidity, atemp và windspeed có phân phối khá tự nhiên
par(mfrow=c(4,2))
par(mar = rep(2, 4))
hist(data$season)
hist(data$weather)
hist(data$humidity)
hist(data$holiday)
hist(data$workingday)
hist(data$temp)
hist(data$atemp)
hist(data$windspeed)
data$season <- data$season %>% as.factor()
data$weather <- data$weather %>% as.factor()
data$holiday <- data$holiday %>% as.factor()
data$workingday <- data$workingday %>% as.factor()2.4 Giả thuyết
Như vậy chúng ta đã có thể hình dung sơ lược về data. Tiếp đến ta sẽ đi sâu hơn về phân tích data này bằng việc đưa ra các giả thuyết và kiểm định xem giả thuyết đó có đúng hay sai.
2.4.1 Hourly trend
Giả thuyết: việc thuê xe đạp sẽ bị ảnh hưởng bởi vấn đề thời gian
data$hour <- substr(data$datetime,12,13)
data$hour <- as.factor(data$hour)
train <- data[as.integer(substr(data$datetime,9,10))<20,]
test <- data[as.integer(substr(data$datetime,9,10))>19,]
train %>%
ggplot(aes(x = hour, y = count)) +
geom_boxplot(fill = "darkgreen") +
labs(x = "Hour", y = "Count of Users") +
theme_minimal()
Dựa vào biểu đồ ta có thể thấy vài điểm như sau
Thời gian thuê xe cao điểm sẽ rơi vào từ 7 - 9h và 17 - 19h
Thời gian thuê xe trung bình sẽ rơi vào từ 10 - 16h
Ta sẽ chia thêm theo khách hàng Casual và Registered
train %>%
ggplot(aes(x = hour, y = casual)) +
geom_boxplot(fill = "darkgreen") +
labs(x = "Hour", y = "Count of Casual Users") +
theme_minimal()
Đối với khách hàng Casual thì trend rất khác so với tổng số khách hàng
train %>%
ggplot(aes(x = hour, y = registered)) +
geom_boxplot(fill = "darkgreen") +
labs(x = "Hour", y = "Count of Registered Users") +
theme_minimal()
Tuy nhiên đối với khách hàng Registered thì trend lại giống với tổng số khách hàng
Ta có thể nhận thấy khi plot số lượng khách hàng sử dụng xe đạp sẽ bị ảnh hưởng rất nhiều bởi outliers và đa phần đến từ nhóm khách hàng vãng lai (casual customers). Để xử lý vấn đề này được tốt hơn và có cái nhìn tổng quan hơn ta sẽ sử dụng hàm log
train %>%
ggplot(aes(x = hour, y = log(count))) +
geom_boxplot(fill = "darkgreen") +
labs(x = "Hour", y = "Log of Count of Registered Users") +
theme_minimal()
Ta nhận thấy được là Hourly Trend có xuất hiện trong nhu cầu sử dụng xe đạp
2.4.2 Daily Hour
Tương tự đó là Daily Trend ta có thể thấy là nhu cầu của việc sử dụng xe đạp tăng cao hơn đối với khách hàng vãng lai vào dịp cuối tuần
date <- substr(data$datetime,1,10)
days <- weekdays(as.Date(date))
data$day <- days
data %>%
ggplot(aes(x = day, y = casual)) +
geom_boxplot(fill = "darkgreen") +
labs(x = "Hour", y = "Count of Casual Users") +
theme_minimal()
data %>%
ggplot(aes(x = day, y = registered)) +
geom_boxplot(fill = "darkgreen") +
labs(x = "Hour", y = "Count of Registered Users") +
theme_minimal()
2.5 Feature Enginering
Để nâng cao khả năng dự báo của model ta sẽ tạo ra những biến mới dựa trên phương pháp Decision Tree
2.5.1 Biến Hour
train$hour <- train$hour %>% as.integer()
test$hour <- test$hour %>% as.integer()
d <- rpart(registered ~ hour,data = train)
fancyRpartPlot(d)
Dựa trên Decision Tree Plot ta sẽ chia biến Hour thành các 7 bins cho nhóm Registered
data <- rbind(train,test)
data <- data %>%
mutate(dp_reg = case_when(
hour < 8 ~ 1,
hour >= 22 ~ 2,
hour > 9 & hour < 18 ~ 3,
hour == 8 ~ 4,
hour == 9 ~ 5,
hour %in% c(20,21) ~ 6,
hour %in% c(18,19) ~ 7))Ta sẽ làm tương tự như vậy đối với khách hàng Casual Users và Temperature cho cả Casual và Registered
2.5.2 Biến Year
Dựa trên phân tích khám phá dữ liệu chúng ta có thể thấy được rằng nhu cầu về sử dụng xe đạp luôn gia tăng qua từng năm nên ta sẽ chia biến Year thành 8 bins là các quý
data$year <- substr(data$datetime,1,4)
data$year <- data$year %>% as.factor
data$month <- month(data$datetime)
data$year_part[data$year == '2011'] <- 1
data$year_part[data$year == '2011' & data$month > 3] <- 2
data$year_part[data$year == '2011' & data$month > 6] <- 3
data$year_part[data$year == '2011' & data$month > 9] <- 4
data$year_part[data$year == '2012'] <- 5
data$year_part[data$year == '2012' & data$month > 3] <- 6
data$year_part[data$year == '2012' & data$month > 6] <- 7
data$year_part[data$year == '2012' & data$month > 9] <- 8
table(data$year_part)##
## 1 2 3 4 5 6 7 8
## 2067 2183 2192 2203 2176 2182 2208 21682.5.3 Biến Day
data$day_type[data$holiday==0 & data$workingday==0] <- "weekend"
data$day_type[data$holiday==1] <- "holiday"
data$day_type[data$holiday==0 & data$workingday==1] <- "working day"
data$day_type <- data$day_type %>% as.factor()2.5.4 Biến Weekend
data$weekend <- 0
data$weekend[data$day=="Sunday" | data$day=="Saturday" ] <- 13 Xây dựng mô hình dự báo
Mô hình được sử dụng trong dự báo là mô hình Random Forest. Ngoài ra cũng có các mô hình khác có thể ứng dụng tốt như: XGBoost, GBM …
- Đầu tiên ta sẽ dự báo cho khách hàng Registered dựa trên hàm LOG của Registered
train$hour <- train$hour %>% as.factor()
test$hour <- test$hour %>% as.factor()
set.seed(999)
fit1 <- randomForest(logreg ~ hour +workingday + day+ holiday + day_type +
temp_reg + humidity +atemp + windspeed + season +
weather + dp_reg + temp_reg + weekend + year + year_part,
data = train,
importance = TRUE,
ntree = 250)
pred1 <- predict(fit1,test)
test$logreg <- pred1- Tương tự với dự báo cho khách hàng Casual dựa trên hàm LOG của Casual
set.seed(999)
fit2 <- randomForest(logcas ~ hour + day_type + day + humidity + atemp +
temp_cas + windspeed + season + weather + holiday +
workingday + dp_cas + temp_cas +weekend + year + year_part,
data = train,
importance = TRUE,
ntree = 250)
pred2 <- predict(fit2,test)
test$logcas <- pred2- Sau đó sẽ biến đổi dạng Log của biến predict đầu ra để đưa ra kết quả cuối cùng
test$registered <- exp(test$logreg) - 1
test$casual <- exp(test$logcas) - 1
test$count <- test$casual + test$registered
result <-data.frame(datetime=test$datetime, count=test$count)