gap.boxplot与boxplot绘图差异原因及带分位数间隙箱线图实现求助

问题根源拆解 #

你观察得非常准确！gap.boxplot在处理公式形式输入（比如mpg~cyl）时，异常值的分组分配逻辑和base R的boxplot不一致，这就是图形差异的核心原因：

• base的boxplot会精准把每个异常值映射到它所属的原始分组（你输出的group长度是2，对应2个异常值各自的分组）；
• 而gap.boxplot在处理公式时，错误地用bxgap$group <- at把异常值分组直接设为所有分组的索引（你的输出里group长度是3，等于分组总数），导致异常值被分配到错误的分组，最终图形错位。

至于向量/列表形式下两者结果一致，是因为这种输入没有复杂的分组映射逻辑，gap.boxplot不需要处理公式解析，直接按独立组处理，逻辑和base保持一致。

解决你的核心需求：带间隙+5%/95%分位数须的箱线图 #

这里给你两种可行方案，按需选择：

方案1：修改gap.boxplot源码修正分组逻辑

直接复制gap.boxplot的源码，把异常值分组的赋值逻辑改成继承base boxplot的正确分组，就能让公式输入的结果和base一致，再自定义须的分位数：

library(plotrix)

# 重写gap.boxplot，修正分组逻辑
my_gap.boxplot <- function (x, gap = 0.05, width = gap, varwidth = FALSE, boxwex = 0.8, 
          staplewex = 0.5, outwex = 0.5, notch = FALSE, outline = TRUE, 
          names, plot = TRUE, border = par("fg"), col = NULL, log = "", 
          pars = list(boxwex = boxwex, staplewex = staplewex, outwex = outwex), 
          horizontal = FALSE, add = FALSE, at = NULL, stats = NULL, ...) 
{
  if (inherits(x, "formula")) {
    bx <- boxplot(x, plot = FALSE, ...)
    bxgap <- bx
    # 关键修改：用base boxplot的正确分组替代原逻辑
    bxgap$group <- bx$group 
    if (missing(names)) 
      names <- bx$names
    if (is.null(at)) 
      at <- seq_along(bx$stats[1, ])
  }
  else {
    if (is.list(x)) {
      if (missing(names)) 
        names <- names(x)
      if (is.null(names)) 
        names <- paste("Group", seq_along(x))
      x <- lapply(x, function(x) x[is.finite(x)])
      if (is.null(stats)) {
        stats <- t(sapply(x, function(x) {
          stats <- boxplot.stats(x)$stats
          if (length(stats) < 5) 
            c(rep(stats[1], 5 - length(stats)), stats)
          else stats
        }))
        n <- sapply(x, length)
        conf <- t(sapply(x, function(x) {
          conf <- boxplot.stats(x)$conf
          if (length(conf) < 2) 
            rep(conf[1], 2)
          else conf
        }))
        out <- unlist(lapply(x, function(x) boxplot.stats(x)$out))
        group <- rep(seq_along(x), sapply(x, function(x) length(boxplot.stats(x)$out)))
      }
      else {
        if (nrow(stats) != length(x)) 
          stop("Stats must have one row per group")
        n <- sapply(x, length)
        conf <- matrix(NA, nrow = length(x), ncol = 2)
        out <- NULL
        group <- NULL
      }
      bxgap <- list(stats = stats, n = n, conf = conf, out = out, 
                    group = group, names = names)
      if (is.null(at)) 
        at <- seq_along(x)
    }
    else {
      x <- x[is.finite(x)]
      if (is.null(stats)) {
        bx <- boxplot.stats(x)
        stats <- bx$stats
        n <- length(x)
        conf <- bx$conf
        out <- bx$out
        group <- rep(1, length(out))
      }
      else {
        if (length(stats) != 5) 
          stop("Stats must be a vector of length 5")
        n <- length(x)
        conf <- c(NA, NA)
        out <- NULL
        group <- NULL
      }
      bxgap <- list(stats = matrix(stats, nrow = 5), n = n, 
                    conf = matrix(conf, nrow = 2), out = out, group = group, 
                    names = deparse(substitute(x)))
      if (is.null(at)) 
        at <- 1
    }
  }
  if (plot) {
    if (!add) {
      xlim <- if (horizontal) 
        range(bxgap$stats, bxgap$out, na.rm = TRUE)
      else range(at)
      ylim <- if (horizontal) 
        range(at)
      else range(bxgap$stats, bxgap$out, na.rm = TRUE)
      if (log != "") {
        if (grepl("x", log)) 
          xlim <- xlim[xlim > 0]
        if (grepl("y", log)) 
          ylim <- ylim[ylim > 0]
      }
      plot(0, 0, type = "n", xlim = xlim, ylim = ylim, xlab = "", 
           ylab = "", log = log, axes = FALSE, ...)
      axis(if (horizontal) 
        1
        else 2)
      axis(if (horizontal) 
        2
        else 1, at = at, labels = names)
      box()
    }
    box_col <- if (is.null(col)) 
      rep(par("bg"), length(at))
    else rep(col, length(at))
    if (varwidth) {
      box_widths <- sqrt(bxgap$n)/max(sqrt(bxgap$n)) * boxwex
    }
    else {
      box_widths <- rep(boxwex, length(at))
    }
    for (i in seq_along(at)) {
      ymin <- bxgap$stats[1, i]
      y25 <- bxgap$stats[2, i]
      ymid <- bxgap$stats[3, i]
      y75 <- bxgap$stats[4, i]
      ymax <- bxgap$stats[5, i]
      if (horizontal) {
        rect(ymin, at[i] - box_widths[i]/2, y25, at[i] + box_widths[i]/2, 
             col = box_col[i], border = border)
        rect(y25, at[i] - box_widths[i]/2, y75, at[i] + box_widths[i]/2, 
             col = box_col[i], border = border)
        rect(y75, at[i] - box_widths[i]/2, ymax, at[i] + box_widths[i]/2, 
             col = box_col[i], border = border)
        lines(c(ymid, ymid), c(at[i] - box_widths[i]/2, at[i] + box_widths[i]/2), 
              col = border)
        if (notch) {
          notch_left <- bxgap$conf[1, i]
          notch_right <- bxgap$conf[2, i]
          lines(c(notch_left, y25, y25, notch_left), c(at[i] - box_widths[i]/2, 
                                                        at[i] - box_widths[i]/2, at[i] + box_widths[i]/2, 
                                                        at[i] + box_widths[i]/2), col = border)
          lines(c(notch_right, y75, y75, notch_right), c(at[i] - box_widths[i]/2, 
                                                         at[i] - box_widths[i]/2, at[i] + box_widths[i]/2, 
                                                         at[i] + box_widths[i]/2), col = border)
        }
        if (outline && length(bxgap$out) > 0) {
          out_x <- bxgap$out[bxgap$group == i]
          out_y <- rep(at[i], length(out_x))
          if (length(out_x) > 0) {
            gap_low <- ymin - gap * (ymax - ymin)
            gap_high <- ymax + gap * (ymax - ymin)
            out_x_low <- out_x[out_x < gap_low]
            out_x_high <- out_x[out_x > gap_high]
            if (length(out_x_low) > 0) {
              lines(c(gap_low, gap_low - width * (ymax - ymin)), 
                    c(at[i], at[i]), col = border)
              points(out_x_low, out_y[out_x < gap_low], pch = 19, 
                     cex = outwex, col = border)
            }
            if (length(out_x_high) > 0) {
              lines(c(gap_high, gap_high + width * (ymax - ymin)), 
                    c(at[i], at[i]), col = border)
              points(out_x_high, out_y[out_x > gap_high], pch = 19, 
                     cex = outwex, col = border)
            }
          }
        }
      }
      else {
        rect(at[i] - box_widths[i]/2, ymin, at[i] + box_widths[i]/2, 
             y25, col = box_col[i], border = border)
        rect(at[i] - box_widths[i]/2, y25, at[i] + box_widths[i]/2, 
             y75, col = box_col[i], border = border)
        rect(at[i] - box_widths[i]/2, y75, at[i] + box_widths[i]/2, 
             ymax, col = box_col[i], border = border)
        lines(c(at[i] - box_widths[i]/2, at[i] + box_widths[i]/2), 
              c(ymid, ymid), col = border)
        if (notch) {
          notch_bottom <- bxgap$conf[1, i]
          notch_top <- bxgap$conf[2, i]
          lines(c(at[i] - box_widths[i]/2, at[i] - box_widths[i]/2, 
                  at[i] + box_widths[i]/2, at[i] + box_widths[i]/2), c(notch_bottom, 
                                                                       y25, y25, notch_bottom), col = border)
          lines(c(at[i] - box_widths[i]/2, at[i] - box_widths[i]/2, 
                  at[i] + box_widths[i]/2, at[i] + box_widths[i]/2), c(notch_top, 
                                                                       y75, y75, notch_top), col = border)
        }
        if (outline && length(bxgap$out) > 0) {
          out_y <- bxgap$out[bxgap$group == i]
          out_x <- rep(at[i], length(out_y))
          if (length(out_y) > 0) {
            gap_low <- ymin - gap * (ymax - ymin)
            gap_high <- ymax + gap * (ymax - ymin)
            out_y_low <- out_y[out_y < gap_low]
            out_y_high <- out_y[out_y > gap_high]
            if (length(out_y_low) > 0) {
              lines(c(at[i], at[i]), c(gap_low, gap_low - width * (ymax - ymin)), 
                    col = border)
              points(out_x[out_y < gap_low], out_y_low, pch = 19, 
                     cex = outwex, col = border)
            }
            if (length(out_y_high) > 0) {
              lines(c(at[i], at[i]), c(gap_high, gap_high + width * (ymax - ymin)), 
                    col = border)
              points(out_x[out_y > gap_high], out_y_high, pch = 19, 
                     cex = outwex, col = border)
            }
          }
        }
      }
    }
    invisible(bxgap)
  }
  else {
    invisible(bxgap)
  }
}

# 自定义统计量函数：返回5%、25%、中位数、75%、95%分位数
custom_stats <- function(x) {
  quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95), na.rm = TRUE)
}

# 用修改后的函数绘制目标图形
my_gap.boxplot(mpg~cyl, data=mtcars, stats=tapply(mtcars$mpg, mtcars$cyl, custom_stats))

方案2：手动整理数据为列表形式

避开公式输入的坑，直接把数据按分组整理成列表，再传入自定义统计量和异常值分组：

library(plotrix)

# 按cyl分组整理数据
grouped_mpg <- list(
  "4" = mtcars$mpg[mtcars$cyl == 4],
  "6" = mtcars$mpg[mtcars$cyl == 6],
  "8" = mtcars$mpg[mtcars$cyl == 8]
)

# 计算每个组的5%/95%分位数统计量
stats_matrix <- t(sapply(grouped_mpg, function(x) {
  quantile(x, c(0.05, 0.25, 0.5, 0.75, 0.95), na.rm = TRUE)
}))

# 提取异常值和对应的分组
outliers <- c()
out_groups <- c()
for (i in 1:length(grouped_mpg)) {
  vals <- grouped_mpg[[i]]
  low <- stats_matrix[i,1]
  high <- stats_matrix[i,5]
  outs <- vals[vals < low | vals > high]
  outliers <- c(outliers, outs)
  out_groups <- c(out_groups, rep(i, length(outs)))
}

# 绘制带间隙的箱线图
gap.boxplot(
  x = grouped_mpg,
  stats = stats_matrix,
  out = outliers,
  group = out_groups,
  gap = 0.5,
  main = "Boxplot with 5%/95% Whiskers and Outlier Gaps"
)

内容的提问来源于stack exchange，提问作者David