1#' Compute normal data ellipses 2#' 3#' The method for calculating the ellipses has been modified from 4#' `car::dataEllipse` (Fox and Weisberg, 2011) 5#' 6#' @references John Fox and Sanford Weisberg (2011). An \R Companion to 7#' Applied Regression, Second Edition. Thousand Oaks CA: Sage. URL: 8#' \url{https://socialsciences.mcmaster.ca/jfox/Books/Companion/} 9#' @param level The level at which to draw an ellipse, 10#' or, if `type="euclid"`, the radius of the circle to be drawn. 11#' @param type The type of ellipse. 12#' The default `"t"` assumes a multivariate t-distribution, and 13#' `"norm"` assumes a multivariate normal distribution. 14#' `"euclid"` draws a circle with the radius equal to `level`, 15#' representing the euclidean distance from the center. 16#' This ellipse probably won't appear circular unless `coord_fixed()` is applied. 17#' @param segments The number of segments to be used in drawing the ellipse. 18#' @inheritParams layer 19#' @inheritParams geom_point 20#' @export 21#' @examples 22#' ggplot(faithful, aes(waiting, eruptions)) + 23#' geom_point() + 24#' stat_ellipse() 25#' 26#' ggplot(faithful, aes(waiting, eruptions, color = eruptions > 3)) + 27#' geom_point() + 28#' stat_ellipse() 29#' 30#' ggplot(faithful, aes(waiting, eruptions, color = eruptions > 3)) + 31#' geom_point() + 32#' stat_ellipse(type = "norm", linetype = 2) + 33#' stat_ellipse(type = "t") 34#' 35#' ggplot(faithful, aes(waiting, eruptions, color = eruptions > 3)) + 36#' geom_point() + 37#' stat_ellipse(type = "norm", linetype = 2) + 38#' stat_ellipse(type = "euclid", level = 3) + 39#' coord_fixed() 40#' 41#' ggplot(faithful, aes(waiting, eruptions, fill = eruptions > 3)) + 42#' stat_ellipse(geom = "polygon") 43stat_ellipse <- function(mapping = NULL, data = NULL, 44 geom = "path", position = "identity", 45 ..., 46 type = "t", 47 level = 0.95, 48 segments = 51, 49 na.rm = FALSE, 50 show.legend = NA, 51 inherit.aes = TRUE) { 52 layer( 53 data = data, 54 mapping = mapping, 55 stat = StatEllipse, 56 geom = geom, 57 position = position, 58 show.legend = show.legend, 59 inherit.aes = inherit.aes, 60 params = list( 61 type = type, 62 level = level, 63 segments = segments, 64 na.rm = na.rm, 65 ... 66 ) 67 ) 68} 69 70#' @rdname ggplot2-ggproto 71#' @format NULL 72#' @usage NULL 73#' @export 74StatEllipse <- ggproto("StatEllipse", Stat, 75 required_aes = c("x", "y"), 76 77 compute_group = function(data, scales, type = "t", level = 0.95, 78 segments = 51, na.rm = FALSE) { 79 calculate_ellipse(data = data, vars = c("x", "y"), type = type, 80 level = level, segments = segments) 81 } 82) 83 84calculate_ellipse <- function(data, vars, type, level, segments){ 85 dfn <- 2 86 dfd <- nrow(data) - 1 87 88 if (!type %in% c("t", "norm", "euclid")) { 89 message("Unrecognized ellipse type") 90 ellipse <- rbind(as.numeric(c(NA, NA))) 91 } else if (dfd < 3) { 92 message("Too few points to calculate an ellipse") 93 ellipse <- rbind(as.numeric(c(NA, NA))) 94 } else { 95 if (type == "t") { 96 v <- MASS::cov.trob(data[,vars]) 97 } else if (type == "norm") { 98 v <- stats::cov.wt(data[,vars]) 99 } else if (type == "euclid") { 100 v <- stats::cov.wt(data[,vars]) 101 v$cov <- diag(rep(min(diag(v$cov)), 2)) 102 } 103 shape <- v$cov 104 center <- v$center 105 chol_decomp <- chol(shape) 106 if (type == "euclid") { 107 radius <- level/max(chol_decomp) 108 } else { 109 radius <- sqrt(dfn * stats::qf(level, dfn, dfd)) 110 } 111 angles <- (0:segments) * 2 * pi/segments 112 unit.circle <- cbind(cos(angles), sin(angles)) 113 ellipse <- t(center + radius * t(unit.circle %*% chol_decomp)) 114 } 115 116 colnames(ellipse) <- vars 117 mat_2_df(ellipse) 118} 119