Walmart

---

Подгрузка пакетов

library(readr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(Metrics)
library(naniar)

Функционалы для временных рядов

Напомним, что наши прогнозы выглядят вот так: \[ \hat{y_i} = \hat{w_0} + \hat{w_1} \cdot x_i^{(1)} + \ldots + \hat{w_k} \cdot x_i^{k} \] Функционалы, которые принято использовать, когда вы работаете с временными рядами, выглядят так:

\[ MAE = \frac{1}{n} \sum_{i=1}^{n} |y_i - \hat y_i| \]

\[ MAPE = \frac{1}{n} \sum_{i=1}^{n} |\frac{y_i - \hat y_i}{y_i}| \]

Данные

Нам предоставлены исторические данные о продажах 45 магазинов Walmart, расположенных в разных регионах. В каждом магазине есть несколько отделов, и нам поручено прогнозировать продажи в каждом магазине для каждого отдела.

Зачем нам прогнозировать продажи?

• Нельзя привозить в магазин слишком мало товара, потребителям его может не хватить. Мало того, что они не принесут нам денег, так ещё и станут менее лояльными: “Не поедем в этот магазин. Там вечно ничего нет. Погнали лучше в Магнолию сходим. Там всё есть.”
• Нельзя привозить в магазин слишком много товара. Его нужно хранить. Это лишние расходы. Более того, товар может протухнуть. Придётся его списывать. Это тоже довольно неприятно.

walmart <- read_csv('data/walmart.csv')

Нам доступны следующие переменные:

• Weekly_Sales – объём продаж в данную неделю в данном отделе
• Store – номер магазина;
• Type – тип магазина;
• Size – размер магазина;
• Dept – номер отдела;
• Date – дата;
• IsHoliday – является ли неделя праздничной;
• Temperature – средняя температура в регионе в градусах по Фаренгейту;
• Fuel_Price – стоимость топлива в регионе;
• MarkDown1-5 – данные, связанные с рекламными уценками, которые запускает Walmart. Данные уценки доступны только после ноября 2011 года и доступны не для всех магазинов. Данные анонимизированы. Непонятно на какие именно товары производилась уценка и в каких количествах. Компании часто анонимизируют данные, когда выкладывают их в открытый доступ.
• CPI – индекс потребительских цен;
• Unemployment – уровень безработицы.
• next_Weekly_Sales – объём продаж в следующую неделю в данном отделе (target)

glimpse(walmart)

## Rows: 418,239
## Columns: 17
## $ Store             <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ Date              <date> 2010-02-05, 2010-02-12, 2010-02-19, 2010-02-26, 20…
## $ Dept              <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ Weekly_Sales      <dbl> 24924.50, 46039.49, 41595.55, 19403.54, 21827.90, 2…
## $ IsHoliday         <lgl> FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FAL…
## $ Type              <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "…
## $ Size              <dbl> 151315, 151315, 151315, 151315, 151315, 151315, 151…
## $ Temperature       <dbl> 42.31, 38.51, 39.93, 46.63, 46.50, 57.79, 54.58, 51…
## $ Fuel_Price        <dbl> 2.572, 2.548, 2.514, 2.561, 2.625, 2.667, 2.720, 2.…
## $ MarkDown1         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ MarkDown2         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ MarkDown3         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ MarkDown4         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ MarkDown5         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ CPI               <dbl> 211.0964, 211.2422, 211.2891, 211.3196, 211.3501, 2…
## $ Unemployment      <dbl> 8.106, 8.106, 8.106, 8.106, 8.106, 8.106, 8.106, 8.…
## $ next_Weekly_Sales <dbl> 46039.49, 41595.55, 19403.54, 21827.90, 21043.39, 2…

Изучение данных

Пропущенные значения

Давайте сначала разберемся с пропущенными значениями. Набор данных достаточно большой, поэтому не будет его рисовать. Посмотрим на статистики.

miss_var_summary(walmart)

## # A tibble: 17 x 3
##    variable          n_miss pct_miss
##    <chr>              <int>    <dbl>
##  1 MarkDown2         310095     74.1
##  2 MarkDown4         285924     68.4
##  3 MarkDown3         283671     67.8
##  4 MarkDown1         270755     64.7
##  5 MarkDown5         270057     64.6
##  6 Store                  0      0  
##  7 Date                   0      0  
##  8 Dept                   0      0  
##  9 Weekly_Sales           0      0  
## 10 IsHoliday              0      0  
## 11 Type                   0      0  
## 12 Size                   0      0  
## 13 Temperature            0      0  
## 14 Fuel_Price             0      0  
## 15 CPI                    0      0  
## 16 Unemployment           0      0  
## 17 next_Weekly_Sales      0      0

Видим, что пропущенные значения есть только в переменных Markdown1 - Markdown5. Их достаточно много, поэтому просто заменить средним нет смысла. Давайте выбросим их из нашего датасета.

walmart <- walmart %>%
  select(-starts_with('MarkDown'))

Теперь в нашем наборе данных нет пропущенных значений.

n_miss(walmart)

## [1] 0

Временной ряд продаж

Давайте немного изучим наши данные перед тем как строить модель. Построим временной ряд продаж.

walmart %>%
  group_by(Date) %>%
  summarise(Mean_Weekly_Sales = mean(Weekly_Sales)) %>%
  ggplot(aes(x = Date, y = Mean_Weekly_Sales)) +
    geom_line(color = 'darkblue') + 
  labs(x = 'Дата',
       y = 'Среднее количество проданных товаров')

Статистики

summary(walmart)

##      Store           Date                 Dept        Weekly_Sales   
##  Min.   : 1.0   Min.   :2010-02-05   Min.   : 1.00   Min.   : -4989  
##  1st Qu.:11.0   1st Qu.:2010-10-08   1st Qu.:18.00   1st Qu.:  2090  
##  Median :22.0   Median :2011-06-17   Median :37.00   Median :  7626  
##  Mean   :22.2   Mean   :2011-06-14   Mean   :44.25   Mean   : 16000  
##  3rd Qu.:33.0   3rd Qu.:2012-02-17   3rd Qu.:74.00   3rd Qu.: 20227  
##  Max.   :45.0   Max.   :2012-10-19   Max.   :99.00   Max.   :693099  
##  IsHoliday           Type                Size         Temperature    
##  Mode :logical   Length:418239      Min.   : 34875   Min.   : -2.06  
##  FALSE:388604    Class :character   1st Qu.: 93638   1st Qu.: 46.66  
##  TRUE :29635     Mode  :character   Median :140167   Median : 62.10  
##                                     Mean   :136750   Mean   : 60.08  
##                                     3rd Qu.:202505   3rd Qu.: 74.36  
##                                     Max.   :219622   Max.   :100.14  
##    Fuel_Price         CPI         Unemployment    next_Weekly_Sales
##  Min.   :2.472   Min.   :126.1   Min.   : 3.879   Min.   : -4989   
##  1st Qu.:2.925   1st Qu.:132.0   1st Qu.: 6.891   1st Qu.:  2090   
##  Median :3.445   Median :182.3   Median : 7.866   Median :  7624   
##  Mean   :3.358   Mean   :171.2   Mean   : 7.968   Mean   : 15990   
##  3rd Qu.:3.735   3rd Qu.:212.4   3rd Qu.: 8.572   3rd Qu.: 20222   
##  Max.   :4.468   Max.   :227.2   Max.   :14.313   Max.   :693099

Видим, что переменная Weekly_Sales и next_Weekly_Sales принимают отрицательные значения. Такого быть не должно, так как это выручка магазина. Давайте выкинем такие значения.

walmart <- walmart %>%
  filter(Weekly_Sales > 0) %>%
  filter(next_Weekly_Sales > 0)

Строим модели

Перед тем как строить модели давайте разделим наш датасет на train и test. Так как данные имеют временную структуру, не будем перемешивать наши данные. Посмотрим за сколько недель у нас есть данных.

walmart$Date %>% unique() %>% length()

## [1] 142

Видим, что у нас есть 142 недели в нашем распоряжении. Давайте возьмем первые 110 недель на обучение, а 32 оставшиеся на тестирование.

weeks <- walmart$Date %>% unique() %>% sort()
weeks[110]

## [1] "2012-03-09"

В трейн у нас попадут недели до 2012-03-09, а тест после этой недели, включая эту самую неделю.

train <- walmart %>%
  filter(Date < "2012-03-09")

test <- walmart %>%
  filter(Date >= "2012-03-09")

Давайте построим сначала некоторые простые модели. Их результат мы можем взять за базовый (baseline).

1. Среднее значение по всем отделам всех магазинов

Возьмем константный алгоритм, который будет всегда предсказывать значение равное среднему значению по всем отделам всех магазинов.

train_1 <- train
test_1 <- test

MEAN <- mean(train_1$Weekly_Sales)

train_1$predict1 <- MEAN
test_1$predict1 <- MEAN

Подсчитаем чему равно MAPE и MAE.

mae(train_1$next_Weekly_Sales, train_1$predict1)

## [1] 15228.08

mae(test_1$next_Weekly_Sales, test_1$predict1)

## [1] 15134.17

mape(train_1$next_Weekly_Sales, train_1$predict1)

## [1] 211.0388

mape(test_1$next_Weekly_Sales, test_1$predict1)

## [1] 848.3803

Видим большое значение MAPE. Скорее всего MAE тоже большое. Давайте посмотрим на среднее значение нашего таргета.

mean(train_1$next_Weekly_Sales)

## [1] 16118.9

Видим, что MAE и правда большое.

2. Среднее значение для каждого отдела в каждом магазине

Магазины и отделы в них бывают разные. Попробуем находить среднее в каждом отделе для каждого магазина.

MEAN_by_Store_Dept <- train %>%
  group_by(Store, Dept) %>%
  summarise(predict2 = mean(Weekly_Sales))

train_1 <- train_1 %>%
  left_join(MEAN_by_Store_Dept, by = c('Store', 'Dept'))

test_1 <- test_1 %>%
  left_join(MEAN_by_Store_Dept, by = c('Store', 'Dept'))

Появляются NA, потому что появились новые отделы в магазинах. Пусть если у нас нет прогноза, то мы возьмем прогноз из прошлой модели.

test_1 <- test_1 %>%
  mutate(predict2 = ifelse(is.na(predict2), predict1, predict2))

mae(train_1$next_Weekly_Sales, train_1$predict2)

## [1] 2781.925

mae(test_1$next_Weekly_Sales, test_1$predict2)

## [1] 2606.157

mape(train_1$next_Weekly_Sales, train_1$predict2)

## [1] 36.15821

mape(test_1$next_Weekly_Sales, test_1$predict2)

## [1] 107.8383

Уже получается лучше.

3. Значение предыдущей недели

Очень хорошим baseline’ом является следующая простая модель. Давайте привозить на следующую неделю столько товара, сколько мы продали на этой неделе.

train_1$predict3 <- train_1$Weekly_Sales
test_1$predict3 <- test_1$Weekly_Sales

mae(train_1$next_Weekly_Sales, train_1$predict3)

## [1] 2254.023

mae(test_1$next_Weekly_Sales, test_1$predict3)

## [1] 1648.434

mape(train_1$next_Weekly_Sales, train_1$predict3)

## [1] 1.176847

mape(test_1$next_Weekly_Sales, test_1$predict3)

## [1] 2.129825

Мы получили достаточно хорошие результаты. Давайте ориентироваться на них.

4. Обычная линейная регрессия

Давайте попробуем построить следующую очень простую регрессию:

\[next\_Weekly\_Sales = w_0 + w_1 Weekly\_Sales\].

Она похожа на то, что мы делали ранее.

train_2 <- train
test_2 <- test

model4 <- lm(next_Weekly_Sales ~ Weekly_Sales , data = train_2)

train_2$predict4 <- predict(model4, train_2)
test_2$predict4 <- predict(model4, test_2)

mae(train_2$next_Weekly_Sales, train_2$predict4)

## [1] 2523.658

mae(test_2$next_Weekly_Sales, test_2$predict4)

## [1] 1941.544

mape(train_2$next_Weekly_Sales, train_2$predict4)

## [1] 13.88451

mape(test_2$next_Weekly_Sales, test_2$predict4)

## [1] 53.67912

Результаты получились не такими хорошими как было раньше. Давайте попробуем нарисовать график рессеяния и понять почему так вышло.

train_2 %>%
  ggplot(aes(Weekly_Sales, next_Weekly_Sales)) + 
    geom_point() + 
    geom_smooth(method = 'lm')

Видно, что облако вытянуто неоднозначно. Есть большие магазины, а есть маленькие. В такой модели мы совсем не учитываем это и мешаем все в одну кучу.

5. Линейная регрессия с магазинами и отделами

Попробуем добавить в регрессию номер магазина и номер отдела. Оставлять переменные такими какие они есть не очень корректно, так как мы говорили, что Label Encoding работает не очень хорошо. Давайте закодируем наши магазины и отделы с помощью One Hot Encoding.

# install.packages("fastDummies")
library(fastDummies)

walmart2 <- walmart %>%
  select(Date, Store, Dept, Weekly_Sales, next_Weekly_Sales)

walmart2 <- dummy_cols(walmart2, select_columns = c('Store', 'Dept'),
                remove_first_dummy = T)

train_3 <- walmart2 %>%
  filter(Date < "2012-03-09")

test_3 <- walmart2 %>%
  filter(Date >= "2012-03-09")

Обучим модель.

model5 <- lm(next_Weekly_Sales ~ . - 1, data = train_3)

train_3$predict <- predict(model5, train_3)
test_3$predict <- predict(model5, test_3)

mae(train_3$next_Weekly_Sales, train_3$predict)

## [1] 2726.129

mae(test_3$next_Weekly_Sales, test_3$predict)

## [1] 2186.507

mape(train_3$next_Weekly_Sales, train_3$predict)

## [1] 17.34136

mape(test_3$next_Weekly_Sales, test_3$predict)

## [1] 55.55517

6. Логарифм таргета

Линейные модели любят, когда таргет да и вообще все переменнын имели нормальное распределение.

ggplot(train, aes(next_Weekly_Sales)) +
  geom_histogram()

Один из легких способов получить нормальное распределение это взять логарифм.

ggplot(train, aes(log(next_Weekly_Sales))) +
  geom_histogram()

Попробуем обучить модель с логарифмом таргета. При прогнозе не забываем возращать значения обратно с помощью экспоненты.

train_3 <- walmart2 %>%
  filter(Date < "2012-03-09") %>%
  mutate(next_Weekly_Sales = log(next_Weekly_Sales),
         Weekly_Sales = log(Weekly_Sales)) %>%
  select(-Date, Store, Dept)

test_3 <- walmart2 %>%
  filter(Date >= "2012-03-09") %>%
  mutate(next_Weekly_Sales = log(next_Weekly_Sales),
         Weekly_Sales = log(Weekly_Sales)) %>%
  select(-Date, Store, Dept)

model5 <- lm(next_Weekly_Sales ~ . - 1, data = train_3)

train_3$predict <- exp(predict(model5, train_3))
test_3$predict <- exp(predict(model5, test_3))

mae(exp(train_3$next_Weekly_Sales), train_3$predict)

## [1] 2304.386

mae(exp(test_3$next_Weekly_Sales), test_3$predict)

## [1] 1741.404

mape(exp(train_3$next_Weekly_Sales), train_3$predict)

## [1] 1.221557

mape(exp(test_3$next_Weekly_Sales), test_3$predict)

## [1] 2.329645

7. Добавляем остальные признаки

walmart3 <- walmart 

walmart3 <- dummy_cols(walmart3, select_columns = c('Store', 'Dept'),
                remove_first_dummy = T)

train_4 <- walmart3 %>%
  filter(Date < "2012-03-09") %>%
  mutate(next_Weekly_Sales = log(next_Weekly_Sales),
         Weekly_Sales = log(Weekly_Sales)) %>%
  select(-Date, Store, Dept)

test_4 <- walmart3 %>%
  filter(Date >= "2012-03-09") %>%
  mutate(next_Weekly_Sales = log(next_Weekly_Sales),
         Weekly_Sales = log(Weekly_Sales)) %>%
  select(-Date, Store, Dept)

model6 <- lm(next_Weekly_Sales ~ . - 1, data = train_4)

train_4$predict <- exp(predict(model6, train_4))
test_4$predict <- exp(predict(model6, test_4))

mae(exp(train_4$next_Weekly_Sales), train_4$predict)

## [1] 2298.427

mae(exp(test_4$next_Weekly_Sales), test_4$predict)

## [1] 1733.385

mape(exp(train_4$next_Weekly_Sales), train_4$predict)

## [1] 1.22473

mape(exp(test_4$next_Weekly_Sales), test_4$predict)

## [1] 2.342098