介绍
看到Y叔为ggord做的添加置信椭圆的geom_ord_ellipse.R（用法见上一篇文章《画个小圈圈》），决定学习一点ggplot图形的语言，对于初学者最好的方法就是照葫芦画瓢，而Y叔的代码自然是最好的模板。我对Y叔的代码进行了逐行的分析，希望以后有需要可以套用。
以下为geom_ord_ellipse.R代码。这个图层的代码其实很短，很简洁，但是如果想要透彻理解还是需要下些功夫的。
##' add confidence ellipse to ordinary plot produced by ggord
##'
##'
##' @title geom_ord_ellipse
##' @param mapping aes mapping
##' @param ellipse_pro confidence value for the ellipse
##' @param fill color to fill the ellipse, NA by default
##' @param ... additional parameters
##' @return ggplot layer
##' @importFrom ggplot2 aes_
##' @importFrom ggplot2 layer
##' @importFrom utils modifyList
##' @export
##' @author Guangchuang Yu
##' @references \url{http://lchblogs.netlify.com/post/2017-12-22-r-addconfellipselda/}
geom_ord_ellipse <- function(mapping = NULL, ellipse_pro = 0.97, fill = NA, ...) {
default_aes <- aes_(color = ~Groups, group = ~Groups)
if (is.null(mapping)) {
mapping <- default_aes
} else {
mapping <- modifyList(default_aes, mapping)
}
layer(
geom = "polygon",
stat = StatOrdEllipse,
mapping = mapping,
position = 'identity',
data = NULL,
params = list(
ellipse_pro = ellipse_pro,
fill = fill,
...
)
)
}
##' @importFrom ggplot2 ggproto
##' @importFrom ggplot2 Stat
##' @importFrom plyr ddply
##' @importFrom grDevices chull
StatOrdEllipse <- ggproto("StatOrdEllipse", Stat,
compute_group = function(self, data, scales, params, ellipse_pro) {
names(data)[1:2] <- c('one', 'two')
theta <- c(seq(-pi, pi, length = 50), seq(pi, -pi, length = 50))
circle <- cbind(cos(theta), sin(theta))
ell <- ddply(data, .(group), function(x) {
if(nrow(x) <= 2) {
return(NULL)
}
sigma <- var(cbind(x$one, x$two))
mu <- c(mean(x$one), mean(x$two))
ed <- sqrt(qchisq(ellipse_pro, df = 2))
data.frame(sweep(circle %*% chol(sigma) * ed, 2, mu, FUN = '+'))
})
names(ell)[2:3] <- c('one', 'two')
ell <- ddply(ell, .(group), function(x) x[chull(x$one, x$two), ])
names(ell) <- c('Groups', 'x', 'y')
return(ell)
},
required_aes = c("x", "y", "group")
)
## . function was from plyr package
. <- function (..., .env = parent.frame()) {
structure(as.list(match.call()[-1]), env = .env, class = "quoted")
}roxygen文档注释
##' add confidence ellipse to ordinary plot produced by ggord
##'
##'
##' @title geom_ord_ellipse
##' @param mapping aes mapping
...
##' @return ggplot layer
##' @importFrom ggplot2 aes_
...
##' @export
##' @author Guangchuang Yu
##' @references \url{http://lchblogs.netlify.com/post/2017-12-22-r-addconfellipselda/}roxygen会根据这一部分对单一函数生成帮助文档，也就是我们在R命令行中输入?FunctionName看到的帮助信息。
其基础格式是（Y叔使用了##'，我觉得好像跟#'没有什么区别？）：
#' @param 函数参数（对应Arguments） 函数的介绍（对应Description）上面代码的注释大多可以顾名思义。比较有意思的是#' @export这行了，roxygen会把这个函数放在NAMESPACE文件中，这样用户便可以调用这个函数。我是可以调用yyplot::geom_ord_ellipse()这个函数的。但yyplot:StatOrdEllipse()这个函数是yyplot的内部函数（注意代码中此函数没有#' @export注释）。如果我调用就会报错：
yyplot::StatOrdEllipse()更多关于roxygen的介绍可以参考这篇文章或官方文档。
ggproto—ggplot2的语言
ggproto是ggplot2模块化、面向对象（Object Oriented）化的核心部分。基础的格式是：
ggproto(`_class` = NULL, `_inherit` = NULL, ...)ggproto是一个很庞大的系统，我目前理解还不是很深入，提供一些参考资料：
官方详细文档
ggproto教程
在理解Y叔这个脚本中我们需要使用最基础的两个模块Geom（创建图层），Stat（数据处理）。
StatOrdEllipse内部函数—ggplot2中的数据处理
我们在作图之前基本都是要对输入数据进行一些数据预处理，比如在做线箱图的时候需要计算中位数、IQR等。在这个脚本中，我们需要做的是计算置信区间椭圆，这一步是通过ggplot::Stat实现的。
StatOrdEllipse <- ggproto("StatOrdEllipse", Stat,
compute_group = function(self, data, scales, params, ellipse_pro) {
## 此处省略...
## 解析见后文
return(ell)
},
required_aes = c("x", "y", "group")
)_class：这个类的名字为StatOrdEllipse_inherit：继承Stat类compute_group：核心处理数据部分，对每一组进行处理，模板为compute_group(self, data, scales, ...)，在这里，它主要接受置信区间（ellipse_pro）参数，返回值为计算好的置信区间轮廓上的点坐标。另外，我认为这里params并不必要。required_aes：创建图层所需要的mapping参数
geom_ord_ellipse函数—创建ggplot2图层
下面便是重头戏，使用上面的Stat来创建一个Geom图层。其实这就是一个普通的函数，只是为了返回一个图层layer。
geom_ord_ellipse <- function(mapping = NULL, ellipse_pro = 0.97, fill = NA, ...) {
default_aes <- aes_(color = ~Groups, group = ~Groups)
if (is.null(mapping)) {
mapping <- default_aes
} else {
mapping <- modifyList(default_aes, mapping)
}
layer(
geom = "polygon",
stat = StatOrdEllipse,
mapping = mapping,
position = 'identity',
data = NULL,
params = list(
ellipse_pro = ellipse_pro,
fill = fill,
...
)
)
}函数的输入值：
mapping：使用过ggplot2都不应该陌生，就是我们写的aes(x=, y=)ellipse_pro：置信区间的概率fill：置信区间内部的颜色，NA会产生透明的椭圆
输入预处理：
这里对输入的mapping做了判断，如果没有输入，就会使用默认值，如果有输入，则替换默认值。
default_aes <- aes_(color = ~Groups, group = ~Groups)：定义了默认的aesthetics，使用aes_时，变量要用双引号引用或使用~（~Groups）modifyList：根据一个list修改另一个list中的值
返回的ggplot::layer：
这个函数的模板是：
layer(geom = NULL, stat = NULL, data = NULL, mapping = NULL,
position = NULL, params = list(), inherit.aes = TRUE,
check.aes = TRUE, check.param = TRUE, subset = NULL, show.legend = NA)对应到Y叔的函数：
geom = "polygon"：一个多边形（椭圆）stat = StatOrdEllipse：使用我们定义的Statmapping = mapping：数据与图形的对应，如x=, y=position = 'identity'：位置的定义data = NULL：从上一图层继承params = list(...)：geom和stat的参数
到这里，对于一些基本的图层，我觉得完全可以套用Y叔的这个模板。在宏基因组公众号中曾经有人问，这个ggord包中能不能把那些向量去掉，或者加入少部分变量的向量，我想用这个模板完全可以实现，只需要用GeomCurve来做个图层就可以了。
计算置信区间椭圆
names(data)[1:2] <- c('one', 'two')
theta <- c(seq(-pi, pi, length = 50), seq(pi, -pi, length = 50))
circle <- cbind(cos(theta), sin(theta))
ell <- ddply(data, .(group), function(x) {
if(nrow(x) <= 2) {
return(NULL)
}
sigma <- var(cbind(x$one, x$two))
mu <- c(mean(x$one), mean(x$two))
ed <- sqrt(qchisq(ellipse_pro, df = 2))
data.frame(sweep(circle %*% chol(sigma) * ed, 2, mu, FUN = '+'))
})
names(ell)[2:3] <- c('one', 'two')
ell <- ddply(ell, .(group), function(x) x[chull(x$one, x$two), ])
names(ell) <- c('Groups', 'x', 'y')这个代码是从ggord的源码改过来的。对于排序图来说，基本最后都会降维到一个低维的空间（2维），方便展示。所以，我们最后需要处理的问题就是根据2维上的散点，计算出这些点分布的可能范围（在2维正态分布的假设下，对协方差使用卡方检验）。其中涉及到我们要把数据分成组（不同椭圆，不同颜色标记），然后对每一组求出上述的范围。这个操作使用ddply实现的：
分而治之
ell <- ddply(data, .(group), function(x) {
## ...
## 解析见后文
})在此，把data根据group拆分成组，然后每一组套用function，最后再把结果组合（rbind）在一起。以下我们介绍每一组是如何处理的。
计算置信椭圆
从统计上来讲这个置信椭圆是这样做的：
做出一个单位圆（半径为单位1）
theta <- c(seq(-pi, pi, length = 50), seq(pi, -pi, length = 50)) ## should be: ## theta <- seq(-pi, pi, length = 50) circle <- cbind(cos(theta), sin(theta))
我们知道单位圆的参数方程为，在这里ggord的作者其实重复了两圈，具体的原因我也没有搞清楚，我认为没有必要。
套用公式把这个圆转换成椭圆
sigma <- var(cbind(x$one, x$two))
mu <- c(mean(x$one), mean(x$two))
ed <- sqrt(qchisq(ellipse_pro, df = 2))
data.frame(sweep(circle %*% chol(sigma) * ed, 2, mu, FUN = '+'))这里的转化实际是：
其中
至于为什么这么做，涉及一些线代推导，可以参见下面两篇文章（特别是第一篇）：
https://jellymatter.wordpress.com/2011/03/31/drawing-confidence-ellipses-and-ellipsoids/#more-507
http://www.visiondummy.com/2014/04/draw-error-ellipse-representing-covariance-matrix/
有点多余？
names(ell)[2:3] <- c('one', 'two')
ell <- ddply(ell, .(group), function(x) x[chull(x$one, x$two), ])这部分的意思是对之前算出的那些椭圆上的点找出对应的凸多边形。我认为对于上述方法得到的点已经在一个凸多边形上了（椭圆）。我没有太明白这一步的必要性，去掉之后并没有发现影响。
作者李陈浩，没错就是那个写了个原型，请我来封装成图层的人，我在他代码（其实来源于ggord包）的基础上，封装了图层，并且修改增强为适应其它ordination plot的通用图层。他对我的代码进行了逐行解析，是学习写ggplot2的好材料，当然这里其实只是写了Stat而已，也就是我之前说的数据处理，画图的Geom使用的是内置的polygon，写图层的能力还应该包括写Geom，以后有机会再跟大家介绍。