1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
|
#' Plot data ellipses.
#'
#' @param level The confidence level at which to draw an ellipse (default is 0.95),
#' or, if \code{type="euclid"}, the radius of the circle to be drawn.
#' @param type The type of ellipse.
#' The default \code{"t"} assumes a multivariate t-distribution, and
#' \code{"norm"} assumes a multivariate normal distribution.
#' \code{"euclid"} draws a circle with the radius equal to \code{level},
#' representing the euclidian distance from the center.
#' This ellipse probably won't appear circular unless \code{coord_fixed()} is applied.
#' @param segments The number of segments to be used in drawing the ellipse.
#' @param na.rm If \code{FALSE} (the default), removes missing values with
#' a warning. If \code{TRUE} silently removes missing values.
#' @inheritParams stat_identity
#'
#' @details The method for calculating the ellipses has been modified from car::ellipse (Fox and Weisberg, 2011)
#'
#' @references
#' John Fox and Sanford Weisberg (2011). An {R} Companion to Applied Regression, Second Edition. Thousand Oaks CA: Sage. URL: http://socserv.socsci.mcmaster.ca/jfox/Books/Companion
#'
#' @export
#' @importFrom MASS cov.trob
#'
#' @examples
#' ggplot(faithful, aes(waiting, eruptions))+
#' geom_point()+
#' stat_ellipse()
#'
#' ggplot(faithful, aes(waiting, eruptions, color = eruptions > 3))+
#' geom_point()+
#' stat_ellipse()
#'
#' ggplot(faithful, aes(waiting, eruptions, color = eruptions > 3))+
#' geom_point()+
#' stat_ellipse(type = "norm", linetype = 2)+
#' stat_ellipse(type = "t")
#'
#' ggplot(faithful, aes(waiting, eruptions, color = eruptions > 3))+
#' geom_point()+
#' stat_ellipse(type = "norm", linetype = 2)+
#' stat_ellipse(type = "euclid", level = 3)+
#' coord_fixed()
#'
#' ggplot(faithful, aes(waiting, eruptions, color = eruptions > 3))+
#' stat_ellipse(geom = "polygon")
stat_ellipse <- function(mapping = NULL, data = NULL, geom = "path", position = "identity", type = "t", level = 0.95, segments = 51, na.rm = FALSE, ...) {
StatEllipse$new(mapping = mapping, data = data, geom = geom, position = position, type = type, level = level, segments = segments, na.rm = na.rm, ...)
}
StatEllipse <- proto(Stat, {
objname <- "ellipse"
required_aes <- c("x", "y")
default_geom <- function(.) GeomPath
calculate_groups <- function(., data, scales, ...){
.super$calculate_groups(., data, scales,...)
}
calculate <- function(., data, scales, type = "t", level = 0.95, segments = 51, na.rm = FALSE, ...){
data <- remove_missing(data, na.rm, vars = c("x","y"), name = "stat_ellipse", finite = TRUE)
ellipse <- calculate_ellipse(data=data, vars= c("x","y"), type=type, level=level, segments=segments)
return(ellipse)
}
})
calculate_ellipse <- function(data, vars, type, level, segments){
dfn <- 2
dfd <- nrow(data) - 1
if (!type %in% c("t", "norm", "euclid")){
message("Unrecognized ellipse type")
ellipse <- rbind(as.numeric(c(NA, NA)))
} else if (dfd < 3){
message("Too few points to calculate an ellipse")
ellipse <- rbind(as.numeric(c(NA, NA)))
} else {
if (type == "t"){
v <- cov.trob(data[,vars])
} else if (type == "norm"){
v <- cov.wt(data[,vars])
} else if (type == "euclid"){
v <- cov.wt(data[,vars])
v$cov <- diag(rep(min(diag(v$cov)), 2))
}
shape <- v$cov
center <- v$center
chol_decomp <- chol(shape)
if (type == "euclid"){
radius <- level/max(chol_decomp)
} else {
radius <- sqrt(dfn * 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))
}
ellipse <- as.data.frame(ellipse)
colnames(ellipse) <- vars
return(ellipse)
}
|
