R/ellipse.R

In abnormally-distributed/cvreg: Cross Validation and Robust Estimation Utilities

#' Create coordinates for a tolerance ellipse
#'
#' @description This takes a pair of variables and creates the coordinates needed to plot
#' a tolerance ellipse.
#'
#' @param data a data frame.
#' @param vars a character vector naming the two variables of interest.
#' @param cov.type a covariance function, preferably from the cvreg, ICS, ICSNP, MNM, or rrcov packages.
#' @param cov.args a list with named arguments to pass to cov.type.
#' @param segments The number of segments to be used in drawing the ellipse.
#' @param level The coverage level at which to draw an ellipse. The default is 0.95 for a 95 pct coverage of the data points. Must be greater than 0.50.
#' @export
#' @return a data frame consisting of coordinates.
#'
create_ellipse <- function(data, vars, cov.type = cov.ogk, cov.args = list(), level = 0.95, segments = 180){
  dfn <- 2
  dfd <- nrow(data) - 1
  if (dfd < 3) {
    message("Too few points to calculate an ellipse")
    ellipse <- rbind(as.numeric(c(NA, NA)))
  }
  else {
    data <- data[,vars]
    cov.args[["x"]] <- data
    v <- do.call(cov.type, cov.args)
    if (!is.null(v$cov)){
      shape <- v$cov
    }
    else if (!is.null(v$scatter)){
      shape <- v$scatter
    }
    else{
      return(cat(crayon::red("Please use a covariance method that supplies a list with an element named either 'cov' or 'scatter'")))
    }

    if (!is.null(v$center)){
      center <- v$center
    }
    else if (!is.null(v$location)){
      center <- v$location
    }
    else{
      return(cat(crayon::red("Please use a covariance method that supplies a list with an element named either 'location' or 'center'")))
    }
    chol_decomp <- chol(shape)
    radius <- sqrt(dfn * stats::qf(level, dfn, dfd))
    angles <- (0:segments) * 2 * pi/segments
    unit.circle <- cbind(cos(angles), sin(angles))
    ellipse <- t(center + radius * t(unit.circle %*% chol_decomp))
  }

  colnames(ellipse) <- vars
  as.data.frame(ellipse)
}

#' Compute robust data ellipses
#'
#' @description This is a modification of  \link[ggplot2]{stat_ellipse} to utilize a variety
#' of robust estimators. This can be very useful for spotting bivariate outliers.
#'
#' @param cov.type a covariance function, preferably from the cvreg, ICS, ICSNP, MNM, or rrcov packages.
#' @param cov.args a list with named arguments to pass to cov.type.
#' @param level The coverage level at which to draw an ellipse. The default is 0.95 for a 95 pct coverage of the data points. Must be greater than 0.50.
#' @param segments The number of segments to be used in drawing the ellipse.
#' @inheritParams ggplot2::layer
#' @inheritParams ggplot2::geom_point
#' @export
#' @examples
#' ggplot(faithful, aes(waiting, eruptions, color = eruptions > 3)) +
#'   geom_point() +
#'   stat_ellipse2(type = "mcd", linetype = 2) +
#'   stat_ellipse2(type = "mgv")
#'
#'
stat_ellipse2 <- function(mapping = NULL, data = NULL,
                          geom = "path", position = "identity",
                          cov.type = cov.ogk,
                          cov.args = list(),
                          level = 0.95,
                          segments = 90,
                          na.rm = TRUE,
                          show.legend = NA,
                          inherit.aes = TRUE,
                          ...) {

  layer(
    data = data,
    mapping = mapping,
    stat = StatEllipse2,
    geom = geom,
    position = position,
    show.legend = show.legend,
    inherit.aes = inherit.aes,
    params = list(
      cov.type = cov.type,
      cov.args = cov.args,
      level = level,
      segments = segments,
      na.rm = na.rm,
      ...
    )
  )
}

#' @rdname ggplot2-ggproto
#' @title StatEllipse2
#' @format NULL
#' @usage NULL
#' @export
StatEllipse2 <- ggproto("StatEllipse2", Stat, required_aes = c("x", "y"),

                        compute_group = function(data, scales, cov.type, cov.args, level = 0.95, segments = 90, na.rm = T) {

                          create_ellipse(data = data, vars = c("x", "y"), cov.type = cov.type,
                                         cov.args = cov.args, level = level, segments = segments)
                        }
)