4

根据下面的帮助，我尝试了这个脚本来绘制带有凸壳的 PCA，但没有成功，知道如何解决它吗？

library(ggbiplot)
library(plyr)

data <-read.csv("C:/Users/AAA.csv")
my.pca <- prcomp(data[,1:9] , scale. = TRUE)


find_hull <- function(my.pca) my.pca[chull(my.pca$x[,1], my.pca$x[,2]), ]
hulls <- ddply(my.pca , "Group", find_hull)

ggbiplot(my.pca, obs.scale = 1, var.scale = 1,groups = data$Group) + 
  scale_color_discrete(name = '') + geom_polygon(data=hulls, alpha=.2) + 
  theme_bw() + theme(legend.direction = 'horizontal', legend.position = 'top')

谢谢。

下面的脚本用省略号绘制 PCA（不推荐使用来自https://github.com/vqv/ggbiplot的稍微修改的示例，因为“opts”已被弃用）

library(ggbiplot)
data(wine)
wine.pca <- prcomp(wine, scale. = TRUE)
g <- ggbiplot(wine.pca, obs.scale = 1, var.scale = 1, 
              groups = wine.class, ellipse = TRUE, circle = TRUE)
g <- g + scale_color_discrete(name = '')
g <- g + theme(legend.direction = 'horizontal', legend.position = 'top')
print(g)

删除椭圆很容易，但我试图用凸包替换它们但没有成功，知道该怎么做吗？

谢谢

4

2 回答 2

4

是的，我们可以为 ggplot 设计一个新的几何图形，然后将它与 ggbiplot 一起使用。这是一个可以做凸包的新geom：

library(ggplot2)
StatBag <- ggproto("Statbag", Stat,
                   compute_group = function(data, scales, prop = 0.5) {

                     #################################
                     #################################
                     # originally from aplpack package, plotting functions removed
                     plothulls_ <- function(x, y, fraction, n.hull = 1,
                                            col.hull, lty.hull, lwd.hull, density=0, ...){
                       # function for data peeling:
                       # x,y : data
                       # fraction.in.inner.hull : max percentage of points within the hull to be drawn
                       # n.hull : number of hulls to be plotted (if there is no fractiion argument)
                       # col.hull, lty.hull, lwd.hull : style of hull line
                       # plotting bits have been removed, BM 160321
                       # pw 130524
                       if(ncol(x) == 2){ y <- x[,2]; x <- x[,1] }
                       n <- length(x)
                       if(!missing(fraction)) { # find special hull
                         n.hull <- 1
                         if(missing(col.hull)) col.hull <- 1
                         if(missing(lty.hull)) lty.hull <- 1
                         if(missing(lwd.hull)) lwd.hull <- 1
                         x.old <- x; y.old <- y
                         idx <- chull(x,y); x.hull <- x[idx]; y.hull <- y[idx]
                         for( i in 1:(length(x)/3)){
                           x <- x[-idx]; y <- y[-idx]
                           if( (length(x)/n) < fraction ){
                             return(cbind(x.hull,y.hull))
                           }
                           idx <- chull(x,y); x.hull <- x[idx]; y.hull <- y[idx];
                         }
                       }
                       if(missing(col.hull)) col.hull <- 1:n.hull
                       if(length(col.hull)) col.hull <- rep(col.hull,n.hull)
                       if(missing(lty.hull)) lty.hull <- 1:n.hull
                       if(length(lty.hull)) lty.hull <- rep(lty.hull,n.hull)
                       if(missing(lwd.hull)) lwd.hull <- 1
                       if(length(lwd.hull)) lwd.hull <- rep(lwd.hull,n.hull)
                       result <- NULL
                       for( i in 1:n.hull){
                         idx <- chull(x,y); x.hull <- x[idx]; y.hull <- y[idx]
                         result <- c(result, list( cbind(x.hull,y.hull) ))
                         x <- x[-idx]; y <- y[-idx]
                         if(0 == length(x)) return(result)
                       }
                       result
                     } # end of definition of plothulls
                     #################################


                     # prepare data to go into function below
                     the_matrix <- matrix(data = c(data$x, data$y), ncol = 2)

                     # get data out of function as df with names
                     setNames(data.frame(plothulls_(the_matrix, fraction = prop)), nm = c("x", "y"))
                     # how can we get the hull and loop vertices passed on also?
                   },

                   required_aes = c("x", "y")
)

#' @inheritParams ggplot2::stat_identity
#' @param prop Proportion of all the points to be included in the bag (default is 0.5)
stat_bag <- function(mapping = NULL, data = NULL, geom = "polygon",
                     position = "identity", na.rm = FALSE, show.legend = NA, 
                     inherit.aes = TRUE, prop = 0.5, alpha = 0.3, ...) {
  layer(
    stat = StatBag, data = data, mapping = mapping, geom = geom, 
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(na.rm = na.rm, prop = prop, alpha = alpha, ...)
  )
}


geom_bag <- function(mapping = NULL, data = NULL,
                     stat = "identity", position = "identity",
                     prop = 0.5, 
                     alpha = 0.3,
                     ...,
                     na.rm = FALSE,
                     show.legend = NA,
                     inherit.aes = TRUE) {
  layer(
    data = data,
    mapping = mapping,
    stat = StatBag,
    geom = GeomBag,
    position = position,
    show.legend = show.legend,
    inherit.aes = inherit.aes,
    params = list(
      na.rm = na.rm,
      alpha = alpha,
      prop = prop,
      ...
    )
  )
}

#' @rdname ggplot2-ggproto
#' @format NULL
#' @usage NULL
#' @export
GeomBag <- ggproto("GeomBag", Geom,
                   draw_group = function(data, panel_scales, coord) {
                     n <- nrow(data)
                     if (n == 1) return(zeroGrob())

                     munched <- coord_munch(coord, data, panel_scales)
                     # Sort by group to make sure that colors, fill, etc. come in same order
                     munched <- munched[order(munched$group), ]

                     # For gpar(), there is one entry per polygon (not one entry per point).
                     # We'll pull the first value from each group, and assume all these values
                     # are the same within each group.
                     first_idx <- !duplicated(munched$group)
                     first_rows <- munched[first_idx, ]

                     ggplot2:::ggname("geom_bag",
                                      grid:::polygonGrob(munched$x, munched$y, default.units = "native",
                                                         id = munched$group,
                                                         gp = grid::gpar(
                                                           col = first_rows$colour,
                                                           fill = alpha(first_rows$fill, first_rows$alpha),
                                                           lwd = first_rows$size * .pt,
                                                           lty = first_rows$linetype
                                                         )
                                      )
                     )


                   },

                   default_aes = aes(colour = "NA", fill = "grey20", size = 0.5, linetype = 1,
                                     alpha = NA, prop = 0.5),

                   handle_na = function(data, params) {
                     data
                   },

                   required_aes = c("x", "y"),

                   draw_key = draw_key_polygon
)

在这里它与 ggbiplot 一起使用，我们设置prop为 1 表示我们要绘制一个包含所有点的多边形：

library(ggbiplot)
data(wine)
wine.pca <- prcomp(wine, scale. = TRUE)
g <- ggbiplot(wine.pca, obs.scale = 1, var.scale = 1, 
              groups = wine.class, ellipse = FALSE, circle = TRUE)
g <- g + scale_color_discrete(name = '')
g <- g + theme(legend.direction = 'horizontal', legend.position = 'top')
g + geom_bag(aes(group = wine.class, fill = wine.class), prop = 1)

在此处输入图像描述

于 2016-03-22T11:29:40.673 回答
3

我们也可以使用 ggbiplot 和一个名为ggpubr的更新的 pkg 来做到这一点：

library(ggpubr)
library(ggbiplot)
data(wine)
wine.pca <- prcomp(wine, scale. = TRUE)

ggbiplot(
  wine.pca,
  obs.scale = 1,
  var.scale = 1,
  groups = wine.class,
  ellipse = FALSE,
  circle = TRUE
) + 
  stat_chull(aes(color = wine.class, 
                 fill = wine.class), 
             alpha = 0.1, 
             geom = "polygon") +
  scale_colour_brewer(palette = "Set1",
                      name = '',
                      guide = 'none') +
  scale_fill_brewer(palette = "Set1",
                      name = '')  +
  theme_minimal()

在此处输入图像描述

我使用scale_colour_brewer和scale_fill_brewer来控制船体和点的颜色，并压制其中一个图例。

为了在多个图中保持相同的颜色，我认为将类别转换为有序因子并确保所有绘制的数据集中都存在因子的每个级别都应该这样做。

于 2018-05-15T06:05:54.287 回答