トップ 履歴 一覧 Farm ソース 検索 ヘルプ RSS ログイン

anomalizeパッケージによる異常値検出と補正


キーワード

最終更新時間:2019年10月09日 01時56分46秒
アフィリエイト・広告について
プライバシーポリシー

この記事はCleaning Anomalies to Reduce Forecast Error by 9% with anomalizeを翻訳しつつ、コードを一部修正し、実際に実行した結果で書き換えたものです。なお、記事の中核となる異常値を補正し、予測値を返す関数が記事の通りだと動作しなかったため、パッケージのvignetteをもとに対処しました。また、ここで言う異常値は、Outlier, 外れ値と訳されることもあります。

はじめに

このチュートリアルでは、私達がどのようにanomalizeパッケージの clean_anomalies() 関数を使い、予測精度を約9%向上させたかを紹介します。

  • 使用するRパッケージ
    • anomalize: 時系列データの異常値検出を行う

私達は、モデル作成・予測を行う前にデータから異常値を検出、補正することで、予測精度を向上させました。これは、anomalizeパッケージの clean_anomalies() 関数を一連のモデル作成・予測ワークフローに組み込んだユースケースとして実現しました。

ワークフロー

私達は、モデル作成の前に、時系列データから異常値を補正するために、以下のワークフローで取り組みました。

  1. 異常値を特定する: time_decompose() 関数で時系列データを分解し、(季節成分、トレンドを除いた) 残差に anomalize() 関数を適用しました。
  2. 異常値を補正する: clean_anomalies() 関数を使い、異常値をトレンドと季節成分から推定される値で置換しました。
  3. モデルを作成し予測する: 学習用データセットからモデルを作成し、検証用データセットを予測しました。異常値補正をした場合、しなかった場合で精度を比較しました。

ステップ1: ライブラリのロード

以下のライブラリを使用します。(訳注: anomalize, tidyquant, timetkの各パッケージはCRANから install.packages() 関数でインストールできます)

library(tidyverse) # データ操作と可視化のためのライブラリ
library(tidyquant) # ビジネス向けのggplot2テーマ集
library(anomalize) # 時系列データからの異常値検出、補正
library(timetk) # 時系列データの機械学習のためのライブラリ
library(knitr) # kable() 関数のため

ステップ2: データの取得

このチュートリアルでは、anomalizeパッケージが提供する tidyverse_cran_downloads データセットを使います。これは、2017年1月1日から2018年3月1日までの期間の、"tidy" なデータ処理を実現するパッケージ群のダウンロード数です。ここでは、ダウンロード数に異常値が含まれるいくつかのパッケージのうち、lubridateパッケージを取り上げます (他のパッケージでチュートリアルを進めても構いません)。

(その前に、異常値検出を行う一連のフローが取り上げられています)

tidyverse_cran_downloads %>%
    time_decompose(count) %>%
    anomalize(remainder) %>%
    time_recompose() %>%
    plot_anomalies(ncol = 3, alpha_dots = 0.3)

まず、filter() 関数でlubridateパッケージのログだけを抽出します。

lubridate_tbl <- tidyverse_cran_downloads %>%
  ungroup() %>%
  filter(package == "lubridate")

以下が実験条件を可視化したものです。学習用データは2018年1月1日までのデータとします。

ggplot(lubridate_tbl, aes(x = date, y = count)) + geom_line() +
geom_vline(xintercept = as.Date("2018-01-01"), color = "red", size = 2) +
annotate("text", x = as.Date("2017-07-01"), y = 15000, label = "学習用データ") +
annotate("text", x = as.Date("2018-02-10"), y = 15000, label = "検証用データ", color = "red") +
theme_tq() + ggtitle(label = "lubridateパッケージのダウンロード数",
                     subtitle = "学習用データ: 2017年12月31日まで 検証用データ: 2018年1月1日以降")

ステップ3: 異常値検出と補正のためのワークフロー

以下のようなワークフローで、異常値を検出、補正します。

  1. counts 列を time_decompose() 関数で分解します。この結果、季節成分、トレンド、残差が得られます。
  2. 負の値 (もしあれば) を修正します。
  3. 残差に対して anomalize() 関数を適用します。結果として、正常値の上限、下限、個別のデータが異常値であるかどうかが返ってきます。
  4. clean_anomalies() 関数を使い、異常値をトレンドと季節成分から補間した値で置換した結果をデータセットに追加します。
lubridate_anomalized_tbl <- lubridate_tbl %>%
    # 1. ダウンロード数 (count) を季節成分、トレンド、残差に分解する
    time_decompose(count) %>%
    # 2. 負の値の補正 (負の値をゼロに置換する)
    mutate(observed = ifelse(observed < 0, 0, observed)) %>%
    # 3. 異常値の検出
    anomalize(remainder) %>%
    # 4. 異常値の修正
    clean_anomalies()

# 観測値と補正後の値を比較する
lubridate_anomalized_tbl %>% 
  filter(anomaly == "Yes") %>%
  select(date, anomaly, observed, observed_cleaned) %>%
  head() %>% 
  kable()

以下の図は、元データと異常値を補正したデータを比較したものです。以降で、異常値を補正した場合、しなかった場合で回帰モデルの精度にどのような影響が出るか検証します。

ggplot(lubridate_anomalized_tbl, aes(x = date, y = observed)) + geom_line() +
geom_line(aes(x = date, y = observed_cleaned), color = "#009999") +
geom_point(data = lubridate_anomalized_tbl %>% filter(anomaly == "Yes"),aes(x = date, y = observed), size = 3, shape = 1) +
theme_tq() + ggtitle(label = "異常値が高レバレッジ値を作り出す", subtitle = "高レバレッジ値が回帰モデルを上方 / 下方に変化させる")

ステップ4: lubridateパッケージのダウンロード数を予測する

はじめに、forecast_mae() 関数を作成します。この関数は、元データと補正後のデータの両方を受け取り、実際のダウンロード数と予測したダウンロード数を返します。詳細は下部の付録を参照してください。

  ステップ4.1: 異常値補正をしないモデル (高レバレッジ値あり)

まず、異常値を補正しない状態でモデルを作成、予測します。手順は以下のような流れです。

  • forecast_mae() 関数は、元データを受け取り、実測値と予測値を返します。
  • 目的変数の平方根を取ります。
  • モデルは、sqrt(元データ) ~ 行番号 + 年 + 四半期 + 月 + 曜日 という線形モデルです。
  • 予測値を2乗したものを出力します。負の値は0に置換します。
lubridate_forecast_with_anomalies_tbl <- lubridate_anomalized_tbl %>%
    # 関数の定義については付録を参照
    forecast_mae(col_train = observed, col_test = observed, trans = "sqrt")

  実測値と予測値の比較

予測値を実測値に重ねて表示しています。

ggplot(lubridate_anomalized_tbl, aes(x = date, y = observed)) + geom_line() +
geom_line(aes(x = date, y = observed_cleaned)) +
geom_point(data = lubridate_anomalized_tbl %>% filter(anomaly == "Yes"), aes(x = date, y = observed), size = 3, shape = 1) +
geom_line(data = lubridate_forecast_with_anomalies_tbl, aes(x = date, y = pred), color = "red", size = 1.5) +
theme_tq() + ggtitle(label = "異常値補正なし", subtitle = "高レバレッジ値が予測値を押し上げている")

予測値が、異常値 (高レバレッジ値) によって、全体に高めにずれていることがわかります。

ggplot(lubridate_anomalized_tbl, aes(x = date, y = observed)) + geom_line() +
geom_line(data = lubridate_forecast_with_anomalies_tbl, aes(x = date, y = pred), color = "red", size = 1) +
geom_smooth(data = lubridate_anomalized_tbl, aes(x = date, y = observed), span = 0.5, color = "gray20") +
geom_smooth(data = lubridate_forecast_with_anomalies_tbl, aes(x = date, y = pred), span = 0.5, color = "red", size = 1.5) +
xlim(as.Date("2018-01-01"), as.Date("2018-03-01")) +
theme_tq() + ggtitle(label = "予測値が高レバレッジ値により上振れしている", subtitle = "異常値補正をしていないデータによる予測")

  誤差の算出

平均絶対誤差 (MAE) は約1570でした。これは、毎日のダウンロード数に約1570の誤差が含まれることを意味します。

lubridate_forecast_with_anomalies_tbl$mae
[1] 1569.613

  ステップ4.2: 異常値補正をしたモデル

次に、clean_anomalies() 関数で異常値を補正したデータを使い、モデルを作成して精度評価します。

  • forecast_mae() 関数は、補正済データを受け取り、実測値と予測値を返します。
  • 目的変数の平方根を取って、モデルを構築します。
  • モデルは、sqrt(補正済データ) ~ 行番号 + 年 + 四半期 + 月 + 曜日 という線形モデルです。
  • 予測値を2乗したものを出力します。負の値は0に置換します。
lubridate_forecast_without_anomalies_tbl <- lubridate_anomalized_tbl %>%
    # 関数の定義については付録を参照
    forecast_mae(col_train = observed_cleaned, col_test = observed, trans = "sqrt")

  実測値と予測値の比較

予測値を実測値に重ねて表示しています。補正済データは黄色で図示しています。

ggplot(lubridate_tbl %>% filter(date < ymd("2018-01-01")), aes(x = date, y = count)) + geom_line(color = "#EEC900") +
geom_line(data = lubridate_tbl %>% filter(date >= ymd("2018-01-01")), aes(x = date, y = count)) +
geom_line(data = lubridate_forecast_without_anomalies_tbl, aes(x = date, y = pred), color = "red", size = 1.5) +
theme_tq() + ggtitle(label = "異常値補正あり", subtitle = "異常値をclean_anomalies()関数で除去した結果")

予測区間を拡大してみると、異常値を補正したモデルでは、トレンドをよく表現できているように見えます。

ggplot(lubridate_anomalized_tbl, aes(x = date, y = observed)) + geom_line() +
geom_line(data = lubridate_forecast_without_anomalies_tbl, aes(x = date, y = pred), color = "red", size = 1) +
geom_smooth(data = lubridate_anomalized_tbl, aes(x = date, y = observed), span = 0.5, color = "gray20") +
geom_smooth(data = lubridate_forecast_without_anomalies_tbl, aes(x = date, y = pred), span = 0.5, color = "red", size = 1.5) +
xlim(as.Date("2018-01-01"), as.Date("2018-03-01")) +
theme_tq() + ggtitle(label = "未知のデータに対する予測値が良くフィットしている", subtitle = "異常値をclean_anomalies()関数で補正したデータによる予測")

  誤差の算出

平均絶対誤差 (MAE) は約1435でした。これは、毎日のダウンロード数に約1435の誤差が含まれることを意味します。

lubridate_forecast_without_anomalies_tbl$mae
[1] 1434.631

9.4%の精度向上

anomalizeパッケージの clean_anomalies() 関数を使うことで、MAEで評価した際に、予測における誤差が9.4%低減できました。

(lubridate_forecast_without_anomalies_tbl$mae - lubridate_forecast_with_anomalies_tbl$mae)
/ lubridate_forecast_without_anomalies_tbl$mae
[1] -0.09408838

まとめ

異常値 (高レバレッジ値) を検出、除去して予測モデルを作成することは、モデルの精度を向上させました。anomalizeパッケージの clean_anomalies() 関数は、異常値検出、除去の仕組みを、分析ワークフローに容易に組み込むことができます。詳細は、anomalizeパッケージのドキュメントを参照してください。

付録: forecast_mae() 関数について

訳注: オリジナルの記事では forecast_download() 関数が紹介されていますが、書いてある通りにしても動作しません。いろいろと調査をした結果、anomalizeパッケージのvignetteをもとに改変すると、記事が意図しているであろう動作を実現できました。

forecast_mae() 関数は、以下のような動作をします。

  • データを学習用と検証用に分割します。分割する区切りは、sep オプションで指定します。
  • trans オプションで指定した方法でデータを変換します。対数変換 (log1p()) か平方根 (sqrt()) または変換しない (none) が指定できます。
  • 日次の時系列データを指定し、予測モデルを作成します。col_train オプションで対象列を指定します。このチュートリアルでは、元データと補正済データをそれぞれ指定し、比較しました。
  • 線形回帰モデルを作成し、検証用データの予測を行います。col_test オプションで対象列を指定します。
  • 実測値と予測値を比較し、MAE (Mean absolute error) を算出します。
forecast_mae <- function(data, col_train, col_test, trans = c("none", "log1p", "sqrt"), sep = "2018-01-01") {
  predict_expr <- enquo(col_train)
  actual_expr <- enquo(col_test)
  
  # 学習用データと検証用データの分割
  train_tbl <- data %>% filter(date < ymd(sep))
  test_tbl  <- data %>% filter(date >= ymd(sep))
  
  # データ変換が指定されている場合の変換式を作成
  pred_form <- quo_name(predict_expr)
  if (trans != "none") pred_form <- str_glue("{trans}({pred_form})")
  
  # 学習用データで線形モデルを作成
  model_formula <- as.formula(paste0(pred_form, " ~ index.num + year + quarter + month.lbl + day + wday.lbl"))
  
  model_glm <- train_tbl %>%
    tk_augment_timeseries_signature() %>%
    glm(model_formula, data = .)
  
  # 検証用データの予測
  # 警告メッセージの回避
  suppressWarnings({
    prediction <- predict(model_glm, newdata = test_tbl %>% tk_augment_timeseries_signature())
    if (trans == "log1p") prediction <- expm1(prediction)
    if (trans == "sqrt")  prediction <- ifelse(prediction < 0, 0, prediction)^2
    actual     <- test_tbl %>% pull(!! actual_expr)
  })
  
  # MAEの算出
  mae <- mean(abs(prediction - actual))
  
  # グラフ描画用にMAE, 日時, 予測値を返す
  res <- list(mae = mae, date = test_tbl$date, pred = prediction)
  return(res)
}


関連ページ: R言語を学ぶための参考書籍リスト
カテゴリ: [R, データ分析, 時系列分析]