Egzamin: Warsztaty Machine Learning w R

Wprowadzenie

W ramach egzaminu wykonano przykładową analizę dostarczonego zbioru danych dotyczącego czasowego wypożaczania rowerów w pewnym mieście w celu uzyskania prognozy godzinowego popytu na tę usługę.

Analizę podzielono na kilka etapów:

Eksploracyjna analiza danych (EDA) w celu uporządkowania zmiennych
Konstrukcja wybranych modeli wraz z walidacją krzyżową i tuningiem:
- regresja liniowa,
- uogólniony model addytywny GAM,
- las losowy,
- XGBoost,
Wykonanie prognozy na skonstruowanych modelach i porównanie wyników.

Opis bazy danych

Nazwa zmiennej	Opis
datetime	hourly date + timestamp
season	1 = spring,
	2 = summer,
	3 = fall,
	4 = winter
holiday	whether the day is considered a holiday
workingday	whether the day is neither a weekend nor holiday
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	temperature in Celsius
atemp	“feels like” temperature in Celsius
humidity	relative humidity
windspeed	wind speed
casual	number of non-registered user rentals initiated
registered	number of registered user rentals initiated
count	number of total rentals

Eksploracyjna analiza i przygotowanie danych

Wstępne rozpoznanie danych

train <- read.csv("E:/R2/Egzamin/train.csv")
head(train)

Zbiór train.csv posiada 12 zmiennych oraz 10886 obserwacji. Wśród tych zmiennych znajdują się dane:

numeryczne (temp, atemp, humidity, windspeed, casual, registered, count),
kategoryczne (season, weather)
w tym binarne (holiday, workingday)
oraz data (datetime).

Usunięcie niepotrzebnych kolumn

Zacznijmy od selekcji niepotrzebnych kolumn, czyli zmiennych casual i registered dających razem wartość zmiennej objaśnianej count.

train <- dplyr::select(train, -c(casual, registered))

Wyciągnięcie danych z kolumny datetime

Rozłóżmy zmienną datetime na lata, miesiące, dni tygodnia oraz godziny.

train %>% 
  mutate(datetime = fastPOSIXct(datetime, "GMT")) %>% 
  mutate(hour = hour(datetime),
         month = month(datetime),
         year = year(datetime),
         weekday = wday(datetime)) -> train
head(train)

Usunięcie braków danych

summary(train)

##     datetime                       season         holiday       
##  Min.   :2011-01-01 00:00:00   Min.   :1.000   Min.   :0.00000  
##  1st Qu.:2011-07-02 07:15:00   1st Qu.:2.000   1st Qu.:0.00000  
##  Median :2012-01-01 20:30:00   Median :3.000   Median :0.00000  
##  Mean   :2011-12-27 05:56:22   Mean   :2.507   Mean   :0.02857  
##  3rd Qu.:2012-07-01 12:45:00   3rd Qu.:4.000   3rd Qu.:0.00000  
##  Max.   :2012-12-19 23:00:00   Max.   :4.000   Max.   :1.00000  
##    workingday        weather           temp           atemp      
##  Min.   :0.0000   Min.   :1.000   Min.   : 0.82   Min.   : 0.76  
##  1st Qu.:0.0000   1st Qu.:1.000   1st Qu.:13.94   1st Qu.:16.66  
##  Median :1.0000   Median :1.000   Median :20.50   Median :24.24  
##  Mean   :0.6809   Mean   :1.418   Mean   :20.23   Mean   :23.66  
##  3rd Qu.:1.0000   3rd Qu.:2.000   3rd Qu.:26.24   3rd Qu.:31.06  
##  Max.   :1.0000   Max.   :4.000   Max.   :41.00   Max.   :45.45  
##     humidity        windspeed          count            hour      
##  Min.   :  0.00   Min.   : 0.000   Min.   :  1.0   Min.   : 0.00  
##  1st Qu.: 47.00   1st Qu.: 7.002   1st Qu.: 42.0   1st Qu.: 6.00  
##  Median : 62.00   Median :12.998   Median :145.0   Median :12.00  
##  Mean   : 61.89   Mean   :12.799   Mean   :191.6   Mean   :11.54  
##  3rd Qu.: 77.00   3rd Qu.:16.998   3rd Qu.:284.0   3rd Qu.:18.00  
##  Max.   :100.00   Max.   :56.997   Max.   :977.0   Max.   :23.00  
##      month             year         weekday     
##  Min.   : 1.000   Min.   :2011   Min.   :1.000  
##  1st Qu.: 4.000   1st Qu.:2011   1st Qu.:2.000  
##  Median : 7.000   Median :2012   Median :4.000  
##  Mean   : 6.521   Mean   :2012   Mean   :3.999  
##  3rd Qu.:10.000   3rd Qu.:2012   3rd Qu.:6.000  
##  Max.   :12.000   Max.   :2012   Max.   :7.000

sum(is.na(train))

## [1] 0

W analizowanym zbiorze nie ma braków danych - nie jest konieczna imputacja.

Rozkłady zmiennych

Zróbmy szybkie histogramy wszystkich zmiennych:

train %>%
  dplyr::select(-c(datetime)) %>% 
  gather() %>% 
  ggplot(aes(value)) +
  facet_wrap(~ key, scales = "free") +
  geom_histogram() +
  theme_minimal()

lillie.test(train$count)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  train$count
## D = 0.14639, p-value < 2.2e-16

lillie.test(train$windspeed)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  train$windspeed
## D = 0.082833, p-value < 2.2e-16

Widzimy, że zmienna objaśniana count posiada rozkład prawoskośny i zgodnie z testem Lillieforsa nie jest on zgodny z rozkładem normalnym (podobnie jak rozkłady zmiennych objaśniających, np. windspeed).

Korelacje między zmiennymi

Sprawdźmy występowanie zależności między zmiennymi.

ggcorrplot(cor(train[,2:14]), hc.order = TRUE, type = "lower",
           lab = TRUE)

Między dwoma parami zmiennych występują bardzo wysokie korelacje:

temp i atemp (r=0.98)
season i month (r=0.97)

ggplot(train, aes(x=temp, y=atemp)) +
  geom_point() +
  geom_smooth(method="lm") +
  theme_minimal()

Prosta analiza danych na podstawie wykresów

Sporządźmy prostą analizę danych w oparciu o wykresy, aby zobaczyć jak się one zachowują względem siebie.

train$season <- factor(train$season, levels=c(1,2,3,4), 
                       labels=c("Wiosna", "Lato", "Jesień", "Zima"))
train$weather <- factor(train$weather, levels=c(1,2,3,4), 
                        labels=c("Bardzo dobra", "Dobra", "Średnia", "Zła"))
train$weekday <- factor(train$weekday, levels=c(1,2,3,4,5,6,7), 
                        labels=c("Poniedziałek", "Wtorek", "Środa", "Czwartek", "Piątek", "Sobota", "Niedziela"))
train %>% 
  group_by(season, hour) %>% 
  summarise(count = mean(count)) -> season_count

train %>% 
  group_by(weather, hour) %>% 
  summarise(count = mean(count)) -> weather_count

train %>% 
  group_by(weekday, hour) %>% 
  summarise(count = mean(count)) -> weekday_count

ggplot(train, aes(x=hour, y=count, color=season)) +
 geom_point(data = season_count, aes(group = season)) +
 geom_line(data = season_count, aes(group = season)) +
  scale_color_discrete(name="Pora roku")+
  xlab("Godzina") +
  ylab("Ilość wypożyczeń") +
 ggtitle("Wypożyczenia rowerów z podziałem na pory roku") + 
  theme_minimal()

ggplot(train, aes(x=hour, y=count, color=weather)) +
  geom_point(data = weather_count, aes(group = weather)) +
  geom_line(data = weather_count, aes(group = weather)) +
  scale_color_discrete(name="Pogoda")+
  xlab("Godzina") +
  ylab("Ilość wypożyczeń") +
  ggtitle("Wypożyczenia rowerów z podziałem na sytuację pogodową") + 
  theme_minimal()

ggplot(train, aes(x=hour, y=count, color=weekday)) +
  geom_point(data = weekday_count, aes(group = weekday)) +
  geom_line(data = weekday_count, aes(group = weekday)) +
  scale_color_discrete(name="Dzień tygodnia") +
  xlab("Godzina") +
  ylab("Ilość wypożyczeń") +
  ggtitle("Wypożyczenia rowerów z podziałem na dzień tygodnia") + 
  theme_minimal()

Większość spostrzeżeń ujawniona na podstawie powyższych wykresów jest intuicyjna:

Najwięcej rowerów wypożyczanych jest w godzinach porannych (między 7:00 a 9:00) oraz popołudniowych (między 16:00 a 18:00). Odpowiada to obserwowanym szczytom ruchu komunikacyjnego i związane jest z dojazdem do i z pracy.
W dniach wolnych od pracy (Sobota, Niedziela) trend ilości wypożyczeń wygląda zupełnie inaczej, wzrasta on w ciągu dnia co zapewne jest spowodowane przejażdżkami rekreacyjnymi.
Najwięcej rowerów wypożyczanych jest podczas dobrej i bardzo dobrej pogody, a w czasie złej pogody bardzo rzadko.
Najmniej rowerów wypożyczanych jest wiosną, a w pozostałych trzech porach roku trend jest podobny - mniej intuicyjna obserwacja.

Przejdźmy do modelowania.

Regresja liniowa

Przygotowanie danych

Wykonajmy podział próby na zbiór treningowy i walidacyjny, tak aby w tym drugim znajdowało się 20% najnowszych danych. Usuńmy też zmienną datetime.

train1 <- head(train[,-1], 8709)
valid1 <- tail(train[-1], 2177)

Konfiguracja scenariusza kroswalidacji

cross.walid <- trainControl(method = "repeatedcv",
                            number = 10,
                            repeats = 10, 
                            returnResamp = 'all')

Trening modelu

model_lm <- train(count ~ .,
               data = train1,
               trControl = cross.walid,
               method = 'lm') 
model_lm$results

Wyniki nie przedstawiają się za dobrze, uzyskaliśmy wartość R² na poziomie 0.4055689.

Kolejne próby pokazały, że dodawanie interakcji lub usuwanie zmiennych silnie ze sobą skorelowanych nie wpłynęło na poprawę skuteczności modelu. Spróbujmy natomiast wytrenować model na zlogarytmowanej zmiennej count.

model_lm <- train(log(count) ~ .,
               data = train1,
               trControl = cross.walid,
               method = 'lm') 
model_lm$results

Nieznacznie poprawiło to nasze wyniki (w odniesieniu do wartości R²), ale jednoznacznie wskazuje to, że nasze dane nie są odpowiednie do regresji liniowej. Niemniej jednak, sprawdźmy prognozę.

Prognoza na zbiorze walidacyjnym

pred_lm <- predict(model_lm, valid1)

ggplot(data=valid1, aes(x=count, y=exp(pred_lm)))+
  geom_point() +
  geom_smooth(method="lm", formula = y~x-1) +
  stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               label.x.npc = "right", label.y.npc = 0,
               formula = y~x-1, parse = TRUE, size = 5) +
  xlab("Ilość wypożyczeń (Obserwacja)") +
  ylab("Ilość wypożyczeń (Prognoza)") +
  ggtitle("Wykres rozrzutu prognozy od obserwacji dla regresji liniowej") +
  theme_minimal()

Widzimy, że rozrzut wartości jest bardzo duży, a wartości zmiennej count są średnio niemal dwukrotnie niedoszacowane. Błędy na zbiorze walidacyjnym wynoszą kolejno:

MAE = 143.959048
MAPE = 1.0582885
RMSE = 206.0077843

Model GAM

Przygotowanie danych

Już wiemy, że dla zlogarytmowanej zmiennej objaśnianej count otrzymamy lepsze wyniki, a także unikniemy wartości ujemnych w prognozie. Postąpmy podobnie jak w przypadku regresji liniowej i przeprowadźmy trening modelu na zmiennej log(count).

Konfiguracja scenariusza kroswalidacji

cross.walid <- trainControl(method = "repeatedcv",
                            number = 10,
                            repeats = 10, 
                            returnResamp = 'all')

Trening modelu

model_gam <- train(log(count) ~ .,
               data = train1,
               trControl = cross.walid,
               method = 'gam') 
model_gam$results

Widzimy, że model GAM lepiej wpasowuje się w nasze dane niż zwykła regresja liniowa.

Prognoza na zbiorze walidacyjnym

pred_gam <- predict(model_gam, valid1)

ggplot(data=valid1, aes(x=count, y=exp(pred_gam)))+
  geom_point() +
  geom_smooth(method="lm", formula = y~x-1) +
  stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               label.x.npc = "right", label.y.npc = 0,
               formula = y~x-1, parse = TRUE, size = 5) +
  xlab("Ilość wypożyczeń (Obserwacja)") +
  ylab("Ilość wypożyczeń (Prognoza)") +
  ggtitle("Wykres rozrzutu prognozy od obserwacji dla modelu GAM") +
  theme_minimal()

Porównanie prognoz i obserwacji zmiennej count na zbiorze walidacyjnym pokazuje, że model lepiej radzi sobie z przewidywaniem ilości wypożyczeń rowerów. Wciąż jednak rozrzut wartości jest dosyć duży. Błędy na zbiorze walidacyjnym wynoszą:

MAE = 89.4395411
MAPE = 0.4755623
RMSE = 128.8581528

Las losowy

Przygotowanie danych

Wykonajmy przykładową binaryzację zmiennych season, weather oraz weekday na zbiorze treningowym i walidacyjnym. Niestety przeprowadzone próby pokazały, że nie przyniosło to dużych korzyści w otrzymanych modelach, a jedynie wydłużyło czas obliczeń.

train1b <- mlr::createDummyFeatures(train1)
valid1b <- mlr::createDummyFeatures(valid1)

Konfiguracja scenariusza kroswalidacji

cross.walid <- trainControl(method = "cv",
                            number = 10,
                            allowParallel = T,
                            returnResamp = 'all')

Trening modelu z tuningiem

Spróbujmy znaleźć optymalną liczbę kolumn.

tunegrid <- expand.grid(.mtry = c(1:24))

cluster <- makeCluster(detectCores() - 1)
registerDoParallel(cluster)

model_rf <- train(count ~ ., 
                data = train1b, 
                method = "rf", 
                replace = T,
                ntree = 200,
                metric = 'RMSE',
                tuneGrid = tunegrid, 
                trControl = cross.walid)
model_rf

## Random Forest 
## 
## 8709 samples
##   24 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 7837, 7839, 7840, 7838, 7839, 7838, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE       Rsquared   MAE      
##    1    139.25526  0.5162304  106.44310
##    2    107.99003  0.6590996   79.02187
##    3     87.77159  0.7556748   61.19685
##    4     75.22961  0.8143507   50.77985
##    5     67.01472  0.8505674   44.34606
##    6     61.03256  0.8753839   39.95735
##    7     56.06076  0.8941737   36.39289
##    8     52.30063  0.9073482   33.76010
##    9     48.75027  0.9191907   31.43853
##   10     46.07337  0.9273778   29.70209
##   11     44.36187  0.9321519   28.48868
##   12     42.63571  0.9370466   27.26125
##   13     41.25952  0.9406867   26.30997
##   14     40.27014  0.9430265   25.56375
##   15     39.25838  0.9456893   24.94056
##   16     39.06184  0.9460205   24.74073
##   17     38.54119  0.9471750   24.33353
##   18     38.43057  0.9473751   24.09964
##   19     38.20620  0.9478116   23.86731
##   20     38.01083  0.9482645   23.78760
##   21     37.97553  0.9482635   23.68383
##   22     38.20068  0.9475364   23.73693
##   23     38.13572  0.9476584   23.55963
##   24     38.37799  0.9469454   23.63398
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 21.

stopCluster(cluster)
registerDoSEQ() 

plot(model_rf)

model_rf$finalModel

## 
## Call:
##  randomForest(x = x, y = y, ntree = 200, mtry = param$mtry, replace = ..1) 
##                Type of random forest: regression
##                      Number of trees: 200
## No. of variables tried at each split: 21
## 
##           Mean of squared residuals: 1419.971
##                     % Var explained: 94.87

Osiągamy wartości R² przy optymalnej ilości kolumn równej 21. Sprawdźmy jak to wygląda na wykresie prognoz.

Prognoza na zbiorze walidacyjnym

pred_rf <- predict(model_rf, valid1b)

ggplot(data=valid1b, aes(x=count, y=pred_rf))+
  geom_point() +
  geom_smooth(method="lm", formula = y~x-1) +
  stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               label.x.npc = "right", label.y.npc = 0,
               formula = y~x-1, parse = TRUE, size = 5) +
  xlab("Ilość wypożyczeń (Obserwacja)") +
  ylab("Ilość wypożyczeń (Prognoza)") +
  ggtitle("Wykres rozrzutu prognozy od obserwacji dla lasu losowego") +
  theme_minimal()

Na wykresie rozrzutu widoczne jest o wiele lepsze dopasowanie prognoz do obserwacji w lasie losowym w stosunku do poprzednich metod. Błędy modelu na zbiorze walidacyjnym wynoszą:

MAE = 49.1839593
MAPE = 0.2864133
RMSE = 76.2186923

XGBoost

Przygotowanie danych

Skorzystajmy z tych samych zbiorów danych co przy lesie losowym. Rozdzielmy dane na zmienne objaśniające i objaśnianą, zamieńmy na obiekt xgb.DMatrix.

train1b_X <- as.matrix(train1b[, -7])
train1b_Y <- as.numeric(train1b[,7])
dtrain <- xgb.DMatrix(data = train1b_X, label = train1b_Y)

valid1b_X <- as.matrix(valid1b[, -7])
valid1b_Y <- as.numeric(valid1b[,7])
dvalid <- xgb.DMatrix(data = valid1b_X, label = valid1b_Y)

Trening modelu

model_xgb = xgb.train(data = dtrain, nround = 150, max_depth = 5, eta = 0.1, subsample = 0.9)

xgb.importance(feature_names = colnames(train1b_X), model_xgb) %>% xgb.plot.importance()

Widzimy, że zdecydowanie najważniejszą zmienną objaśniającą jest zmienna hour i to ona w największym stopniu wpływa na zmienną objaśnianą. Ważny jest też dzień tygodnia, czy warunki pogodowe (temperatura, wilgotność).

Tuning hiperparametrów (random search)

cluster <- makeCluster(detectCores() - 1)
registerDoParallel(cluster)

cross.walid <- trainControl(method = "cv",
                        number = 10,
                        search = "random",
                        allowParallel = TRUE)

model_xgb <- caret::train(x = train1b_X,
                          y = train1b_Y,
                          trControl = cross.walid,
                          tuneLength = 50,
                          method = "xgbTree")

stopCluster(cluster)
registerDoSEQ()

model_xgb$finalModel

## ##### xgb.Booster
## raw: 3.6 Mb 
## call:
##   xgboost::xgb.train(params = list(eta = param$eta, max_depth = param$max_depth, 
##     gamma = param$gamma, colsample_bytree = param$colsample_bytree, 
##     min_child_weight = param$min_child_weight, subsample = param$subsample), 
##     data = x, nrounds = param$nrounds, objective = "reg:linear")
## params (as set within xgb.train):
##   eta = "0.0417642481841613", max_depth = "8", gamma = "5.21795803215355", colsample_bytree = "0.618658572621644", min_child_weight = "20", subsample = "0.94847874797415", objective = "reg:linear", silent = "1"
## xgb.attributes:
##   niter
## callbacks:
##   cb.print.evaluation(period = print_every_n)
## # of features: 24 
## niter: 632
## nfeatures : 24 
## xNames : holiday workingday temp atemp humidity windspeed hour month year season.Wiosna season.Lato season.Jesień season.Zima weather.Bardzo.dobra weather.Dobra weather.Średnia weather.Zła weekday.Poniedziałek weekday.Wtorek weekday.Środa weekday.Czwartek weekday.Piątek weekday.Sobota weekday.Niedziela 
## problemType : Regression 
## tuneValue :
##    nrounds max_depth        eta    gamma colsample_bytree min_child_weight
## 6     632         8 0.04176425 5.217958        0.6186586               20
##   subsample
## 6 0.9484787
## obsLevels : NA 
## param :
##  list()

Prognoza na zbiorze walidacyjnym

pred_xgb <- predict(model_xgb, valid1b_X)

ggplot(data=valid1b, aes(x=count, y=pred_xgb))+
  geom_point() +
  geom_smooth(method="lm", formula = y~x-1) +
  stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
               label.x.npc = "right", label.y.npc = 0,
               formula = y~x-1, parse = TRUE, size = 5) +
  xlab("Ilość wypożyczeń (Obserwacja)") +
  ylab("Ilość wypożyczeń (Prognoza)") +
  ggtitle("Wykres rozrzutu prognozy od obserwacji dla XGBoost") +
  theme_minimal()

Na pewno zauważalne jest bardzo dobre dopasowanie prognoz. W stosunku do lasu losowego nastąpiło zmniejszenie błędu RMSE. Wartości błędów na zbiorze walidacyjnym wynoszą tutaj:

MAE = 44.9797052
MAPE = 0.2919877
RMSE = 66.0639642

Prognoza na zbiorze testowym i porównanie wyników

Wczytanie i przygotowanie zbioru testowego

test <- read.csv("E:/R2/Egzamin/test.csv")

test %>% 
  mutate(datetime = fastPOSIXct(datetime, "GMT")) %>% 
  mutate(hour = hour(datetime),
         month = month(datetime),
         year = year(datetime),
         weekday = wday(datetime)) -> test

test$season <- factor(test$season, levels=c(1,2,3,4), 
                       labels=c("Wiosna", "Lato", "Jesień", "Zima"))
test$weather <- factor(test$weather, levels=c(1,2,3,4), 
                        labels=c("Bardzo dobra", "Dobra", "Średnia", "Zła"))
test$weekday <- factor(test$weekday, levels=c(1,2,3,4,5,6,7), 
                        labels=c("Poniedziałek", "Wtorek", "Środa", "Czwartek", "Piątek", "Sobota", "Niedziela"))

testb <- mlr::createDummyFeatures(test)

Policzenie prognozy dla każdego modelu

test$count_lm<- exp(predict(model_lm, test))
test$count_gam<- exp(predict(model_gam, test))
test$count_rf<- predict(model_rf, testb)
test$count_xgb<- predict(model_xgb, testb)

Wykres porównawczy

Połączmy zestaw danych train oraz test i pokażmy na wykresie trendu wybrany okres czasu (tutaj od 15 marca 2011 do 5 kwietnia 2011), aby zobaczyć, która prognoza wygląda najbardziej wiarygodnie.

Najpierw zobaczmy pojedyncze wyniki, a następnie zbiorcze (żeby ułatwić interpretację wykresu z wieloma liniami).

train <- read.csv("E:/R2/Egzamin/train.csv")
train <- dplyr::select(train, -c(casual, registered))
train %>% 
  mutate(datetime = fastPOSIXct(datetime, "GMT")) %>% 
  mutate(hour = hour(datetime),
         month = month(datetime),
         year = year(datetime),
         weekday = wday(datetime)) -> train

train$season <- factor(train$season, levels=c(1,2,3,4), 
                       labels=c("Wiosna", "Lato", "Jesień", "Zima"))
train$weather <- factor(train$weather, levels=c(1,2,3,4), 
                        labels=c("Bardzo dobra", "Dobra", "Średnia", "Zła"))
train$weekday <- factor(train$weekday, levels=c(1,2,3,4,5,6,7), 
                        labels=c("Poniedziałek", "Wtorek", "Środa", "Czwartek", "Piątek", "Sobota", "Niedziela"))

data <- merge(train, test, all = T)

ggplot(data = data) +
  geom_path(aes(x = datetime, y = count), color = "grey", size = 1.5) + 
  geom_path(aes(x = datetime, y = count_lm), color = "red", size = 1.5, alpha=0.8) +
  scale_x_datetime(limits = ymd_h(c("2011-03-15 00", "2011-04-05 23"))) +
  scale_y_continuous(limits=c(0,400)) +
    ggtitle("Wykres trendu ilości wypożyczeń rowerów z prognozą regresji liniowej") +
  theme_minimal()

ggplot(data = data) +
  geom_path(aes(x = datetime, y = count), color = "grey", size = 1.5) + 
  geom_path(aes(x = datetime, y = count_gam), color = "blue", size = 1.5, alpha=0.8) +
  scale_x_datetime(limits = ymd_h(c("2011-03-15 00", "2011-04-05 23"))) +
  scale_y_continuous(limits=c(0,400)) +
  ggtitle("Wykres trendu ilości wypożyczeń rowerów z prognozą modelu GAM") +
  theme_minimal()

ggplot(data = data) +
  geom_path(aes(x = datetime, y = count), color = "grey", size = 1.5) + 
  geom_path(aes(x = datetime, y = count_rf), color = "green", size = 1.5) +
  scale_x_datetime(limits = ymd_h(c("2011-03-15 00", "2011-04-05 23"))) +
  scale_y_continuous(limits=c(0,400)) +
  ggtitle("Wykres trendu ilości wypożyczeń rowerów z prognozą lasu losowego") +
  theme_minimal()

ggplot(data = data) +
  geom_path(aes(x = datetime, y = count), color = "grey", size = 1.5) + 
  geom_path(aes(x = datetime, y = count_xgb), color= "orange", size = 1.5, alpha=0.8) +
  scale_x_datetime(limits = ymd_h(c("2011-03-15 00", "2011-04-05 23"))) +
    ggtitle("Wykres trendu ilości wypożyczeń rowerów z prognozą modelu XGBoost") +
  scale_y_continuous(limits=c(0,400)) +
  theme_minimal()

ggplot(data = data) +
  geom_path(aes(x = datetime, y = count), color = "grey", size = 1.5) + 
  geom_path(aes(x = datetime, y = count_rf), color = "green", size = 1.5, alpha=0.8) +
  geom_path(aes(x = datetime, y = count_xgb), color= "orange", size = 1.5, alpha=0.8) +
  geom_path(aes(x = datetime, y = count_lm), color = "red", size = 1.5, alpha=0.8) +
  geom_path(aes(x = datetime, y = count_gam), color = "blue", size = 1.5, alpha=0.8) +
  scale_x_datetime(limits = ymd_h(c("2011-03-15 00", "2011-04-05 23"))) +
  scale_y_continuous(limits=c(0,400)) +
  theme_minimal()

Kolorem szarym zaznaczono wartości zmiennej count w zbiorze treningowym. Legenda dla prognoz zbioru testowego:

regresja liniowa
model GAM
las losowy
XGBoost

Można zauważyć, że najlepiej w trend danych wpasowują się wyniki lasu losowego oraz XGBoost, przy czym w drugim przypadku występują braki danych. Model regresji liniowej oraz GAM zaniżają wartości zmiennej count (co było również widoczne na wykresach rozrzutu). Model GAM dodatkowo ma problem z przewidywaniem wartości minimalnych ilości wypożyczeń rowerów.

Podsumujmy tabelą błędów otrzymanych na zbiorze walidacyjnym dla poszczególnych modeli.

stats_lm <- as.vector(c(MAE(valid1$count, exp(pred_lm)), MAPE(valid1$count, exp(pred_lm)), RMSE(valid1$count, exp(pred_lm))))
stats_gam <- as.vector(c(MAE(valid1$count, exp(pred_gam)), MAPE(valid1$count, exp(pred_gam)), RMSE(valid1$count, exp(pred_gam))))
stats_rf <- as.vector(c(MAE(valid1b$count, pred_rf), MAPE(valid1b$count, pred_rf), RMSE(valid1b$count, pred_rf)))
stats_xgb <- as.vector(c(MAE(valid1b$count, pred_xgb), MAPE(valid1b$count, pred_xgb), RMSE(valid1b$count, pred_xgb)))

stats <- rbind(stats_lm, stats_gam, stats_rf, stats_xgb)
rownames(stats) <- c("LM", "GAM", "RF", "XGB")
colnames(stats) <- c("MAE", "MAPE", "RMSE")
head(stats)

##           MAE      MAPE      RMSE
## LM  143.95905 1.0582885 206.00778
## GAM  89.43954 0.4755623 128.85815
## RF   49.18396 0.2864133  76.21869
## XGB  44.97971 0.2919877  66.06396

Spośród omawianych, zdecydowanie najlepiej wypadł las losowy oraz XGBoost. W wypadku tego drugiego otrzymano nawet nieco mniejsze wartości błędów i lepsze dopasowanie wartości prognoz do obserwacji.

Egzamin: Warsztaty Machine Learning w R

Adriana Szulecka, szulecka@agh.edu.pl

24 luty, 2019

Wprowadzenie

Opis bazy danych

Eksploracyjna analiza i przygotowanie danych

Wstępne rozpoznanie danych

Usunięcie niepotrzebnych kolumn

Wyciągnięcie danych z kolumny datetime

Usunięcie braków danych

Rozkłady zmiennych

Korelacje między zmiennymi

Prosta analiza danych na podstawie wykresów

Regresja liniowa

Przygotowanie danych

Konfiguracja scenariusza kroswalidacji

Trening modelu

Prognoza na zbiorze walidacyjnym

Model GAM

Przygotowanie danych

Konfiguracja scenariusza kroswalidacji

Trening modelu

Prognoza na zbiorze walidacyjnym

Las losowy

Przygotowanie danych

Konfiguracja scenariusza kroswalidacji

Trening modelu z tuningiem

Prognoza na zbiorze walidacyjnym

XGBoost

Przygotowanie danych

Trening modelu

Tuning hiperparametrów (random search)

Prognoza na zbiorze walidacyjnym

Prognoza na zbiorze testowym i porównanie wyników

Wczytanie i przygotowanie zbioru testowego

Policzenie prognozy dla każdego modelu

Wykres porównawczy