用R语言绘制y轴不连续的柱形图

26,671

R作图包plotrix提供了不连续y轴（或者称断裂y轴）图形的绘制，原barplot函数的beside参数都不能用，图形也不怎么如意：

  1. library(plotrix)
  2. par(mar=c(3,3,1,1))
  3. par(mgp=c(2,0.5,0))
  4. y1 <- c(75, 130, 4, 3, 5, 10, 100, 1, 150, 110)
  5. y2 <- c(60, 120, 3, 8, 6, 12, 100, 2, 180, 90)
  6. plotrix::gap.barplot(rbind(y1,y2), gap=c(15,50), beside=TRUE, ylab="Level", xlab="Sample")
  7. ## Warning: "beside" is not a graphical parameter
  8. ## Warning: "beside" is not a graphical parameter
  9. ## Warning: "beside" is not a graphical parameter
  10. ## Warning: "beside" is not a graphical parameter
  11. ## Warning: "beside" is not a graphical parameter

用R语言绘制y轴不连续的柱形图-图片1

下面是使用自编函数gap.barplot（代码在后面）绘制的图形。函数可以手动设置断点，也可以由函数自动计算。断点位置的符号表示提供了平行线和zigzag两种，并且可设置背景颜色、大小、线型、平行线旋转角度等。参数使用方法请参看函数说明。

  1. datax <- na.omit(airquality)[,1:4]
  2. cols <- terrain.colors(ncol(datax) - 1)
  3. layout(matrix(1:4, ncol=2))
  4. set.seed(0)
  5. for (ndx in 1:4){
  6.     dt <- datax[sample(rownames(datax), 10), ]
  7.     dt <- cbind(dt, dt[, -1]* 0.1)
  8.     par(mar=c(1, 3, 0.5, 0.5))
  9.     brkt <- sample(c('normal', 'zigzag'), 1)
  10.     gap.barplot(dt, y.cols=2:4, sd.cols=5:7, col=cols, brk.type=brkt, max.fold=5, ratio=2, cex.error=0.3)
  11. }

用R语言绘制y轴不连续的柱形图-图片2

  1. #' 使用R基本绘图函数绘制y轴不连续的柱形图
  2. #'
  3. #' 绘制y轴不连续的柱形图，具有误差线添加功能。断点位置通过btm和top参数设置，如果不设置，函数可自动计算合适的断点位置。
  4. #' @title gap.barplot function
  5. #' @param df 长格式的data.frame，即数据框中每一列为一组绘图数据。
  6. #' @param y.cols 用做柱形图y值的数据列（序号或名称），一列为一组。
  7. #' @param sd.cols 与y值列顺序对应的误差值的数据列（序号或名称）。
  8. #' @param btm 低位断点。如果btm和top均不设置，程序将自动计算和设置断点位置。
  9. #' @param top 高位断点。
  10. #' @param min.range 自动计算断点的阈值：最大值与最小值的最小比值
  11. #' @param max.fold 自动计算断点时最大值与下方数据最大值的最大倍数比
  12. #' @param ratio 断裂后上部与下部y轴长度的比例。
  13. #' @param gap.width y轴断裂位置的相对物理宽度（非坐标轴实际刻度）
  14. #' @param brk.type 断点类型，可设为normal或zigzag
  15. #' @param brk.bg 断点处的背景颜色
  16. #' @param brk.srt 断点标记线旋转角度
  17. #' @param brk.size 断点标记线的大小（长度）
  18. #' @param brk.col 断点标记线的颜色
  19. #' @param brk.lwd 断点标记线的线宽
  20. #' @param cex.error 误差线相对长度，默认为1
  21. #' @param ... 其他传递给R基本绘图函数barplot的参数
  22. #' @return 返回barplot的原始返回值，即柱形图的x坐标
  23. #' @examples
  24. #' datax <- na.omit(airquality)[,1:4]
  25. #' cols <- cm.colors(ncol(datax))
  26. #' layout(matrix(1:6, ncol=2))
  27. #' set.seed(0)
  28. #' for (ndx in 1:6){
  29. #'     dt <- datax[sample(rownames(datax), 10), ]
  30. #'     par(mar=c(0.5,2,0.5,0.5))
  31. #'     brkt <- sample(c('normal', 'zigzag'), 1)
  32. #'     gap.barplot(dt, col=cols, brk.type=brkt, max.fold=5, ratio=2)
  33. #' }
  34. #' @author ZG Zhao
  35. #' @export
  36. gap.barplot <- function(df, y.cols=1:ncol(df), sd.cols=NULL, btm=NULL, top=NULL, min.range=10, max.fold=5, ratio=1, gap.width=1,
  37.                         brk.type='normal', brk.bg='white', brk.srt=135, brk.size=1, brk.col='black', brk.lwd=1, cex.error=1, ...){
  38.     if (missing(df)) stop('No data provided.')
  39.     if (is.numeric(y.cols)) ycol <- y.cols else ycol <- colnames(df)==y.cols
  40.     if (!is.null(sd.cols))
  41.         if (is.numeric(sd.cols)) scol <- sd.cols else scol <- colnames(df)==sd.cols
  42.     ## Arrange data
  43.     opts <- options()
  44.     options(warn=-1)
  45.     y <- t(df[, ycol])
  46.     colnames(y) <- NULL
  47.     if(missing(sd.cols)) sdx <- 0 else sdx <- t(df[, scol])
  48.     sdu <- y + sdx
  49.     sdd <- y - sdx
  50.     ylim <- c(0, max(sdu) * 1.05)
  51.     ## 如果没有设置btm或top，自动计算
  52.     if (is.null(btm) | is.null(top)){
  53.         autox <- .auto.breaks(dt=sdu, min.range=min.range, max.fold=max.fold)
  54.         if (autox$flag){
  55.             btm <- autox$btm
  56.             top <- autox$top
  57.         } else {
  58.             xx <- barplot(y, beside=TRUE, ylim=ylim, ...)
  59.             if (!missing(sd.cols)) errorbar(xx, y, sdu - y, horiz=FALSE, cex=cex.error)
  60.             box()
  61.             return(invisible(xx))
  62.         }
  63.     }
  64.     ## Set up virtual y limits
  65.     halflen <- btm - ylim[1]
  66.     xlen <- halflen * 0.1 * gap.width
  67.     v_tps1 <- btm + xlen        # virtual top positions
  68.     v_tps2 <- v_tps1 + halflen * ratio
  69.     v_ylim <- c(ylim[1], v_tps2)
  70.     r_tps1 <- top               # real top positions
  71.     r_tps2 <- ylim[2]
  72.     ## Rescale data
  73.     lmx <- summary(lm(c(v_tps1, v_tps2)~c(r_tps1, r_tps2)))
  74.     lmx <- lmx$coefficients
  75.     sel1 <- y > top
  76.     sel2 <- y >=btm & y <=top
  77.     y[sel1] <- y[sel1] * lmx[2] + lmx[1]
  78.     y[sel2] <- btm + xlen/2
  79.     sel1 <- sdd > top
  80.     sel2 <- sdd >=btm & sdd <=top
  81.     sdd[sel1] <- sdd[sel1] * lmx[2] + lmx[1]
  82.     sdd[sel2] <- btm + xlen/2
  83.     sel1 <- sdu > top
  84.     sel2 <- sdu >=btm & sdu <=top
  85.     sdu[sel1] <- sdu[sel1] * lmx[2] + lmx[1]
  86.     sdu[sel2] <- btm + xlen/2
  87.     ## bar plot
  88.     xx <- barplot(y, beside=TRUE, ylim=v_ylim, axes = FALSE, names.arg=NULL, ...)
  89.     ## error bars
  90.     if(!missing(sd.cols)) errorbar(xx, y, sdu - y, horiz=FALSE, cex=cex.error)
  91.     ## Real ticks and labels    
  92.     brks1 <- pretty(seq(0, btm, length=10), n=4)
  93.     brks1 <- brks1[brks1 >= 0 & brks1 < btm]
  94.     brks2 <- pretty(seq(top, r_tps2, length=10), n=4)
  95.     brks2 <- brks2[brks2 > top & brks2 <= r_tps2]
  96.     labx <- c(brks1, brks2)
  97.     ## Virtual ticks
  98.     brks <- c(brks1, brks2 * lmx[2] + lmx[1])
  99.     axis(2, at=brks, labels=labx)
  100.     box()
  101.     ## break marks
  102.     pos <- par("usr")
  103.     xyratio <- (pos[2] - pos[1])/(pos[4] - pos[3])
  104.     xlen <- (pos[2] - pos[1])/50 * brk.size
  105.     px1 <- pos[1] - xlen
  106.     px2 <- pos[1] + xlen
  107.     px3 <- pos[2] - xlen
  108.     px4 <- pos[2] + xlen
  109.     py1 <- btm
  110.     py2 <- v_tps1
  111.     rect(px1, py1, px4, py2, col=brk.bg, xpd=TRUE, border=brk.bg)
  112.     x1 <- c(px1, px1, px3, px3)
  113.     x2 <- c(px2, px2, px4, px4)
  114.     y1 <- c(py1, py2, py1, py2)
  115.     y2 <- c(py1, py2, py1, py2)
  116.     px <- .xy.adjust(x1, x2, y1, y2, xlen, xyratio, angle=brk.srt*pi/90)
  117.     if (brk.type=='zigzag'){
  118.         x1 <- c(x1, px1, px3)
  119.         x2 <- c(x2, px2, px4)
  120.         if (brk.srt > 90){
  121.             y1 <- c(y1, py2, py2)
  122.             y2 <- c(y2, py1, py1)
  123.         } else {
  124.             y1 <- c(y1, py1, py1)
  125.             y2 <- c(y2, py2, py2)
  126.         }
  127.     }
  128.     if (brk.type=='zigzag') {
  129.         px$x1 <- c(pos[1], px2, px1, pos[2], px4, px3)
  130.         px$x2 <- c(px2, px1, pos[1], px4, px3, pos[2])
  131.         mm <- (v_tps1 - btm)/3
  132.         px$y1 <- rep(c(v_tps1, v_tps1 - mm, v_tps1 - 2 * mm), 2)
  133.         px$y2 <- rep(c(v_tps1 - mm, v_tps1 - 2 * mm, btm), 2)
  134.     }
  135.     par(xpd=TRUE)
  136.     segments(px$x1, px$y1, px$x2, px$y2, lty=1, col=brk.col, lwd=brk.lwd)
  137.     options(opts)
  138.     par(xpd=FALSE)
  139.     invisible(xx)
  140. }
  141. ## 绘制误差线的函数
  142. errorbar <- function(x, y, sd.lwr, sd.upr, horiz=FALSE, cex=1, ...)
  143. {
  144.     if(missing(sd.lwr) & missing(sd.upr)) return(NULL)
  145.     if(missing(sd.upr)) sd.upr <- sd.lwr
  146.     if(missing(sd.lwr)) sd.lwr <- sd.upr
  147.     if(!horiz){
  148.         arrows(x, y, y1=y-sd.lwr, length=0.1*cex, angle=90, ...)
  149.         arrows(x, y, y1=y+sd.upr, length=0.1*cex, angle=90, ...)
  150.     } else{
  151.         arrows(y, x, x1=y-sd.lwr, length=0.1*cex, angle=90, ...)
  152.         arrows(y, x, x1=y+sd.upr, length=0.1*cex, angle=90, ...)
  153.     }
  154. }
  155. .xy.adjust <- function(x1, x2, y1, y2, xlen, xyratio, angle){
  156.     xx1 <- x1 - xlen * cos(angle)
  157.     yy1 <- y1 + xlen * sin(angle)/xyratio
  158.     xx2 <- x2 + xlen * cos(angle)
  159.     yy2 <- y2 - xlen * sin(angle)/xyratio
  160.     return(list(x1=xx1, x2=xx2, y1=yy1, y2=yy2))
  161. }
  162. ## 自动计算断点位置的函数
  163. .auto.breaks <- function(dt, min.range, max.fold){
  164.     datax <- sort(as.vector(dt))
  165.     flags <- FALSE
  166.     btm <- top <- NULL
  167.     if (max(datax)/min(datax) < min.range) return(list(flag=flags, btm=btm, top=top))
  168.     m <- max(datax)
  169.     btm <- datax[2]
  170.     i <- 3
  171.     while(m/datax[i] > max.fold){
  172.         btm <- datax[i]
  173.         flags <- TRUE
  174.         i <- i + 1
  175.     }
  176.     if (flags) {
  177.         btm <- btm + 0.05 * btm
  178.         x <- 2
  179.         top <- datax[i] * (x - 1)/x
  180.         while (top < btm) {
  181.             x <- x + 1
  182.             top <- datax[i] * (x - 1)/x
  183.             if (x > 100) {
  184.                 flags <- FALSE
  185.                 break
  186.             }
  187.         }
  188.     }
  189.     return(list(flag=flags, btm=btm, top=top))
  190. }

原文来自：http://blog.csdn.net/u014801157/article/details/24372371

评论  2  访客  1
    • 安在便好 0
      1F
      回复

      这个函数怎么可以加横坐标各组名称

    发表评论

    匿名网友

    拖动滑块以完成验证
    加载失败