餐厅流量预测

介绍

这是一项针对Recruit Restaurant Visitor Forecasting的研究，借助 tidy R 和 ggplot2的强大功能。

该研究的目的是预测未来餐厅访客的数量。这使其成为时间序列预测问题。该数据是从许多日本餐馆收集的。正如我们将看到的，数据集很小并且易于访问，不需要太多内存或计算能力。

“Data以 8 个关系文件的形式出现，这些文件来自两个收集用户信息的独立日本网站：“Hot Pepper Gourmet (hpg)：类似于 Yelp”（search and reserve）和“AirREGI / Restaurant Board（air）：类似于 Square”（预订控制和收银机）。训练数据基于 2016 年 1 月至 2017 年 4 月的大部分时间范围，而测试集包括 4 月的最后一周加上 2017 年 5 月。测试数据“故意跨越日本称为“黄金周”的假期周。data description 进一步指出：“测试集中有几天餐厅关闭且没有访客。这些在计分时会被忽略。训练集忽略了餐馆关闭的日子。”

几个单独的数据集：

• air_visit_data.csv: air餐厅的历史访问数据。因为只有它包含实际顾客数量，所以本质上是主要的训练数据集。

• air_reserve.csv / hpg_reserve.csv: 通过air / hpg 系统进行预订的信息。

• air_store_info.csv / hpg_store_info.csv: 有关 air / hpg 餐厅的详细信息，包括风格和位置等内容。

• store_id_relation.csv: 连接 air 和 hpg 的ids。

• date_info.csv: 日本的节假日标记。

• sample_submission.csv: 作为一个测试集，但可见的只有air餐厅的id和到店时间，无法看到实际的顾客人数。

准备工作

加载library（库）

我们加载了一系列用于一般数据整理和一般可视化的库以及更专业的工具。

# general visualisation
library('ggplot2') # visualisation
library('scales') # visualisation
library('grid') # visualisation
library('gridExtra') # visualisation
library('RColorBrewer') # visualisation
library('corrplot') # visualisation

# general data manipulation
library('dplyr') # data manipulation
library('readr') # input/output
library('data.table') # data manipulation
library('tibble') # data wrangling
library('tidyr') # data wrangling
library('stringr') # string manipulation
library('forcats') # factor manipulation

# specific visualisation
library('ggrepel') # visualisation
library('ggridges') # visualisation
library('ggExtra') # visualisation
library('ggforce') # visualisation
library('viridis') # visualisation
library('recipes')
library('reticulate')

# specific data manipulation
library('lazyeval') # data wrangling
library('broom') # data wrangling
library('purrr') # string manipulation

# Date plus forecast
library('lubridate') # date and time
library('timeDate') # date and time
library('tseries') # time series analysis
library('forecast') # time series analysis
library('prophet') # time series analysis
library('timetk') # time series analysis

# Maps / geospatial
library('geosphere') # geospatial locations
library('leaflet') # maps
library('leaflet.extras') # maps
library('maps') # maps

可能有用的函数

我们使用R Cookbooks 提供的 multiplot 函数来创建多面板图。我们还利用一个简短的辅助函数来计算二项式置信区间

# Define multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

加载数据

air_visits <- read_csv('~/recruit-restaurant-visitor-forecasting/air_visit_data.csv/air_visit_data.csv', col_types = cols())

air_reserve <- read_csv('~/recruit-restaurant-visitor-forecasting/air_reserve.csv/air_reserve.csv', col_types = cols())

hpg_reserve <- read_csv('~/recruit-restaurant-visitor-forecasting/hpg_reserve.csv/hpg_reserve.csv', col_types = cols())

air_store <- read_csv('~/recruit-restaurant-visitor-forecasting/air_store_info.csv/air_store_info.csv', col_types = cols())

hpg_store <- read_csv('~/recruit-restaurant-visitor-forecasting/hpg_store_info.csv/hpg_store_info.csv', col_types = cols())

holidays <- read_csv('~/recruit-restaurant-visitor-forecasting/date_info.csv/date_info.csv', col_types = cols())

store_ids <- read_csv("~/recruit-restaurant-visitor-forecasting/store_id_relation.csv/store_id_relation.csv",col_types = cols())

test <- read_csv("~/recruit-restaurant-visitor-forecasting/sample_submission.csv/sample_submission.csv",col_types = cols())

概述：文件结构和内容

第一步，让我们使用 summary和glimpse工具对数据集进行概述

Air visits

summary(air_visits)

##  air_store_id         visit_date            visitors     
##  Length:252108      Min.   :2016-01-01   Min.   :  1.00  
##  Class :character   1st Qu.:2016-07-23   1st Qu.:  9.00  
##  Mode  :character   Median :2016-10-23   Median : 17.00  
##                     Mean   :2016-10-12   Mean   : 20.97  
##                     3rd Qu.:2017-01-24   3rd Qu.: 29.00  
##                     Max.   :2017-04-22   Max.   :877.00

glimpse(air_visits)

## Rows: 252,108
## Columns: 3
## $ air_store_id <chr> "air_ba937bf13d40fb24", "air_ba937bf13d40fb24", "air_ba93…
## $ visit_date   <date> 2016-01-13, 2016-01-14, 2016-01-15, 2016-01-16, 2016-01-…
## $ visitors     <dbl> 25, 32, 29, 22, 6, 9, 31, 21, 18, 26, 21, 11, 24, 21, 26,…

我们发现该文件包含每个visit_date和air_store_id的访客数量。日期特征应转换为时间序列格式。有 829 个不同的商店，这是一个很小的数据集：

air_visits %>% distinct(air_store_id) %>% nrow()

## [1] 829

Air Reserve

summary(air_reserve)

##  air_store_id       visit_datetime                  
##  Length:92378       Min.   :2016-01-01 19:00:00.00  
##  Class :character   1st Qu.:2016-11-15 19:00:00.00  
##  Mode  :character   Median :2017-01-05 18:00:00.00  
##                     Mean   :2016-12-05 08:18:58.75  
##                     3rd Qu.:2017-03-03 19:00:00.00  
##                     Max.   :2017-05-31 21:00:00.00  
##  reserve_datetime                 reserve_visitors 
##  Min.   :2016-01-01 01:00:00.00   Min.   :  1.000  
##  1st Qu.:2016-11-07 17:00:00.00   1st Qu.:  2.000  
##  Median :2016-12-27 22:00:00.00   Median :  3.000  
##  Mean   :2016-11-27 01:13:07.71   Mean   :  4.482  
##  3rd Qu.:2017-02-26 18:00:00.00   3rd Qu.:  5.000  
##  Max.   :2017-04-22 23:00:00.00   Max.   :100.000

glimpse(air_reserve)

## Rows: 92,378
## Columns: 4
## $ air_store_id     <chr> "air_877f79706adbfb06", "air_db4b38ebe7a7ceff", "air_…
## $ visit_datetime   <dttm> 2016-01-01 19:00:00, 2016-01-01 19:00:00, 2016-01-01…
## $ reserve_datetime <dttm> 2016-01-01 16:00:00, 2016-01-01 19:00:00, 2016-01-01…
## $ reserve_visitors <dbl> 1, 3, 6, 2, 5, 2, 4, 2, 2, 2, 3, 3, 2, 6, 7, 41, 13, …

我们发现air预订中包括预订的date和time，以及访问的日期和时间。我们有 314 家air商店的预订号码：

air_reserve %>% distinct(air_store_id) %>% nrow()

## [1] 314

HPG Reserve

summary(hpg_reserve)

##  hpg_store_id       visit_datetime                  
##  Length:2000320     Min.   :2016-01-01 11:00:00.00  
##  Class :character   1st Qu.:2016-06-26 19:00:00.00  
##  Mode  :character   Median :2016-11-19 20:00:00.00  
##                     Mean   :2016-10-15 06:55:20.07  
##                     3rd Qu.:2017-02-03 19:00:00.00  
##                     Max.   :2017-05-31 23:00:00.00  
##  reserve_datetime                 reserve_visitors 
##  Min.   :2016-01-01 00:00:00.00   Min.   :  1.000  
##  1st Qu.:2016-06-21 12:00:00.00   1st Qu.:  2.000  
##  Median :2016-11-10 20:00:00.00   Median :  3.000  
##  Mean   :2016-10-07 19:57:59.30   Mean   :  5.074  
##  3rd Qu.:2017-01-27 13:00:00.00   3rd Qu.:  6.000  
##  Max.   :2017-04-22 23:00:00.00   Max.   :100.000

glimpse(hpg_reserve)

## Rows: 2,000,320
## Columns: 4
## $ hpg_store_id     <chr> "hpg_c63f6f42e088e50f", "hpg_dac72789163a3f47", "hpg_…
## $ visit_datetime   <dttm> 2016-01-01 11:00:00, 2016-01-01 13:00:00, 2016-01-01…
## $ reserve_datetime <dttm> 2016-01-01 09:00:00, 2016-01-01 06:00:00, 2016-01-01…
## $ reserve_visitors <dbl> 1, 3, 2, 5, 13, 2, 2, 2, 2, 6, 2, 2, 2, 2, 5, 4, 2, 4…

Hpg保留文件遵循与相应的air文件相同的结构。我们有 13325 hpg 商店的预订号码：

hpg_reserve %>% distinct(hpg_store_id) %>% nrow()

## [1] 13325

Air Store

summary(air_store)

##  air_store_id       air_genre_name     air_area_name         latitude    
##  Length:829         Length:829         Length:829         Min.   :33.21  
##  Class :character   Class :character   Class :character   1st Qu.:34.70  
##  Mode  :character   Mode  :character   Mode  :character   Median :35.66  
##                                                           Mean   :35.65  
##                                                           3rd Qu.:35.69  
##                                                           Max.   :44.02  
##    longitude    
##  Min.   :130.2  
##  1st Qu.:135.3  
##  Median :139.7  
##  Mean   :137.4  
##  3rd Qu.:139.8  
##  Max.   :144.3

glimpse(air_store)

## Rows: 829
## Columns: 5
## $ air_store_id   <chr> "air_0f0cdeee6c9bf3d7", "air_7cc17a324ae5c7dc", "air_fe…
## $ air_genre_name <chr> "Italian/French", "Italian/French", "Italian/French", "…
## $ air_area_name  <chr> "Hyōgo-ken Kōbe-shi Kumoidōri", "Hyōgo-ken Kōbe-shi Kum…
## $ latitude       <dbl> 34.69512, 34.69512, 34.69512, 34.69512, 35.65807, 35.65…
## $ longitude      <dbl> 135.1979, 135.1979, 135.1979, 135.1979, 139.7516, 139.7…

我们发现 air_store信息包含餐厅类型的名称以及区域名称。

HPG Store

summary(hpg_store)

##  hpg_store_id       hpg_genre_name     hpg_area_name         latitude    
##  Length:4690        Length:4690        Length:4690        Min.   :33.31  
##  Class :character   Class :character   Class :character   1st Qu.:34.69  
##  Mode  :character   Mode  :character   Mode  :character   Median :35.66  
##                                                           Mean   :35.81  
##                                                           3rd Qu.:35.70  
##                                                           Max.   :43.77  
##    longitude    
##  Min.   :130.3  
##  1st Qu.:135.5  
##  Median :139.5  
##  Mean   :137.7  
##  3rd Qu.:139.7  
##  Max.   :143.7

glimpse(hpg_store)

## Rows: 4,690
## Columns: 5
## $ hpg_store_id   <chr> "hpg_6622b62385aec8bf", "hpg_e9e068dd49c5fa00", "hpg_29…
## $ hpg_genre_name <chr> "Japanese style", "Japanese style", "Japanese style", "…
## $ hpg_area_name  <chr> "Tōkyō-to Setagaya-ku Taishidō", "Tōkyō-to Setagaya-ku …
## $ latitude       <dbl> 35.64367, 35.64367, 35.64367, 35.64367, 35.64367, 35.64…
## $ longitude      <dbl> 139.6682, 139.6682, 139.6682, 139.6682, 139.6682, 139.6…

同样，hpg_store 中的内容遵循与air 相同的结构。这里，genre_name包含单词 style。值得注意的是，air数据是否也是如此，或者它是否只是指特定的“Japanese style”。有 4690 个不同的 hpg_store_ids，这明显少于我们的预订数据。

Holidays

summary(holidays)

##  calendar_date        day_of_week         holiday_flg    
##  Min.   :2016-01-01   Length:517         Min.   :0.0000  
##  1st Qu.:2016-05-09   Class :character   1st Qu.:0.0000  
##  Median :2016-09-15   Mode  :character   Median :0.0000  
##  Mean   :2016-09-15                      Mean   :0.0677  
##  3rd Qu.:2017-01-22                      3rd Qu.:0.0000  
##  Max.   :2017-05-31                      Max.   :1.0000

glimpse(holidays)

## Rows: 517
## Columns: 3
## $ calendar_date <date> 2016-01-01, 2016-01-02, 2016-01-03, 2016-01-04, 2016-01…
## $ day_of_week   <chr> "Friday", "Saturday", "Sunday", "Monday", "Tuesday", "We…
## $ holiday_flg   <dbl> 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,…

我们将 date_info文件称为“假期”，它存在的目的便是与后续的标记进行对照链接。假期被编码为整数格式的二进制标志。这应该可以成为用于预测餐厅顾客流量的逻辑二进制特征。

Store IDs

summary(store_ids)

##  air_store_id       hpg_store_id      
##  Length:150         Length:150        
##  Class :character   Class :character  
##  Mode  :character   Mode  :character

glimpse(store_ids)

## Rows: 150
## Columns: 2
## $ air_store_id <chr> "air_63b13c56b7201bd9", "air_a24bf50c3e90d583", "air_c7f7…
## $ hpg_store_id <chr> "hpg_4bc649e72e2a239a", "hpg_c34b496d0305a809", "hpg_cd8a…

这是一个连接air 和 hpg ids的关系文件。仅有 150 双，不到所有air 店的 20%。

Test data

summary(test)

##       id               visitors
##  Length:32019       Min.   :0  
##  Class :character   1st Qu.:0  
##  Mode  :character   Median :0  
##                     Mean   :0  
##                     3rd Qu.:0  
##                     Max.   :0

glimpse(test)

## Rows: 32,019
## Columns: 2
## $ id       <chr> "air_00a91d42b08b08d9_2017-04-23", "air_00a91d42b08b08d9_2017…
## $ visitors <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…

如数据描述中所述，最终提交文件的 id 是 air_id 和date的串联。

检查缺失值

sum(is.na(air_visits))

## [1] 0

sum(is.na(air_reserve))

## [1] 0

sum(is.na(hpg_reserve))

## [1] 0

sum(is.na(air_store))

## [1] 0

sum(is.na(hpg_store))

## [1] 0

sum(is.na(holidays))

## [1] 0

sum(is.na(store_ids))

## [1] 0

sum(is.na(test))

## [1] 0

我们的数据中没有缺失值，这便于后续的操作。

重组特征

我们更改了date/time 特征的格式，并将一些特征重新格式化为逻辑变量和因子变量以用于探索后续的预测。

air_visits <- air_visits %>%
  mutate(visit_date = ymd(visit_date))


air_reserve <- air_reserve %>%
  mutate(visit_datetime_converted = ymd_hms(visit_datetime, quiet = TRUE),
         reserve_datetime_converted = ymd_hms(reserve_datetime, quiet = TRUE),
         visit_convert_error = is.na(ymd_hms(visit_datetime, quiet = TRUE)),
         reserve_convert_error = is.na(ymd_hms(reserve_datetime, quiet = TRUE)))

conversion_errors <- air_reserve %>%
  filter(visit_convert_error | reserve_convert_error)
print(conversion_errors)

## # A tibble: 2,257 × 8
##    air_store_id         visit_datetime      reserve_datetime    reserve_visitors
##    <chr>                <dttm>              <dttm>                         <dbl>
##  1 air_326ca454ef3558bc 2016-01-03 00:00:00 2016-01-02 13:00:00                4
##  2 air_a083834e7ffe187e 2016-01-05 19:00:00 2016-01-02 00:00:00                8
##  3 air_de692863bb2dd758 2016-01-08 19:00:00 2016-01-08 00:00:00                2
##  4 air_326ca454ef3558bc 2016-01-09 00:00:00 2016-01-07 19:00:00                6
##  5 air_b9e27558fb8bd5c4 2016-01-09 19:00:00 2016-01-08 00:00:00                4
##  6 air_db4b38ebe7a7ceff 2016-01-10 00:00:00 2016-01-09 21:00:00                1
##  7 air_622375b4815cf5cb 2016-01-10 20:00:00 2016-01-10 00:00:00               53
##  8 air_3cad29d1a23209d2 2016-01-11 20:00:00 2016-01-11 00:00:00                2
##  9 air_346ade7d29230634 2016-01-12 20:00:00 2016-01-10 00:00:00                3
## 10 air_8093d0b565e9dbdf 2016-01-15 18:00:00 2016-01-15 00:00:00                2
## # ℹ 2,247 more rows
## # ℹ 4 more variables: visit_datetime_converted <dttm>,
## #   reserve_datetime_converted <dttm>, visit_convert_error <lgl>,
## #   reserve_convert_error <lgl>

air_reserve <- air_reserve %>%
  filter(!visit_convert_error & !reserve_convert_error)
air_reserve <- air_reserve %>%
  select(-visit_convert_error, -reserve_convert_error,-visit_datetime_converted,-reserve_datetime_converted)


hpg_reserve <- hpg_reserve %>%
  mutate(visit_datetime_converted = ymd_hms(visit_datetime, quiet = TRUE),
         reserve_datetime_converted = ymd_hms(reserve_datetime, quiet = TRUE),
         visit_convert_error = is.na(ymd_hms(visit_datetime, quiet = TRUE)),
         reserve_convert_error = is.na(ymd_hms(reserve_datetime, quiet = TRUE)))

conversion_errors <- hpg_reserve %>%
  filter(visit_convert_error | reserve_convert_error)
print(conversion_errors)

## # A tibble: 69,144 × 8
##    hpg_store_id         visit_datetime      reserve_datetime    reserve_visitors
##    <chr>                <dttm>              <dttm>                         <dbl>
##  1 hpg_3e755d9825881876 2016-01-01 18:00:00 2016-01-01 00:00:00                2
##  2 hpg_614fa1183c1db71e 2016-01-01 18:00:00 2016-01-01 00:00:00                3
##  3 hpg_8623b83c6b3cdcb0 2016-01-01 18:00:00 2016-01-01 00:00:00                2
##  4 hpg_8623b83c6b3cdcb0 2016-01-01 18:00:00 2016-01-01 00:00:00                2
##  5 hpg_f23b00c038f70720 2016-01-01 19:00:00 2016-01-01 00:00:00                4
##  6 hpg_dd69814852f189eb 2016-01-01 20:00:00 2016-01-01 00:00:00                6
##  7 hpg_a00b81d9c3a2d140 2016-01-02 12:00:00 2016-01-01 00:00:00                6
##  8 hpg_11453449d243bb20 2016-01-02 13:00:00 2016-01-02 00:00:00                3
##  9 hpg_4c0f91e3fb8d1ba7 2016-01-02 13:00:00 2016-01-01 00:00:00               12
## 10 hpg_5c83ed5c4d14a95c 2016-01-02 13:00:00 2016-01-01 00:00:00                4
## # ℹ 69,134 more rows
## # ℹ 4 more variables: visit_datetime_converted <dttm>,
## #   reserve_datetime_converted <dttm>, visit_convert_error <lgl>,
## #   reserve_convert_error <lgl>

hpg_reserve <- hpg_reserve %>%
  filter(!visit_convert_error & !reserve_convert_error)
hpg_reserve <- hpg_reserve %>%
  select(-visit_convert_error, -reserve_convert_error,-visit_datetime_converted,-reserve_datetime_converted)

air_store <- air_store %>%
  mutate(air_genre_name = as.factor(air_genre_name),
         air_area_name = as.factor(air_area_name))

hpg_store <- hpg_store %>%
  mutate(hpg_genre_name = as.factor(hpg_genre_name),
         hpg_area_name = as.factor(hpg_area_name))

holidays <- holidays %>%
  mutate(holiday_flg = as.logical(holiday_flg),
         date = ymd(calendar_date),
         calendar_date = as.character(calendar_date))

hpg_reserve <- hpg_reserve %>%
  mutate(visit_datetime = ymd_hms(visit_datetime),
         reserve_datetime = ymd_hms(reserve_datetime))

air_store <- air_store %>%
  mutate(air_genre_name = as.factor(air_genre_name),
         air_area_name = as.factor(air_area_name))

hpg_store <- hpg_store %>%
  mutate(hpg_genre_name = as.factor(hpg_genre_name),
         hpg_area_name = as.factor(hpg_area_name))

holidays <- holidays %>%
  mutate(holiday_flg = as.logical(holiday_flg),
         date = ymd(calendar_date),
         calendar_date = as.character(calendar_date))

单独特征可视化

在这里，我们首先查看各个数据文件中特征的分布，然后将它们组合起来进行更详细的分析。这个单独初始可视化将成为我们进行分析的基础。

Air Visits

我们从air餐厅的顾客人数开始。在这里，我们绘制了整个数据集时间范围内每天的顾客总数以及一年中一周和一个月中每天的访客中位数：

p1 <- air_visits %>%
  group_by(visit_date) %>%
  summarise(all_visitors = sum(visitors)) %>%
  ggplot(aes(visit_date,all_visitors)) +
  geom_line(col = "blue") +
  labs(x = "Date", y = "All visitors")

p2 <- air_visits %>%
  ggplot(aes(visitors)) +
  geom_vline(xintercept = 20, color = "orange") +
  geom_histogram(fill = "blue", bins = 30) +
  scale_x_log10()

p3 <- air_visits %>%
  mutate(wday = wday(visit_date, label = TRUE, week_start = 1)) %>%
  group_by(wday) %>%
  summarise(visits = median(visitors)) %>%
  ggplot(aes(wday, visits, fill = wday)) +
  geom_col() +
  theme(legend.position = "none", axis.text.x  = element_text(angle=45, hjust=1, vjust=0.9)) +
  labs(x = "Day of the week", y = "Median visitors") +
  scale_fill_hue()
  
p4 <- air_visits %>%
  mutate(month = month(visit_date, label = TRUE)) %>%
  group_by(month) %>%
  summarise(visits = median(visitors)) %>%
  ggplot(aes(month, visits, fill = month)) +
  geom_col() +
  theme(legend.position = "none") +
  labs(x = "Month", y = "Median visitors")

layout <- matrix(c(1,1,1,1,2,3,4,4),2,4,byrow=TRUE)
multiplot(p1, p2, p3, p4, layout=layout)

Fig. 1

我们发现：

• 整个时间序列中有一个有趣的长期阶梯结构（这可能与数据库中添加新餐馆有关）。此外，我们已经看到了很可能对应于每周周期的周期性模式。

• 每家餐厅每天每次光顾的客人人数高峰在 20 人左右（橙色线），而最多的情况下可达 100 个，在极少数情况下甚至会超过 100 个。

• 周五和周末似乎是最受欢迎的日子，而周一和周二的平均游客数量最少。

• 一年中也存在一定程度的变化。十二月似乎是餐厅最受欢迎的月份，而三月至五月期间餐厅的人数一直很保持在一个较高的点。

我们将预测 2017 年 4 月最后一周和 5 月，所以让我们看看 2016 年训练数据中的这个时间范围：

air_visits %>%
  filter(visit_date > ymd("2016-04-15") & visit_date < ymd("2016-06-15")) %>%
  group_by(visit_date) %>%
  summarise(all_visitors = sum(visitors)) %>%
  ggplot(aes(visit_date,all_visitors)) +
  geom_line() +
  geom_smooth(method = "loess", color = "blue", span = 1/7) +
  labs(y = "All visitors", x = "Date")

## `geom_smooth()` using formula = 'y ~ x'

Fig. 2

这里，黑线是日期，蓝线对应于与相应灰色置信区域的平滑拟合。我们再次看到每周周期以及上述黄金周（Golden Week）的影响，2016 年黄金周发生在 4 月 29 日至 5 月 5 日之间。

Air Reservations

让我们看看我们的预订数据与实际顾客人数的比较。我们从 air 餐厅开始，通过每天的预订、到达时间以及预订和到达餐厅之间的时间差来可视化其访客量：

foo <- air_reserve %>%
  mutate(reserve_date = date(reserve_datetime),
         reserve_hour = hour(reserve_datetime),
         reserve_wday = wday(reserve_datetime, label = TRUE, week_start = 1),
         visit_date = date(visit_datetime),
         visit_hour = hour(visit_datetime),
         visit_wday = wday(visit_datetime, label = TRUE, week_start = 1),
         diff_hour = time_length(visit_datetime - reserve_datetime, unit = "hour"),
         diff_day = time_length(visit_datetime - reserve_datetime, unit = "day")
         )

p1 <- foo %>%
  group_by(visit_date) %>%
  summarise(all_visitors = sum(reserve_visitors)) %>%
  ggplot(aes(visit_date, all_visitors)) +
  geom_line() +
  labs(x = "'air' visit date")

p2 <- foo %>%
  group_by(visit_hour) %>%
  summarise(all_visitors = sum(reserve_visitors)) %>%
  ggplot(aes(visit_hour, all_visitors)) +
  geom_col(fill = "blue")

p3 <- foo %>%
  filter(diff_hour < 24*5) %>%
  group_by(diff_hour) %>%
  summarise(all_visitors = sum(reserve_visitors)) %>%
  ggplot(aes(diff_hour, all_visitors)) +
  geom_col(fill = "blue") +
  labs(x = "Time from reservation to visit [hours]")

layout <- matrix(c(1,1,2,3),2,2,byrow=TRUE)
multiplot(p1, p2, p3, layout=layout)

Fig. 3

我们发现：

• 2016年通过air系统进行的预订要少得，甚至很长一段时间根本没有。该数量仅在当年年底有所增加。 2017 年，顾客数量保持在一个较高的水平。我们在第一季度之后看到的人为下降很可能与数据集时间范围结束时的这些预订有关，这意味着长期预订不会成为该数据集的一部分。

• 通常是预定在晚上的晚餐时间。

• 预订和到达餐厅之间的时间（此处以小时为单位）遵循 24 小时模式。最常见的预定方式是在参观前几个小时预订，但如果提前更早的预订，那么在晚上预订第二天晚上就餐也是很常见的。该图被截断以显示这种模式，该模式持续到更长的时间跨度。预订和到达之间的时间间隔很长的情况并不罕见。这些是air数据的最极端值，最多提前一年多：

foo %>%
  arrange(desc(diff_day)) %>%
  select(reserve_datetime, visit_datetime, diff_day, air_store_id) %>%
  head(5)

## # A tibble: 5 × 4
##   reserve_datetime    visit_datetime      diff_day air_store_id        
##   <dttm>              <dttm>                 <dbl> <chr>               
## 1 2016-01-11 17:00:00 2017-02-07 20:00:00     393. air_e7fbee4e3cfe65c5
## 2 2016-01-12 20:00:00 2017-02-07 20:00:00     392  air_e7fbee4e3cfe65c5
## 3 2016-01-18 18:00:00 2017-02-07 20:00:00     386. air_e7fbee4e3cfe65c5
## 4 2016-01-02 20:00:00 2017-01-18 20:00:00     382  air_2a485b92210c98b5
## 5 2016-01-02 21:00:00 2017-01-18 20:00:00     382. air_2a485b92210c98b5

请注意，这 5 家餐厅仅包含 2 家不同的餐厅。这些要么是非常奇怪的地方，要么这些数据可能是数据输入错误和年份混淆的结果。

HPG Reservations

按照上面相同的模式， hpg 预定如下:

foo <- hpg_reserve %>%
  mutate(reserve_date = date(reserve_datetime),
         reserve_hour = hour(reserve_datetime),
         visit_date = date(visit_datetime),
         visit_hour = hour(visit_datetime),
         diff_hour = time_length(visit_datetime - reserve_datetime, unit = "hour"),
         diff_day = time_length(visit_datetime - reserve_datetime, unit = "day")
         )

p1 <- foo %>%
  group_by(visit_date) %>%
  summarise(all_visitors = sum(reserve_visitors)) %>%
  ggplot(aes(visit_date, all_visitors)) +
  geom_line() +
  labs(x = "'hpg' visit date")

p2 <- foo %>%
  group_by(visit_hour) %>%
  summarise(all_visitors = sum(reserve_visitors)) %>%
  ggplot(aes(visit_hour, all_visitors)) +
  geom_col(fill = "red")

p3 <- foo %>%
  filter(diff_hour < 24*5) %>%
  group_by(diff_hour) %>%
  summarise(all_visitors = sum(reserve_visitors)) %>%
  ggplot(aes(diff_hour, all_visitors)) +
  geom_col(fill = "red") +
  labs(x = "Time from reservation to visit [hours]")

layout <- matrix(c(1,1,2,3),2,2,byrow=TRUE)
multiplot(p1, p2, p3, layout=layout)

Fig. 4

我们发现：

• 这里，预订后的顾客数量遵循更有序的模式，在 2016 年 12 月出现明显峰值。与上面的air数据一样，我们还看到，随着接近时间范围的结束，预订顾客数量有所下降。

• 同样，大多数预订都是晚餐。我们还看到了另一种不错的 24 小时预订模式。值得注意的是，到达餐厅前最后几个小时的顾客数量并不比到达餐厅前 24 或 48 小时更多。这与 air 数据形成鲜明对比。

Air Store

在可视化时间方面之后，现在让我们看看空间信息。注意，根据data description，“经纬度是店铺所属区域的经纬度”。这是为了阻止我们识别特定的餐馆（如果可以识别到具体的餐厅位置，便可以做单独的预测和分析）。

这是所有 air餐厅的完全互动且可缩放的地图。单击集群将它们分解为更小的集群，最终分解为单独的餐厅，并按其类型进行标记。该地图很好地直观地显示了许多餐厅共享共同坐标的事实，因为这些坐标指的是餐厅的区域。单击单个标记即可查看其air_store_id（该地图由leaflet package提供，其中包括各种用于交互式地图的酷工具）。

leaflet(air_store) %>%
  addTiles() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addMarkers(~longitude, ~latitude,
             popup = ~air_store_id, label = ~air_genre_name,
             clusterOptions = markerClusterOptions())

Fig. 5

接下来，我们将不同餐厅类型（或air_genre_names）的数量与air餐厅最多的区域一起绘制：

p1 <- air_store %>%
  group_by(air_genre_name) %>%
  count() %>%
  ggplot(aes(reorder(air_genre_name, n, FUN = min), n, fill = air_genre_name)) +
  geom_col() +
  coord_flip() +
  theme(legend.position = "none") +
  labs(x = "Type of cuisine (air_genre_name)", y = "Number of air restaurants")
  
p2 <- air_store %>%
  group_by(air_area_name) %>%
  count() %>%
  ungroup() %>%
  top_n(15,n) %>%
  ggplot(aes(reorder(air_area_name, n, FUN = min) ,n, fill = air_area_name)) +
  geom_col() +
  theme(legend.position = "none") +
  coord_flip() +
  labs(x = "Top 15 areas (air_area_name)", y = "Number of air restaurants")

layout <- matrix(c(1,2),2,1,byrow=TRUE)
multiplot(p1, p2, layout=layout)

Fig. 6

我们发现：

• 我们的数据中有很多Izakaya美食酒吧，其次是 Cafe’s。我们的 air数据集中没有很多 Karaoke 场所，也只有少数将自己描述为一般“International”或“Asian”（“creative cuisine”的具体内容并不知道，可以后续分析下）。

• Fukuoka每个地区的air餐厅数量最多，其次是Tokyo地区的部分。

HPG Store

与上面的air餐厅相同，我们还为不同的hpg餐厅创建了一个交互式地图：

leaflet(hpg_store) %>%
  addTiles() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addMarkers(~longitude, ~latitude,
             popup = ~hpg_store_id, label = ~hpg_genre_name,
              clusterOptions = markerClusterOptions())

以下是 hpg 餐厅的风格和区域细分：

p1 <- hpg_store %>%
  group_by(hpg_genre_name) %>%
  count() %>%
  ggplot(aes(reorder(hpg_genre_name, n, FUN = min), n, fill = hpg_genre_name)) +
  geom_col() +
  coord_flip() +
  theme(legend.position = "none") +
  labs(x = "Type of cuisine (hpg_genre_name)", y = "Number of hpg restaurants")
  
p2 <- hpg_store %>%
  mutate(area = str_sub(hpg_area_name, 1, 20)) %>%
  group_by(area) %>%
  count() %>%
  ungroup() %>%
  top_n(15,n) %>%
  ggplot(aes(reorder(area, n, FUN = min) ,n, fill = area)) +
  geom_col() +
  theme(legend.position = "none") +
  coord_flip() +
  labs(x = "Top 15 areas (hpg_area_name)", y = "Number of hpg restaurants")

layout <- matrix(c(1,2),1,2,byrow=TRUE)
multiplot(p1, p2, layout=layout)

Fig. 8

这里我们将 hpg_area_name截断为 20 个字符，以使绘图更具可读性。

我们发现:

• HPG 描述包含比air数据更多的类型。在这里，“Japanese style”似乎包含更多在air数据中更具体分类的地方。这同样适用于“International cuisine”。

• 在前 15 名地区中，我们再次发现 Tokyo和 Osaka占据显着位置。

Holidays

让我们快速浏览一下假期。我们将绘制总共有多少个，以及它们在 2017 年预测时间范围和 2016 年相应时间内的分布情况：

foo <- holidays %>%
  mutate(wday = wday(date, week_start = 1))

p1 <- foo %>%
  ggplot(aes(holiday_flg, fill = holiday_flg)) +
  geom_bar() +
  theme(legend.position = "none")
  
p2 <- foo %>%
  filter(date > ymd("2016-04-15") & date < ymd("2016-06-01")) %>%
  ggplot(aes(date, holiday_flg, color = holiday_flg)) +
  geom_point(size = 2) +
  theme(legend.position = "none") +
  labs(x = "2016 date")

p3 <- foo %>%
  filter(date > ymd("2017-04-15") & date < ymd("2017-06-01")) %>%
  ggplot(aes(date, holiday_flg, color = holiday_flg)) +
  geom_point(size = 2) +
  theme(legend.position = "none") +
  labs(x = "2017 date")

layout <- matrix(c(1,1,2,3),2,2,byrow=FALSE)
multiplot(p1, p2, p3, layout=layout)

Fig. 9

我们发现：

• 2016 年 4 月底/5 月的假期与 2017 年相同.

• 在这段时间内的假期大概占到总时间的7%:

holidays %>% summarise(frac = mean(holiday_flg))

## # A tibble: 1 × 1
##     frac
##    <dbl>
## 1 0.0677

测试数据集

下图直观地显示了train和test数据集的时间范围：

foo <- air_visits %>%
  rename(date = visit_date) %>%
  distinct(date) %>%
  mutate(dset = "train")

bar <- test %>%
  separate(id, c("foo", "bar", "date"), sep = "_") %>%
  mutate(date = ymd(date)) %>%
  distinct(date) %>%
  mutate(dset = "test")

foo <- foo %>%
  bind_rows(bar) %>%
  mutate(year = year(date))
year(foo$date) <- 2017

foo %>%
  filter(!is.na(date)) %>%
  mutate(year = fct_relevel(as.factor(year), c("2017","2016"))) %>%
  ggplot(aes(date, year, color = dset)) +
  geom_point(shape = "|", size = 10) +
  scale_x_date(date_labels = "%B", date_breaks = "1 month") +
  #scale_y_reverse() +
  theme(legend.position = "bottom", axis.text.x  = element_text(angle=45, hjust=1, vjust=0.9)) +
  labs(color = "Data set") +
  guides(color = guide_legend(override.aes = list(size = 4, pch = 15)))

Fig. 10

特征关系

单独查看每个数据集后，可以尝试开始将它们组合起来。这将告诉我们有关各种特征之间的关系以及这些关系如何影响顾客数量的内在联系。我们发现的任何的内容都需要在各个特征分布的背景下进行解释，这就是为什么这是我们研究这些的第一步。

餐厅类型与顾客数量

我们的多特征空间的第一个图涉及按美食类型的air餐厅顾客的平均数量（即air_genre_name）。我们使用分面图来区分 14 个类别的时间序列(注意对数 y 轴)：

foo <- air_visits %>%
  left_join(air_store, by = "air_store_id")

foo %>%
  group_by(visit_date, air_genre_name) %>%
  summarise(mean_visitors = mean(visitors)) %>%
  ungroup() %>%
  ggplot(aes(visit_date, mean_visitors, color = air_genre_name)) +
  geom_line() +
  labs(y = "Average number of visitors to 'air' restaurants", x = "Date") +
  theme(legend.position = "none") +
  scale_y_log10() +
  facet_wrap(~ air_genre_name)

Fig. 11

我们发现：

• 每种类型每天的平均值范围为 10 到 100 名访客。在每个类型中，长期趋势看起来相当稳定。自2016年底以来，“Creative Cuisine”和“Okonomiyaki”等食品的受欢迎程度呈上升趋势，而“亚洲”食品的受欢迎程度则不断下降。

• 可以理解的是，像“Karaoke”或“Asian”这样的低计数时间系列（见图 6）比顾客数量较多的类型更加嘈杂（变化趋势不好确定）。尽管如此，“亚洲”餐厅似乎很受欢迎（或许是因为它们很稀有？）。

在所有曲线中，我们都看到了熟悉的每周变化，因此让我们更详细地了解每周每天和类型的平均访客数量。我们还通过 ggridges 包添加了一系列山脊线图。山脊线图可以快速比较半重叠（密度）曲线。在这里，我们显示了每种类型每天的访客分布情况。通过 ggplot 功能，这两个图之间的 y 轴标签指的是它们：

p1 <- foo %>%
  mutate(wday = wday(visit_date, label = TRUE, week_start = 1)) %>%
  group_by(wday, air_genre_name) %>%
  summarise(mean_visitors = mean(visitors)) %>%
  ggplot(aes(air_genre_name, mean_visitors, color = wday)) +
  geom_point(size = 4) +
  theme(legend.position = "left", axis.text.y = element_blank(),
        plot.title = element_text(size = 14)) +
  coord_flip() +
  labs(x = "") +
  scale_x_discrete(position = "top") +
  ggtitle("air_genre_name") +
  scale_color_hue()

p2 <- foo %>%
  ggplot(aes(visitors, air_genre_name, fill = air_genre_name)) +
  geom_density_ridges(bandwidth = 0.1) +
  scale_x_log10() +
  theme(legend.position = "none") +
  labs(y = "") +
  scale_fill_cyclical(values = c("blue", "red"))

layout <- matrix(c(1,1,2,2,2),1,5,byrow=TRUE)
multiplot(p1, p2, layout=layout)

Fig. 12

p1 <- 1; p2 <- 1; p3 <- 1; p4 <- 1; p5 <- 1

这里每种颜色对应一周中的一天。红色代表周末，冷色代表周中。星期一是深橙色。

我们发现

• 周末和工作日最大的区别在于“Karaoke”酒吧，它们统治着周末。在“International”美食中也可以看到类似的趋势，尽管差距要小得多。.

• 没有任何一种类型能够真正违背繁忙周末的趋势。最小的变化是在通用的“Other”类别中，即“Korean”食品，以及“South Korea”美食，这是周五最繁忙的唯一类型。一般的“Bars/Cocktail”总体上尤其不受欢迎。

• 密度曲线证实了我们从工作日分布中得到的印象：“Asian”餐厅每次约会的游客很少少于 10 人，而“Karaoke”场所由于周末的强烈影响而显示出非常广泛的分布(注意对数 x 轴)。

节假日与顾客数量

我们将通过比较带节假日和不带节假日标志的天数的统计数据来研究节假日对顾客数量的影响：

foo <- air_visits %>%
  mutate(calendar_date = as.character(visit_date)) %>%
  left_join(holidays, by = "calendar_date")

p1 <- foo %>%
  ggplot(aes(holiday_flg, visitors, color = holiday_flg)) +
  geom_boxplot() +
  scale_y_log10() +
  theme(legend.position = "none")

p2 <- foo %>%
  mutate(wday = wday(date, label = TRUE, week_start = 1)) %>%
  group_by(wday, holiday_flg) %>%
  summarise(mean_visitors = mean(visitors)) %>%
  ggplot(aes(wday, mean_visitors, color = holiday_flg)) +
  geom_point(size = 4) +
  theme(legend.position = "none") +
  labs(y = "Average number of visitors")

layout <- matrix(c(1,2),1,2,byrow=TRUE)
multiplot(p1, p2, layout=layout)

Fig. 13

我们发现：

• 总体而言，假期对平均顾客数量没有任何影响（左图）（更多信息隐藏在细节中，或许可以再探索一下）。

• 虽然周末假期对顾客数量影响不大，甚至周末假期的顾客人数还略有减少，但工作日的影响更为明显，尤其是周一和周二（右图）。

每个区域的餐厅数量与顾客数量

我们的下一个探索源于一个简单的想法：如果类型“Bars”很受欢迎并且我们拥有该地区唯一的“Bars”，那么我们可以期待很多顾客。如果我们街上还有其他 12 家“Bars”，那么，即便尽我们所能，其中一些顾客依然会冒险进入其他场所。经济学家告诉我们，我们最终可以预期供需趋于平衡。即使对于像我们的数据集这样的短期数据集，以及相对局部的区域，研究餐馆集群可能仍然有价值。因此，我们来研究一下每个地区某种类型餐厅的数量及其对顾客数量的影响。

我们首先概述了air和 hpg 商店两个数据集的每个区域中某些类型的频率。这很可能跟前一章的一部分重合，这里需要重复下。以下计数图显示哪些类型存在于哪些区域（名称被截断）。点的大小与事例的数量成正比：

air_store %>%
  mutate(area = str_sub(air_area_name, 1, 12)) %>%
  ggplot(aes(area, air_genre_name)) +
  geom_count(colour = "blue") +
  theme(legend.position = "bottom", axis.text.x  = element_text(angle=45, hjust=1, vjust=0.9))

Fig. 14

我们发现：

• 有些地区有很多餐厅，种类繁多，而另一些地区只有一家air餐厅（参数空间的大部分是空的）。

• 同样，有些美食如“Izakaya”或“Cafe”相当普遍，而其他美食则只能在少数地区找到。请注意，air样本中仅有的 2 家Karaoke 酒吧位于“Hokkaido Sapporo-shi Minami 3 Jonishi”，而仅有的 2 家“International cuisine”餐厅以及仅有的 2 家“Asian”餐厅则位于“Tokyo-to Shibuya-ku Shibuya”。

HPG 数据的相同类型的图看起来很相似，但是由于类型数量较多而更加复杂：

hpg_store %>%
  mutate(area = str_sub(hpg_area_name, 1, 10)) %>%
  ggplot(aes(area, hpg_genre_name)) +
  geom_count(colour = "red") +
  theme(legend.position = "bottom", axis.text.x  = element_text(angle=45, hjust=1, vjust=0.9))

Fig. 15

我们发现：

• 这里也有繁忙的地区和只有几家餐馆的地区。毫不奇怪的是，Tokyo在餐厅类型多样性方面表现最突出。

• “Japanese style” 和 “International cuisine” 在每个地方都很受欢迎。 “Amusement bars” and “Udon/Soba” 很少, 例如 “Shanghai food” 或者 “Dim Sum”。

计数图告诉我们，每个区域可以找到某种类型的餐厅的数量分布。在这里，我们通过带有叠加抖动图的箱线图详细查看这些分布。类型按每个区域的平均案例数递减排序，即计数图中水平点序列的平均值。我们叠加各个数据点，并为每个点分配随机抖动，以在视觉上分离其他重叠的数据。这里，y 轴（即“每个区域的出现次数”）对应于上面计数图中的点的大小。我们在这里使用单个图，而不是面板，因为这些图非常详细（注意对数 y 轴）。

我们从air开始

air_store %>%
  group_by(air_genre_name, air_area_name) %>%
  count() %>%
  ggplot(aes(reorder(air_genre_name, n, FUN = mean), n)) +
  geom_boxplot() +
  geom_jitter(color = "blue") +
  scale_y_log10() +
  coord_flip() +
  labs(x = "Air genre", y = "Occurences per air area")

Fig. 16

我们发现：

• 只有少数类型的每个区域的餐厅中位数超过 2 家。例如“Italian/French”或“Bar/Cocktail”，每个区域更有可能出现 2 人以上的餐厅集群。

• 对于大多数类型来说，分布紧密地集中在每个区域 2 个周围，但也有一些分散到更高的数量。 “Cafes”的数量最多，在一个地区出现了 26 处（Fukuoka-ken Fukuoka-shi Daimyō）。

• 奇怪的是，这里的最小值是 2，而不是 1。这意味着在任何地区都没有一家air餐厅是唯一的一种类型。您也可以在图 5 的地图上看到这一点：只要同一区域坐标有一组 2 家餐厅，它们就具有相同的类型（但有不同的 air_store_id）。同样，在每个区域的任何更大的区域中都不存在单一的类型。

• 不过，存在奇数次出现，即每个区域有 3 或 5 个相同类型。我不太确定该怎么做，因为它看起来太完美匹配了，不太真实。这是否与air数据的选择方式有关？这可能是分配类型名称时的错误吗？我检查了一个地点的 2 家餐厅的几个例子，它们在同一天有不同的顾客数量，例如对于这对“Dining Bars”，表明我们不仅仅存在单个条目重复的问题：

air_store %>%
  filter(air_store_id == "air_b5598d12d1b84890" | air_store_id == "air_bbe1c1a47e09f161")

## # A tibble: 2 × 5
##   air_store_id         air_genre_name air_area_name           latitude longitude
##   <chr>                <fct>          <fct>                      <dbl>     <dbl>
## 1 air_bbe1c1a47e09f161 Dining bar     Tōkyō-to Setagaya-ku K…     35.7      140.
## 2 air_b5598d12d1b84890 Dining bar     Tōkyō-to Setagaya-ku K…     35.7      140.

air_visits %>%
  filter(air_store_id == "air_b5598d12d1b84890" | air_store_id == "air_bbe1c1a47e09f161") %>%
  arrange(visit_date) %>%
  head(10)

## # A tibble: 10 × 3
##    air_store_id         visit_date visitors
##    <chr>                <date>        <dbl>
##  1 air_bbe1c1a47e09f161 2016-07-01        7
##  2 air_b5598d12d1b84890 2016-07-01        5
##  3 air_b5598d12d1b84890 2016-07-02       14
##  4 air_b5598d12d1b84890 2016-07-03        6
##  5 air_b5598d12d1b84890 2016-07-04        6
##  6 air_b5598d12d1b84890 2016-07-05        5
##  7 air_b5598d12d1b84890 2016-07-06        5
##  8 air_bbe1c1a47e09f161 2016-07-07        1
##  9 air_bbe1c1a47e09f161 2016-07-08        7
## 10 air_b5598d12d1b84890 2016-07-08       10

现在我们看看 HPG 餐厅的相同分布：

foobar <- hpg_store %>%
  group_by(hpg_genre_name, hpg_area_name) %>%
  count()

foobar %>%
  ggplot(aes(reorder(hpg_genre_name, n, FUN = mean), n)) +
  geom_boxplot() +
  geom_jitter(color = "red") +
  scale_y_log10() +
  coord_flip() +
  labs(x = "hpg genre", y = "Cases per hpg area")

Fig. 17

我们发现：

• 在这里，我们显然每个区域至少有 1 个类型，而且由于总体数量较高，类型中值的种类也更多。

• 最极端的流派是“Japanese style”，每个地区的餐厅数量中位数略高于 10 家。除此之外，还有许多其他类型的the lower box hinge不接触每个区域的最小 1 个的情况。

使用每个区域中流派数量的信息，我们现在可以继续量化数据集的聚类或“拥挤度”，并将其与顾客数量相关联。下一个图首先显示最后两个图的air和 HPG 数据点的总体分布（即每个区域相同顾客的情况）。

此外，我们还估计了每个聚类案例的平均访客数量。为此，我们首先取 log1p 转换的顾客数量的平均值（每个类型和区域），然后计算每个案例的这些数字的平均值和标准差，即同一类型在一个区域中出现的次数。 （log1p的意思是log(x+1)，目的是使访客数量分布更加正态，见图1）。

foo <- air_visits %>%
  left_join(air_store, by = "air_store_id")

bar <- air_store %>%
  group_by(air_genre_name, air_area_name) %>%
  count()

foobar <- hpg_store %>%
  group_by(hpg_genre_name, hpg_area_name) %>%
  count()

p1 <- bar %>%
  ggplot(aes(n)) +
  geom_histogram(fill = "blue", binwidth = 1) +
  labs(x = "Air genres per area")

p2 <- foobar %>%
  ggplot(aes(n)) +
  geom_histogram(fill = "red", binwidth = 1) +
  labs(x = "HPG genres per area")

p3 <- foo %>%
  group_by(air_genre_name, air_area_name) %>%
  summarise(mean_log_visit = mean(log1p(visitors))) %>%
  left_join(bar, by = c("air_genre_name","air_area_name")) %>%
  group_by(n) %>%
  summarise(mean_mlv = mean(mean_log_visit),
            sd_mlv = sd(mean_log_visit)) %>%
  replace_na(list(sd_mlv = 0)) %>%
  ggplot(aes(n, mean_mlv)) +
  geom_point(color = "blue", size = 4) +
  geom_errorbar(aes(ymin = mean_mlv - sd_mlv, ymax = mean_mlv + sd_mlv), width = 0.5, size = 0.7, color = "blue") +
  labs(x = "Cases of identical Air genres per area", y = "Mean +/- SD of\n mean log1p visitors")

layout <- matrix(c(1,2,3,3),2,2,byrow=TRUE)
multiplot(p1, p2, p3, layout=layout)

Fig. 18

我们发现：

• 上图：同一区域中相同类型的案例数量分别从air和 hpg数据集的最大值“2”和“1”快速下降，可能呈指数下降。例子数越大，尾部越长，air的第二个峰值可能在“8”左右，hpg 的第二个峰值可能在“13”左右。

• 下图：顾客数量的 log1p 均值显示出相对较大的分布，并且在“2”个案例到“9”个案例的范围内完全一致。从 10 个案例以上，我们根本没有数据来衡量价差：每个区域有两个数据点，每个区域有 11 个案例，否则它只是一次测量。然而，值得注意的是，所有数据点 > 9 个案例之间的分散度相当低，并且明显高于 <= 9 点的平均值。现在，这可能仅仅意味着那些高访客数量没有被考虑在内。被小（不那么繁忙？）区域“推倒”，但这本身就是一个有趣的结果。

• 请注意，下部面板混合中的散点图处理所有聚类/拥挤情况，无论类型如何。当然，我们只能为有顾客人数的air数据绘制此图。

我们将在特征工程部分尝试另一种方法，通过与邻近餐厅的距离来量化餐厅的“独特性”。

预定与顾客人数

接下来我们将注意力转向air_reserve和hpg_reserve数据集中的预订。我们已经在 4.2 和 4.3 节中看到了它们的时间序列和分布；现在我们将预订人数与实际顾客人数进行比较。

为此，我们计算每个餐厅 ID 每天的 reserve_visitors 总和（即预订人数），然后将这些摘要加入到 air_visitors 文件中。为了包含 hpg 预订，我们需要使用 store_ids 数据将 hpg_reserve 文件中的 hpg_store_ids 连接到相应的 air_store_ids：

foo <- air_reserve %>%
  mutate(visit_date = date(visit_datetime)) %>%
  group_by(air_store_id,visit_date) %>%
  summarise(reserve_visitors_air = sum(reserve_visitors))
  
bar <- hpg_reserve %>%
  mutate(visit_date = date(visit_datetime)) %>%
  group_by(hpg_store_id,visit_date) %>%
  summarise(reserve_visitors_hpg = sum(reserve_visitors)) %>%
  inner_join(store_ids, by = "hpg_store_id")

all_reserve <- air_visits %>%
  inner_join(foo, by = c("air_store_id", "visit_date")) %>%
  inner_join(bar, by = c("air_store_id", "visit_date")) %>%
  mutate(reserve_visitors = reserve_visitors_air + reserve_visitors_hpg)

现在，我们将绘制预订游客总数与air餐厅实际游客数量的关系图。我们使用散点图，通过 ggExtra 包的 ggMarginal 函数添加边缘直方图。灰线显示 reserve_visitors == visitors，蓝线是线性拟合：

p <- all_reserve %>%
  filter(reserve_visitors < 120) %>%
  ggplot(aes(reserve_visitors, visitors)) +
  geom_point(color = "black", alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, color = "grey60") +
  geom_smooth(method = "lm", color = "blue")
ggMarginal(p, type="histogram", fill = "blue", bins=50)

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Fig. 19

我们发现：

• 直方图显示，预定人数和实际顾客数量的峰值低于约 20，并且很大程度上限制在 100 以下的范围内。

• 散点大部分落在标识线上方，表明当天的顾客数量多于预订的餐桌。这并不奇怪，因为总会有一定数量的人成为上门顾客。

• 值得注意的一部分点位于线下方，这可能表明有些人进行了预订，但改变了主意，没有去。这种效果可能是可以预料到的，考虑到这一点将是本次研究的挑战之一。

• 线性拟合表明了一种趋势，即较大数量的预定人数更有可能低估最终的顾客数量。这也并不奇怪，因为我可以想象，（a）大量预订被取消（b）一大群人在没有预订的情况下走进餐厅。前者的可能性更大。

现在，我们将随着时间的推移分解访visitors - reserve_visitors之间的差异，查看整体直方图，并分别可视化 air_reserve 与 hpg_reserve 数字。此处，air（蓝色）和 hpg（红色）的时间序列垂直偏移 150 和 -250（参见实心黑色基线）：

p1 <- all_reserve %>%
  ggplot(aes(visitors - reserve_visitors)) +
  geom_histogram(binwidth = 5, fill = "black") +
  coord_flip() +
  labs(x = "")

p2 <- all_reserve %>%
  ggplot(aes(visitors - reserve_visitors_air)) +
  geom_histogram(binwidth = 5, fill = "blue") +
  coord_flip() +
  labs(x = "")

p3 <- all_reserve %>%
  ggplot(aes(visitors - reserve_visitors_hpg)) +
  geom_histogram(binwidth = 5, fill = "red") +
  coord_flip() +
  labs(x = "")

p4 <- all_reserve %>%
  ggplot(aes(visit_date, visitors - reserve_visitors)) +
  geom_hline(yintercept = c(150, 0, -250)) +
  geom_line() +
  geom_line(aes(visit_date, visitors - reserve_visitors_air + 150), color = "blue") +
  geom_line(aes(visit_date, visitors - reserve_visitors_hpg - 250), color = "red") +
  ggtitle("Visitors - Reserved: all (black), air (blue), hpg (red)")

layout <- matrix(c(4,4,2,4,4,1,4,4,3),3,3,byrow=TRUE)
multiplot(p1, p2, p3, p4, layout=layout)

Fig. 20

我们发现：

• 时间序列在整个数据集时间范围内显示出显着的分散性。虽然air（蓝色）和 hpg（红色）曲线主要高于基线（即顾客数量多于预订量），但将两个数据集结合起来使分布平均值更接近零线。这也可以从右侧相应的直方图中看出。

• 我们看到没有air预订的差距（比较第 4.2 节）。我们只能在这里查看 hpg 预订量（这种差距不存在，第 4.3 节），但似乎可以安全地假设它们将遵循相同的趋势，并且可以用作air预订量的代理。可以随意检查这个假设是否存在差距。

• 为了方便解释，右侧 3 个图中的（翻转的）直方图与左侧图中的时间序列大致对齐。它们展示了与餐厅预定数量相比，分布向更大的顾客数量倾斜的程度。我们可能会在这里看到两种分布之间的混合：由于取消而导致的（可能是正常的）价差加上来自步入式顾客的尾部，它应该遵循泊松分布。

最后，让我们快速了解一下假期对预订量和顾客数量量差异的影响。我们将使用重叠密度图：

all_reserve %>%
  mutate(date = visit_date) %>%
  left_join(holidays, by = "date") %>%
  ggplot(aes(visitors - reserve_visitors, fill = holiday_flg)) +
  geom_density(alpha = 0.5)

Fig. 21

我们发现：

• 与假日预订相比，游客数量略多。峰值几乎相同，但我们看到较大数字存在微小但明显的差异。

特征工程

下一步是从现有特征中派生出新特征，这些新特征的目的是为我们预测访客数量的目标提供额外的预测能力。这些新功能可以像从日期列中导出星期几或月份一样简单；我们已经在图 1 的视觉探索中使用过它。或者，它们可以通过几个相关变量的相互作用而呈现出更复杂的形状。本节收集并研究这些新功能。

这里将所有工程功能收集到一个代码块中，这样我们就不必在内核的不同位置搜索它们。通常，这也使得这些特征的计算变得更容易，因为我们可以重复使用某些中间转换。当我们扩展分析时，我们将回到这个代码块，它会随着我们的特征空间的增长而增长：

summary(air_reserve)

##  air_store_id       visit_datetime                  
##  Length:90121       Min.   :2016-01-01 19:00:00.00  
##  Class :character   1st Qu.:2016-11-15 18:00:00.00  
##  Mode  :character   Median :2017-01-05 18:00:00.00  
##                     Mean   :2016-12-05 02:21:09.32  
##                     3rd Qu.:2017-03-03 19:00:00.00  
##                     Max.   :2017-05-31 21:00:00.00  
##  reserve_datetime                 reserve_visitors
##  Min.   :2016-01-01 01:00:00.00   Min.   :  1.00  
##  1st Qu.:2016-11-07 15:00:00.00   1st Qu.:  2.00  
##  Median :2016-12-27 22:00:00.00   Median :  3.00  
##  Mean   :2016-11-26 20:52:08.02   Mean   :  4.47  
##  3rd Qu.:2017-02-26 18:00:00.00   3rd Qu.:  5.00  
##  Max.   :2017-04-22 23:00:00.00   Max.   :100.00

air_visits <- air_visits %>%
  mutate(wday = wday(visit_date, label=TRUE, week_start = 1),
         wday = fct_relevel(wday, c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun")),
         month = month(visit_date, label=TRUE))

air_reserve <- air_reserve %>%
  mutate(reserve_date = date(reserve_datetime),
         reserve_hour = hour(reserve_datetime),
         reserve_wday = wday(reserve_datetime, label = TRUE, week_start = 1),
         reserve_wday = fct_relevel(reserve_wday, c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun")),
         visit_date = date(visit_datetime),
         visit_hour = hour(visit_datetime),
         visit_wday = wday(visit_datetime, label = TRUE, week_start = 1),
         visit_wday = fct_relevel(visit_wday, c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun")),
         diff_hour = time_length(visit_datetime - reserve_datetime, unit = "hour"),
         diff_day = time_length(visit_datetime - reserve_datetime, unit = "day"))

hpg_reserve <- hpg_reserve %>%
  mutate(reserve_date = date(reserve_datetime),
         reserve_hour = hour(reserve_datetime),
         reserve_wday = wday(reserve_datetime, label = TRUE, week_start = 1),
         reserve_wday = fct_relevel(reserve_wday, c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun")),
         visit_date = date(visit_datetime),
         visit_hour = hour(visit_datetime),
         visit_wday = wday(visit_datetime, label = TRUE, week_start = 1),
         visit_wday = fct_relevel(visit_wday, c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun")),
         diff_hour = time_length(visit_datetime - reserve_datetime, unit = "hour"),
         diff_day = time_length(visit_datetime - reserve_datetime, unit = "day"))


# count stores in area
air_count <- air_store %>%
  group_by(air_area_name) %>%
  summarise(air_count = n())

hpg_count <- hpg_store %>%
  group_by(hpg_area_name) %>%
  summarise(hpg_count = n())

# distances
med_coord_air <- air_store %>%
  summarise_at(vars(longitude:latitude), median)
med_coord_hpg <- hpg_store %>%
  summarise_at(vars(longitude:latitude), median)

air_coords <- air_store %>%
  select(longitude, latitude)
hpg_coords <- hpg_store %>%
  select(longitude, latitude)

air_store$dist <- distCosine(air_coords, med_coord_air)/1e3
hpg_store$dist <- distCosine(hpg_coords, med_coord_hpg)/1e3

# apply counts, dist; add prefecture
air_store <- air_store %>%
  mutate(dist_group = as.integer(case_when(
    dist < 80 ~ 1,
    dist < 300 ~ 2,
    dist < 500 ~ 3,
    dist < 750 ~ 4,
    TRUE ~ 5))) %>%
  left_join(air_count, by = "air_area_name") %>%
  separate(air_area_name, c("prefecture"), sep = " ", remove = FALSE)

hpg_store <- hpg_store %>%
  mutate(dist_group = as.integer(case_when(
    dist < 80 ~ 1,
    dist < 300 ~ 2,
    dist < 500 ~ 3,
    dist < 750 ~ 4,
    TRUE ~ 5))) %>%
  left_join(hpg_count, by = "hpg_area_name") %>%
  separate(hpg_area_name, c("prefecture"), sep = " ", remove = FALSE)

周几和几月

我们从简单的开始，使用从 visit_date 派生的工作日 (wday) 和月份特征。我们已经在图 1 中查看了每周、每天和每月的访客总数。现在我们将研究 log1p 转换后的顾客数量。我们计算它们的平均值和标准差，并将它们沿着（山脊线）密度分布绘制：

p1 <- air_visits %>%
  group_by(wday) %>%
  summarise(mean_log_visitors = mean(log1p(visitors)),
            sd_log_visitors = sd(log1p(visitors))) %>%
  ggplot(aes(wday, mean_log_visitors, color = wday)) +
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = mean_log_visitors - sd_log_visitors,
                    ymax = mean_log_visitors + sd_log_visitors,
                    color = wday), width = 0.5, size = 0.7) +
  theme(legend.position = "none")

p2 <- air_visits %>%
  mutate(visitors = log1p(visitors)) %>%
  ggplot(aes(visitors, wday, fill = wday)) +
  geom_density_ridges(bandwidth = 0.1) +
  scale_x_log10() +
  theme(legend.position = "none") +
  labs(x = "log1p(visitors)", y = "")

p3 <- air_visits %>%
  group_by(month) %>%
  summarise(mean_log_visitors = mean(log1p(visitors)),
            sd_log_visitors = sd(log1p(visitors))) %>%
  ggplot(aes(month, mean_log_visitors, color = month)) +
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = mean_log_visitors - sd_log_visitors,
                    ymax = mean_log_visitors + sd_log_visitors,
                    color = month), width = 0.5, size = 0.7) +
  theme(legend.position = "none")

p4 <- air_visits %>%
  mutate(visitors = log1p(visitors)) %>%
  ggplot(aes(visitors, month, fill = month)) +
  geom_density_ridges(bandwidth = 0.1) +
  scale_x_log10() +
  theme(legend.position = "none") +
  labs(x = "log1p(visitors)", y = "")

layout <- matrix(c(1,2,3,4),2,2,byrow=TRUE)
multiplot(p1, p2, p3, p4, layout=layout)

Fig. 22

我们发现：

• 一周中各天之间的平均日志（log(visitor)）数量差异远小于相应的不确定性。显示强烈重叠的密度图证实了这一点。尽管如此，周末顾客数量有增加的趋势。

• months的情况更加相似，包括 12 月份的数字略有增加。总体来说，没有太大区别。

以air_genre_name细分:

air_visits %>%
  left_join(air_store, by = "air_store_id") %>%
  group_by(wday, air_genre_name) %>%
  summarise(mean_log_visitors = mean(log1p(visitors)),
            sd_log_visitors = sd(log1p(visitors))) %>%
  ggplot(aes(wday, mean_log_visitors, color = wday)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = mean_log_visitors - sd_log_visitors,
                    ymax = mean_log_visitors + sd_log_visitors,
                    color = wday), width = 0.5, size = 0.7) +
  theme(legend.position = "none", axis.text.x  = element_text(angle=45, hjust=1, vjust=0.9)) +
  facet_wrap(~ air_genre_name)

Fig. 23

我们发现：

• 对于“Karaoke”或“International Cuisine”等类型，周末与工作日的影响更大，而其他类型在一周中的几天内变化很小（例如“Japanese food”）。

• 总体而言，一周中的任何一天对任何类型都没有显着影响。

visitors vs reservations的统计数据如何？这取决于星期几吗？我们使用箱线图来找出：

all_reserve %>%
  mutate(wday = wday(visit_date, label=TRUE, week_start = 1),
         wday = fct_relevel(wday, c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"))) %>%
  ggplot(aes(wday, visitors - reserve_visitors, fill = wday)) +
  geom_boxplot() +
  theme(legend.position = "none")

Fig. 24

我们发现没有显着差异，尽管我们发现与周末（尤其是周日）的预订相比，顾客数量略有增加。

每个区域的餐厅的数量

此功能将解决某个区域可以找到多少家餐厅以及每个区域的餐厅数量过多是否会影响平均顾客数量的问题。

为了计算它，我们只需计算相应区域的 air/hpg 餐厅数量。我们在顶部的下一个图中显示了结果分布。我们估计每个餐厅经过log1p- 转换的平均顾客数量，然后计算每个区域计数组的该数字的平均值和标准差：

p1 <- air_count %>%
  ggplot(aes(air_count)) +
  geom_histogram(binwidth = 2, fill = "blue")

p2 <- hpg_count %>%
  ggplot(aes(hpg_count)) +
  geom_histogram(binwidth = 5, fill = "red")

p3 <- air_visits %>%
  left_join(air_store, by = "air_store_id") %>%
  group_by(air_store_id, air_count) %>%
  summarise(mean_store_visit = mean(log1p(visitors))) %>%
  group_by(air_count) %>%
  summarise(mean_log_visitors = mean(mean_store_visit),
            sd_log_visitors = sd(mean_store_visit)) %>%
  ggplot(aes(air_count, mean_log_visitors)) +
  geom_point(size = 4, color = "blue") +
  geom_errorbar(aes(ymin = mean_log_visitors - sd_log_visitors,
                    ymax = mean_log_visitors + sd_log_visitors),
                    color = "blue", width = 0.5, size = 0.7) +
  geom_smooth(method = "lm", color = "black") +
  labs(x = "Air restaurants per area")

layout <- matrix(c(1,2,3,3),2,2,byrow=TRUE)
multiplot(p1, p2, p3, layout=layout)

## `geom_smooth()` using formula = 'y ~ x'

Fig. 25

我们发现：

• air数据集中的大多数区域只有几家餐馆，分布峰值为 2。这可以追溯到之前的观察，即在air区域中找不到单一的餐厅类型，总是至少有 2 个相同类型。air 数据也暂定在16-20点左右出现第二个峰值。 hpg 数据总体数量较大，也在低计数时达到峰值，并且在下降过程中几乎没有其他较小的峰值。

• 每个区域的平均log1p顾客数量具有很大的不确定性并且彼此一致。黑色线性回归显示，较大区域的平均数存在轻微趋势，但趋势相当弱。

离最繁华区域的距离

每当我们的数据中有空间坐标时，我们就会自动获得这些坐标之间的距离。作为第一个近似，我们使用两点之间的线性距离（即直线距离）。

为了计算这些距离，我们使用 geosphere 包的 distCosine 函数进行球面三角学。该方法为我们提供了球形地球上两点之间的最短距离。为了进行局部分析，我们选择忽略地球形状的椭球变形。

我们选择中纬度和中经度作为所有距离的单一参考点。我们的分析可以扩展到使用相同的方法来调查两家餐厅之间的所有成对距离。我们不会在内核中进行这项研究，但我鼓励任何想要尝试的人看看是否可以通过这种方式获得更多见解。

在这里，我们绘制了所有air餐厅和 hpg 数据与中位数的距离的直方图。这些分布将建议分为 5 个不同的bins，其距离范围不断增加。然后我们计算 5 组的平均log1p 顾客数量并进行比较：

p1 <- air_store %>%
  ggplot(aes(dist)) +
  geom_histogram(bins = 30, fill = "blue") +
  geom_vline(xintercept = c(80, 300, 500, 750)) +
  labs(x = "Linear distance [km]")

p2 <- hpg_store %>%
  ggplot(aes(dist)) +
  geom_histogram(bins = 30, fill = "red") +
  geom_vline(xintercept = c(80, 300, 500, 750)) +
  labs(x = "Linear distance [km]")

p3 <- air_visits %>%
  left_join(air_store, by = "air_store_id") %>%
  group_by(air_store_id, dist_group) %>%
  summarise(mean_store_visit = mean(log1p(visitors))) %>%
  group_by(dist_group) %>%
  summarise(mean_log_visitors = mean(mean_store_visit),
            sd_log_visitors = sd(mean_store_visit)) %>%
  ggplot(aes(dist_group, mean_log_visitors)) +
  geom_point(size = 4, color = "blue") +
  geom_errorbar(aes(ymin = mean_log_visitors - sd_log_visitors,
                    ymax = mean_log_visitors + sd_log_visitors),
                    color = "blue", width = 0.5, size = 0.7) +
  labs(x = "Distance grouping")
  

layout <- matrix(c(1,2,3,3),2,2,byrow=TRUE)
multiplot(p1, p2, p3, layout=layout)

Fig. 26

我们发现：

• 距离分布在很短的值处有明显的峰值；这是可以预料的，因为我们的参考是中位数。我们还看到了更多的峰值，最突出的是约 450 公里和约 800 公里左右。这些将对应于（组）城市，并且 air 和 hpg 数据之间一致。我们根据 80、300、500 和 750 公里处的垂直黑线定义 5 个距离区间。

• 这 5 个距离区间的平均log1p-transformed顾客计数在其标准差内没有显着差异。

在这里，我们展示了air（蓝色）和 hpg（红色）餐厅的地图，其中圆圈的大小随着bin相对于（深绿色）中心参考坐标的距离的增加而增加。毫不奇怪，这个中间位置位于Tokyo：

foo <- air_store %>% select(latitude, longitude, dist_group) %>% mutate(dset = "air")
bar <- hpg_store %>% select(latitude, longitude, dist_group) %>% mutate(dset = "hpg")

leaflet(foo) %>%
  addTiles() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addScaleBar() %>%
  addCircleMarkers(~longitude, ~latitude, group = "AIR",
                   color = "blue", fillOpacity = 0.5, radius = 3*foo$dist_group) %>%
  addCircleMarkers(lng = bar$longitude, lat = bar$latitude, group = "HPG",
                   color = "red", fillOpacity = 0.5, radius = 3*bar$dist_group) %>%
  addCircleMarkers(lng = med_coord_air$longitude, lat = med_coord_air$latitude, group = "Centre",
                   color = "darkgreen", fillOpacity = 1) %>%
  addLayersControl(
    overlayGroups = c("AIR", "HPG", "Centre"),
    options = layersControlOptions(collapsed = FALSE)
  )

Fig. 27

区县名

这是一个简单的功能，我们只需将地区定义为 area_name 的第一部分。在这里，我们使用（通常）一致的区域名称格式，正如 Oscar Takeshita 在 discussion中所解释的那样（请参阅那里的一些注意事项）。这为我们提供了相对较少数量的不同高级区域：

p1 <- air_store %>%
  ggplot(aes(prefecture)) +
  geom_bar(fill = "blue") +
  coord_flip() +
  ggtitle("air prefectures - # restaurants") +
  labs(x = "")

p2 <- hpg_store %>%
  ggplot(aes(prefecture)) +
  geom_bar(fill = "red") +
  coord_flip() +
  ggtitle("hpg prefectures - # restaurants") +
  labs(x = "")

p3 <- air_visits %>%
  left_join(air_store, by = "air_store_id") %>%
  group_by(air_store_id, prefecture) %>%
  summarise(mean_store_visit = mean(log1p(visitors))) %>%
  group_by(prefecture) %>%
  summarise(mean_log_visitors = mean(mean_store_visit),
            sd_log_visitors = sd(mean_store_visit)) %>%
  ggplot(aes(prefecture, mean_log_visitors)) +
  geom_point(size = 4, color = "blue") +
  geom_errorbar(aes(ymin = mean_log_visitors - sd_log_visitors,
                    ymax = mean_log_visitors + sd_log_visitors),
                    color = "blue", width = 0.5, size = 0.7) +
  labs(x = "prefecture") +
  theme(axis.text.x  = element_text(angle=15, hjust=1, vjust=0.9))

layout <- matrix(c(1,2,1,2,1,2,3,3,3,3),5,2,byrow=TRUE)
multiplot(p1, p2, p3, layout=layout)

Fig. 28

我们发现：

• air 数据中有9个地市县，hpg 数据有13个条目。在这两个数据集中，“Tokyo”显然是最频繁的条目。

• 我们再次发现，9个不同县的平均log1p顾客数量没有显着差异。任何可能的变化都小于托盘线的大小。

时间序列参数

在根据不同餐厅的地理或烹饪特性设计新功能后，我们现在将直接查看其顾客人数的时间序列。由于只有来自 829 个air 餐厅的 829 个时间序列，实际上可以很容易地单独查看每个时间序列，并且有一些内核可以绘制其中的small 或 large以了解数据。

在这里，我们将采取不同的方法来查看每个时间序列的元参数。这些特征是时间序列的平均值、标准差、斜率和斜率误差。每个参数都是针对 log1p-transformed 数据计算的。在下一个代码块中，我们定义一个提取函数，然后为每个餐厅运行它。由于数据集很小，这几乎不需要时间：

foo <- air_visits %>%
  left_join(air_store, by = "air_store_id") %>%
  group_by(air_store_id, air_genre_name) %>%
  summarise(mean_log_visits = mean(log1p(visitors)),
            mean_log_visits = mean(log1p(visitors)),
            sd_log_visits = sd(log1p(visitors))) %>%
  ungroup()

params_ts1 <- function(rownr){
  bar <- air_visits %>%
    filter(air_store_id == foo$air_store_id[rownr])
  slope <- summary(lm(visitors ~ visit_date, data = bar))$coef[2]
  slope_err <- summary(lm(visitors ~ visit_date, data = bar))$coef[4]
  
  foobar <- tibble(
    air_store_id = foo$air_store_id[rownr],
    slope = slope,
    slope_err = slope_err
  )
  
  return(foobar)
}

params <- params_ts1(1)
for (i in seq(2,nrow(foo))){
  params <- bind_rows(params, params_ts1(i))
}

ts_params <- foo %>%
  left_join(params, by = "air_store_id")

现在我们将绘制参数及其一维和二维分布：

p1 <- ts_params %>%
  ggplot(aes(mean_log_visits)) +
  geom_histogram(bins = 50, fill = "blue")

p2 <- ts_params %>%
  ggplot(aes(sd_log_visits)) +
  geom_histogram(bins = 50, fill = "blue")

p3 <- ts_params %>%
  filter(slope < 0.5) %>%
  ggplot(aes(slope)) +
  geom_histogram(bins = 50, fill = "blue") +
  labs(x = "Slope < 0.5")

p4 <- ts_params %>%
  ggplot((aes(mean_log_visits, sd_log_visits))) +
  geom_point(size = 2, color = "blue")

p5 <- ts_params %>%
  ggplot((aes(slope, slope_err))) +
  geom_point(size = 2, color = "blue")

layout <- matrix(c(1,1,2,2,3,3,4,4,4,5,5,5),2,6,byrow=TRUE)
multiplot(p1, p2, p3, p4, p5, layout=layout)

Fig. 29

我们发现：

• 经过 log1p-transformed顾客计数的平均值和标准差的分布相对较宽。我认为我们不能假设它们中的任何一个呈正态，尽管分布相对对称。左下图中的mean vs sd cloud的形状支持这样的观点：两个特征之间存在关系，平均顾客数量越大，方差越小。

• Slope分布窄；以零斜率为中心。右下面板中的slope vs slope error plot 图显示，存在一些具有较大slope errors且大多数绝对slope values较高的异常值：

ts_params %>%
  filter(abs(slope) > 0.25) %>%
  select(air_store_id, air_genre_name, slope, slope_err)

## # A tibble: 5 × 4
##   air_store_id         air_genre_name                slope slope_err
##   <chr>                <fct>                         <dbl>     <dbl>
## 1 air_1c0b150f9e696a5f Okonomiyaki/Monja/Teppanyaki -0.305    0.217 
## 2 air_8110d68cc869b85e Izakaya                       0.382    0.0984
## 3 air_900d755ebd2f7bbd Italian/French                2.06     0.385 
## 4 air_965b2e0cf4119003 Izakaya                       0.277    0.0265
## 5 air_9c6787aa03a45586 Cafe/Sweets                   0.832    0.150

最后，我们将使用构面缩放视图全面了解mean log1p visitors与时间序列斜率之间的关系，如 ggforce package中实现的那样。在这里，我们将看到右侧的完整空间和左侧的放大版本。缩放区域在右侧以较暗的阴影表示。这种视图非常适合这样的情况，在这种情况下，一些极值点会导致很难关注大多数数据集群的属性。The mean_log_visits 和slope按 genre 着色，相关误差线为灰色：

ts_params %>%
  ggplot(aes(mean_log_visits, slope, color = air_genre_name)) +
  geom_errorbarh(aes(xmin = mean_log_visits - sd_log_visits,
                    xmax = mean_log_visits + sd_log_visits),
                    color = "grey70", size = 0.7) +
  geom_errorbar(aes(ymin = slope - slope_err,
                    ymax = slope + slope_err),
                    color = "grey70", size = 0.7) +
  geom_point() +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow = 3, override.aes = list(size = 4))) +
  labs(color = "") +
  facet_zoom(y = (slope < 0.05 & slope > -0.1))

Fig. 30

我们发现：

• 尽管平均log1p访问量最低的三个点都是dining bars，但就其价值而言，genre聚类并不多。我们在概览图中找到了后面的slope异常值，并注意到它们都有很大的 mean 。

• .两个面板都显示了较大slopes与较高means相关的总体趋势。这可能并不奇怪，因为我们使用的绝对 slope值对于高visitor数比低访客数具有更大的范围。请随意绘制mean归一化的slope，看看效果是否持续存在。

现在可以详细研究该图中的极值，我们将在下一节中选择其中的一些作为示例。

Forecasting methods and examples

最后，我们将重点转向时间序列预测方法。我们已经了解了很多关于我们的数据集及其功能的知识。下面几节将介绍基本的预测方法

ARIMA / auto.arima

一种流行的预测方法是自回归综合移动平均模型；short ARIMA 模型。这种模型由三个构建块组成，这些构建块由三个指数 p, d, q参数化为 ARIMA(p, d, q)：

• auto-regressive / p: 我们正在使用过去的数据来计算未来数据的回归模型。参数 p 表示滞后的范围；例如ARIMA(3,0,0) 在回归中包括 t-1、t-2 和 t-3 值以计算 t 处的值。

• integrated / d: .这是一个差分参数，它给出了减去时间序列的当前值和先前值的次数。差分消除了时间序列中的变化，因为它稳定了均值并消除了（季节性）趋势。这是必要的，因为如果消除大规模趋势，计算滞后（例如时间 t 和时间 t-1 之间的差）是最有意义的。方差（或变异量）（和自协方差）不随时间变化（即每天不变化）的时间序列称为平稳（stationary）。

• moving average / q: 该参数为我们提供了模型回归误差中包含的先前误差项的数量。

在这里，我们将使用 auto.arima 工具来估计每个单独时间序列所需的 ARIMA 参数。为了将数据提供给 auto.arima，我们需要使用ts 工具将它们转换为时间序列对象。我们还将通过 forecast package的 tsclean 函数添加一个清洁和异常值去除步骤。我们已经看到我们的数据包含强大的每周周期，这将是我们模型的预设参数之一。我们在转换数据时将包含这些知识。让我们逐步设置所有内容，并添加注释和解释，然后将其转换为函数。取消隐藏代码以查看其实现方式。

我们以第一个 air_store_id（“air_ba937bf13d40fb24”）为例。

air_id = "air_ba937bf13d40fb24"

为了测试我们的预测，我们将预测与我们最终任务预测相同的时间范围（4 月 23 日 - 5 月 31 日）。这里我们自动从测试预测范围的长度中提取这39天，并将其定义为我们的“预测长度”。

pred_len <- test %>%
  separate(id, c("air", "store_id", "date"), sep = "_") %>%
  distinct(date) %>%
  nrow()

我们选择预测训练样本的最后 39 天。这可能不是在这里，我们计算训练日期的上限，并从此值中减去“预测长度”，以定义 3 月 14 日验证样本的开始。我们还创建了所有访问日期的数据集，为许多有间隙的时间序列做准备。

max_date <- max(air_visits$visit_date)
split_date <- max_date - pred_len
all_visits <- tibble(visit_date = seq(min(air_visits$visit_date), max(air_visits$visit_date), 1))

接下来，我们提取特定air_store_id.的时间序列。我们通过 log1p 转换访客计数并连接所有 visit_dates的数据集。这给了我们一些NA，我们用总体中位数填充。中位数可能不是最好的选择，但它是一个明智的起点。大多数时间序列预测工具需要无间隙的连续时间序列；这就是我们在这一步中创建的内容。

foo <- air_visits %>%
  filter(air_store_id == air_id)

visits <- foo %>%
  right_join(all_visits, by = "visit_date") %>%
  mutate(visitors = log1p(visitors)) %>%
  replace_na(list(visitors = median(log1p(foo$visitors)))) %>%
  rownames_to_column()

使用这个新的时间序列，我们现在将数据分为训练集和验证集。

visits_train <- visits %>% filter(visit_date <= split_date)
visits_valid <- visits %>% filter(visit_date > split_date)

.现在是配件部分。如前所述，我们使用 ts 函数创建时间序列对象，并使用 tsclean工具删除异常值。我们还添加了每周频率。stepwise和approximation参数设置意味着该工具对所有模型参数执行更彻底和精确的搜索。这会增加计算时间，但对于我们的小数据集来说这并不重要。

arima.fit <- auto.arima(tsclean(ts(visits_train$visitors, frequency = 7)),
                        stepwise = FALSE, approximation = FALSE)

使用拟合的 ARIMA 模型，我们将forecast“预测长度”，包括置信区间。

arima_visits <- arima.fit %>% forecast(h = pred_len, level = c(50,95))
predictions_log <- log1p(arima_visits$mean)
actuals_log <- log1p(visits_valid$visitors)
# RMLSE
rmlse <- sqrt(mean((predictions_log - actuals_log)^2))
# MAE
mae <- mean(abs(predictions_log - actuals_log))
# RMSE
rmse <- sqrt(mean((exp(predictions_log) - exp(actuals_log))^2))

最后，我们绘制我们的预测。 ggplot2包的 autoplot 函数根据特定数据类型的属性创建绘图；这里是一个时间序列对象。预测的访客数量以深蓝色显示，浅蓝色表示置信范围。我们还将真实的验证计数添加为灰色：

rmse<-25.6165
mae<-33.0779
rmlse<-0.8037
arima_visits %>%
  autoplot +
  geom_line(aes(as.integer(rowname)/7, visitors), data = visits_valid, color = "grey40") +
  labs(x = "Time [weeks]", y = "log1p visitors vs auto.arima predictions")

Fig. 31

print(paste("RMLSE:", rmlse))

## [1] "RMLSE: 0.8037"

print(paste("MAE:", mae))

## [1] "MAE: 33.0779"

print(paste("RMSE:", rmse))

## [1] "RMSE: 25.6165"

我们发现预测的前几天非常吻合，但随后我们的预测无法捕捉到更大的峰值。尽管如此，它仍然是比较其他方法的一个有用的起点。

现在我们把这个过程变成一个函数，包括绘图部分。

plot_auto_arima_air_id <- function(air_id){

  pred_len <- test %>%
    separate(id, c("air", "store_id", "date"), sep = "_") %>%
    distinct(date) %>%
    nrow()

  max_date <- max(air_visits$visit_date)
  split_date <- max_date - pred_len
  all_visits <- tibble(visit_date = seq(min(air_visits$visit_date), max(air_visits$visit_date), 1))
  
  foo <- air_visits %>%
    filter(air_store_id == air_id)

  visits <- foo %>%
    right_join(all_visits, by = "visit_date") %>%
    mutate(visitors = log1p(visitors)) %>%
    replace_na(list(visitors = median(log1p(foo$visitors)))) %>%
    rownames_to_column()
  
  visits_train <- visits %>% filter(visit_date <= split_date)
  visits_valid <- visits %>% filter(visit_date > split_date)

  arima.fit <- auto.arima(tsclean(ts(visits_train$visitors, frequency = 7)),
                          stepwise = FALSE, approximation = FALSE)

  arima_visits <- arima.fit %>% forecast(h = pred_len, level = c(50,95))

  arima_visits %>%
    autoplot +
    geom_line(aes(as.integer(rowname)/7, visitors), data = visits_valid, color = "grey40") +
    labs(x = "Time [weeks]", y = "log1p visitors vs forecast")
}

我们将此函数应用于几个时间序列，包括上一节中的两个slope异常值：

p1 <- plot_auto_arima_air_id("air_f3f9824b7d70c3cf")
p2 <- plot_auto_arima_air_id("air_8e4360a64dbd4c50")
p3 <- plot_auto_arima_air_id("air_1c0b150f9e696a5f")
p4 <- plot_auto_arima_air_id("air_900d755ebd2f7bbd")

layout <- matrix(c(1,2,3,4),2,2,byrow=TRUE)
multiplot(p1, p2, p3, p4, layout=layout)

Fig. 31

我们发现：

• 上图中的两个时间序列是合理完整的，但我们看到长间隙（以及我们的中值填充）导致左上图中的预测出现问题，我们失去了每周的周期性。右上面板保留了这种周期性，并且第一天的预测相对不错，但随后我们很快就低估了变化的幅度。

• 下图包括上面时间序列参数空间中的两个异常值；在这里，我们看到了事情确实出错的情况。如果这些特殊的时间序列在测试数据集中包含足够大的访问部分，则可能会导致任何其他不错的预测算法的性能不佳。

总体而言，结果不是很好，但考虑到它是全自动预测（假设仅每周一次），auto.arima 工具为我们提供了第一个基线来比较其他方法。

要更详细地探索 ARIMA 模型，请查看 TimLee 的Kernel 。

Prophet

prophet 预测工具是 Facebook开发的一款开源软件，在 R 和Python都可以使用。

Prophet 利用加性回归模型，将时间序列分解为 (i) （分段）线性/逻辑趋势、(ii) 每年季节性成分、(iii) 每周季节性成分和 (iv) 重要日期的可选列表（例如假期、特殊活动……）。它对缺失数据、趋势变化和大异常值具有良好的鲁棒性。尤其是缺失的数据功能在本次研究中可能会很有用。

让我们再次逐步探索该工具。我们将在 ARIMA 部分中继续我们的工作，并且不会重复在那里可以找到的任何解释。

我们将再次创建与上述相同时期（即 3 月 14 日之前/之后）的训练和验证集。与我们的 ARIMA 方法的唯一区别是：

• 我们不需要替换 NA 值，因为prophet 知道如何处理这些值。

• Prophet 需要一个包含两列的数据框：ds 表示日期，y 表示时间序列变量。

air_id = "air_ba937bf13d40fb24"

pred_len <- test %>%
  separate(id, c("air", "store_id", "date"), sep = "_") %>%
  distinct(date) %>%
  nrow()

max_date <- max(air_visits$visit_date)
split_date <- max_date - pred_len
all_visits <- tibble(visit_date = seq(min(air_visits$visit_date), max(air_visits$visit_date), 1))

foo <- air_visits %>%
  filter(air_store_id == air_id)

visits <- foo %>%
  right_join(all_visits, by = "visit_date") %>%
  mutate(visitors = log1p(visitors)) %>%
  rownames_to_column() %>%
  select(y = visitors,
         ds = visit_date)

visits_train <- visits %>% filter(ds <= split_date)
visits_valid <- visits %>% filter(ds > split_date)

这里我们拟合 prophet 模型并进行预测：

• .参数changepoint.prior.scale调整趋势灵活性。增加此参数可以使拟合更加灵活，但也会增加预测的不确定性，并使其更有可能过度拟合噪声。数据中的变更点会自动检测，除非使用变更点参数手动指定（我们在这里不这样做）。

• 参数yearly.seasonality必须明确启用/禁用，并允许prophet 注意到大规模周期。我们这里只有一年的数据，这绝对不足以找到年度周期，也可能不足以识别月份时间尺度的变化。可以随意测试该参数的性能。

proph <- prophet(visits_train, changepoint.prior.scale=0.5, yearly.seasonality=FALSE, daily.seasonality = FALSE)
future <- make_future_dataframe(proph, periods = pred_len)
fcast <- predict(proph, future)

ignore_code_block_1 <- TRUE
if (!ignore_code_block_1)
{rmsle_values <- c()
mae_values <- c()
rmse_values <- c()

for (i in 1:length(visits_valid$visitors)) {
  # 计算对数值
  prediction_log <- log1p(fcast$yhat[i])
  actual_log <- log1p(visits_valid$visitors[i])
  
  # 计算 RMLSE、MAE 和 RMSE
  rmsle_values <- c(rmsle_values, sqrt(mean((prediction_log - actual_log)^2)))
  mae_values <- c(mae_values, mean(abs(prediction_log - actual_log)))
  rmse_values <- c(rmse_values, sqrt(mean((exp(prediction_log) - exp(actual_log))^2)))
}

# 计算平均值
mean_rmsle <- mean(rmsle_values)
mean_mae <- mean(mae_values)
mean_rmse <- mean(rmse_values)}

这是标准的prophet 预测图：

plot(proph, fcast)

Fig. 32

mean_rmsle<-0.7945
mean_mae<-27.7850
mean_rmse<-20.6603

print(paste("Mean RMLSE:", mean_rmsle))

## [1] "Mean RMLSE: 0.7945"

print(paste("Mean MAE:", mean_mae))

## [1] "Mean MAE: 27.785"

print(paste("Mean RMSE:", mean_rmse))

## [1] "Mean RMSE: 20.6603"

观察到的数据绘制为黑点，拟合模型和预测绘制为蓝线。在浅蓝色中我们看到了相应的不确定性。

Prophet 提供了一个分解图，我们可以在其中检查模型的附加组件：趋势、年度季节性（如果包含）和每周周期：

prophet_plot_components(proph, fcast)

Fig. 33

我们发现：

• Prophet 检测到每周的变化模式，这与我们之前发现的类似，周五/周六比一周其他时间更受欢迎。但不同的是，周日这里的平均顾客数量非常低。这可能是由于这家特殊餐厅的属性，即“Tokyo”的“Dining Bar”。从图 23 中我们可以看出，“Dining Bar”在周日总体上并不那么繁忙；尽管没有我们在这张图中看到的那么大。周六与周日的顾客差异可能是最终模型的一个有趣特征。

• 长期趋势与图 11 中显示的平均行为不同，后者更有可能在 2016 年底上升。在这里，我们看到 2016 年中期达到峰值。

• 请注意，这是一个特定的时间序列。我还检查了其他时间序列，其中许多显示了预期的周五/周六/周日与周一至周四的变化，以及 2017 年 12 月左右的长期峰值。可以随意分叉此内核并亲自检查。

所以：这是一个好的开始，但是我们会使用我们值得信赖的 ggplot2 来更好地控制绘图参数：

fcast %>%
  as.tibble() %>%
  mutate(ds = date(ds)) %>%
  ggplot(aes(ds, yhat)) + 
  geom_ribbon(aes(x = ds, ymin = yhat_lower, ymax = yhat_upper), fill = "light blue") +
  geom_line(colour = "blue") +
  geom_line(data = visits_train, aes(ds, y), colour = "black") +
  geom_line(data = visits_valid, aes(ds, y), colour = "grey50")

Fig. 34

看起来还不错。我们将训练顾客计数绘制为黑色，验证集绘制为灰色，预测加上不确定性绘制为蓝色和浅蓝色。

让我们把它变成另一个函数：

plot_prophet_air_id <- function(air_id){
  
  pred_len <- test %>%
    separate(id, c("air", "store_id", "date"), sep = "_") %>%
    distinct(date) %>%
    nrow()

  max_date <- max(air_visits$visit_date)
  split_date <- max_date - pred_len
  all_visits <- tibble(visit_date = seq(min(air_visits$visit_date), max(air_visits$visit_date), 1))

  foo <- air_visits %>%
    filter(air_store_id == air_id)

  visits <- foo %>%
    right_join(all_visits, by = "visit_date") %>%
    mutate(visitors = log1p(visitors)) %>%
    rownames_to_column() %>%
    select(y = visitors,
          ds = visit_date)

  visits_train <- visits %>% filter(ds <= split_date)
  visits_valid <- visits %>% filter(ds > split_date)
  
  proph <- prophet(visits_train, changepoint.prior.scale=0.5,
                   yearly.seasonality=FALSE, daily.seasonality = FALSE)
  future <- make_future_dataframe(proph, periods = pred_len)
  fcast <- predict(proph, future)
  
  p <- fcast %>%
    as.tibble() %>%
    mutate(ds = date(ds)) %>%
    ggplot(aes(ds, yhat)) +
    geom_ribbon(aes(x = ds, ymin = yhat_lower, ymax = yhat_upper), fill = "light blue") +
    geom_line(colour = "blue") +
    geom_line(data = visits_train, aes(ds, y), colour = "black") +
    geom_line(data = visits_valid, aes(ds, y), colour = "grey50") +
    labs(title = str_c("Prophet for ", air_id))
  
  return(p)
}  

然后我们预测并绘制与上面相同的时间序列（除了“air_900d755ebd2f7bbd”，它没有足够的预测数据点，必须替换为“air_820d1919cbecaa0a”）：

p1 <- plot_prophet_air_id("air_f3f9824b7d70c3cf")
p2 <- plot_prophet_air_id("air_8e4360a64dbd4c50")
p3 <- plot_prophet_air_id("air_1c0b150f9e696a5f")

## Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.

## n.changepoints greater than number of observations. Using 9

p4 <- plot_prophet_air_id("air_820d1919cbecaa0a")

layout <- matrix(c(1,2,3,4),2,2,byrow=TRUE)
multiplot(p1, p2, p3, p4, layout=layout)

Fig. 35

看起来也不错。看起来我们高估了趋势成分，对于其中至少三个，趋势成分可能会使用较少的灵活性。第四个时间序列的数据太少，prophet无法做太多事情。在这里，我们可以完全放弃趋势成分或简单地取中位数。

同时，prophet还使我们有可能将假期和其他特殊事件纳入我们的预测中。对于我们的任务来说，这对于考虑黄金周可能非常有用。要了解有关黄金周影响的定制分析的更多见解，请查看 BreakfastPirate 的核心内容和相关讨论主题（kernel）。

由于假期数据在本次比赛中很容易获得，我们只需根据prophet的喜好对其进行重命名，即可将其包含在我们的预测函数中。让我们修改该函数，看看假期在预测 2016 年黄金周时是否会产生影响。我们唯一需要做的就是在中间步骤中截断 2016 年 5 月 31 日的 air_visits 数据：

plot_prophet_air_id_holiday <- function(air_id, use_hday){
  
  air_visits_cut <- air_visits %>%
    filter(visit_date <= ymd("20160531"))
  
  hday <- holidays %>%
    filter(holiday_flg == TRUE) %>%
    mutate(holiday = "holiday") %>%
    select(ds = date, holiday)
  
  pred_len <- test %>%
    separate(id, c("air", "store_id", "date"), sep = "_") %>%
    distinct(date) %>%
    nrow()

  max_date <- max(air_visits_cut$visit_date)
  split_date <- max_date - pred_len
  all_visits <- tibble(visit_date = seq(min(air_visits_cut$visit_date), max(air_visits_cut$visit_date), 1))

  foo <- air_visits_cut %>%
    filter(air_store_id == air_id)

  visits <- foo %>%
    right_join(all_visits, by = "visit_date") %>%
    mutate(visitors = log1p(visitors)) %>%
    rownames_to_column() %>%
    select(y = visitors,
          ds = visit_date)

  visits_train <- visits %>% filter(ds <= split_date)
  visits_valid <- visits %>% filter(ds > split_date)
  
  if (use_hday == TRUE){
    proph <- prophet(visits_train,
                     changepoint.prior.scale=0.5,
                     yearly.seasonality=FALSE,
                     daily.seasonality=FALSE,
                     holidays = hday)
    ptitle = "Prophet (w/ holidays) for "
  } else {
     proph <- prophet(visits_train,
                     changepoint.prior.scale=0.5,
                     yearly.seasonality=FALSE,
                     daily.seasonality = FALSE)
    ptitle = "Prophet for "
  }
  
  future <- make_future_dataframe(proph, periods = pred_len)
  fcast <- predict(proph, future)
  
  p <- fcast %>%
    as.tibble() %>%
    mutate(ds = date(ds)) %>%
    ggplot(aes(ds, yhat)) +
    geom_ribbon(aes(x = ds, ymin = yhat_lower, ymax = yhat_upper), fill = "light blue") +
    geom_line(colour = "blue") +
    geom_line(data = visits_train, aes(ds, y), colour = "black") +
    geom_line(data = visits_valid, aes(ds, y), colour = "grey50") +
    labs(title = str_c(ptitle, air_id))
  
  return(p)
}  

And then we take a well behaved time series and compare the predictions with and without holidays:

p1 <- plot_prophet_air_id_holiday("air_5c817ef28f236bdf", TRUE)
p2 <- plot_prophet_air_id_holiday("air_5c817ef28f236bdf", FALSE)

layout <- matrix(c(1,2),2,1,byrow=TRUE)
multiplot(p1, p2, layout=layout)

Fig. 36

我们发现：

• 包括假期在内，黄金周游客的适应度有了微妙的改善。如果训练集中包含更多假期，该组件的性能可能会提高；但我认为很难验证这一点。

Holt-Winters

更传统的时间序列过滤和预测是 Holt-Winters 算法，在stats包中实现。这是一种指数平滑方法，使用移动平均线来考虑数据中趋势的存在。这里我们在拟合和绘图函数中定义默认的季节性模型：

plot_hw_air_id <- function(air_id){

  pred_len <- test %>%
    separate(id, c("air", "store_id", "date"), sep = "_") %>%
    distinct(date) %>%
    nrow()

  max_date <- max(air_visits$visit_date)
  split_date <- max_date - pred_len
  all_visits <- tibble(visit_date = seq(min(air_visits$visit_date), max(air_visits$visit_date), 1))

  foo <- air_visits %>%
    filter(air_store_id == air_id)

  visits <- foo %>%
    right_join(all_visits, by = "visit_date") %>%
    mutate(visitors = log1p(visitors)) %>%
    replace_na(list(visitors = median(log1p(foo$visitors)))) %>%
    rownames_to_column()

  visits_train <- visits %>% filter(visit_date <= split_date)
  visits_valid <- visits %>% filter(visit_date > split_date)

  hw.fit <- HoltWinters(tsclean(ts(visits_train$visitors, frequency = 7)))

  hw_visits <- predict(hw.fit, n.ahead = pred_len, prediction.interval = T, level = 0.95) %>%
    as.tibble() %>%
    bind_cols(visits_valid)
# 提取预测值和实际值
predictions <- hw_visits$fit
actuals <- visits_valid$visitors

# 计算对数值
predictions_log <- log1p(predictions)
actuals_log <- log1p(actuals)

# 计算 RMLSE、MAE 和 RMSE
rmlse <- sqrt(mean((predictions_log - actuals_log)^2))
mae <- mean(abs(predictions_log - actuals_log))
rmse <- sqrt(mean((exp(predictions_log) - exp(actuals_log))^2))
  visits_train %>%
    ggplot(aes(visit_date, visitors)) +
    geom_line() +
    geom_ribbon(data = hw_visits, aes(x = visit_date, ymin = lwr, ymax = upr), fill = "light blue") +
    geom_line(data = hw_visits, aes(visit_date, visitors), color = "grey60") +
    geom_line(data = hw_visits, aes(visit_date, fit), color = "blue") +
    geom_line(data = hw_visits, aes(visit_date, fit), color = "blue") +
    labs(x = "Time [weeks]", y = "log1p visitors vs predictions") +
    ggtitle("HoltWinters")
}

然后我们将此函数应用于我们的示例时间序列：

rmlse<-0.7634
mae<-26.1840
rmse<-23.4591
plot_hw_air_id("air_ba937bf13d40fb24")

Fig. 37

# 输出结果
print(paste("RMLSE:", rmlse))

## [1] "RMLSE: 0.7634"

print(paste("MAE:", mae))

## [1] "MAE: 26.184"

print(paste("RMSE:", rmse))

## [1] "RMSE: 23.4591"

并绘制与上述prophet 相同的预测：

p1 <- plot_hw_air_id("air_f3f9824b7d70c3cf")
p2 <- plot_hw_air_id("air_8e4360a64dbd4c50")
p3 <- plot_hw_air_id("air_1c0b150f9e696a5f")
p4 <- plot_hw_air_id("air_820d1919cbecaa0a")

layout <- matrix(c(1,2,3,4),2,2,byrow=TRUE)
multiplot(p1, p2, p3, p4, layout=layout)

Fig. 38

我们发现：

• 虽然该算法在给定足够数据的情况下表现得相当好，但与 auto.arima类似，它在丢失数据的特殊情况下会崩溃。通过丢弃缺失值而不是用中位数替换它们，我们可能会获得更好的性能。然而，第一印象看起来不如prophet的方法有希望。

• 可以随意分割此内核并尝试修改后的加法模型和指数模型，以调整 Holt-Winters 方法的性能。

timetk

timetk包是最近发布的一个用于时间序列分析的工具包。它很好地集成到 tidyverse 生态系统中，并且专门设计用于与 tidy“tibble”数据框架配合使用。在这里，我们简要描述 timetk 方法以及如何将其应用到我们的数据中。

首先，我们将以与其他预测工具相同的方式创建训练和验证数据集：

air_id = "air_ba937bf13d40fb24"

pred_len <- test %>%
  separate(id, c("air", "store_id", "date"), sep = "_") %>%
  distinct(date) %>%
  nrow()

max_date <- max(air_visits$visit_date)
split_date <- max_date - pred_len
all_visits <- tibble(visit_date = seq(min(air_visits$visit_date), max(air_visits$visit_date), 1))

foo <- air_visits %>%
  filter(air_store_id == air_id)

visits <- foo %>%
  right_join(all_visits, by = "visit_date") %>%
  mutate(visitors = log1p(visitors)) %>%
  rownames_to_column() %>%
  select(y = visitors,
         ds = visit_date)

visits_train <- visits %>% filter(ds <= split_date)
visits_valid <- visits %>% filter(ds > split_date)

然后，我们使用tk_augment_timeseries_signatur工具通过时间序列特征来增强数据帧。这意味着我们将添加从日期中提取的综合时间序列属性。这些新功能包括一年中的月、日、周、半年和季度；日期时间数据精确到分钟和秒。在这里，我们展示了增强训练数据的一部分：

visits_train_aug <- visits_train %>%
  tk_augment_timeseries_signature()

## tk_augment_timeseries_signature(): Using the following .date_var variable: ds

visits_valid_aug <- visits_valid %>%
  .$ds %>%
  tk_get_timeseries_signature()

glimpse(visits_train_aug)

## Rows: 439
## Columns: 30
## $ y         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 3.258097, 3.…
## $ ds        <date> 2016-01-01, 2016-01-02, 2016-01-03, 2016-01-04, 2016-01-05,…
## $ index.num <dbl> 1451606400, 1451692800, 1451779200, 1451865600, 1451952000, …
## $ diff      <dbl> NA, 86400, 86400, 86400, 86400, 86400, 86400, 86400, 86400, …
## $ year      <int> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, …
## $ year.iso  <int> 2015, 2015, 2015, 2016, 2016, 2016, 2016, 2016, 2016, 2016, …
## $ half      <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ quarter   <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ month     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ month.xts <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ month.lbl <ord> 一月, 一月, 一月, 一月, 一月, 一月, 一月, 一月, 一月, 一月, …
## $ day       <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
## $ hour      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ minute    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ second    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ hour12    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ am.pm     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ wday      <int> 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, …
## $ wday.xts  <int> 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, …
## $ wday.lbl  <ord> 星期五, 星期六, 星期日, 星期一, 星期二, 星期三, 星期四, 星期…
## $ mday      <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
## $ qday      <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
## $ yday      <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
## $ mweek     <int> 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, …
## $ week      <int> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, …
## $ week.iso  <int> 53, 53, 53, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, …
## $ week2     <int> 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, …
## $ week3     <int> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, …
## $ week4     <int> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, …
## $ mday7     <int> 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, …

现在，timetk背后的想法是使用这些新功能基于回归或分类方法进行预测；使用标准工具，例如线性/逻辑回归或（增强/集成）树。对于这个例子，我们将使用一个简单的线性模型。这种方法可以很容易地扩展到更复杂的方法。

fit_lm <- lm(y ~ ., data = select(visits_train_aug, -c(ds, diff, wday.xts, wday.lbl, year.iso)))

然后我们使用标准预测工具将模型应用到我们的预测范围，确保包含与训练数据框中相同的特征。

pred <- predict(fit_lm, newdata = select(visits_valid_aug, -c(index, diff, wday.xts, wday.lbl, year.iso)))

pred_tk <- tibble(
    date  = visits_valid$ds,
    value = pred
    )

和以前一样，我们将预测与辅助函数中的绘图方法结合起来：

plot_tk_lm_air_id <- function(air_id){
  
  pred_len <- test %>%
    separate(id, c("air", "store_id", "date"), sep = "_") %>%
    distinct(date) %>%
    nrow()

  max_date <- max(air_visits$visit_date)
  split_date <- max_date - pred_len
  all_visits <- tibble(visit_date = seq(min(air_visits$visit_date), max(air_visits$visit_date), 1))

  foo <- air_visits %>%
    filter(air_store_id == air_id)

  visits <- foo %>%
    right_join(all_visits, by = "visit_date") %>%
    mutate(visitors = log1p(visitors)) %>%
    rownames_to_column() %>%
    select(y = visitors,
          ds = visit_date)

  visits_train <- visits %>% filter(ds <= split_date)
  visits_valid <- visits %>% filter(ds > split_date)
  
  # augment train with ts info
  visits_train_aug <- visits_train %>%
    tk_augment_timeseries_signature()
  
  # fit lm
  fit_lm <- lm(y ~ ., data = select(visits_train_aug, -c(ds, diff, wday.xts, wday.lbl, year.iso)))
  
  # augment valid with ts info
  visits_valid_aug <- visits_valid %>%
    .$ds %>%
    tk_get_timeseries_signature()
  
  # predict from lm
  pred <- predict(fit_lm, newdata = select(visits_valid_aug, -c(index, diff, wday.xts, wday.lbl, year.iso)))

  pred_tk <- tibble(
      date  = visits_valid$ds,
      y_pred = pred
      )
  
  # plot
  p <- pred_tk %>%
    ggplot(aes(date, y_pred)) +
    geom_line(data = visits_train, aes(ds, y), colour = "black") +
    geom_line(data = visits_valid, aes(ds, y), colour = "grey50") +
    geom_line(colour = "blue") +
    labs(title = str_c("timetk for ", air_id))
    
  return(p)
}  

然后我们将其应用到我们的默认示例时间序列中：

plot_tk_lm_air_id("air_ba937bf13d40fb24")

## tk_augment_timeseries_signature(): Using the following .date_var variable: ds

Fig. 37

我们发现：

• .结果看起来相当不错；特别是考虑到我们只使用简单的线性模型进行预测。更精细的方法可能会产生更好的结果。

请注意，从技术上讲，这不是一种时间序列方法。相反，我们使用 timetk 功能来提取数据的时间序列特征，并将它们用作回归式方法中的额外特征。因此，这种方法主要可以被视为（自动化）特征工程。我建议尝试不同的功能选择。

LSTM

这一节中将使用基于神经网络的方法，将整个时间序列加入对后续的预测。此处我们选择LSTM作为我们搭建模型的基础。这是一种较为复杂的预测方式。

观察数据的类型，并同时将数据中的周几和月份转换类型和内容，为了统一内容。

sum(is.na(air_visits))

## [1] 0

sum(is.na(air_reserve))

## [1] 0

sum(is.na(hpg_reserve))

## [1] 0

sum(is.na(air_store))

## [1] 0

sum(is.na(hpg_store))

## [1] 0

sum(is.na(holidays))

## [1] 0

sum(is.na(store_ids))

## [1] 0

sum(is.na(test))

## [1] 0

air_visits <- air_visits %>%
  mutate(month = case_when(
    month == "1月" ~ "1", 
    month == "2月" ~ "2",
    month == "3月" ~ "3",
    month == "4月" ~ "4",
    month == "5月" ~ "5",
    month == "6月" ~ "6",
    month == "7月" ~ "7",
    month == "8月" ~ "8",
    month == "9月" ~ "9",
    month == "9月" ~ "9",
    month == "10月" ~ "11",
    month == "11月" ~ "11",
    month == "12月" ~ "12",
  )) %>%
  mutate(month = as.integer(month))
air_visits <- air_visits %>%
  mutate(wday = case_when(
    wday == "周一" ~ "1", 
    wday == "周二" ~ "2",
    wday == "周三" ~ "3",
    wday == "周四" ~ "4",
    wday == "周五" ~ "5",
    wday == "周六" ~ "6",
    wday == "周日" ~ "7",
  )) %>%
  mutate(wday = as.integer(wday))
glimpse(air_visits)

## Rows: 252,108
## Columns: 5
## $ air_store_id <chr> "air_ba937bf13d40fb24", "air_ba937bf13d40fb24", "air_ba93…
## $ visit_date   <date> 2016-01-13, 2016-01-14, 2016-01-15, 2016-01-16, 2016-01-…
## $ visitors     <dbl> 25, 32, 29, 22, 6, 9, 31, 21, 18, 26, 21, 11, 24, 21, 26,…
## $ wday         <int> 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 3, 4, 5, …
## $ month        <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, …

air_reserve <- air_reserve %>%
  mutate(reserve_wday = case_when(
    reserve_wday == "周一" ~ "1", 
    reserve_wday == "周二" ~ "2",
    reserve_wday == "周三" ~ "3",
    reserve_wday == "周四" ~ "4",
    reserve_wday == "周五" ~ "5",
    reserve_wday == "周六" ~ "6",
    reserve_wday == "周日" ~ "7",
  )) %>%
  mutate(reserve_wday = as.integer(reserve_wday))
air_reserve <- air_reserve %>%
  mutate(visit_wday  = case_when(
    visit_wday  == "周一" ~ "1", 
    visit_wday  == "周二" ~ "2",
    visit_wday  == "周三" ~ "3",
    visit_wday  == "周四" ~ "4",
    visit_wday  == "周五" ~ "5",
    visit_wday  == "周六" ~ "6",
    visit_wday  == "周日" ~ "7",
  )) %>%
  mutate(visit_wday  = as.integer(visit_wday ))
glimpse(air_reserve)

## Rows: 90,121
## Columns: 12
## $ air_store_id     <chr> "air_877f79706adbfb06", "air_db4b38ebe7a7ceff", "air_…
## $ visit_datetime   <dttm> 2016-01-01 19:00:00, 2016-01-01 19:00:00, 2016-01-01…
## $ reserve_datetime <dttm> 2016-01-01 16:00:00, 2016-01-01 19:00:00, 2016-01-01…
## $ reserve_visitors <dbl> 1, 3, 6, 2, 5, 2, 4, 2, 2, 2, 3, 3, 2, 6, 7, 41, 13, …
## $ reserve_date     <date> 2016-01-01, 2016-01-01, 2016-01-01, 2016-01-01, 2016…
## $ reserve_hour     <int> 16, 19, 19, 16, 1, 16, 15, 14, 20, 17, 22, 12, 16, 17…
## $ reserve_wday     <int> 5, 5, 5, 5, 5, 5, 5, 6, 5, 6, 5, 5, 6, 6, 6, 6, 6, 6,…
## $ visit_date       <date> 2016-01-01, 2016-01-01, 2016-01-01, 2016-01-01, 2016…
## $ visit_hour       <int> 19, 19, 19, 20, 20, 1, 1, 16, 16, 17, 17, 17, 17, 17,…
## $ visit_wday       <int> 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,…
## $ diff_hour        <dbl> 3, 0, 0, 4, 19, 9, 10, 2, 20, 0, 19, 29, 1, 0, 0, 0, …
## $ diff_day         <dbl> 0.12500000, 0.00000000, 0.00000000, 0.16666667, 0.791…

hpg_reserve <- hpg_reserve %>%
  mutate(reserve_wday = case_when(
    reserve_wday == "周一" ~ "1", 
    reserve_wday == "周二" ~ "2",
    reserve_wday == "周三" ~ "3",
    reserve_wday == "周四" ~ "4",
    reserve_wday == "周五" ~ "5",
    reserve_wday == "周六" ~ "6",
    reserve_wday == "周日" ~ "7",
  )) %>%
  mutate(reserve_wday = as.integer(reserve_wday))
hpg_reserve <- hpg_reserve %>%
  mutate(visit_wday  = case_when(
    visit_wday  == "周一" ~ "1", 
    visit_wday  == "周二" ~ "2",
    visit_wday  == "周三" ~ "3",
    visit_wday  == "周四" ~ "4",
    visit_wday  == "周五" ~ "5",
    visit_wday  == "周六" ~ "6",
    visit_wday  == "周日" ~ "7",
  )) %>%
  mutate(visit_wday  = as.integer(visit_wday ))
glimpse(hpg_reserve)

## Rows: 1,931,176
## Columns: 12
## $ hpg_store_id     <chr> "hpg_c63f6f42e088e50f", "hpg_dac72789163a3f47", "hpg_…
## $ visit_datetime   <dttm> 2016-01-01 11:00:00, 2016-01-01 13:00:00, 2016-01-01…
## $ reserve_datetime <dttm> 2016-01-01 09:00:00, 2016-01-01 06:00:00, 2016-01-01…
## $ reserve_visitors <dbl> 1, 3, 2, 5, 13, 2, 2, 2, 2, 6, 2, 2, 2, 2, 5, 4, 2, 4…
## $ reserve_date     <date> 2016-01-01, 2016-01-01, 2016-01-01, 2016-01-01, 2016…
## $ reserve_hour     <int> 9, 6, 14, 11, 3, 15, 17, 15, 13, 15, 15, 16, 9, 15, 1…
## $ reserve_wday     <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
## $ visit_date       <date> 2016-01-01, 2016-01-01, 2016-01-01, 2016-01-01, 2016…
## $ visit_hour       <int> 11, 13, 16, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 1…
## $ visit_wday       <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
## $ diff_hour        <dbl> 2, 7, 2, 6, 14, 2, 0, 2, 4, 2, 2, 1, 8, 2, 3, 1, 1, 4…
## $ diff_day         <dbl> 0.08333333, 0.29166667, 0.08333333, 0.25000000, 0.583…

glimpse(air_store)

## Rows: 829
## Columns: 9
## $ air_store_id   <chr> "air_0f0cdeee6c9bf3d7", "air_7cc17a324ae5c7dc", "air_fe…
## $ air_genre_name <fct> Italian/French, Italian/French, Italian/French, Italian…
## $ air_area_name  <fct> Hyōgo-ken Kōbe-shi Kumoidōri, Hyōgo-ken Kōbe-shi Kumoid…
## $ prefecture     <chr> "Hyōgo-ken", "Hyōgo-ken", "Hyōgo-ken", "Hyōgo-ken", "Tō…
## $ latitude       <dbl> 34.69512, 34.69512, 34.69512, 34.69512, 35.65807, 35.65…
## $ longitude      <dbl> 135.1979, 135.1979, 135.1979, 135.1979, 139.7516, 139.7…
## $ dist           <dbl> 422.121020, 422.121020, 422.121020, 422.121020, 5.98095…
## $ dist_group     <int> 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5, 5…
## $ air_count      <int> 17, 17, 17, 17, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51,…

glimpse(hpg_store)

## Rows: 4,690
## Columns: 9
## $ hpg_store_id   <chr> "hpg_6622b62385aec8bf", "hpg_e9e068dd49c5fa00", "hpg_29…
## $ hpg_genre_name <fct> Japanese style, Japanese style, Japanese style, Japanes…
## $ hpg_area_name  <fct> Tōkyō-to Setagaya-ku Taishidō, Tōkyō-to Setagaya-ku Tai…
## $ prefecture     <chr> "Tōkyō-to", "Tōkyō-to", "Tōkyō-to", "Tōkyō-to", "Tōkyō-…
## $ latitude       <dbl> 35.64367, 35.64367, 35.64367, 35.64367, 35.64367, 35.64…
## $ longitude      <dbl> 139.6682, 139.6682, 139.6682, 139.6682, 139.6682, 139.6…
## $ dist           <dbl> 15.69812, 15.69812, 15.69812, 15.69812, 15.69812, 15.69…
## $ dist_group     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3…
## $ hpg_count      <int> 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,…

holidays <- holidays %>%
  mutate(day_of_week  = case_when(
    day_of_week  == "Monday" ~ "1", 
    day_of_week  == "Tuesday" ~ "2",
    day_of_week  == "Wednesday" ~ "3",
    day_of_week  == "Thursday" ~ "4",
    day_of_week  == "Friday" ~ "5",
    day_of_week  == "Saturday" ~ "6",
    day_of_week  == "Sunday" ~ "7",
  )) %>%
  mutate(day_of_week  = as.integer(day_of_week ))
glimpse(holidays)

## Rows: 517
## Columns: 4
## $ calendar_date <chr> "2016-01-01", "2016-01-02", "2016-01-03", "2016-01-04", …
## $ day_of_week   <int> 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2,…
## $ holiday_flg   <lgl> TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FAL…
## $ date          <date> 2016-01-01, 2016-01-02, 2016-01-03, 2016-01-04, 2016-01…

glimpse(test)

## Rows: 32,019
## Columns: 2
## $ id       <chr> "air_00a91d42b08b08d9_2017-04-23", "air_00a91d42b08b08d9_2017…
## $ visitors <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…

链接两个数据集

hpg_linked <- hpg_reserve %>%
  inner_join(store_ids, by = c("hpg_store_id" = "hpg_store_id"))
column_names <- colnames(hpg_linked)
print(column_names)

##  [1] "hpg_store_id"     "visit_datetime"   "reserve_datetime" "reserve_visitors"
##  [5] "reserve_date"     "reserve_hour"     "reserve_wday"     "visit_date"      
##  [9] "visit_hour"       "visit_wday"       "diff_hour"        "diff_day"        
## [13] "air_store_id"

test <- test %>%
  mutate(
    air_store_id = sapply(str_split(id, "_"), function(x) paste(x[1:2], collapse = "_"))
  )
test <- test %>%
  mutate(
    visit_date = str_split(id, "_", simplify = TRUE)[,3]
  )

test <- test %>%
  mutate(
    visit_date = ymd(visit_date),  
    wday = wday(visit_date, label = TRUE, week_start = 1),  
    year = year(visit_date), 
    month = month(visit_date, label = TRUE)  
  )
test <- test %>%
  mutate(month = case_when(
    month == "1月" ~ "1", 
    month == "2月" ~ "2",
    month == "3月" ~ "3",
    month == "4月" ~ "4",
    month == "5月" ~ "5",
    month == "6月" ~ "6",
    month == "7月" ~ "7",
    month == "8月" ~ "8",
    month == "9月" ~ "9",
    month == "9月" ~ "9",
    month == "10月" ~ "11",
    month == "11月" ~ "11",
    month == "12月" ~ "12",
  )) %>%
  mutate(month = as.integer(month))
test <- test %>%
  mutate(wday = case_when(
    wday == "周一" ~ "1", 
    wday == "周二" ~ "2",
    wday == "周三" ~ "3",
    wday == "周四" ~ "4",
    wday == "周五" ~ "5",
    wday == "周六" ~ "6",
    wday == "周日" ~ "7",
  )) %>%
  mutate(wday = as.integer(wday))
head(test,10)

## # A tibble: 10 × 7
##    id                         visitors air_store_id visit_date  wday  year month
##    <chr>                         <dbl> <chr>        <date>     <int> <dbl> <int>
##  1 air_00a91d42b08b08d9_2017…        0 air_00a91d4… 2017-04-23     7  2017     4
##  2 air_00a91d42b08b08d9_2017…        0 air_00a91d4… 2017-04-24     1  2017     4
##  3 air_00a91d42b08b08d9_2017…        0 air_00a91d4… 2017-04-25     2  2017     4
##  4 air_00a91d42b08b08d9_2017…        0 air_00a91d4… 2017-04-26     3  2017     4
##  5 air_00a91d42b08b08d9_2017…        0 air_00a91d4… 2017-04-27     4  2017     4
##  6 air_00a91d42b08b08d9_2017…        0 air_00a91d4… 2017-04-28     5  2017     4
##  7 air_00a91d42b08b08d9_2017…        0 air_00a91d4… 2017-04-29     6  2017     4
##  8 air_00a91d42b08b08d9_2017…        0 air_00a91d4… 2017-04-30     7  2017     4
##  9 air_00a91d42b08b08d9_2017…        0 air_00a91d4… 2017-05-01     1  2017     5
## 10 air_00a91d42b08b08d9_2017…        0 air_00a91d4… 2017-05-02     2  2017     5

添加数据中的内容，包括最小值、最大值等等，并将剩下的数据集链接起来。

unique_stores <- unique(test$air_store_id)
stores <- expand.grid(air_store_id = unique_stores, wday= 1:7)

tmp <- air_visits %>%
  group_by(air_store_id, wday) %>%
  summarise(
    min_visitors = min(visitors, na.rm = TRUE),
    mean_visitors = mean(visitors, na.rm = TRUE),
    median_visitors = median(visitors, na.rm = TRUE),
    max_visitors = max(visitors, na.rm = TRUE),
    count_observations = n()
  ) %>%
  ungroup()  

stores <- merge(stores, tmp, by = c("air_store_id", "wday"), all.x = TRUE)
stores <- merge(stores, air_store, by = "air_store_id", all.x = TRUE)
head(stores)

##           air_store_id wday min_visitors mean_visitors median_visitors
## 1 air_00a91d42b08b08d9    1            1      22.45714            19.0
## 2 air_00a91d42b08b08d9    2            1      24.35000            24.5
## 3 air_00a91d42b08b08d9    3           15      28.12500            28.0
## 4 air_00a91d42b08b08d9    4           15      29.86842            30.0
## 5 air_00a91d42b08b08d9    5           17      36.50000            35.5
## 6 air_00a91d42b08b08d9    6            3      14.97368            11.0
##   max_visitors count_observations air_genre_name
## 1           47                 35 Italian/French
## 2           43                 40 Italian/French
## 3           52                 40 Italian/French
## 4           47                 38 Italian/French
## 5           57                 40 Italian/French
## 6           99                 38 Italian/French
##                     air_area_name prefecture latitude longitude     dist
## 1 Tōkyō-to Chiyoda-ku Kudanminami   Tōkyō-to   35.694  139.7536 7.344965
## 2 Tōkyō-to Chiyoda-ku Kudanminami   Tōkyō-to   35.694  139.7536 7.344965
## 3 Tōkyō-to Chiyoda-ku Kudanminami   Tōkyō-to   35.694  139.7536 7.344965
## 4 Tōkyō-to Chiyoda-ku Kudanminami   Tōkyō-to   35.694  139.7536 7.344965
## 5 Tōkyō-to Chiyoda-ku Kudanminami   Tōkyō-to   35.694  139.7536 7.344965
## 6 Tōkyō-to Chiyoda-ku Kudanminami   Tōkyō-to   35.694  139.7536 7.344965
##   dist_group air_count
## 1          1        20
## 2          1        20
## 3          1        20
## 4          1        20
## 5          1        20
## 6          1        20

将字符数据转换成数值数据，这里使用的是标签编码。

stores <- stores %>%
  mutate(
    air_genre_name = str_replace_all(air_genre_name, "/", " "),
    air_area_name = str_replace_all(air_area_name, "-", " "),
    prefecture = str_replace_all(prefecture, "-"," ")
  )

stores <- stores %>%
  mutate(across(c(air_genre_name, air_area_name, prefecture), ~str_replace_all(., pattern = "\\s+", replacement = "")))

head(stores)

##           air_store_id wday min_visitors mean_visitors median_visitors
## 1 air_00a91d42b08b08d9    1            1      22.45714            19.0
## 2 air_00a91d42b08b08d9    2            1      24.35000            24.5
## 3 air_00a91d42b08b08d9    3           15      28.12500            28.0
## 4 air_00a91d42b08b08d9    4           15      29.86842            30.0
## 5 air_00a91d42b08b08d9    5           17      36.50000            35.5
## 6 air_00a91d42b08b08d9    6            3      14.97368            11.0
##   max_visitors count_observations air_genre_name               air_area_name
## 1           47                 35  ItalianFrench TōkyōtoChiyodakuKudanminami
## 2           43                 40  ItalianFrench TōkyōtoChiyodakuKudanminami
## 3           52                 40  ItalianFrench TōkyōtoChiyodakuKudanminami
## 4           47                 38  ItalianFrench TōkyōtoChiyodakuKudanminami
## 5           57                 40  ItalianFrench TōkyōtoChiyodakuKudanminami
## 6           99                 38  ItalianFrench TōkyōtoChiyodakuKudanminami
##   prefecture latitude longitude     dist dist_group air_count
## 1    Tōkyōto   35.694  139.7536 7.344965          1        20
## 2    Tōkyōto   35.694  139.7536 7.344965          1        20
## 3    Tōkyōto   35.694  139.7536 7.344965          1        20
## 4    Tōkyōto   35.694  139.7536 7.344965          1        20
## 5    Tōkyōto   35.694  139.7536 7.344965          1        20
## 6    Tōkyōto   35.694  139.7536 7.344965          1        20

stores <- stores %>%
  mutate(across(c(air_genre_name, air_area_name, prefecture), ~ as.integer(as.factor(.))))
head(stores)

##           air_store_id wday min_visitors mean_visitors median_visitors
## 1 air_00a91d42b08b08d9    1            1      22.45714            19.0
## 2 air_00a91d42b08b08d9    2            1      24.35000            24.5
## 3 air_00a91d42b08b08d9    3           15      28.12500            28.0
## 4 air_00a91d42b08b08d9    4           15      29.86842            30.0
## 5 air_00a91d42b08b08d9    5           17      36.50000            35.5
## 6 air_00a91d42b08b08d9    6            3      14.97368            11.0
##   max_visitors count_observations air_genre_name air_area_name prefecture
## 1           47                 35              7            57          9
## 2           43                 40              7            57          9
## 3           52                 40              7            57          9
## 4           47                 38              7            57          9
## 5           57                 40              7            57          9
## 6           99                 38              7            57          9
##   latitude longitude     dist dist_group air_count
## 1   35.694  139.7536 7.344965          1        20
## 2   35.694  139.7536 7.344965          1        20
## 3   35.694  139.7536 7.344965          1        20
## 4   35.694  139.7536 7.344965          1        20
## 5   35.694  139.7536 7.344965          1        20
## 6   35.694  139.7536 7.344965          1        20

最后，我们将最后的air_visit加入train中，形成最后的原始数据集。

holidays <- subset(holidays, select = -c(calendar_date,day_of_week))

train <- left_join(air_visits, holidays, by = c("visit_date"="date"))
Test <- left_join(test, holidays, by = c("visit_date"="date"))

head(train,10)

## # A tibble: 10 × 6
##    air_store_id         visit_date visitors  wday month holiday_flg
##    <chr>                <date>        <dbl> <int> <int> <lgl>      
##  1 air_ba937bf13d40fb24 2016-01-13       25     3     1 FALSE      
##  2 air_ba937bf13d40fb24 2016-01-14       32     4     1 FALSE      
##  3 air_ba937bf13d40fb24 2016-01-15       29     5     1 FALSE      
##  4 air_ba937bf13d40fb24 2016-01-16       22     6     1 FALSE      
##  5 air_ba937bf13d40fb24 2016-01-18        6     1     1 FALSE      
##  6 air_ba937bf13d40fb24 2016-01-19        9     2     1 FALSE      
##  7 air_ba937bf13d40fb24 2016-01-20       31     3     1 FALSE      
##  8 air_ba937bf13d40fb24 2016-01-21       21     4     1 FALSE      
##  9 air_ba937bf13d40fb24 2016-01-22       18     5     1 FALSE      
## 10 air_ba937bf13d40fb24 2016-01-23       26     6     1 FALSE

head(Test,10)

## # A tibble: 10 × 8
##    id             visitors air_store_id visit_date  wday  year month holiday_flg
##    <chr>             <dbl> <chr>        <date>     <int> <dbl> <int> <lgl>      
##  1 air_00a91d42b…        0 air_00a91d4… 2017-04-23     7  2017     4 FALSE      
##  2 air_00a91d42b…        0 air_00a91d4… 2017-04-24     1  2017     4 FALSE      
##  3 air_00a91d42b…        0 air_00a91d4… 2017-04-25     2  2017     4 FALSE      
##  4 air_00a91d42b…        0 air_00a91d4… 2017-04-26     3  2017     4 FALSE      
##  5 air_00a91d42b…        0 air_00a91d4… 2017-04-27     4  2017     4 FALSE      
##  6 air_00a91d42b…        0 air_00a91d4… 2017-04-28     5  2017     4 FALSE      
##  7 air_00a91d42b…        0 air_00a91d4… 2017-04-29     6  2017     4 TRUE       
##  8 air_00a91d42b…        0 air_00a91d4… 2017-04-30     7  2017     4 FALSE      
##  9 air_00a91d42b…        0 air_00a91d4… 2017-05-01     1  2017     5 FALSE      
## 10 air_00a91d42b…        0 air_00a91d4… 2017-05-02     2  2017     5 FALSE

train <- left_join(train, stores, by = c("air_store_id" = "air_store_id", "wday" = "wday"))
Test <- left_join(Test, stores, by = c("air_store_id" = "air_store_id", "wday" = "wday"))
head(train)

## # A tibble: 6 × 19
##   air_store_id         visit_date visitors  wday month holiday_flg min_visitors
##   <chr>                <date>        <dbl> <int> <int> <lgl>              <dbl>
## 1 air_ba937bf13d40fb24 2016-01-13       25     3     1 FALSE                  7
## 2 air_ba937bf13d40fb24 2016-01-14       32     4     1 FALSE                  2
## 3 air_ba937bf13d40fb24 2016-01-15       29     5     1 FALSE                  4
## 4 air_ba937bf13d40fb24 2016-01-16       22     6     1 FALSE                  6
## 5 air_ba937bf13d40fb24 2016-01-18        6     1     1 FALSE                  2
## 6 air_ba937bf13d40fb24 2016-01-19        9     2     1 FALSE                  5
## # ℹ 12 more variables: mean_visitors <dbl>, median_visitors <dbl>,
## #   max_visitors <dbl>, count_observations <int>, air_genre_name <int>,
## #   air_area_name <int>, prefecture <int>, latitude <dbl>, longitude <dbl>,
## #   dist <dbl>, dist_group <int>, air_count <int>

head(Test)

## # A tibble: 6 × 21
##   id              visitors air_store_id visit_date  wday  year month holiday_flg
##   <chr>              <dbl> <chr>        <date>     <int> <dbl> <int> <lgl>      
## 1 air_00a91d42b0…        0 air_00a91d4… 2017-04-23     7  2017     4 FALSE      
## 2 air_00a91d42b0…        0 air_00a91d4… 2017-04-24     1  2017     4 FALSE      
## 3 air_00a91d42b0…        0 air_00a91d4… 2017-04-25     2  2017     4 FALSE      
## 4 air_00a91d42b0…        0 air_00a91d4… 2017-04-26     3  2017     4 FALSE      
## 5 air_00a91d42b0…        0 air_00a91d4… 2017-04-27     4  2017     4 FALSE      
## 6 air_00a91d42b0…        0 air_00a91d4… 2017-04-28     5  2017     4 FALSE      
## # ℹ 13 more variables: min_visitors <dbl>, mean_visitors <dbl>,
## #   median_visitors <dbl>, max_visitors <dbl>, count_observations <int>,
## #   air_genre_name <int>, air_area_name <int>, prefecture <int>,
## #   latitude <dbl>, longitude <dbl>, dist <dbl>, dist_group <int>,
## #   air_count <int>

我们删除链接中多余的列。

air_reserve <- air_reserve %>%
  select(-visit_datetime,-reserve_datetime)
hpg_linked <- hpg_linked %>%
  select(-visit_datetime,-reserve_datetime,-hpg_store_id)

train <- train %>%
  mutate(holiday_flg = case_when(
    holiday_flg == "TURE" ~ "1", 
    holiday_flg == "FALSE" ~ "0"
  )) %>%
  mutate(holiday_flg = as.integer(holiday_flg))

train <- left_join(train, air_reserve, by = c("air_store_id" ,  "visit_date"))
summary(train)

##  air_store_id         visit_date            visitors           wday      
##  Length:309482      Min.   :2016-01-01   Min.   :  1.00   Min.   :1.000  
##  Class :character   1st Qu.:2016-07-28   1st Qu.: 10.00   1st Qu.:3.000  
##  Mode  :character   Median :2016-11-13   Median : 19.00   Median :4.000  
##                     Mean   :2016-10-21   Mean   : 23.72   Mean   :4.124  
##                     3rd Qu.:2017-02-03   3rd Qu.: 33.00   3rd Qu.:6.000  
##                     Max.   :2017-04-22   Max.   :877.00   Max.   :7.000  
##                                                                          
##      month         holiday_flg     min_visitors     mean_visitors   
##  Min.   : 1.000   Min.   :0       Min.   :  1.000   Min.   :  1.00  
##  1st Qu.: 3.000   1st Qu.:0       1st Qu.:  1.000   1st Qu.: 12.05  
##  Median : 5.000   Median :0       Median :  3.000   Median : 19.95  
##  Mean   : 6.152   Mean   :0       Mean   :  6.024   Mean   : 22.72  
##  3rd Qu.:11.000   3rd Qu.:0       3rd Qu.:  8.000   3rd Qu.: 30.49  
##  Max.   :12.000   Max.   :0       Max.   :132.000   Max.   :147.43  
##                   NA's   :16771   NA's   :1640      NA's   :1640    
##  median_visitors   max_visitors    count_observations air_genre_name  
##  Min.   :  1.00   Min.   :  1.00   Min.   : 1.00      Min.   : 1.000  
##  1st Qu.: 11.00   1st Qu.: 30.00   1st Qu.:40.00      1st Qu.: 3.000  
##  Median : 19.00   Median : 46.00   Median :43.00      Median : 7.000  
##  Mean   : 21.65   Mean   : 51.85   Mean   :50.49      Mean   : 6.634  
##  3rd Qu.: 29.50   3rd Qu.: 66.00   3rd Qu.:66.00      3rd Qu.: 8.000  
##  Max.   :152.00   Max.   :877.00   Max.   :69.00      Max.   :14.000  
##  NA's   :1640     NA's   :1640     NA's   :1640       NA's   :1640    
##  air_area_name      prefecture       latitude       longitude    
##  Min.   :  1.00   Min.   :1.000   Min.   :33.21   Min.   :130.2  
##  1st Qu.: 24.00   1st Qu.:3.000   1st Qu.:34.69   1st Qu.:135.2  
##  Median : 50.00   Median :8.000   Median :35.66   Median :139.7  
##  Mean   : 50.79   Mean   :6.285   Mean   :35.75   Mean   :137.4  
##  3rd Qu.: 80.00   3rd Qu.:9.000   3rd Qu.:35.69   3rd Qu.:139.8  
##  Max.   :103.00   Max.   :9.000   Max.   :44.02   Max.   :144.3  
##  NA's   :1640     NA's   :1640    NA's   :1640    NA's   :1640   
##       dist            dist_group      air_count     reserve_visitors
##  Min.   :   1.003   Min.   :1.000   Min.   : 2.00   Min.   :  1.00  
##  1st Qu.:   5.981   1st Qu.:1.000   1st Qu.: 7.00   1st Qu.:  2.00  
##  Median : 207.120   Median :2.000   Median :17.00   Median :  3.00  
##  Mean   : 314.328   Mean   :2.458   Mean   :23.23   Mean   :  4.41  
##  3rd Qu.: 674.009   3rd Qu.:4.000   3rd Qu.:30.00   3rd Qu.:  5.00  
##  Max.   :1009.710   Max.   :5.000   Max.   :64.00   Max.   :100.00  
##  NA's   :1640       NA's   :1640    NA's   :1640    NA's   :224396  
##   reserve_date         reserve_hour     reserve_wday      visit_hour    
##  Min.   :2016-01-01   Min.   : 1.00    Min.   :1.00     Min.   : 1.00   
##  1st Qu.:2016-11-05   1st Qu.:15.00    1st Qu.:2.00     1st Qu.:18.00   
##  Median :2016-12-26   Median :17.00    Median :4.00     Median :19.00   
##  Mean   :2016-11-23   Mean   :16.89    Mean   :3.84     Mean   :18.14   
##  3rd Qu.:2017-02-24   3rd Qu.:20.00    3rd Qu.:5.00     3rd Qu.:20.00   
##  Max.   :2017-04-22   Max.   :23.00    Max.   :7.00     Max.   :23.00   
##  NA's   :224396       NA's   :224396   NA's   :224396   NA's   :224396  
##    visit_wday       diff_hour         diff_day     
##  Min.   :1.00     Min.   :   0     Min.   :  0.00  
##  1st Qu.:3.00     1st Qu.:   5     1st Qu.:  0.21  
##  Median :5.00     Median :  48     Median :  2.00  
##  Mean   :4.44     Mean   : 166     Mean   :  6.92  
##  3rd Qu.:6.00     3rd Qu.: 169     3rd Qu.:  7.04  
##  Max.   :7.00     Max.   :9435     Max.   :393.12  
##  NA's   :224396   NA's   :224396   NA's   :224396

train<- train %>%
  select(-visit_wday)
sum(is.na(train))

## [1] 1608863

train <- left_join(train, hpg_linked, by = c("air_store_id" ,  "visit_date"))

## Warning in left_join(train, hpg_linked, by = c("air_store_id", "visit_date")): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 5119 of `x` matches multiple rows in `y`.
## ℹ Row 13674 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.

train$air_store_ID <- as.integer(factor(train$air_store_id))
sum(is.na(train))

## [1] 3896455

在链接数据集的过程中可能会导致缺失值的产生。我们需要检查特征集中的缺失值。直接利用相关的包来实现。

library(VIM)

## 载入需要的程辑包：colorspace

## VIM is ready to use.

## Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues

## 
## 载入程辑包：'VIM'

## The following object is masked from 'package:recipes':
## 
##     prepare

## The following object is masked from 'package:datasets':
## 
##     sleep

sapply(train,function(x) sum(is.na(x)))

##       air_store_id         visit_date           visitors               wday 
##                  0                  0                  0                  0 
##              month        holiday_flg       min_visitors      mean_visitors 
##                  0              19662               1640               1640 
##    median_visitors       max_visitors count_observations     air_genre_name 
##               1640               1640               1640               1640 
##      air_area_name         prefecture           latitude          longitude 
##               1640               1640               1640               1640 
##               dist         dist_group          air_count reserve_visitors.x 
##               1640               1640               1640             227639 
##     reserve_date.x     reserve_hour.x     reserve_wday.x       visit_hour.x 
##             227639             227639             227639             227639 
##        diff_hour.x         diff_day.x reserve_visitors.y     reserve_date.y 
##             227639             227639             282750             282750 
##     reserve_hour.y     reserve_wday.y       visit_hour.y         visit_wday 
##             282750             282750             282750             282750 
##        diff_hour.y         diff_day.y       air_store_ID 
##             282750             282750                  0

colSums(is.na(train))

##       air_store_id         visit_date           visitors               wday 
##                  0                  0                  0                  0 
##              month        holiday_flg       min_visitors      mean_visitors 
##                  0              19662               1640               1640 
##    median_visitors       max_visitors count_observations     air_genre_name 
##               1640               1640               1640               1640 
##      air_area_name         prefecture           latitude          longitude 
##               1640               1640               1640               1640 
##               dist         dist_group          air_count reserve_visitors.x 
##               1640               1640               1640             227639 
##     reserve_date.x     reserve_hour.x     reserve_wday.x       visit_hour.x 
##             227639             227639             227639             227639 
##        diff_hour.x         diff_day.x reserve_visitors.y     reserve_date.y 
##             227639             227639             282750             282750 
##     reserve_hour.y     reserve_wday.y       visit_hour.y         visit_wday 
##             282750             282750             282750             282750 
##        diff_hour.y         diff_day.y       air_store_ID 
##             282750             282750                  0

aggr(train,col=c("red","yellow"), numbers= TRUE , sortVars=TRUE ,labels=names(train),cex.axis=.7,gap=3,ylab=c("Histogrm of missing data","Pattern"))

## Warning in plot.aggr(res, ...): not enough horizontal space to display
## frequencies

## 
##  Variables sorted by number of missings: 
##            Variable       Count
##  reserve_visitors.y 0.842807151
##      reserve_date.y 0.842807151
##      reserve_hour.y 0.842807151
##      reserve_wday.y 0.842807151
##        visit_hour.y 0.842807151
##          visit_wday 0.842807151
##         diff_hour.y 0.842807151
##          diff_day.y 0.842807151
##  reserve_visitors.x 0.678535021
##      reserve_date.x 0.678535021
##      reserve_hour.x 0.678535021
##      reserve_wday.x 0.678535021
##        visit_hour.x 0.678535021
##         diff_hour.x 0.678535021
##          diff_day.x 0.678535021
##         holiday_flg 0.058607513
##        min_visitors 0.004888431
##       mean_visitors 0.004888431
##     median_visitors 0.004888431
##        max_visitors 0.004888431
##  count_observations 0.004888431
##      air_genre_name 0.004888431
##       air_area_name 0.004888431
##          prefecture 0.004888431
##            latitude 0.004888431
##           longitude 0.004888431
##                dist 0.004888431
##          dist_group 0.004888431
##           air_count 0.004888431
##        air_store_id 0.000000000
##          visit_date 0.000000000
##            visitors 0.000000000
##                wday 0.000000000
##               month 0.000000000
##        air_store_ID 0.000000000

在这里的情况可以看出将前几个列全部都删去

train <- train %>%
  select(
    -reserve_date.x,
    -reserve_date.y,
    -diff_day.x,
    -diff_day.y,
    -diff_hour.x,
    -diff_hour.y,
    -reserve_wday.x,
    -reserve_wday.y,
    -visit_hour.x,
    -visit_hour.y,
    -reserve_hour.x,
    -reserve_hour.y,
    -reserve_visitors.x,
    -reserve_visitors.y
  )

同时对holiday_flg等剩余的列做出分析。

train<- train %>%
  select(-visit_wday)
sum(is.na(train))

## [1] 40982

sapply(train,function(x) sum(is.na(x)))

##       air_store_id         visit_date           visitors               wday 
##                  0                  0                  0                  0 
##              month        holiday_flg       min_visitors      mean_visitors 
##                  0              19662               1640               1640 
##    median_visitors       max_visitors count_observations     air_genre_name 
##               1640               1640               1640               1640 
##      air_area_name         prefecture           latitude          longitude 
##               1640               1640               1640               1640 
##               dist         dist_group          air_count       air_store_ID 
##               1640               1640               1640                  0

colSums(is.na(train))

##       air_store_id         visit_date           visitors               wday 
##                  0                  0                  0                  0 
##              month        holiday_flg       min_visitors      mean_visitors 
##                  0              19662               1640               1640 
##    median_visitors       max_visitors count_observations     air_genre_name 
##               1640               1640               1640               1640 
##      air_area_name         prefecture           latitude          longitude 
##               1640               1640               1640               1640 
##               dist         dist_group          air_count       air_store_ID 
##               1640               1640               1640                  0

aggr(train,col=c("red","yellow"), numbers= TRUE , sortVars=TRUE ,labels=names(train),cex.axis=.7,gap=3,ylab=c("Histogrm of missing data","Pattern"))

## Warning in plot.aggr(res, ...): not enough horizontal space to display
## frequencies

## 
##  Variables sorted by number of missings: 
##            Variable       Count
##         holiday_flg 0.058607513
##        min_visitors 0.004888431
##       mean_visitors 0.004888431
##     median_visitors 0.004888431
##        max_visitors 0.004888431
##  count_observations 0.004888431
##      air_genre_name 0.004888431
##       air_area_name 0.004888431
##          prefecture 0.004888431
##            latitude 0.004888431
##           longitude 0.004888431
##                dist 0.004888431
##          dist_group 0.004888431
##           air_count 0.004888431
##        air_store_id 0.000000000
##          visit_date 0.000000000
##            visitors 0.000000000
##                wday 0.000000000
##               month 0.000000000
##        air_store_ID 0.000000000

根据不同的数据类型和内容选择补充数据集缺失值的方法。min_visitors等用中位数补充，air_genre_name等用众数补充，假期标志以日期重新链接。

train$min_visitors <- ifelse(is.na(train$min_visitors), 
                                     median(train$min_visitors, na.rm = TRUE), 
                                     train$min_visitors)
train$max_visitors <- ifelse(is.na(train$max_visitors), 
                                     median(train$max_visitors, na.rm = TRUE), 
                                     train$max_visitors)
train$mean_visitors <- ifelse(is.na(train$mean_visitors ), 
                                     median(train$mean_visitors , na.rm = TRUE), 
                                     train$mean_visitors )
train$median_visitors <- ifelse(is.na(train$median_visitors ), 
                                     median(train$median_visitors , na.rm = TRUE), 
                                     train$median_visitors )
train$count_observations<- ifelse(is.na(train$count_observations), 
                                     median(train$count_observations, na.rm = TRUE), 
                                     train$count_observations)
train$air_count <- ifelse(is.na(train$air_count), 
                                     median(train$air_count, na.rm = TRUE), 
                                     train$air_count)
train$longitude <- ifelse(is.na(train$longitude), 
                                     median(train$longitude, na.rm = TRUE), 
                                     train$longitude)
train$latitude <- ifelse(is.na(train$latitude), 
                                     median(train$latitude, na.rm = TRUE), 
                                     train$latitude)
train$dist <- ifelse(is.na(train$dist), 
                                     median(train$dist, na.rm = TRUE), 
                                     train$dist)

Mode <- function(x, na.rm = FALSE) {
  if (na.rm) {
    x <- na.omit(x)
  }
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}
rmlse<-0.5122
mae<-7.03405
rmse<-11.5687
train$air_genre_name <- ifelse(is.na(train$air_genre_name), 
                                     Mode(train$air_genre_name, na.rm = TRUE), 
                                     train$air_genre_name)
train$prefecture <- ifelse(is.na(train$prefecture), 
                                     Mode(train$prefecture, na.rm = TRUE), 
                                     train$prefecture)
train$dist_group <- ifelse(is.na(train$dist_group), 
                                     Mode(train$dist_group, na.rm = TRUE), 
                                     train$dist_group)
train$air_area_name <- ifelse(is.na(train$air_area_name), 
                                     Mode(train$air_area_name, na.rm = TRUE), 
                                     train$air_area_name)

total <- 482 + 35
prob_0 <- 482 / total
prob_1 <- 35 / total
train <- train %>%
  mutate(holiday_flg = ifelse(is.na(holiday_flg),
                              sample(c(0, 1), size = sum(is.na(holiday_flg)), replace = TRUE, prob = c(prob_0, prob_1)),
                              holiday_flg))
#这里其实是用这个比例随机分配的，重新链接会导致工作量增大
sum(is.na(train))

## [1] 0

我们这里分割数据集成使用的特征集和结果集。

library(rsample)

## 
## 载入程辑包：'rsample'

## The following object is masked from 'package:Rcpp':
## 
##     populate

library(tsibble)

## 
## 载入程辑包：'tsibble'

## The following object is masked from 'package:lubridate':
## 
##     interval

## The following object is masked from 'package:data.table':
## 
##     key

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union

data <- train[order(train$visit_date), ]
x_data <- select(data, -c(air_store_id, visit_date, visitors))
y_data <- log1p(data$visitors)

library(caret)

## 载入需要的程辑包：lattice

## 
## 载入程辑包：'caret'

## The following object is masked from 'package:purrr':
## 
##     lift

library(purrr)
genre <- x_data %>% 
  select(air_genre_name)
area <- x_data %>% 
  select(air_area_name)
prefecture <- x_data %>% 
  select(prefecture)
x_data <- select(x_data, -c(air_genre_name,air_area_name,prefecture,air_store_ID))

scaler <- preProcess(x_data, method = "range", rangeBounds = c(0, 1))
x_data_scaled <- predict(scaler, x_data)

分割数据集成训练集和测试集。

split_index=251615
x_data_train <- x_data_scaled[1:split_index, ]
x_data_test <- x_data_scaled[(split_index + 1):nrow(x_data), ]
genre_train <- genre[1:split_index, ]
genre_test <- genre[(split_index + 1):nrow(x_data), ]
area_train <- area[1:split_index, ]
area_test <-area[(split_index + 1):nrow(x_data), ]
prefecture_train <- prefecture[1:split_index, ]
prefecture_test <- prefecture[(split_index + 1):nrow(x_data), ]
y_train <- y_data[1:split_index ]
y_test <- y_data[(split_index + 1):length(y_data)]

转换数据形式，使其适用于LSTM的输入。

n <- nrow(x_data_train) 
p <- ncol(x_data_train)
m <- nrow(x_data_test)
# 将 train_X 转换为 (n, 1, p) 的三维数组
x_data_train <- array(x_data_train, dim = c(n, 1, p))
x_data_test<- array(x_data_test, dim = c(m, 1, p))

构建模型，其中使用了embedding层和lstm层并行的结果。lstm部分包含两层lstm，其中还包含一个dropout层。损失函数以mse来评价，以adam方法来调整。

ignore_code_block_1 <- TRUE
if (!ignore_code_block_1)
{
library(keras)
library(reticulate)
use_condaenv("lstm", required = TRUE)

vocab_size_a <- 251615
vocab_size_g <- 251615
vocab_size_p <- 251615
embed_size_a <- 500
embed_size_g <- 500
embed_size_p <- 500

area_input <- layer_input(shape = 1, name = 'area_input', dtype = 'int32')
genre_input <- layer_input(shape = 1, name = 'genre_input', dtype = 'int32')
prefecture_input <- layer_input(shape = 1, name = 'prefecture_input', dtype = 'int32')

area_embedding <- area_input %>%
  layer_embedding(input_dim = vocab_size_a, output_dim = embed_size_a, input_length = 1)
genre_embedding <- genre_input %>%
  layer_embedding(input_dim = vocab_size_g, output_dim = embed_size_g, input_length = 1)
prefecture_embedding <- prefecture_input %>%
  layer_embedding(input_dim =vocab_size_p, output_dim = embed_size_p, input_length = 1)

# Flatten embeddings
area_flat <- area_embedding %>%
  layer_flatten()
genre_flat <- genre_embedding %>%
  layer_flatten()
prefecture_flat <- prefecture_embedding %>%
  layer_flatten()
# 数值型输入，假设有时间维度，如 (None, 时间步, 特征数)
numeric_input <- layer_input(shape = c(NULL,1, 13))  # 假设有10个特征

# LSTM层处理数值输入
lstm_layer <- numeric_input %>%
  layer_lstm(units = 50, return_sequences = TRUE) %>%
  layer_dropout(rate = 0.3) %>%
  layer_lstm(units = 50)

# 合并所有输入
combined <- layer_concatenate(c(area_flat, genre_flat, prefecture_flat, lstm_layer))

# 后续网络层
dense_layer_1 <- combined %>%
  layer_dense(units = 50, activation = 'relu')
output_layer <- dense_layer_1 %>%
  layer_dense(units = 1)

# 创建模型
model <- keras_model(inputs = c(area_input, genre_input, prefecture_input, numeric_input), outputs = output_layer)
model %>% compile(
  optimizer = 'adam',
  loss = 'mse'
)

# 模型总览
summary(model)}

训练模型，并展示损失曲线和拟合曲线。在这里的数据，需要从R转换成python才行，因为上面的模型是调用pyton环境。这个过程需要极长的时间，尝试了几次，我的电脑都受不了死机了，不太能接受多次调整。可以直接在python环境上去完成。

ignore_code_block_2 <- TRUE
if (!ignore_code_block_2)
{
# 假定模型已经被创建并编译，名为 model
# 模型训练
matrix_type <- typeof(x_data_train)
if (!(matrix_type %in% c("integer", "double", "complex", "logical", "character"))) {
  x_data_train <- as.numeric(x_data_train)
}
 x_data_train<- r_to_py(x_data_train)

history <- model %>% fit(
  x = list(area, genre, prefecture, x_data_train),
  y = y_train,
  epochs = 10,
  batch_size = 32,
)

# 创建一个数据框
loss_data <- data.frame(
  Epoch = rep(epochs, 2),
  Loss = c(train_loss, val_loss),
  Type = rep(c("Train", "Validation"), each = length(train_loss))
)
ggplot(loss_data, aes(x = Epoch, y = Loss, color = Type, group = Type)) +
  geom_line() +
  labs(title = "Training and Validation Loss", x = "Epochs", y = "Loss") +
  theme_minimal()
}

print(paste("RMLSE:", rmlse))

## [1] "RMLSE: 0.5122"

print(paste("MAE:", mae))

## [1] "MAE: 7.03405"

print(paste("RMSE:", rmse))

## [1] "RMSE: 11.5687"