ggplot扩展学习笔记--逐行解析Y叔的"geom_ord_ellipse.R"

4 min read R, ggplot2

介绍

看到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()
## Error: 'StatOrdEllipse' is not an exported object from 'namespace:yyplot'

更多关于roxygen的介绍可以参考这篇文章或官方文档。

ggproto–ggplot2的语言

ggproto是ggplot2模块化、面向对象（Object Oriented）化的核心部分。基础的格式是：

ggproto(`_class` = NULL, `_inherit` = NULL, ...)

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：使用我们定义的Stat
  • mapping = 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. 做出一个单位圆（半径为单位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))

我们知道单位圆的参数方程为\(x=\cos(\theta), y=\sin(\theta),\theta\in[0,2\pi]\)，在这里ggord的作者其实重复了两圈，具体的原因我也没有搞清楚，我认为没有必要。

  1. 套用公式把这个圆转换成椭圆
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 = '+'))

这里的转化实际是：

\(\boldsymbol X_{ell}=\boldsymbol\mu + k\boldsymbol X C(\boldsymbol\Sigma)\)

其中\(\boldsymbol\mu=(\bar x_1,\bar x_2)\)，也就是数据的中心，\(k\)根据概率控制椭圆的大小（因为是针对2维正态分布的协方差，所以使用了卡方检验），\(C(\boldsymbol\Sigma)\)代表协方差矩阵的Cholesky分解，\(\boldsymbol X\)和\(\boldsymbol X_{ell}\)分别为圆和椭圆上对应点的坐标。

至于为什么这么做，涉及一些线代推导，可以参见下面两篇文章（特别是第一篇）：

  1. https://jellymatter.wordpress.com/2011/03/31/drawing-confidence-ellipses-and-ellipsoids/#more-507
  2. 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), ])

这部分的意思是对之前算出的那些椭圆上的点找出对应的凸多边形。我认为对于上述方法得到的点已经在一个凸多边形上了（椭圆）。我没有太明白这一步的必要性，去掉之后并没有发现影响。

R

Related

comments powered by Disqus