R/smooth.R

In tf: S3 Classes and Methods for Tidy Functional Data

#' @title Simple smoothing of `tf` objects
#'
#' @description Apply running means or medians, `lowess` or Savitzky-Golay
#'   filtering to smooth functional data. This does nothing for `tfb`-objects,
#'   which should be smoothed by using a smaller basis / stronger penalty.
#'
#' @details `tf_smooth.tfd` overrides/automatically sets some defaults of the
#'   used methods:
#'
#'   - **`lowess`** uses a span parameter of `f` = 0.15 (instead of 0.75)
#'   by default.
#'   - **`rollmean`/`median`** use a window size of `k` = $<$number of
#'   grid points$>$/20 (i.e., the nearest odd integer to that) and sets `fill=
#'   "extend"` (i.e., constant extrapolation to replace missing values at the
#'   extremes of the domain) by default. Use `fill= NA` for `zoo`'s default
#'   behavior of shortening the smoothed series.
#'   - **`savgol`** uses a window size of `k` = $<$number of
#'   grid points$>$/10 (i.e., the nearest odd integer to that).
#'
#' @param x a `tf` object containing functional data.
#' @param method one of `"lowess"` (see [stats::lowess()]), `"rollmean"`,
#'   `"rollmedian"` (see [zoo::rollmean()]) or `"savgol"` (see [pracma::savgol()]).
#' @param verbose give lots of diagnostic messages? Defaults to `TRUE`.
#' @param ... arguments for the respective `method`. See details.
#' @returns a smoothed version of the input. For some methods/options, the
#'   smoothed functions may be shorter than the original ones (at both ends).
#' @export
#' @family tidyfun nonparametric smoothers
tf_smooth <- function(x, ...) {
  UseMethod("tf_smooth")
}

#' @export
#' @rdname tf_smooth
tf_smooth.tfb <- function(x, verbose = TRUE, ...) {
  if (verbose) {
    cli::cli_inform(c(
      `!` = "You called {.fn tf_smooth} on a {.cls tfb} object, not on a {.cls tfd} object --
       just use a smaller basis or stronger penalization.",
      i = "Returning unchanged {.cls tfb} object."
    ))
  }
  x
}

#' @importFrom pracma savgol
#' @rdname tf_smooth
#' @export
#' @examples
#' library(zoo)
#' library(pracma)
#' f <- tf_sparsify(tf_jiggle(tf_rgp(4, 201, nugget = 0.05)))
#' f_lowess <- tf_smooth(f, "lowess")
#' # these methods ignore the distances between arg-values:
#' f_mean <- tf_smooth(f, "rollmean")
#' f_median <- tf_smooth(f, "rollmedian", k = 31)
#' f_sg <- tf_smooth(f, "savgol", fl = 31)
#' layout(t(1:4))
#' plot(f, points = FALSE, main = "original")
#' plot(f_lowess,
#'   points = FALSE, col = "blue", main = "lowess (default,\n span 0.9 in red)"
#' )
#' lines(tf_smooth(f, "lowess", f = 0.9), col = "red", alpha = 0.2)
#' plot(f_mean,
#'   points = FALSE, col = "blue", main = "rolling means &\n medians (red)"
#' )
#' lines(f_median, col = "red", alpha = 0.2) # note constant extrapolation at both ends!
#' plot(f, points = FALSE, main = "original and\n savgol (red)")
#' lines(f_sg, col = "red")
tf_smooth.tfd <- function(
  x,
  method = c("lowess", "rollmean", "rollmedian", "savgol"),
  verbose = TRUE,
  ...
) {
  method <- match.arg(method)
  smoother <- get(method, mode = "function")
  dots <- list(...)
  nas <- is.na(x)
  x_evals <- tf_evaluations(x)[!nas]
  # nocov start
  if (method %in% c("savgol", "rollmean", "rollmedian")) {
    if (verbose && !is_equidist(x)) {
      cli::cli_inform(c(
        x = "Non-equidistant arg-values in {.arg x} ignored by {.val {method}}."
      ))
    }
    if (startsWith(method, "rollm")) {
      if (is.null(dots$k)) {
        dots$k <- ceiling(0.05 * min(tf_count(x)))
        dots$k <- dots$k + !(dots$k %% 2) # make uneven
        if (verbose)
          cli::cli_inform(
            "Using {.code k = {dots$k}} observations for rolling data window."
          )
      }
      if (is.null(dots$fill)) {
        if (verbose)
          cli::cli_inform(
            "Setting {.code fill = 'extend'} for start/end values."
          )
        dots$fill <- "extend"
      }
    } else if (is.null(dots$fl)) {
      dots$fl <- ceiling(0.15 * min(tf_count(x)))
      dots$fl <- dots$fl + !(dots$fl %% 2) # make uneven
      if (verbose)
        cli::cli_inform(
          "Using {.code fl = {dots$fl}} observations for rolling data window."
        )
    }
    smoothed <- map(
      x_evals,
      \(x) do.call(smoother, append(list(x), dots))
    )
  }
  # nocov end
  if (method == "lowess") {
    if (is.null(dots$f)) {
      dots$f <- 0.15
      if (verbose)
        cli::cli_inform(
          "Using {.code f = {dots$f}} as smoother span for {.fn lowess}."
        )
    }
    smoothed <- map(
      x_evals,
      \(x) do.call(smoother, append(list(x), dots))$y
    )
  }
  x_smoothed <- vector(length(x), mode = "list")
  x_smoothed[!nas] <- smoothed
  tfd(
    x_smoothed,
    tf_arg(x),
    evaluator = !!attr(x, "evaluator_name"),
    domain = tf_domain(x)
  )
}

#' @export
tf_smooth.default <- function(x, ...) .NotYetImplemented()