Github Code: Jun4871 Github

개요

캐글의 자전거 수요량 데이터를 사용하여 Training Data 와 Test Data 로 나누어 수요량을 예측해보고자 한다. 평가는 RMSLE(Root Mean Square Logarithmic Error)로 하며 이는 RMSE에 각각 log를 씌운것으로 오차에 대해서 과대평가된 항목보다 과소평가된 항목에 패널티를 더 주게 된다. 식은 다음과 같다.

\[ \sqrt{\frac{1}{n} \sum(\log(a_i + 1) - \log(a_i+1))2 } \]

라이브러리 활성화

데이터 분석에 필요한 라이브러리 활성화를 활성화 시켜주도록 한다. 사용할 라이브러리는 다음과 같다.

library(ggplot2) # 시각화 
library(tidyverse) # 전처리
library(lubridate) # 날짜 
library(stringr)
library(caret)
library(readr) # 데이터 로딩
library(gridExtra) # 차트 분할
library(xgboost)
library(Metrics)

데이터 셋

예측을 위해서는 모델을 만들어야 하는데, 이 때 모델을 만들기 위한 데이터가 필요하다. 이것을 트레이닝 세트라고 하고, 또한 모델이 신규 데이터에 대해서도 정확한지 테스트할 데이터도 필요한데 이것을 테스트 세트라고 한다. 트레이닝 세트는 모델을 구축하는 알고리즘(회귀, 의사결정나무 등)에 사용할 데이터로, 알고리즘이 출력변수를 가장 잘 예측할 수 있는 올바른 인자를 설정하도록 한다. 테스트 세트는 모델의 결과가 정확한지 검증하는 결과를 만드는 데 필요하다. read_csv() 함수를 사용하면 날짜형 데이터로 바로 변환이 된다.

# train_set 및 test_set 할당
train_set <-read_csv("train.csv")
test_set <- read_csv("test.csv")
SampleSubmission <- read_csv("sampleSubmission.csv")

탐색적 데이터 분석

Bike Sharing Demand 의 데이터는 다음과 같다. 이중에서 casual, registered는 test_set에 없어서 사용하지 않을 것이고, 날짜와 날씨 데이터로 이루어져 있다.

datetime: 년-월-일 시간 데이터
season: 1 = 봄, 2 = 여름, 3 = 가을, 4 = 겨울
holiday: 공휴일 또는 주말
workingday: 공휴일, 주말을 제외한 평일
weather
- 1: 매우 맑음(Clear, Few clouds, Partly cloudy, Partly cloudy)
- 2: 맑음(Mist + Cloudy, Mist + Broken clouds, Mist + Few clouds, Mist)
- 3: 나쁨(Light Snow, Light Rain + Thunderstorm + Scattered clouds, Light Rain + Scattered clouds)
- 4: 매우 나쁨(Heavy Rain + Ice Pallets + Thunderstorm + Mist, Snow + Fog)
temp: 기온
atemp: 체감온도 정도로 보자
humidity: 상대 습도
windspeed: 바람의 세기
casual: 미등록 사용자 렌탈 수
registered: 등록된 사용자 렌탈수
count: 렌탈한 총 합

dim(train_set)
str(train_set) 
summary(train_set)
sum(is.na(train_set))

dim(test_set)
str(test_set)
summary(test_set)
sum(is.na(test_set))

트레인 셋 가공

모델링 전에 데이터를 가공하는 작업으로, 컬럼을 추가 또는 삭제하는 등 알맞게 데이터를 정제하는 과정이다. 데이터를 불러오고 전처리를 해주도록 하자. xgboost를 사용하기 위해선, 각각의 데이터는 모두 숫자형으로 표현 되어야 한다. 날짜 데이터는 년, 월, 시간, 요일로 따로 데이터를 불러 오도록 하자.

train_set <- train_set %>% 
  select(-casual, -registered) %>% 
  mutate(
    year = year(datetime),
    month = month(datetime),
    hour = hour(datetime),
    wday = wday(datetime))


test_set <- test_set %>% 
  mutate(
    year = year(datetime),
    month = month(datetime),
    hour = hour(datetime),
    wday = wday(datetime))

시각화

이렇게 데이터를 파악 해봤으니 이를 시각화 해보도록 하자. 시각화의 편의성을 위해서 숫자형 데이터는 모두 라벨링을 해주었다. 대신 train_set_vis라고 train_set를 따로 지정해주어 혼선이 없도록 해주었다.

train_set_vis <- train_set
dim(train_set_vis)

## [1] 10886    14

train_set_vis$season  <- factor(train_set_vis$season, labels = c("Spring", "Summer", "Fall", "Winter"))
train_set_vis$weather <- factor(train_set_vis$weather, labels = c("Good", "Normal", "Bad", "Very Bad"))
train_set_vis$holiday <- factor(train_set_vis$holiday)
train_set_vis$workingday <- factor(train_set_vis$workingday)
train_set_vis$year <- factor(train_set_vis$year)
train_set_vis$month <- factor(train_set_vis$month)
train_set_vis$wday <- factor(train_set_vis$wday, labels = c("Sun","Mon", "Tue","Wed","Thu","Fir","Sat"))

각 변수별 count와의 관계

다음 그림은 각 변수가 count에 어떤 관계를 보이는지 알려주고 있다. 분석을 하기전에 우리가 예측을 해야하는 것은 count이기 때문에 count를 기준으로 각 변수의 관계를 그렸으며, 한눈에 보기위해 grid 함수를 사용하였다. 자세한 내용은 다음 참고를 보자.

참고: map in tydiverse (lapply 보다 빠르게 병렬 처리를 해줄수 있다.) 참고: grid.arrange in r (ggplot을 grid화 시켜서 표현을 할 수 있다.)

# count 칼럼이 10번째인데 이것을 제외시켰음. which()를 사용하지 않았을 때는 논리연산으로 수행되어 TRUE FALSE 로 구분하고 FALSE 로 뜬다.
non_hour_list <- (colnames(train_set_vis) != "count") %>% # which() 를 사용해 "count"인 인덱스를 제외.
  which()  

# 병렬처리로 non_hour_list에 function을 적용한다. fuction은 첫번째, train_set_vis의 컬럼을 df_list에 
lst <- map(non_hour_list, function(i) {
  df_list <- colnames(train_set_vis)[i]

  train_set_vis %>%
    select(df_list, count) %>%
    rename(aa = df_list) %>%
    ggplot(aes(aa,count)) +
    geom_point(alpha=.2,color = "#008ABC") +
    labs(title = paste0(df_list," vs count"), x = df_list, y = "",color=df_list) +
    theme_bw() +
    theme(legend.position = "bottom")
})

grid.arrange(grobs=lst, ncol=2)

  # df_list <- colnames(train_set_vis)[1]
  # 
  # train_set_vis %>%
  #   select(df_list, count) %>%
  #   rename(aa = df_list) %>%
  #   ggplot(aes(aa,count)) +
  #   geom_point(alpha=.2,color = "#008ABC") +
  #   labs(title = paste0(df_list," vs count"), x = df_list, y = "",color=df_list) +
  #   theme_bw() +
  #   theme(legend.position = "bottom")
  # 
  # 
  #   train_set_vis %>%
  #   ggplot(aes(year,count)) +
  #   geom_point(alpha=.2,color = "#008ABC") +
  #   labs(title = paste0(df_list," vs count"), x = df_list, y = "",color=df_list) +
  #   theme_bw() +
  #   theme(legend.position = "bottom")

hour와 다른 변수들 간의 관계

우리가 count를 예측하기 위해 다른 변수들과의 관계를 보니, hour과 count가 가장 영향을 끼치는 것으로 확인되었다. 지난 번 블로그에서 서울시 유동인구 데이터 분석을 해보니, 이러한 결과를 예상했었으나, 그 때와는 다르게 여기서는 시간에 적당한 데이터가 포함되어 있는 것을 확인할 수 있었다. 따라서 hour와 다른 변수들은 어떤 관계를 가지고 있는지 보고자 했다.

season: 봄, 여름, 가을, 겨울 중 봄, 가을이 가장 높게 나올줄 알았으나, 봄이 가장 낮게 나왔다.
holiday: 0과 1이 뚜렸하게 구분되어 있다.
workingday: 역시 0과 1이 뚜렸하게 구분되어 있다.
weather: 날씨 역시 계절과 같이 뚜렸하게 차이를 보일 것으로 예측했었다. 특히 날씨가 좋을 때 자전거 수요량이 늘어났음을 확인할 수 있었다.
year: 11년도에 비해 12년도에 더 높아져 있었다. 이 추세로 13년도가 더 높아질 것을 예상했으나 성급한 추측인 것 같다.
month: 월별로는 뚜렸한 특성을 찾을 수 없었다.
wday: 평일과 주말의 구분이 확연하게 확인된다. 주말은 나들이 하기 좋은시간인 14쯤이 피크였고, 평일은 오전 8시경과 오후 16경이 피크였다.

factor_list <- sapply(train_set_vis, is.factor) %>% 
  which()


lst <- lapply(factor_list, function(x) {
  df_list <- colnames(train_set_vis)[x]
  
  train_set_vis %>% 
    rename(aa = df_list) %>% 
    group_by(aa, hour) %>% 
    summarise(count = sum(count)) %>% 
    ggplot(aes( x = hour, y = count, group = aa, colour = aa)) +
    labs(title = paste0("count by ", df_list), color = df_list) +
    theme_bw() +
    geom_line()
})





grid.arrange(grobs=lst, ncol=2)

XGboost

램덤포레스트, 의사결정 모델 등등 여러모델이 있지만 이번에는 xgboost를 사용해볼 것이다. xgboost는 데이터 프레임이 아니라 matrix형식이 되어야 한다.

count to log & train/test set 분리

RMSLE 값을 최소화하기 위해서 Count 변수에 Log를 적용할 것이다. 날짜 데이터는 년, 월, 시간, 요일 등으로 구분하였으니 모델링 작업 시 제외하자.

train_set$count <- log1p(train_set$count)

x_train <- train_set %>% 
  select(-count, -datetime) %>% 
  as.matrix()

y_train <- train_set$count 

x_test <- test_set %>% 
  select(-datetime) %>% 
  as.matrix()

gridsearch / cross validation

Xgboost 를 사용하는데 할 때, 모델에 다양한 조건부여하는데 있어 Gird search와 Cross validation 을 이용할 것이다. 검증 단계에서만 좋은 결괏값이 나오고 실 성능은 그렇지 않을수도 있기 때문에 Cross validation을 적용하였다.

dtrain <- xgb.DMatrix(x_train, label = y_train)


searchGridSubCol <- expand.grid(subsample = c(0.5, 0.6),
                                colsample_bytree = c(0.5, 0.6),
                                max_depth = c(7:15),
                                min_child = seq(1),
                                eta = c(0.05,0.1,0.15)
) # 2 * 2* 9 * 1 * 3 * 150 = ?

ntrees <- 10

rmseErrorsHyperparameters <- apply(searchGridSubCol, 1, function(parameterList){

  #Extract Parameters to test
  currentSubsampleRate <- parameterList[["subsample"]]
  currentColsampleRate <- parameterList[["colsample_bytree"]]
  currentDepth <- parameterList[["max_depth"]]
  currentEta <- parameterList[["eta"]]
  currentMinChild <- parameterList[["min_child"]]

  xgboostModelCV <- xgb.cv(data =  dtrain, nrounds = ntrees, nfold = 5, showsd = TRUE,
                           metrics = "rmse", verbose = TRUE, "eval_metric" = "rmse",
                           "objective" = "reg:linear", "max.depth" = currentDepth, "eta" = currentEta,
                           "subsample" = currentSubsampleRate, "colsample_bytree" = currentColsampleRate
                           , print_every_n = 10, "min_child_weight" = currentMinChild, booster = "gbtree",
                           early_stopping_rounds = 10)

  xvalidationScores <- as.data.frame(xgboostModelCV$evaluation_log)
  rmse <- tail(xvalidationScores$test_rmse_mean, 1)
  trmse <- tail(xvalidationScores$train_rmse_mean,1)
  output <- return(c(rmse, trmse,currentSubsampleRate, currentColsampleRate, currentDepth, currentEta, currentMinChild))})


output <- as.data.frame(t(rmseErrorsHyperparameters))
varnames <- c("TestRMSE", "TrainRMSE", "SubSampRate", "ColSampRate", "Depth", "eta", "currentMinChild")
names(output) <- varnames

gridsearch / cross validation 결과

tail(output)

##     TestRMSE TrainRMSE SubSampRate ColSampRate Depth  eta currentMinChild
## 103 1.034230 1.0035410         0.5         0.6    14 0.15               1
## 104 1.018344 0.9885732         0.6         0.6    14 0.15               1
## 105 1.100907 1.0602544         0.5         0.5    15 0.15               1
## 106 1.079219 1.0424332         0.6         0.5    15 0.15               1
## 107 1.020200 0.9935810         0.5         0.6    15 0.15               1
## 108 1.044239 1.0089566         0.6         0.6    15 0.15               1

xgboost 모델

model <- xgb.train(
  data = dtrain, 
  max_depth = 10,
  nround = 150,
  eta = 0.15,
  subsample = 0.6,
  colsample_bytree = 0.6,
  min_child_wight = 1
)

xgb.importance(feature_names = colnames(x_train), model) %>% 
  xgb.plot.importance()

pred <- predict(model, x_test) %>% 
  expm1()


solution = data.frame(datetime = test_set$datetime, count = pred)
write.csv(solution, "solution.csv", row.names = FALSE)

출처 : https://unfinishedgod.github.io/docs/kaggle/bike_sharing_demand/Bike_Sharing_Demand.html

자전거 대여 수요예측

Ryu

2020-03-23

개요