量化Covid19对SFO机场客运航空交通的影响

专题介绍：R是一种广泛用于数据分析和统计计算的强大语言，于上世纪90年代开始发展起来。得益于全世界众多爱好者的无尽努力，大家继而开发出了一种基于R但优于R基本文本编辑器的R Studio（用户的界面体验更好）。也正是由于全世界越来越多的数据科学社区和用户对R包的慷慨贡献，让R语言在全球范围内越来越流行。其中一些R包，例如MASS，SparkR， ggplot2，使数据操作,可视化和计算功能越来越强大。R是用于统计分析、绘图的语言和操作环境。R是属于GNU系统的一个自由、免费、源代码开放的软件，它是一个用于统计计算和统计制图的优秀工具。R作为一种统计分析软件，是集统计分析与图形显示于一体的。它可以运行于UNIX、Windows和Macintosh的操作系统上，而且嵌入了一个非常方便实用的帮助系统，相比于其他统计分析软件，R的学术性开发比较早，适合生物学和医学等学术学科的科研人员使用。

是新朋友吗？记得先点R语言关注我哦～

《R文章翻译》专栏·第2篇

文 | R扫地僧

5155字 | 15分钟阅读

【R语言】已开通R语言社群，五湖四海，天南地北，各行各业，有缘相聚，共享R事，雕刻数据，求解问题，以创价值。喜乐入群者，请加微信号luqin360，或扫描文末二维码，添加为好友，同时附上R-入群。有朋自远方来，不亦乐乎，并诚邀入群，以达相互学习和进步之美好心愿。

R语言出品

题目：Quantify the Covid19 Impact on the SFO Airport Passenger Air Traffic

来源：

https://ramikrispin.github.io/2021/01/covid19-effect/

编译：R扫地僧

不幸的是，Covid19大流行仍对大多数主要行业产生了重大影响。虽然对于某些行业来说，影响是积极的（例如在线零售，互联网和蒸汽供应商等），但对于其他行业则是不利的（例如交通，旅游，娱乐等）。在这两种情况下，我们都可以利用时间序列建模来量化Covid19对行业的影响。

量化Covid19大流行影响（无论是正面还是负面）的一种简单方法将包括以下步骤：

将序列分为Covid前和Covid后
使用Covid前的数据训练时间序列模型。假设没有大流行，这将使我们能够模拟序列的价值
使用训练的模型预测Covid后序列的值
使用预测值与实际值之间的差值（Covid后序列）来量化Covid19对序列的影响

为了演示这种方法，我将使用来自sfo包的sfo_passengers数据集。sfo_passengers数据集提供了2005年7月至2020年9月旧金山国际机场(SFO)空中交通的月度统计数据。关于数据集的更多细节可以在这个帖子和文档中找到。

为了进行分析，我们将使用以下软件包：

sfo-旅客航空交通数据集
dplyr-数据准备
plotly-数据可视化
TSstudio-时间序列分析和预测

library(sfo)
library(dplyr)
library(plotly)
library(TSstudio)

数据
如上所述，sfo_passengers数据集提供2005年以来SFO机场的航空客运量每月统计数据。包括按不同类别(如营运航空公司、地区、终点站等)统计的每月旅客人数。让我们加载数据:

data("sfo_passengers")
str(sfo_passengers)

在将数据转换为时间序列格式之前，我们将activity_period转换为日期格式:

df <- sfo_passengers %>%
  mutate(date = as.Date(paste(substr(sfo_passengers$activity_period, 1,4), 
                              substr(sfo_passengers$activity_period, 5,6), 
                              "01", sep ="/")))

接下来，我们将数据集转换为时间序列格式，将乘客按日期变量分组:

df <- df %>%
  group_by(date) %>%
  summarise(y = sum(passenger_count), .groups = "drop")

head(df)

现在，我们有一个按月的时间序列:

plot_ly(data = df,
        x = ~ date,
        y = ~ y,
        type = "scatter", 
        mode = "line",
        name = "Total Passengers") %>%
  add_segments(x = as.Date("2020-02-01"), 
               xend = as.Date("2020-02-01"),
               y = min(df$y),
               yend = max(df$y) * 1.05,
               line = list(color = "black", dash = "dash"),
               showlegend = FALSE) %>%
  add_annotations(text = "Pre-Covid19",
                  x = as.Date("2018-09-01"),
                  y = max(df$y) * 1.05, 
                  showarrow = FALSE) %>%
  add_annotations(text = "Post-Covid19",
                  x = as.Date("2021-08-01"),
                  y = max(df$y) * 1.05, 
                  showarrow = FALSE) %>%
  layout(title = "Total Number of Air Passengers - SFO Airport",
         yaxis = list(title = "Number of Passengers"),
         xaxis = list(title = "Source: San Francisco Open Data Portal"))

从上图可以看出，2020年3月以来，Covid19的影响十分明显。为量化疫情对旅客人数的影响，我们将该序列分为以下两个序列:

Covid19之前序列-所有的观察结果都是在2020年3月之前
Covid19之后序列-所有观察结果(包括)2020年3月之后

我们将使用Covid19之前序列来训练一个时间序列模型。这将使我们能够预测乘客人数，就像没有大流行一样。

pre_covid <- df %>% 
  dplyr::filter(date < as.Date("2020-03-01")) %>%
  dplyr::arrange(date)

post_covid <- df %>% 
  dplyr::filter(date >= as.Date("2020-03-01")) %>%
  dplyr::arrange(date)

我们将使用Covid前序列训练一个时间序列模型。这将使我们能够预测因没有这种Covid大流行下乘客人数。一旦我们用Covid19前数据预测Covid19后序列的相应观测结果，我们就可以量化Covid19对总乘客人数的影响。

分析数据

在我们预测该序列之前，让我们对该序列进行快速探索性分析，以确定其主要特征。我们将使用TSstudio包可视化pre_covid序列。

注意，这个包还不支持tsibble对象。因此，我们首先将series转换为ts对象:

ts.obj <- ts(pre_covid$y, start = c(2005, 7), frequency = 12)

大流行爆发前的序列：

ts_plot(ts.obj,
        title = "Total Number of Air Passengers - SFO Airport",
        Ytitle = "Number of Passengers",
        slider = TRUE)

像大多数描述月度航空客流量的序列一样，该序列具有很强的月度季节性模式和正向趋势。您还可以注意到，季节性成分的振荡自2017年以来(与前几年相比)变得更大。让我们使用ts_seasonal函数来创建该序列的季节图:

ts_seasonal(ts.obj = ts.obj, type = "all")

从季节图中可以看出，每个月的季节效应都在不断增加，而序列却在逐年增长。

类似地，我们可以使用ACF和PACF函数来回顾与其过去滞后的序列相关性:

ts_cor(ts.obj = ts.obj, lag.max = 36)

正如预期的那样，我们可以看到该序列与第一和季节性滞后之间的强相关性。我们将利用这些信息为季节性数据选择时间序列模型。

预测Covid19前序列

我最喜欢的预测策略之一是将不同时间序列模型之间的赛马和回溯测试结合起来作为一种训练方法。回溯测试是时间序列等价的机器学习交叉验证训练方法。这里的想法很简单—用回测测试每个模型，然后在不同的测试分区上选择平均表现最好的一个。

TSstudio包中的train_model函数使我们能够使用forecast和stats包中的模型无缝地应用这个策略。为了简单起见，我们将使用不同的ETS、Holt-Winters模型和开箱即用的auto.arima和tslm模型。为了回测，我们将把序列分成6个测试分区，每12个月间隔3个月。

methods参数定义要使用的模型，train_method参数定义回测的设置。可以在这里找到关于这个函数的更多细节。

methods <- list(ets1 = list(method = "ets",
                            method_arg = list(opt.crit = "lik"),
                            notes = "ETS model opt.crit=lik"),
                ets2 = list(method = "ets",
                            method_arg = list(opt.crit = "amse"),
                            notes = "ETS model opt.crit=amse"),
                ets3 = list(method = "ets",
                            method_arg = list(opt.crit = "mse"),
                            notes = "ETS model opt.crit=mse"),
                auto_arima = list(method = "auto.arima",
                                  notes = "Auto ARIMA"),
                hw1 = list(method = "HoltWinters",
                           method_arg = NULL,
                           notes = "HoltWinters Model"),
                hw2 = list(method = "HoltWinters",
                           method_arg = list(seasonal = "multiplicative"),
                           notes = "HW with multip. seasonality"),
                tslm = list(method = "tslm",
                            method_arg = list(formula = input ~ trend + season),
                            notes = "tslm with trend and seasonal"))


train_method = list(partitions = 6,
                    sample.out = 12,
                    space = 3)

在定义了方法和train_method参数之后，我们将使用train_model函数来训练模型。请注意，预测地平线被设置为post_covid序列的长度。此外，我们将MAPA设置为误差度量，以评估不同模型在测试分区上的性能:

md <- train_model(input = ts.obj,
                  methods = methods,
                  train_method = train_method,
                  horizon = nrow(post_covid),
                  error = "MAPE")

根据来自train_model函数的排行榜表，在不同测试分区上平均表现最佳的模型是Holt-Winters模型（第一个版本-hw1）。与评估的其他模型相比，该模型平均获得最低的MAPE（2.76％）和RMSE（149046）。此外，该模型对模型预测区间的覆盖率接近完美，对于80％和95％的预测区间，平均覆盖率分别为79.2％和95.8％。我们可以使用plot_error函数查看每个模型在不同分区上的误差分布：

plot_error(md)

plot_model使我们能够在回测的不同测试分区上动画化每个模型的预测值:

plot_model(md)

我们将选择霍尔特-温特斯模型(hw1)来计算covid – 19效应。我们将把选定的预测添加到post_covid数据集:

post_covid$yhat <- as.numeric(md$forecast$hw1$forecast$mean)
post_covid$upper95 <- as.numeric(md$forecast$hw1$forecast$upper[,2])
post_covid$lower95 <- as.numeric(md$forecast$hw1$forecast$lower[,2])

量化Covid19的影响
在我们添加预测值后，很容易计算出Covid19对旧金山国际机场乘客人数的月度影响:

post_covid$passengers_loss <- post_covid$y - post_covid$yhat
post_covid

2020年3月至9月期间，由于大流行，预计减少的乘客总数:

sum(post_covid$passengers_loss)

同样，我们可以可视化新冠病毒对航空客运的影响:

plot_ly() %>%
  add_ribbons(x = post_covid$date,
              ymin = post_covid$y,
              ymax = post_covid$yhat,
              line = list(color = 'rgba(255, 0, 0, 0.05)'),
              fillcolor = 'rgba(255, 0, 0, 0.6)',
              name = "Estimated Loss") %>%
  add_segments(x = as.Date("2020-02-01"), 
               xend = as.Date("2020-02-01"),
               y = min(df$y),
               yend = max(df$y) * 1.05,
               line = list(color = "black", dash = "dash"),
               showlegend = FALSE) %>%
  add_annotations(text = "Pre-Covid19",
                  x = as.Date("2017-09-01"),
                  y = max(df$y) * 1.05, 
                  showarrow = FALSE) %>%
  add_annotations(text = "Post-Covid19",
                  x = as.Date("2020-09-01"),
                  y = max(df$y) * 1.05, 
                  showarrow = FALSE) %>%
  add_annotations(text = paste("Estimated decrease in", "<br>",
                               "passengers volume: ~30M",
                               sep = ""),
                  x = as.Date("2020-05-01"),
                  y = 2 * 10 ^ 6, 
                  arrowhead = 1,
                  ax = -130,
                  ay = -40,
                  showarrow = TRUE) %>%
  add_lines(x = df$date,
            y = df$y,
            line = list(color = "#1f77b4"),
            name = "Actual") %>%
  layout(title = "Covid19 Impact on SFO Air Passenger Traffic",
         yaxis = list(title = "Number of Passengers"),
         xaxis = list(title = "Time Series Model - Holt-Winters",
                      range = c(as.Date("2015-01-01"), as.Date("2021-01-01"))),
         legend = list(x = 0, y = 0.95))

应用领域

一旦我们估计了乘客人数的减少，就可以量化Covid19造成的损失。例如，如果每位旅客平均要支付10美元的机场税，那么在特定时期内，估计的税收损失约为3亿美元。

另外，由于underline预测是点估计。因此，您可以利用预测间隔为航空客运量的下降提供一个范围。

最后但并非最不重要的一点是，您可以使用自上而下的方法并为数据中的某些可用类别进行预测。例如，基于航空公司、地区、国内/国际航班的乘客下降。

参考的信息

更多关于这篇文章中使用的包和工具的细节

数据集

sfo包 – https://github.com/RamiKrispin/sfo
San Francisco Data Portal – https://datasf.org/opendata/
dplyr包 – https://dplyr.tidyverse.org/

时间序列分析

TSstudio包 – https://github.com/RamiKrispin/TSstudio
train_model function documentation – https://ramikrispin.github.io/TSstudio/reference/train_model.html
Great explanation about backtesting from Uber Engineering blog – https://eng.uber.com/omphalos/

数据可视化

plotly包 – https://plotly.com/r/

附录：完整代码

library(sfo)
library(dplyr)
library(plotly)
library(TSstudio)


# 数据集
data("sfo_passengers")
str(sfo_passengers)

df <- sfo_passengers %>%
  mutate(date = as.Date(paste(substr(sfo_passengers$activity_period, 1,4), 
                              substr(sfo_passengers$activity_period, 5,6), 
                              "01", sep ="/"))) 

df <- df %>%
  group_by(date) %>%
  summarise(y = sum(passenger_count), .groups = "drop")

head(df) 

plot_ly(data = df,
        x = ~ date,
        y = ~ y,
        type = "scatter", 
        mode = "line",
        name = "Total Passengers") %>%
  add_segments(x = as.Date("2020-02-01"), 
               xend = as.Date("2020-02-01"),
               y = min(df$y),
               yend = max(df$y) * 1.05,
               line = list(color = "black", dash = "dash"),
               showlegend = FALSE) %>%
  add_annotations(text = "Pre-Covid19",
                  x = as.Date("2018-09-01"),
                  y = max(df$y) * 1.05, 
                  showarrow = FALSE) %>%
  add_annotations(text = "Post-Covid19",
                  x = as.Date("2021-08-01"),
                  y = max(df$y) * 1.05, 
                  showarrow = FALSE) %>%
  layout(title = "Total Number of Air Passengers - SFO Airport",
         yaxis = list(title = "Number of Passengers"),
         xaxis = list(title = "Source: San Francisco Open Data Portal"))

# 数据集划分
pre_covid <- df %>% 
  dplyr::filter(date < as.Date("2020-03-01")) %>%
  dplyr::arrange(date)

post_covid <- df %>% 
  dplyr::filter(date >= as.Date("2020-03-01")) %>%
  dplyr::arrange(date)

ts.obj <- ts(pre_covid$y, start = c(2005, 7), frequency = 12)

ts_plot(ts.obj,
        title = "Total Number of Air Passengers - SFO Airport",
        Ytitle = "Number of Passengers",
        slider = TRUE)

ts_seasonal(ts.obj = ts.obj, type = "all")

ts_cor(ts.obj = ts.obj, lag.max = 36)


methods <- list(ets1 = list(method = "ets",
                            method_arg = list(opt.crit = "lik"),
                            notes = "ETS model opt.crit=lik"),
                ets2 = list(method = "ets",
                            method_arg = list(opt.crit = "amse"),
                            notes = "ETS model opt.crit=amse"),
                ets3 = list(method = "ets",
                            method_arg = list(opt.crit = "mse"),
                            notes = "ETS model opt.crit=mse"),
                auto_arima = list(method = "auto.arima",
                                  notes = "Auto ARIMA"),
                hw1 = list(method = "HoltWinters",
                           method_arg = NULL,
                           notes = "HoltWinters Model"),
                hw2 = list(method = "HoltWinters",
                           method_arg = list(seasonal = "multiplicative"),
                           notes = "HW with multip. seasonality"),
                tslm = list(method = "tslm",
                            method_arg = list(formula = input ~ trend + season),
                            notes = "tslm with trend and seasonal"))


train_method = list(partitions = 6,
                    sample.out = 12,
                    space = 3)

md <- train_model(input = ts.obj,
                  methods = methods,
                  train_method = train_method,
                  horizon = nrow(post_covid),
                  error = "MAPE")

plot_error(md)

plot_model(md)

post_covid$yhat <- as.numeric(md$forecast$hw1$forecast$mean)
post_covid$upper95 <- as.numeric(md$forecast$hw1$forecast$upper[,2])
post_covid$lower95 <- as.numeric(md$forecast$hw1$forecast$lower[,2])

post_covid$passengers_loss <- post_covid$y - post_covid$yhat
post_covid

sum(post_covid$passengers_loss)

plot_ly() %>%
  add_ribbons(x = post_covid$date,
              ymin = post_covid$y,
              ymax = post_covid$yhat,
              line = list(color = 'rgba(255, 0, 0, 0.05)'),
              fillcolor = 'rgba(255, 0, 0, 0.6)',
              name = "Estimated Loss") %>%
  add_segments(x = as.Date("2020-02-01"), 
               xend = as.Date("2020-02-01"),
               y = min(df$y),
               yend = max(df$y) * 1.05,
               line = list(color = "black", dash = "dash"),
               showlegend = FALSE) %>%
  add_annotations(text = "Pre-Covid19",
                  x = as.Date("2017-09-01"),
                  y = max(df$y) * 1.05, 
                  showarrow = FALSE) %>%
  add_annotations(text = "Post-Covid19",
                  x = as.Date("2020-09-01"),
                  y = max(df$y) * 1.05, 
                  showarrow = FALSE) %>%
  add_annotations(text = paste("Estimated decrease in", "<br>",
                               "passengers volume: ~30M",
                               sep = ""),
                  x = as.Date("2020-05-01"),
                  y = 2 * 10 ^ 6, 
                  arrowhead = 1,
                  ax = -130,
                  ay = -40,
                  showarrow = TRUE) %>%
  add_lines(x = df$date,
            y = df$y,
            line = list(color = "#1f77b4"),
            name = "Actual") %>%
  layout(title = "Covid19 Impact on SFO Air Passenger Traffic",
         yaxis = list(title = "Number of Passengers"),
         xaxis = list(title = "Time Series Model - Holt-Winters",
                      range = c(as.Date("2015-01-01"), as.Date("2021-01-01"))),
         legend = list(x = 0, y = 0.95))