R/anomalize.R

In business-science/timetk: A Tool Kit for Working with Time Series

#' Automatic group-wise Anomaly Detection
#'
#' `anomalize()` is used to detect anomalies in time series data,
#' either for a single time series or for multiple time series grouped by a specific column.
#'
#' @param .data A `tibble` or `data.frame` with a time-based column
#' @param .date_var A column containing either date or date-time values
#' @param .value A column containing numeric values
#' @param .frequency Controls the seasonal adjustment (removal of seasonality).
#'  Input can be either "auto", a time-based definition (e.g. "2 weeks"),
#'  or a numeric number of observations per frequency (e.g. 10).
#'  Refer to [tk_get_frequency()].
#' @param .trend Controls the trend component.
#'  For STL, trend controls the sensitivity of the LOESS smoother, which is used to remove the remainder.
#'  Refer to [tk_get_trend()].
#' @param .method The outlier detection method. Default: "stl". Currently
#'  "stl" is the only method. "twitter" is planned.
#' @param .iqr_alpha Controls the width of the "normal" range. Lower values are more conservative
#'  while higher values are less prone to incorrectly classifying "normal" observations.
#' @param .clean_alpha Controls the threshold for cleaning
#' the outliers. The default is 0.75, which means that the anomalies will be
#' cleaned using the 0.75 * lower or upper bound of the recomposed time series,
#' depending on the direction of the anomaly.
#' @param .max_anomalies The maximum percent of anomalies permitted to be identified.
#' @param .message A boolean. If `TRUE`, will output information related to automatic frequency
#' and trend selection (if applicable).
#'
#' @return
#' A `tibble` or `data.frame` with the following columns:
#'  - observed: original data
#'  - seasonal: seasonal component
#'  - seasadaj: seasonal adjusted
#'  - trend: trend component
#'  - remainder: residual component
#'  - anomaly: Yes/No flag for outlier detection
#'  - anomaly score: distance from centerline
#'  - anomaly direction: -1, 0, 1 inidicator for direction of the anomaly
#'  - recomposed_l1: lower level bound of recomposed time series
#'  - recomposed_l2: upper level bound of recomposed time series
#'  - observed_clean: original data with anomalies interpolated
#'
#' @details
#'
#' The `anomalize()` method for anomaly detection that implements a 2-step process to
#' detect outliers in time series.
#'
#' __Step 1: Detrend & Remove Seasonality using STL Decomposition__
#'
#' The decomposition separates the "season" and "trend" components from the "observed" values
#' leaving the "remainder" for anomaly detection.
#'
#' The user can control two parameters: frequency and trend.
#'
#' 1. `.frequency`: Adjusts the "season" component that is removed from the "observed" values.
#' 2. `.trend`: Adjusts the trend window (t.window parameter from [stats::stl()] that is used.
#'
#' The user may supply both `.frequency` and `.trend` as time-based durations (e.g. "6 weeks") or
#' numeric values (e.g. 180) or "auto", which predetermines the frequency and/or trend based on
#' the scale of the time series using the [tk_time_scale_template()].
#'
#' __Step 2: Anomaly Detection__
#'
#' Once "trend" and "season" (seasonality) is removed, anomaly detection is performed on the "remainder".
#' Anomalies are identified, and boundaries (recomposed_l1 and recomposed_l2) are determined.
#'
#' The Anomaly Detection Method uses an inner quartile range (IQR) of +/-25 the median.
#'
#' _IQR Adjustment, alpha parameter_
#'
#' With the default `alpha = 0.05`, the limits are established by expanding
#' the 25/75 baseline by an IQR Factor of 3 (3X).
#' The _IQR Factor = 0.15 / alpha_ (hence 3X with alpha = 0.05):
#'
#' - To increase the IQR Factor controlling the limits, decrease the alpha,
#' which makes it more difficult to be an outlier.
#' - Increase alpha to make it easier to be an outlier.
#'
#'
#' - The IQR outlier detection method is used in `forecast::tsoutliers()`.
#' - A similar outlier detection method is used by Twitter's `AnomalyDetection` package.
#' - Both Twitter and Forecast tsoutliers methods have been implemented in Business Science's `anomalize`
#'  package.
#'
#'
#' @references
#' 1. CLEVELAND, R. B., CLEVELAND, W. S., MCRAE, J. E., AND TERPENNING, I.
#'  STL: A Seasonal-Trend Decomposition Procedure Based on Loess.
#'  Journal of Official Statistics, Vol. 6, No. 1 (1990), pp. 3-73.
#'
#' 2. Owen S. Vallis, Jordan Hochenbaum and Arun Kejariwal (2014).
#'  A Novel Technique for Long-Term Anomaly Detection in the Cloud. Twitter Inc.
#'
#' @examples
#' library(dplyr)
#'
#' walmart_sales_weekly %>%
#'     filter(id %in% c("1_1", "1_3")) %>%
#'     group_by(id) %>%
#'     anomalize(Date, Weekly_Sales)
#'
#' @name anomalize
#' @export
#'
anomalize <- function(
        .data,
        .date_var,
        .value,

        .frequency = "auto",
        .trend = "auto",

        .method = "stl",

        .iqr_alpha = 0.05,
        .clean_alpha = 0.75,
        .max_anomalies = 0.2,

        .message = TRUE
) {

    # Tidyeval Setup
    date_var_expr <- rlang::enquo(.date_var)
    value_expr    <- rlang::enquo(.value)

    # Checks
    if (!is.data.frame(.data)) {
        stop(call. = FALSE, "anomalize(.data) is not a data-frame or tibble. Please supply a data.frame or tibble.")
    }
    if (rlang::quo_is_missing(date_var_expr)) {
        stop(call. = FALSE, "anomalize(.date_var) is missing. Please supply a date or date-time column.")
    }
    if (rlang::quo_is_missing(value_expr)) {
        stop(call. = FALSE, "anomalize(.value) is missing. Please supply a numeric column.")
    }

    if (.method != "stl") {
        stop(call. = FALSE, "anomalize(.method): Only 'stl' is currently implemented.")
    }

    UseMethod("anomalize", .data)

}


#' @export
anomalize.data.frame <- function(
    .data,
    .date_var,
    .value,

    .frequency = "auto",
    .trend = "auto",

    .method = "stl",

    .iqr_alpha = 0.05,
    .clean_alpha = 0.75,
    .max_anomalies = 0.2,

    .message = TRUE
) {

    # STL Decomposition (observed, season, trend, remainder, seasadj)
    ret <- .data %>%
        tk_stl_diagnostics(
            .date_var  = !! rlang::enquo(.date_var),
            .value     = !! rlang::enquo(.value),
            .frequency = .frequency,
            .trend     = .trend,
            .message   = .message
        )

    # Detect Anomalies (remainder_l1, remainder_l2, anomaly)
    ret <- ret %>%
        mutate_anomalies_2(
            target    = remainder,
            alpha     = .iqr_alpha,
            max_anoms = .max_anomalies
        )

    # Recomposition
    ret <- ret %>%
        dplyr::mutate(
            recomposed_l1 = season + trend + remainder_l1,
            recomposed_l2 = season + trend + remainder_l2
        ) %>%
        dplyr::select(-remainder_l1, -remainder_l2)


    # Clean
    ret <- ret %>%
        dplyr::mutate(
            observed_clean = dplyr::case_when(
                anomaly_direction == -1 ~ recomposed_l1 + (1-.clean_alpha)*(recomposed_l2-recomposed_l1)/2,
                anomaly_direction == 1 ~ recomposed_l2 - (1-.clean_alpha)*(recomposed_l2-recomposed_l1)/2,
                TRUE ~ observed
            )
        )


    return(ret)

}

#' @export
anomalize.grouped_df <- function(
        .data,
        .date_var,
        .value,

        .frequency = "auto",
        .trend = "auto",

        .method = "stl",

        .iqr_alpha = 0.05,
        .clean_alpha = 0.75,
        .max_anomalies = 0.2,

        .message = TRUE
) {

    # Tidy Eval Setup
    value_expr    <- rlang::enquo(.value)
    date_var_expr <- rlang::enquo(.date_var)
    group_names   <- dplyr::group_vars(.data)

    # Process groups individually
    .data %>%
        tidyr::nest() %>%
        dplyr::mutate(nested.col = purrr::map(
            .x         = data,
            .f         = function(df) anomalize(
                .data          = df,
                .date_var      = !! date_var_expr,
                .value         = !! value_expr,
                .frequency     = .frequency,
                .trend         = .trend,
                .method        = .method,
                .iqr_alpha     = .iqr_alpha,
                .clean_alpha   = .clean_alpha,
                .max_anomalies = .max_anomalies,
                .message       = .message
            )
        )) %>%
        dplyr::select(-"data") %>%
        tidyr::unnest(cols = nested.col) %>%
        dplyr::group_by_at(.vars = group_names)

}

# UTILITIES ----

mutate_anomalies_2 <- function(data, target, alpha = 0.05, max_anoms = 0.20) {

    # Checks
    if (missing(target)) stop('Error in anomalize(): argument "target" is missing, with no default', call. = FALSE)

    # Setup
    target_expr <- rlang::enquo(target)
    x           <- data %>% dplyr::pull(!! target_expr)

    # Explicitly call functions
    outlier_tbl <- iqr_2(x = x, alpha = alpha, max_anoms = max_anoms, verbose = TRUE)

    ret <- data %>% dplyr::bind_cols(outlier_tbl)

    return(ret)

}


# This is anomalize::iqr()
iqr_2 <- function(x, alpha = 0.05, max_anoms = 0.2, verbose = FALSE) {

    quantile_x <- stats::quantile(x, prob = c(0.25, 0.75), na.rm = TRUE)
    iq_range   <- quantile_x[[2]] - quantile_x[[1]]
    limits     <- quantile_x + (0.15 / alpha) * iq_range * c(-1, 1)

    outlier_idx      <- ((x < limits[1]) | (x > limits[2]))

    centerline <- sum(limits) / 2
    score <- abs(x - centerline)

    outlier_reported <- dplyr::case_when(
        x < limits[1] ~ "Yes",
        x > limits[2] ~ "Yes",
        TRUE ~ "No"
    )

    direction <- dplyr::case_when(
        x < limits[1] ~ -1,
        x > limits[2] ~ 1,
        TRUE ~ 0
    )

    remainder_l1 <- limits[1]
    remainder_l2 <- limits[2]

    return(
        tibble::tibble(
            anomaly = outlier_reported,
            anomaly_direction = direction,
            anomaly_score = score,
            remainder_l1 = remainder_l1,
            remainder_l2 = remainder_l2
        )
    )

}