시계열 이상치 탐색 in R

이상치 탐색

이상치는 시계열 데이터 상의 추세나 계절성에 반하여 나타나는 특별한 데이터를 말하는데 영어로는 outlier 혹은 anomaly라고 한다. 이 이상치는 측정상의 오류나 데이터 자체의 오류일 수도 있지만 특정한 이유로 인해 일시적으로 발생된 데이터일 수도 있다. 시계열 데이터가 아닌 일반 관측치 데이터의 경우는 데이터의 분포에서 IQR 값이 1.5 IQR을 넘어가는 값을 이상치로 보기도 하고 회귀분석을 통해 이상치를 찾아낼 수 있다. 그러나 시계열 데이터는 추세와 계절성이라는 데이터 자체적인 특성이 있기 때문에 일반적 관측치 데이터에서는 측정되지 않는 이상치를 가진다.

이상치는 그 원인을 파악하지 않고 분석에서 제외하거나 다른 값으로 대체하는 것은 피해야한다. 이 이상치를 통해서 시계열 데이터의 숨겨진 특성이나 추가적 정보를 얻을 수 있다.

따라서 어떤 값이 이상치인지를 탐색하는 것이 매우 중요하다. 이 포스트에서는 이상치를 탐색하기 위해 사용하는 방법에 대해 알아본다.

이상치 데이터의 탐지를 위해 월별 취업자수를 사용한다.

library(forecast)

employees.ts[,2] |> autoplot() +
  labs(title = '취업자수 추이', y = '취업자수', x = '기간')

forecast::tsoutlier(), forecast::tsclean()

R에서 사용하는 대표적인 시계열 패키지인 forecast 패키지에는 이상치를 발견하기 위한 tsoutlier()와 이상치를 처리하는 tsclean()을 제공한다. 이 두개의 함수를 사용하기 위해서는 ts 클래스로 설정된 시계열 데이터가 필요하다.

tsoutlier()는 STL 분해를 통해 시계열 데이터를 계절성 데이터, 추세 데이터, 리마인더 데이터의 덧셈 모델로 분해하고 추세와 계절성을 일정 부분 벗어난 데이터를 이상치로 찾아내어 이 데이터의 인덱스와 보정값을 출력한다.

tsoutliers(employees.ts[,2])

## $index
## [1]  2 85 86 88
## 
## $replacements
## [1] 24355.60 26557.37 26409.16 26829.50

tsclean()은 tsoutlier()에서 발견된 이상치가 보정된 시계열 데이터를 생성해준다.

tsclean(employees.ts[,2])

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2013 24287.00 24355.60 24736.00 25322.00 25610.00 25686.00 25681.00 25513.00
## 2014 25050.00 25116.00 25463.00 25985.00 26112.00 26178.00 26280.00 26183.00
## 2015 25392.00 25471.00 25766.00 26153.00 26431.00 26440.00 26538.00 26369.00
## 2016 25646.00 25615.00 25980.00 26325.00 26613.00 26718.00 26765.00 26696.00
## 2017 25878.00 25979.00 26443.00 26744.00 26992.00 27020.00 27078.00 26904.00
## 2018 26213.00 26083.00 26555.00 26868.00 27064.00 27126.00 27083.00 26907.00
## 2019 26232.00 26346.00 26805.00 27038.00 27322.00 27408.00 27383.00 27358.00
## 2020 26557.37 26409.16 26609.00 26829.50 26930.00 27055.00 27106.00 27085.00
##           Sep      Oct      Nov      Dec
## 2013 25701.00 25798.00 25795.00 25248.00
## 2014 26213.00 26247.00 26260.00 25679.00
## 2015 26487.00 26519.00 26469.00 26098.00
## 2016 26697.00 26746.00 26762.00 26347.00
## 2017 27011.00 27026.00 27019.00 26604.00
## 2018 27055.00 27090.00 27184.00 26638.00
## 2019 27404.00 27509.00 27515.00 27154.00
## 2020 27012.00 27088.00 27241.00 26526.00

tsclean(employees.ts[,2]) |> autoplot()

timetk::plot_anomaly_diagnostics(), timetk::tk_anomaly_diagnostics()

timetk 패키지에서는 이상치를 발견해서 플롯으로 알려주는 plot_anomaly_diagnostics()와 이상치 분석에 사용된 데이터를 보여주는 tk_anomaly_diagnostics()를 제공한다. timetk 패키지를 사용하기 때문에 tsibble 클래스로 된 시계열 데이터를 사용해야 한다.

이 함수들도 tsoutlier()와 유사하게 STL 분해 알고리즘을 사용하여 추세, 계절성, 리마인더로 분해하고 이 중 나머지 값이 이상치에 해당하는지를 판단하여 해당 데이터를 이상치로 판단한다.

나머지 값의 이상치를 판단하는 기준은 나머지 값의 분포에서 IQR의 배수를 alpha 매개변수를 통해 찾아내는데 alpha가 기본값인 0.05로 설정되면 3*IQR로 계산되어 이상치를 판단한다. 따라서 이상치 발견을 보다 민감하게 하기 위해서는 alpha 값을 높인다.

library(tsibble)
library(timetk)

employees.tsibble %>%
  plot_anomaly_diagnostics(time, total, .alpha = 0.05, .interactive = FALSE)

employees.tsibble %>%
  plot_anomaly_diagnostics(time, total, .alpha = 0.1, .interactive = FALSE)

tk_anomaly_diagnostics()는 앞서 plot_anomaly_diagnostics()에서 이상치를 찾아내기 위한 기초 데이터를 보여준다. 아래의 데이터를 보면 reminder_l1과 reminder_l2의 값이 나오는데 reminder가 이 두 값의 범위를 벗어나면 이상치로 판단하게 된다.

employees.tsibble |>
  tk_anomaly_diagnostics(time, total) |>
  filter(lubridate::year(time) >= 2020) 

time	observed	season	trend	remainder	seasadj	remainder_l1	remainder_l2	anomaly	recomposed_l1	recomposed_l2
2020-01-01	26800	-666.06647	27019.66	446.404742	27466.07	-189.221	194.8856	Yes	26164.37	26548.48
2020-02-01	26838	-644.60435	26971.94	510.666958	27482.60	-189.221	194.8856	Yes	26138.11	26522.22
2020-03-01	26609	-313.99114	26924.21	-1.221903	26922.99	-189.221	194.8856	No	26421.00	26805.11
2020-04-01	26562	48.45603	26891.38	-377.840141	26513.54	-189.221	194.8856	Yes	26750.62	27134.73
2020-05-01	26930	272.94896	26858.56	-201.504131	26657.05	-189.221	194.8856	Yes	26942.28	27326.39
2020-06-01	27055	303.19452	26842.90	-91.091772	26751.81	-189.221	194.8856	No	26956.87	27340.98
2020-07-01	27106	313.54494	26827.24	-34.784270	26792.46	-189.221	194.8856	No	26951.56	27335.67
2020-08-01	27085	188.46232	26814.32	82.216256	26896.54	-189.221	194.8856	No	26813.56	27197.67
2020-09-01	27012	227.03723	26801.40	-16.440735	26784.96	-189.221	194.8856	No	26839.22	27223.33
2020-10-01	27088	262.62739	26789.93	35.438655	26825.37	-189.221	194.8856	No	26863.34	27247.45
2020-11-01	27241	230.75411	26778.46	231.781476	27010.25	-189.221	194.8856	Yes	26820.00	27204.10
2020-12-01	26526	-222.36354	26767.39	-19.027155	26748.36	-189.221	194.8856	No	26355.81	26739.91

anomalize 패키지

anomalize 패키지는 시계열 이상치의 탐색, 플로팅, 전처리 등을 위한 다양한 함수를 제공한다. 이 패키지는 timetk의 확장판 패키지로 볼 수 있는기 때문에 사용하는 데이터도 tsibble 클래스 데이터이어야 한다.

이 패키지를 사용할 때는 보통 다음의 세 함수를 사용하는 형태로 사용한다.

1. time_decompose(): 시계열을 계절, 추세 및 나머지 구성 요소로 분해한다.
2. anomalize(): 나머지 구성요소에 이상 탐지 방법을 적용한다.
3. time_recompose(): “정상” 데이터와 비정상 데이터를 구분하는 한계를 계산한다.

if (!require(anomalize)) {
  install.packages('anomalize')
  library(anomalize)
}

employees.tsibble |> 
   time_decompose(total, merge = TRUE) %>%
    anomalize(remainder) %>% 
    time_recompose()

## # A time tibble: 96 x 12
## # Index: time
##    time       total employees.edu observed season  trend remainder remainder_l1
##    <date>     <int>         <int>    <dbl>  <dbl>  <dbl>     <dbl>        <dbl>
##  1 2013-01-01 24287          1710    24287 -677.  25047.    -83.0         -270.
##  2 2013-02-01 24215          1681    24215 -649.  25099.   -235.          -270.
##  3 2013-03-01 24736          1716    24736 -296.  25151.   -119.          -270.
##  4 2013-04-01 25322          1745    25322   50.1 25203.     68.9         -270.
##  5 2013-05-01 25610          1774    25610  281.  25255.     74.8         -270.
##  6 2013-06-01 25686          1786    25686  306.  25306.     73.8         -270.
##  7 2013-07-01 25681          1813    25681  321.  25358.      2.50        -270.
##  8 2013-08-01 25513          1811    25513  198.  25409.    -93.8         -270.
##  9 2013-09-01 25701          1794    25701  235.  25460.      6.07        -270.
## 10 2013-10-01 25798          1790    25798  257.  25511.     30.1         -270.
## # ... with 86 more rows, and 4 more variables: remainder_l2 <dbl>,
## #   anomaly <chr>, recomposed_l1 <dbl>, recomposed_l2 <dbl>

이렇게 생성된 데이터는 plot_anomalies()를 통해 이상치를 플로팅 할 수 있다.

employees.tsibble |> 
   time_decompose(total, merge = TRUE) %>%
    anomalize(remainder) %>% 
    time_recompose() |> 
    plot_anomaly_decomposition()