在 R 中使用 Anomalize 检测异常

Anomalize：简要介绍

谈到异常时，我们指的是离群点或特殊事件。在小型数据集中识别这些事件很容易，可以用箱线图等简单的分析图完成。但当切换到大型数据集时，尤其是时间序列场景，情况会变得复杂。时间序列是在一段时间内按固定时间间隔采集的数据，分析后通常会呈现趋势或季节性。在这些情况下识别异常有一定难度。

这时就要用到用于时间序列分析中的异常检测包 anomalize。它是一种以时间为基础、可从单一时间序列扩展到多时间序列的整洁（tidy）异常检测算法。

其实现依托了一些现有的包和方法，可以说是将已有资源与可扩展的方法相结合。

提供了帮助的开源工作包括：

Twitter 的 AnomalyDetection 包：可在 Github 获取（CRAN 上的 anaomalyDetection 是不同的作品）。
Rob Hyndman 的 forecast 包中的 forecast::tsoutliers() 函数。
Javier Lopez-de-lacalle 在 CRAN 上的 tsoutliers 包。

以上这些包和函数被用于集成到一个可扩展的工作流中。

关于 anomalize 的工作流，它分为三个部分：

使用 time_decompose() 进行时间序列分解。
使用 anomalize() 对剩余项进行异常检测。
使用 time_recompose() 对异常的上下界进行重构转换。

异常检测工作流

时间序列分解

第一步是使用 time_decompose() 进行时间序列分解。需要对某个分组执行检测的测量值或数值会被分解为四列：observed（观测值）、season（季节项）、trend（趋势项）和 remainder（剩余项）。默认使用的分解方法是 stl，即利用 Loess 平滑器的季节性分解。

Loess 回归是平滑高波动时间序列最常用的方法之一，它在局部邻域内拟合多次回归。也可以理解为将数据分段，并对每段分别进行回归。这在时间序列中很有用，因为我们知道时间（在此为自变量 X）的取值范围。当时间序列中趋势主导于季节性时，该方法表现良好。

这里的趋势是指跨越多次观测的长期增长；而季节性是指按分钟、小时或每周等周期发生的循环模式。

还有第二种可用于时间序列季节性分解的技术，即基于中位数的 Twitter 方法，AnomalyDetection 包也使用了该方法。它在去除季节性部分上与 STL 类似。不同之处在于趋势的去除：它使用数据的分段中位数（在指定间隔处进行一次或多次中位数划分），而非拟合平滑器。当时间序列中季节性主导于趋势时，该方法表现更好。

接下来看看 time_decompose() 函数的输出，如上所述它会生成 4 列：

observed：实际值。
season：季节性或循环趋势。默认是每周季节性。
trend：长期趋势。默认跨度为 3 个月。
remainder：用于分析离群点。等于 observed 减去 season 和 trend。

time_decompose() 函数包含一个参数 merge，将其设为 TRUE 可在保留原始数据的同时保留生成的各列。

# Using data package provided in the anomalize package and taking single time series of package purrr

purrr_package = tidyverse_cran_downloads%>%
  filter(package == "purrr")%>%
  ungroup()

purrr_anomaly  = purrr_package %>%
  time_decompose(count)

purrr_anomaly%>% glimpse()
技术支持

 ## Observations: 425
## Variables: 5
## $ date      <date> 2017-01-01, 2017-01-02, 2017-01-03, 2017-01-04, 201...
## $ observed  <dbl> 550, 1012, 1515, 1702, 1696, 1613, 860, 796, 2008, 2...
## $ season    <dbl> -2158.8893, 692.6571, 1087.5708, 1052.3294, 939.9377...
## $ trend     <dbl> 1496.712, 1511.009, 1525.307, 1539.604, 1553.901, 15...
## $ remainder <dbl> 1212.1777, -1191.6661, -1097.8773, -889.9334, -797.8...技术支持

在剩余项中检测异常

完成时间序列分解后，若剩余项具备用于异常检测的所需特性，则再次处理并生成三个新列。

remainder_l1：剩余项的下限。
remainder_l2：剩余项的上限。
anomaly：指示该观测是否为异常的列。

异常是会扭曲分布的高杠杆点。anomalize 实现了两种对高杠杆点有抗性的方式：

IQR：四分位距
GESD：广义极端学生化偏差检验

IQR

该方法与 forecast 包中 tsoutliers() 函数所用方法类似。IQR 取一个分布，使用第 25% 与第 75% 四分位范围来确定剩余项的分布。默认将上下限设为内四分位距的 3 倍，任何超出该范围的剩余项都视为异常。

GESD

GESD 会逐步评估异常，移除影响最严重的点，并重新计算检验统计量和临界值。更直观地说，就是在迭代识别异常后重新计算范围。

IQR 与 GESD 各有优缺点。IQR 相对更快，因为没有循环参与；但由于异常会偏移中位数，而 GESD 在流程中去除了该影响，因此 GESD 的准确性往往更高。

purrr_anomaly = purrr_anomaly%>%
  anomalize(remainder)

purrr_anomaly%>% glimpse()
技术支持

## Observations: 425
## Variables: 8
## $ date         <date> 2017-01-01, 2017-01-02, 2017-01-03, 2017-01-04, ...
## $ observed     <dbl> 550, 1012, 1515, 1702, 1696, 1613, 860, 796, 2008...
## $ season       <dbl> -2158.8893, 692.6571, 1087.5708, 1052.3294, 939.9...
## $ trend        <dbl> 1496.712, 1511.009, 1525.307, 1539.604, 1553.901,...
## $ remainder    <dbl> 1212.1777, -1191.6661, -1097.8773, -889.9334, -79...
## $ remainder_l1 <dbl> -4330.511, -4330.511, -4330.511, -4330.511, -4330...
## $ remainder_l2 <dbl> 4400.459, 4400.459, 4400.459, 4400.459, 4400.459,...
## $ anomaly      <chr> "No", "No", "No", "No", "No", "No", "No", "No", "...技术支持

异常上下界转换

工作流的最后一步是使用 time_recompose 在观测值周围创建上下界。它会将 season、trend、remainder_l1 和 remainder_l2 重构为新的界限：

recomposed_l1：围绕观测值的离群下界。
recomposed_l2：围绕观测值的离群上界。

purrr_anomaly = purrr_anomaly%>%
  time_recompose()

purrr_anomaly%>% glimpse()
技术支持

## Observations: 425
## Variables: 10
## $ date          <date> 2017-01-01, 2017-01-02, 2017-01-03, 2017-01-04,...
## $ observed      <dbl> 550, 1012, 1515, 1702, 1696, 1613, 860, 796, 200...
## $ season        <dbl> -2158.8893, 692.6571, 1087.5708, 1052.3294, 939....
## $ trend         <dbl> 1496.712, 1511.009, 1525.307, 1539.604, 1553.901...
## $ remainder     <dbl> 1212.1777, -1191.6661, -1097.8773, -889.9334, -7...
## $ remainder_l1  <dbl> -4330.511, -4330.511, -4330.511, -4330.511, -433...
## $ remainder_l2  <dbl> 4400.459, 4400.459, 4400.459, 4400.459, 4400.459...
## $ anomaly       <chr> "No", "No", "No", "No", "No", "No", "No", "No", ...
## $ recomposed_l1 <dbl> -4992.689, -2126.845, -1717.634, -1738.578, -183...
## $ recomposed_l2 <dbl> 3738.281, 6604.125, 7013.336, 6992.392, 6894.298...技术支持

异常可视化

purrr_anomaly%>%
  plot_anomaly_decomposition()+
  ggtitle("Plotting Anomalies")
技术支持

参数调优

若不调整完全依赖数据的参数，异常检测器的建模就是不完整的。

下面进入参数调整。由于工作流的每个层级各自完成不同任务，因此参数也各不相同。

分解参数的调整

如前所述，趋势与季节性是时间序列分解的基础，因此在 time_decompose 函数中也应在频率与趋势上进行调整。默认值会自动分配：两种方法（STL、Twitter）的频率默认均为 7 天；趋势在 STL 中默认 91 天，在 Twitter 方法中默认 85 天。

您可以根据需要调整其中一个或两个参数，但调整前请仔细观察；盲目修改可能导致分解过程过拟合或欠拟合。

purrr_package %>%
  time_decompose(count, frequency = "auto", trend = "2 weeks")%>%
  anomalize(remainder)%>%
  plot_anomaly_decomposition()+
  ggtitle("Trend = 2 Weeks / Frequency = auto ")
技术支持

异常检测参数的调整

如上所述，这里依据剩余项的取值计算界限来判定离群点。alpha 和 max_anoms 是控制 anomalize() 函数的两个参数。简单说，alpha 控制界带宽度，默认 0.05；减少其值会扩大区间带，使某一点更难被判为异常。

purrr_package%>%
  time_decompose(count)%>%
  anomalize(remainder, alpha = 0.05)%>%
  time_recompose()%>%
  plot_anomalies(time_recompose = T)+
  ggtitle("alpha = 0.05")
技术支持

purrr_package%>%
  time_decompose(count)%>%
  anomalize(remainder, alpha = 0.025)%>%
  time_recompose()%>%
  plot_anomalies(time_recompose = T)+
  ggtitle("alpha = 0.025")
技术支持

max_anom 参数控制可被视为异常的数据百分比。当 alpha 难以调优且您希望聚焦于最强的异常时，该参数很有用。

purrr_package%>%
  time_decompose(count)%>%
  anomalize(remainder, alpha = 0.2, max_anoms = 0.2)%>%
  time_recompose()%>%
  plot_anomalies(time_recompose = T)+
  ggtitle("20% anomaly Allowed")
技术支持

purrr_package%>%
  time_decompose(count)%>%
  anomalize(remainder, alpha = 0.2, max_anoms = 0.05)%>%
  time_recompose()%>%
  plot_anomalies(time_recompose = T)+
  ggtitle("5% anomaly Allowed")
技术支持

参考资料

如果您想进一步学习 R，请参加 DataCamp 的 Introduction to Time Series Analysis 课程。

查看我们的 Time Series Analysis using R: Tutorial。

主题

R 数据科学

主题

技术