diff --git a/DESCRIPTION b/DESCRIPTION index a772116..73740ce 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: gpindex Title: Generalized Price and Quantity Indexes -Version: 0.6.3.9002 +Version: 0.6.3.9004 Authors@R: c( person("Steve", "Martin", role = c("aut", "cre", "cph"), email = "marberts@protonmail.com", @@ -33,4 +33,4 @@ Collate: 'helpers.R' 'means.R' 'weights.R' 'contributions.R' 'price_indexes.R' Config/testthat/edition: 3 VignetteBuilder: quarto Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.2 +RoxygenNote: 7.3.3 diff --git a/NAMESPACE b/NAMESPACE index 39b4301..b298275 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -48,6 +48,7 @@ export(nested_contributions2) export(nested_mean) export(nested_transmute) export(nested_transmute2) +export(outliers) export(paasche_index) export(quantity_index) export(quartile_method) diff --git a/NEWS.md b/NEWS.md index 60c04d7..9a30f2b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,9 @@ - Fixed a bug where `hb_transform()` would give an error with missing values. +- Added a new function `outliers()` to replace existing outlier functions, which +are now deprecated. + ## gpindex 0.6.3 - Bumped minimum version of R to >= 4.1. diff --git a/R/outliers.R b/R/outliers.R index dde13b1..30bd926 100644 --- a/R/outliers.R +++ b/R/outliers.R @@ -2,27 +2,32 @@ #' #' Standard cutoff-based methods for detecting outliers with price relatives. #' -#' Each of these functions constructs an interval of the form \eqn{[b_l(x) - +#' This function constructs an interval of the form \eqn{[b_l(x) - #' c_l \times l(x), b_u(x) + c_u \times u(x)]}{[bl(x) - cl * l(x), bu(x) + cu * #' u(x)]} and assigns a value in `x` as `TRUE` if that value does not -#' belong to the interval, `FALSE` otherwise. The methods differ in how -#' they construct the values \eqn{b_l(x)}{bl(x)}, \eqn{b_u(x)}{bu(x)}, +#' belong to the interval, `FALSE` otherwise. The different methods differ in +#' how they construct the values \eqn{b_l(x)}{bl(x)}, \eqn{b_u(x)}{bu(x)}, #' \eqn{l(x)}, and \eqn{u(x)}. Any missing values in `x` are ignored when #' calculating the cutoffs, but will return `NA`. #' -#' The fixed cutoff method is the simplest, and just uses the interval -#' \eqn{[c_l, c_u]}{[cl, cu]}. -#' #' The quartile method and Tukey algorithm are described in paragraphs 5.113 to #' 5.135 of the CPI manual (2020), as well as by Rais (2008) and Hutton (2008). #' The resistant fences method is an alternative to the quartile method, and is #' described by Rais (2008) and Hutton (2008). The Kimber method is yet another #' alternative. Quantile-based methods often #' identify price relatives as outliers because the distribution is -#' concentrated around 1; setting `a > 0` puts a floor on the minimum +#' concentrated around 1; setting `scale > 0` puts a floor on the minimum #' dispersion between quantiles as a fraction of the median. See the references #' for more details. +# nolint start +#' +#' | | \eqn{b_l(x)}{bl(x)} | \eqn{b_u(x)}{bu(x)} | \eqn{l(x)} | \eqn{u(x)} | +#' | --- | -- | --- | --- | --- | +#' | Quartile | \eqn{Q_2(x)} | \eqn{Q_2(x)} | \eqn{Q_2(x) - Q_1(x)} | \eqn{Q_3(x) - Q_2(x)} | +#' | Resistant fences | \eqn{Q_1(x)} | \eqn{Q_3(x)} | \eqn{Q_3(x) - Q_1(x)} | \eqn{Q_3(x) - Q_1(x)} | +#' | Kimber | \eqn{Q_1(x)} | \eqn{Q_3(x)} | \eqn{Q_2(x) - Q_1(x)} | \eqn{Q_3(x) - Q_2(x)} | #' +# nolint end #' The robust Z-score is the usual method to identify relatives in the #' (asymmetric) tails of the distribution, simply replacing the mean with the #' median, and the standard deviation with the median absolute deviation. @@ -40,12 +45,14 @@ #' #' @param x A numeric vector, usually of price relatives. These can be #' made with, e.g., [back_period()]. -#' @param cu,cl A number giving the upper and lower cutoffs for each -#' element of `x`. -#' @param a A number between 0 and 1 giving the scale factor for the +#' @param upper,lower,cu,cl A number giving the upper and lower cutoffs for +#' each element of `x`. +#' @param method The outlier detection method, one `"quartile"` (the default), +#' `"resistant-fences"`, `"kimber"`, `"robust-z"`, or `"tukey"`. +#' @param scale,a A number between 0 and 1 giving the scale factor for the #' median to establish the minimum dispersion between quartiles for each #' element of `x`. The default does not set a minimum dispersion. -#' @param type See [quantile()]. +#' @param quantile_type,type See [quantile()]. #' #' @returns #' A logical vector, the same length as `x`, that is `TRUE` if the @@ -79,11 +86,11 @@ #' #' x <- rlnorm(10) #' -#' fixed_cutoff(x) -#' robust_z(x) -#' quartile_method(x) -#' resistant_fences(x) # always identifies fewer outliers than above -#' tukey_algorithm(x) +#' outliers(x, method = "robust-z") +#' outliers(x, method = "quartile") +#' # Always identifies fewer outliers than above. +#' outliers(x, method = "resistant-fences") +#' outliers(x, method = "tukey") #' #' log(x) #' hb_transform(x) @@ -91,134 +98,140 @@ #' # Works the same for grouped data. #' #' f <- c("a", "b", "a", "a", "b", "b", "b", "a", "a", "b") -#' grouped(quartile_method)(x, group = f) +#' grouped(outliers)(x, group = f) #' -#' @name outliers #' @export -quartile_method <- function(x, cu = 2.5, cl = cu, a = 0, type = 7) { +outliers <- function( + x, + upper = 2.5, + lower = upper, + method = c( + "quartile", + "resistant-fences", + "kimber", + "robust-z", + "tukey" + ), + scale = 0, + quantile_type = 7 +) { + method <- match.arg(method) x <- as.numeric(x) - cu <- as.numeric(cu) - if (cu < 0) { - stop("'cu' must be greater than 0") + upper <- as.numeric(upper) + if (upper < 0) { + stop("'upper' must be greater than 0") } - cl <- as.numeric(cl) - if (cl < 0) { - stop("'cl' must be greater than 0") + lower <- as.numeric(lower) + if (lower < 0) { + stop("'lower' must be greater than 0") } - a <- as.numeric(a) - if (a < 0 || a > 1) { - stop("'a' must be between 0 and 1") + scale <- as.numeric(scale) + if (scale < 0 || scale > 1) { + stop("'scale' must be between 0 and 1") } + if (method %in% c("quartile", "resistant-fences", "kimber")) { + q <- stats::quantile( + x, + c(0.25, 0.5, 0.75), + names = FALSE, + na.rm = TRUE, + type = quantile_type + ) + if (method == "quartile") { + u <- q[2L] + upper * pmax.int(q[3L] - q[2L], abs(scale * q[2L])) + l <- q[2L] - lower * pmax.int(q[2L] - q[1L], abs(scale * q[2L])) + } else if (method == "resistant-fences") { + iqr <- pmax.int(q[3L] - q[1L], abs(scale * q[2L])) + u <- q[3L] + upper * iqr + l <- q[1L] - lower * iqr + } else { + u <- q[3L] + upper * pmax.int(q[3L] - q[2L], abs(scale * q[2L])) + l <- q[1L] - lower * pmax.int(q[2L] - q[1L], abs(scale * q[2L])) + } + } else if (method == "robust-z") { + med <- stats::quantile( + x, + 0.5, + names = FALSE, + na.rm = TRUE, + type = quantile_type + ) + s <- pmax.int(stats::mad(x, med, na.rm = TRUE), abs(scale * med)) + u <- med + upper * s + l <- med - lower * s + } else if (method == "tukey") { + q <- stats::quantile( + x, + c(0.05, 0.95), + names = FALSE, + na.rm = TRUE, + type = quantile_type + ) + tail <- x < q[1L] | x > q[2L] + ts <- x[x != 1 & !tail] + if (length(ts) == 0L) { + return(tail) + } + # In some versions m is the median. + m <- mean(ts, na.rm = TRUE) + u <- min(m + upper * (mean(ts[ts >= m], na.rm = TRUE) - m), q[2L]) + l <- max(m - lower * (m - mean(ts[ts <= m], na.rm = TRUE)), q[1L]) + } + x > u | x < l +} - q <- stats::quantile( - x, - c(0.25, 0.5, 0.75), - names = FALSE, - na.rm = TRUE, - type = type +#' Quartile method +#' @rdname outliers +#' @export +quartile_method <- function(x, cu = 2.5, cl = cu, a = 0, type = 7) { + warning( + "this function is deprecated and will be removed in a future version;", + " use 'outliers()' instead" ) - x <- x - q[2L] - u <- cu * pmax.int(q[3L] - q[2L], abs(a * q[2L])) - l <- -cl * pmax.int(q[2L] - q[1L], abs(a * q[2L])) - x > u | x < l + outliers(x, cu, cl, method = "quartile", a, type) } #' Resistant fences #' @rdname outliers #' @export resistant_fences <- function(x, cu = 2.5, cl = cu, a = 0, type = 7) { - x <- as.numeric(x) - cu <- as.numeric(cu) - if (cu < 0) { - stop("'cu' must be greater than 0") - } - cl <- as.numeric(cl) - if (cl < 0) { - stop("'cl' must be greater than 0") - } - a <- as.numeric(a) - if (a < 0 || a > 1) { - stop("'a' must be between 0 and 1") - } - - q <- stats::quantile( - x, - c(0.25, 0.5, 0.75), - names = FALSE, - na.rm = TRUE, - type = type + warning( + "this function is deprecated and will be removed in a future version;", + " use 'outliers()' instead" ) - iqr <- pmax.int(q[3L] - q[1L], abs(a * q[2L])) - u <- q[3L] + cu * iqr - l <- q[1L] - cl * iqr - x > u | x < l + outliers(x, cu, cl, method = "resistant-fences", a, type) } #' Kimber method #' @rdname outliers #' @export kimber_method <- function(x, cu = 2.5, cl = cu, a = 0, type = 7) { - x <- as.numeric(x) - cu <- as.numeric(cu) - if (cu < 0) { - stop("'cu' must be greater than 0") - } - cl <- as.numeric(cl) - if (cl < 0) { - stop("'cl' must be greater than 0") - } - a <- as.numeric(a) - if (a < 0 || a > 1) { - stop("'a' must be between 0 and 1") - } - - q <- stats::quantile( - x, - c(0.25, 0.5, 0.75), - names = FALSE, - na.rm = TRUE, - type = type + warning( + "this function is deprecated and will be removed in a future version;", + " use 'outliers()' instead" ) - u <- q[3L] + cu * pmax.int(q[3L] - q[2L], abs(a * q[2L])) - l <- q[1L] - cl * pmax.int(q[2L] - q[1L], abs(a * q[2L])) - x > u | x < l + outliers(x, cu, cl, method = "kimber", a, type) } #' Robust z-score #' @rdname outliers #' @export robust_z <- function(x, cu = 2.5, cl = cu) { - x <- as.numeric(x) - cu <- as.numeric(cu) - if (cu < 0) { - stop("'cu' must be greater than 0") - } - cl <- as.numeric(cl) - if (cl < 0) { - stop("'cl' must be greater than 0") - } - - med <- stats::median(x, na.rm = TRUE) - s <- stats::mad(x, na.rm = TRUE) - x <- x - med - u <- cu * s - l <- -cl * s - x > u | x < l + warning( + "this function is deprecated and will be removed in a future version;", + " use 'outliers()' instead" + ) + outliers(x, cu, cl, method = "robust-z") } #' Fixed cutoff #' @rdname outliers #' @export fixed_cutoff <- function(x, cu = 2.5, cl = 1 / cu) { - x <- as.numeric(x) - cu <- as.numeric(cu) - if (cu < 0) { - stop("'cu' must be greater than 0") - } - cl <- as.numeric(cl) - if (cl < 0) { - stop("'cl' must be greater than 0") - } + warning( + "this function is deprecated and will be removed in a future version;", + " use 'outliers()' instead" + ) x > cu | x < cl } @@ -226,34 +239,11 @@ fixed_cutoff <- function(x, cu = 2.5, cl = 1 / cu) { #' @rdname outliers #' @export tukey_algorithm <- function(x, cu = 2.5, cl = cu, type = 7) { - x <- as.numeric(x) - cu <- as.numeric(cu) - if (cu < 0) { - stop("'cu' must be greater than 0") - } - cl <- as.numeric(cl) - if (cl < 0) { - stop("'cl' must be greater than 0") - } - - q <- stats::quantile( - x, - c(0.05, 0.95), - names = FALSE, - na.rm = TRUE, - type = type + warning( + "this function is deprecated and will be removed in a future version;", + " use 'outliers()' instead" ) - tail <- x < q[1L] | x > q[2L] - ts <- x[x != 1 & !tail] - if (length(ts) == 0L) { - return(tail) - } - # In some versions m is the median. - m <- mean(ts, na.rm = TRUE) - x <- x - m - u <- cu * (mean(ts[ts >= m], na.rm = TRUE) - m) - l <- -cl * (m - mean(ts[ts <= m], na.rm = TRUE)) - x > u | x < l | tail + outliers(x, cu, cl, method = "tukey", quantile_type = type) } #' HB transform diff --git a/man/outliers.Rd b/man/outliers.Rd index 36a9dba..25cf045 100644 --- a/man/outliers.Rd +++ b/man/outliers.Rd @@ -11,6 +11,15 @@ \alias{hb_transform} \title{Outlier detection for price relatives} \usage{ +outliers( + x, + upper = 2.5, + lower = upper, + method = c("quartile", "resistant-fences", "kimber", "robust-z", "tukey"), + scale = 0, + quantile_type = 7 +) + quartile_method(x, cu = 2.5, cl = cu, a = 0, type = 7) resistant_fences(x, cu = 2.5, cl = cu, a = 0, type = 7) @@ -29,14 +38,17 @@ hb_transform(x) \item{x}{A numeric vector, usually of price relatives. These can be made with, e.g., \code{\link[=back_period]{back_period()}}.} -\item{cu, cl}{A number giving the upper and lower cutoffs for each -element of \code{x}.} +\item{upper, lower, cu, cl}{A number giving the upper and lower cutoffs for +each element of \code{x}.} + +\item{method}{The outlier detection method, one \code{"quartile"} (the default), +\code{"resistant-fences"}, \code{"kimber"}, \code{"robust-z"}, or \code{"tukey"}.} -\item{a}{A number between 0 and 1 giving the scale factor for the +\item{scale, a}{A number between 0 and 1 giving the scale factor for the median to establish the minimum dispersion between quartiles for each element of \code{x}. The default does not set a minimum dispersion.} -\item{type}{See \code{\link[=quantile]{quantile()}}.} +\item{quantile_type, type}{See \code{\link[=quantile]{quantile()}}.} } \value{ A logical vector, the same length as \code{x}, that is \code{TRUE} if the @@ -47,26 +59,29 @@ corresponding element of \code{x} is identified as an outlier, Standard cutoff-based methods for detecting outliers with price relatives. } \details{ -Each of these functions constructs an interval of the form \eqn{[b_l(x) - +This function constructs an interval of the form \eqn{[b_l(x) - c_l \times l(x), b_u(x) + c_u \times u(x)]}{[bl(x) - cl * l(x), bu(x) + cu * u(x)]} and assigns a value in \code{x} as \code{TRUE} if that value does not -belong to the interval, \code{FALSE} otherwise. The methods differ in how -they construct the values \eqn{b_l(x)}{bl(x)}, \eqn{b_u(x)}{bu(x)}, +belong to the interval, \code{FALSE} otherwise. The different methods differ in +how they construct the values \eqn{b_l(x)}{bl(x)}, \eqn{b_u(x)}{bu(x)}, \eqn{l(x)}, and \eqn{u(x)}. Any missing values in \code{x} are ignored when calculating the cutoffs, but will return \code{NA}. -The fixed cutoff method is the simplest, and just uses the interval -\eqn{[c_l, c_u]}{[cl, cu]}. - The quartile method and Tukey algorithm are described in paragraphs 5.113 to 5.135 of the CPI manual (2020), as well as by Rais (2008) and Hutton (2008). The resistant fences method is an alternative to the quartile method, and is described by Rais (2008) and Hutton (2008). The Kimber method is yet another alternative. Quantile-based methods often identify price relatives as outliers because the distribution is -concentrated around 1; setting \code{a > 0} puts a floor on the minimum +concentrated around 1; setting \code{scale > 0} puts a floor on the minimum dispersion between quantiles as a fraction of the median. See the references -for more details. +for more details.\tabular{lllll}{ + \tab \eqn{b_l(x)}{bl(x)} \tab \eqn{b_u(x)}{bu(x)} \tab \eqn{l(x)} \tab \eqn{u(x)} \cr + Quartile \tab \eqn{Q_2(x)} \tab \eqn{Q_2(x)} \tab \eqn{Q_2(x) - Q_1(x)} \tab \eqn{Q_3(x) - Q_2(x)} \cr + Resistant fences \tab \eqn{Q_1(x)} \tab \eqn{Q_3(x)} \tab \eqn{Q_3(x) - Q_1(x)} \tab \eqn{Q_3(x) - Q_1(x)} \cr + Kimber \tab \eqn{Q_1(x)} \tab \eqn{Q_3(x)} \tab \eqn{Q_2(x) - Q_1(x)} \tab \eqn{Q_3(x) - Q_2(x)} \cr +} + The robust Z-score is the usual method to identify relatives in the (asymmetric) tails of the distribution, simply replacing the mean with the @@ -88,11 +103,11 @@ set.seed(1234) x <- rlnorm(10) -fixed_cutoff(x) -robust_z(x) -quartile_method(x) -resistant_fences(x) # always identifies fewer outliers than above -tukey_algorithm(x) +outliers(x, method = "robust-z") +outliers(x, method = "quartile") +# Always identifies fewer outliers than above. +outliers(x, method = "resistant-fences") +outliers(x, method = "tukey") log(x) hb_transform(x) @@ -100,7 +115,7 @@ hb_transform(x) # Works the same for grouped data. f <- c("a", "b", "a", "a", "b", "b", "b", "a", "a", "b") -grouped(quartile_method)(x, group = f) +grouped(outliers)(x, group = f) } \references{ diff --git a/tests/Examples/gpindex-Ex.Rout.save b/tests/Examples/gpindex-Ex.Rout.save index 43b21bc..7bfa81d 100644 --- a/tests/Examples/gpindex-Ex.Rout.save +++ b/tests/Examples/gpindex-Ex.Rout.save @@ -745,15 +745,14 @@ Warning in back_period(period) : > > x <- rlnorm(10) > -> fixed_cutoff(x) - [1] TRUE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE -> robust_z(x) +> outliers(x, method = "robust-z") [1] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE -> quartile_method(x) +> outliers(x, method = "quartile") [1] FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE -> resistant_fences(x) # always identifies fewer outliers than above +> # Always identifies fewer outliers than above. +> outliers(x, method = "resistant-fences") [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE -> tukey_algorithm(x) +> outliers(x, method = "tukey") [1] FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE > > log(x) @@ -766,7 +765,7 @@ Warning in back_period(period) : > # Works the same for grouped data. > > f <- c("a", "b", "a", "a", "b", "b", "b", "a", "a", "b") -> grouped(quartile_method)(x, group = f) +> grouped(outliers)(x, group = f) [1] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE > > diff --git a/tests/testthat/test-outliers.R b/tests/testthat/test-outliers.R index 5cdff8e..57f1460 100644 --- a/tests/testthat/test-outliers.R +++ b/tests/testthat/test-outliers.R @@ -3,40 +3,63 @@ set.seed(4321) x <- c(1, 2, 1, 0.5, 1, 10, 1, 0.5, 0.2, 0.05) test_that("outlier methods work", { - expect_equal(fixed_cutoff(x), x > 2.5 | x < 1 / 2.5) expect_equal( - quartile_method(x), + outliers(x, method = "quartile"), x > median(x) + (quantile(x, 0.75) - quantile(x, 0.5)) * 2.5 | x < median(x) - (quantile(x, 0.5) - quantile(x, 0.25)) * 2.5 ) expect_equal( - resistant_fences(x), + outliers(x, method = "resistant-fences"), x > quantile(x, 0.75) + (quantile(x, 0.75) - quantile(x, 0.25)) * 2.5 | x < quantile(x, 0.25) - (quantile(x, 0.75) - quantile(x, 0.25)) * 2.5 ) expect_equal( - kimber_method(x), + outliers(x, method = "kimber"), x > quantile(x, 0.75) + (quantile(x, 0.75) - quantile(x, 0.5)) * 2.5 | x < quantile(x, 0.25) - (quantile(x, 0.5) - quantile(x, 0.25)) * 2.5 ) - expect_true(sum(resistant_fences(x)) <= sum(quartile_method(x))) - expect_equal(robust_z(x), abs(x - median(x)) / mad(x) > 2.5) + expect_true( + sum(outliers(x, method = "resistant-fences")) <= + sum(outliers(x, method = "quartile")) + ) + expect_equal( + outliers(x, method = "robust-z"), + abs(x - median(x)) / mad(x) > 2.5 + ) - expect_equal(tukey_algorithm(integer(0)), logical(0)) - expect_equal(tukey_algorithm(2), FALSE) + expect_equal(outliers(integer(0), method = "tukey"), logical(0)) + expect_equal(outliers(2, method = "tukey"), FALSE) expect_equal( - tukey_algorithm(seq(0.1, 2, by = 0.2)), + outliers(seq(0.1, 2, by = 0.2), method = "tukey"), c(TRUE, rep(FALSE, 8), TRUE) ) - expect_equal(tukey_algorithm(c(NA, 1, 2, 3)), c(NA, TRUE, FALSE, TRUE)) + expect_equal( + outliers(c(NA, 1, 2, 3), method = "tukey"), + c(NA, TRUE, FALSE, TRUE) + ) }) test_that("outliers work with NAs", { - expect_identical(resistant_fences(x), resistant_fences(c(NA, x))[-1]) - expect_identical(quartile_method(x), quartile_method(c(NA, x))[-1]) - expect_identical(robust_z(x), robust_z(c(NA, x))[-1]) - expect_identical(tukey_algorithm(x), tukey_algorithm(c(NA, x))[-1]) - expect_identical(kimber_method(x), kimber_method(c(NA, x))[-1]) + expect_identical( + outliers(x, method = "resistant-fences"), + outliers(c(NA, x), method = "resistant-fences")[-1] + ) + expect_identical( + outliers(x, method = "quartile"), + outliers(c(NA, x), method = "quartile")[-1] + ) + expect_identical( + outliers(x, method = "robust-z"), + outliers(c(NA, x), method = "robust-z")[-1] + ) + expect_identical( + outliers(x, method = "tukey"), + outliers(c(NA, x), method = "tukey")[-1] + ) + expect_identical( + outliers(x, method = "kimber"), + outliers(c(NA, x), method = "kimber")[-1] + ) }) test_that("hb transform works", { @@ -48,7 +71,7 @@ test_that("hb transform works", { }) test_that("recycling gives an error", { - expect_error(quartile_method(x, cl = rep(2.5, 10))) - expect_error(quartile_method(x, cu = rep(2.5, 0))) - expect_error(quartile_method(x, a = rep(0, 11))) + expect_error(outliers(x, lower = rep(2.5, 10))) + expect_error(outliers(x, upper = rep(2.5, 0))) + expect_error(outliers(x, scale = rep(0, 11))) })