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