R 机器学习预测时间序列模型

Jamesjin63

发布于 2022-10-25 06:53:24

99100

代码可运行

文章被收录于专栏：EpiHubEpiHub

运行总次数：0

代码可运行

机器学习在时间序列数据上应用

随着疫情的变化，急性传染病数据经常会随时间变化，我们通过对每天传染病的记录，就形成了时间序列数据，周期可以是天，周，月，年。目前我们经常会用到ARIMA来预测疾病在未来的变化趋势。

image.png

但是随着机器学习的广泛应用，在时间序列上，也可以采用机器学习发方法去预测，结果比传统的ARIMA EST更加快速，简洁，准确。这次将要介绍关于的时间序列预测的Modeltime包，旨在加快模型评估，选择和预测的速度。modeltime通过将tidymodels机器学习软件包生态系统集成到简化的工作流中以进行tidyverse预测来实现此目的。modeltime结合了机器学习模型，经典模型和自动化模型等。

image.png

主要优点：简化数据建模预测流程，包括数据建模，评估，预测及输出

预测的系统工作流程。modeltime_table（），modeltime_calibrate（）和modeltime_refit（）
结合Tidymodels以期加入机器学习算法。如XGBoost，GLMnet，Stan，Random Forest等
改进传统时间序列模型。如arima_boost（），prophet_boost（）

1.数据

我们选取bike_sharing_daily时间序列数据集，其中包括自行车每日的使用数据。这里只需要日期与当日的使用量“date” and “value”。然后可以简单绘制一下。注意这里的时间序列是tibble格式。

install.packages("modeltime")
devtools::install_github("business-science/timetk")
library(tidyverse)
library(tidymodels)
library(modeltime)
library(timetk)   
library(lubridate)

# data
bike_transactions_tbl <- bike_sharing_daily %>%
  select(dteday, cnt) %>%
  set_names(c("date", "value")) 

bike_transactions_tbl

# plot
bike_transactions_tbl %>%
  plot_time_series(date, value, .interactive = FALSE)

2.模型构建

下面继续对数据进行处理，在建模之前，我们需要将数据分成train与test。使用time_series_split()来分割我们的数据，assess = "3 months"来确定后三个月为test数据集，cumulative = TRUE指定前面部分为train。

# Split your time series into training and testing sets
splits <- bike_transactions_tbl %>%
  time_series_split(assess = "3 months", cumulative = TRUE)
  
splits
 <Analysis/Assess/Total>
<641/90/731>
  
# plot
splits %>%
  tk_time_series_cv_plan() %>%
  plot_time_series_cv_plan(date, value, .interactive = FALSE)

image.png

image.png

2.1经典ARIMA与Prophet

因为model time可以简化数据建模流程，包括简化参数的设置，自动模型通常是已经包含了自动化的建模方法。包括“自动ARIMA”和“自动ETS”功能以及“ Prophet”算法。主要包含三个参数设置：

Model Spec: 指定预测模型种类(e.g. arima_reg(), prophet_reg())
Engine: 指定模型
Fit Model: 加载trian数据

结下来我们尝试建立ARIMA模型与Prophet模型

# ARIMA
model_fit_arima <- arima_reg() %>%
  set_engine("auto_arima") %>%
  fit(value ~ date, training(splits))

model_fit_arima 
parsnip model object

Fit time:  426ms 
Series: outcome 
ARIMA(0,1,3) with drift 

Coefficients:
          ma1      ma2      ma3   drift
      -0.6106  -0.1868  -0.0673  9.3169
s.e.   0.0396   0.0466   0.0398  4.6225

sigma^2 estimated as 730568:  log likelihood=-5227.22
AIC=10464.44   AICc=10464.53   BIC=10486.74

# Prophet
model_fit_prophet <- prophet_reg() %>%
  set_engine("prophet", yearly.seasonality = TRUE) %>%
  fit(value ~ date, training(splits))

model_fit_prophet
parsnip model object

Fit time:  455ms 
PROPHET Model
- growth: 'linear'
- n.changepoints: 25
- seasonality.mode: 'additive'
- extra_regressors: 0

2.2机器学习算法

机器学习模型前面设置比自动化经典模型更为复杂。通常在进行机器学习建模之前，对数据进行预处理，称之为workflow一般过程如下：

创建预处理配方 Preprocessing Recipe
创建模型规格 Model Specifications
使用工作流将模型规格和预处理相结合，并拟合模型 Fit Model

首先，我将使用配recipe()创建预处理数据的先前步骤。该过程使用“日期”列创建了我要建模的45个新的列。这些列包含了时间序列的详细信息及傅立叶变化的数据。

recipe_spec <- recipe(value ~ date, training(splits)) %>%
  step_timeseries_signature(date) %>%
  step_rm(contains("am.pm"), contains("hour"), contains("minute"),
          contains("second"), contains("xts")) %>%
  step_fourier(date, period = 365, K = 5) %>%
  step_dummy(all_nominal())

recipe_spec %>% prep() %>% juice()

有了recipe，我们可以建立机器学习模型了。为什么需要recipe是因为在tidymodel里面，设置了建立机器学习模型的一套准则，感兴趣可以去：

机器学习模型

这里我们新建了glmnet与RF模型。主要参数包括：

Start with a workflow()
Add a Model Spec: add_model(model_spec_glmnet)
Add Preprocessing: add_recipe(recipe_spec %>% step_rm(date))
Fit the Workflow: fit(training(splits))

# Elastic NET model
model_spec_glmnet <- linear_reg(penalty = 0.01, mixture = 0.5) %>%
  set_engine("glmnet")

workflow_fit_glmnet <- workflow() %>%
  add_model(model_spec_glmnet) %>%
  add_recipe(recipe_spec %>% step_rm(date)) %>%
  fit(training(splits))

# Random Forest
model_spec_rf <- rand_forest(trees = 500, min_n = 50) %>%
  set_engine("randomForest")

workflow_fit_rf <- workflow() %>%
  add_model(model_spec_rf) %>%
  add_recipe(recipe_spec %>% step_rm(date)) %>%
  fit(training(splits))

2.3机器学习结合传统算法

Prophet Boost算法将Prophet与XGBoost结合使用，从而获得了两全其美的效果（即Prophet Automation + Machine Learning）

首先使用prophet对单个时间序列建模
使用通过预处理配方提供的回归数据（生成的45个新列），并使用XGBoost模型对prophet残差进行回归

model_spec_prophet_boost <- prophet_boost() %>%
  set_engine("prophet_xgboost", yearly.seasonality = TRUE) 

workflow_fit_prophet_boost <- workflow() %>%
  add_model(model_spec_prophet_boost) %>%
  add_recipe(recipe_spec) %>%
  fit(training(splits))
  
## [05:22:38] WARNING: amalgamation/../src/learner.cc:480: 
## Parameters: { validation } might not be used.
## 
##   This may not be accurate due to some parameters are only used in language bindings but
##   passed down to XGBoost core.  Or some parameters are not used but slip through this
##   verification. Please open an issue if you find above cases.

workflow_fit_prophet_boost

3.模型调整

image.png

Modeltime工作流程旨在加速模型评估和选择。现在我们有了几个时间序列模型，让我们对其进行分析，并通过模型时间工作流程预测未来变化趋势。 Modeltime使用ID来定位我们之前建立的模型，以帮助我们识别模型。让我们将模型添加到modeltime_table（）中。 [图片上传中...(image.png-ccae46-1593703307979-0)]

model_table <- modeltime_table(
  model_fit_arima, 
  model_fit_prophet,
  workflow_fit_glmnet,
  workflow_fit_rf,
  workflow_fit_prophet_boost
) 

model_table
## # Modeltime Table
## # A tibble: 5 x 3
##   .model_id .model     .model_desc              
##       <int> <list>     <chr>                    
## 1         1 <fit[+]>   ARIMA(0,1,3) WITH DRIFT  
## 2         2 <fit[+]>   PROPHET                  
## 3         3 <workflow> GLMNET                   
## 4         4 <workflow> RANDOMFOREST             
## 5         5 <workflow> PROPHET W/ XGBOOST ERRORS

模型校准用于量化误差并估计置信区间。我们将使用modeltime_calibrate（）函数对test数据进行模型校准。将生成两个新列（“ .type”和“ .calibration_data”），其中最重要的是“ .calibration_data”。包括l测试集的实际值，拟合值和残差。

calibration_table <- model_table %>%
  modeltime_calibrate(testing(splits))

calibration_table
## # Modeltime Table
## # A tibble: 5 x 5
##   .model_id .model     .model_desc               .type .calibration_data
##       <int> <list>     <chr>                     <chr> <list>           
## 1         1 <fit[+]>   ARIMA(0,1,3) WITH DRIFT   Test  <tibble [90 × 4]>
## 2         2 <fit[+]>   PROPHET                   Test  <tibble [90 × 4]>
## 3         3 <workflow> GLMNET                    Test  <tibble [90 × 4]>
## 4         4 <workflow> RANDOMFOREST              Test  <tibble [90 × 4]>
## 5         5 <workflow> PROPHET W/ XGBOOST ERRORS Test  <tibble [90 × 4]>

4.模型预测

利用校准后的数据，我们可以对test数据进行预测。使用modeltime_forecast（）生成测试集的预测数据。使用plot_modeltime_forecast（）绘制预测结果。

calibration_table %>%
  modeltime_forecast(actual_data = bike_transactions_tbl) %>%
  plot_modeltime_forecast(.interactive = FALSE)

image.png

5.模型评估

接下来，计算预测值与实际值的差异来评估模型好坏。使用modeltime_accuracy（）生成样本外准确性指标数据。使用table_modeltime_accuracy（）生成预测结果表格

calibration_table %>%
  modeltime_accuracy() %>%
  table_modeltime_accuracy(.interactive = FALSE)

image.png

后续还有对最优模型进行预测，使用全部数据作为training，预测未来1年的变化趋势。

参考

Modeltime GitHub Page - Give it a Star if you like it! Timetk Documentation - Data wrangling, visualization, and preprocessing for time series. Tidymodels.org - The tidymodels framework is a collection of packages for modeling and machine learning using tidyverse principles.

本文参与腾讯云自媒体同步曝光计划，分享自作者个人站点/博客。

原始发表：2020-07-02，如有侵权请联系 cloudcommunity@tencent.com 删除

神经网络