如何绘制分位数带（在 R 中）答案

【问题标题】：How to plot quantil band (in R)如何绘制分位数带（在 R 中）
【发布时间】：2012-09-21 10:23:58
【问题描述】：

我有一个 CSV 文件，其中包含我感兴趣的每个（Java GC）事件的行。该对象由亚秒级时间戳（非等距）和一些变量组成。对象如下所示：

gcdata <- read.table("http://bernd.eckenfels.net/view/gc1001.ygc.csv",header=TRUE,sep=",", dec=".")
start = as.POSIXct(strptime("2012-01-01 00:00:00", format="%Y-%m-%d %H:%M:%S"))
gcdata.date = gcdata$Timestamp + start
gcdata = gcdata[,2:7] # remove old date col
gcdata=data.frame(date=gcdata.date,gcdata)
str(gcdata)

结果

'data.frame':   2997 obs. of  7 variables:
 $ date           : POSIXct, format: "2012-01-01 00:00:06" "2012-01-01 00:00:06" "2012-01-01 00:00:18" ...
 $ Distance.s.    : num  0 0.165 11.289 9.029 11.161 ...
 $ YGUsedBefore.K.: int  1610619 20140726 20148325 20213304 20310849 20404772 20561918 21115577 21479211 21544930 ...
 $ YGUsedAfter.K. : int  7990 15589 80568 178113 272036 429182 982841 1346475 1412181 1355412 ...
 $ Promoted.K.    : int  0 0 0 0 8226 937 65429 71166 62548 143638 ...
 $ YGCapacity.K.  : int  22649280 22649280 22649280 22649280 22649280 22649280 22649280 22649280 22649280 22649280 ...
 $ Pause.s.       : num  0.0379 0.022 0.0287 0.0509 0.109 ...

在这种情况下，我关心暂停时间（以秒为单位）。我想绘制一个图表，它将显示每个（挂钟）小时的平均值，基本上是一条线，2% 和 98% 是一条灰色走廊，最大值（每小时内）是一条红线。

我做了一些工作，但是使用q98函数很丑，不得不使用多行语句似乎很浪费，我不知道如何实现q02和q98之间的灰色区域：

q02 <- function(x, ...) {  x <- quantile(x,probs=c(0.2)) }
q98 <- function(x, ...) {  x <- quantile(x,probs=c(0.98)) }
hours = droplevels(cut(gcdata$date, breaks="hours")) # can I have 2 hours?
plot(aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=max),ylim=c(0,2), col="red", ylab="Pause(s)", xlab="Days") # Is always black?
lines(aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=q98),ylim=c(0,2), col="green")
lines(aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=q02),ylim=c(0,2), col="green")
lines(aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=mean),ylim=c(0,2), col="blue")

现在这会产生一个图表，其中黑点为最大值，蓝线为每小时平均值，上下 0,2 + 0,98 绿线。我认为有一条灰色走廊会更好读，也许是一条最大（红色）虚线并以某种方式修复轴标签。有什么建议？（文件在上面）

【问题讨论】：

删除了更好的剪切+粘贴提示。

标签： java r plot quantile

【解决方案1】：

您必须尝试polygon。这段代码很有用：

y98 = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=q98)
y02 = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=q02)
ymax = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=max)
ymin = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=min)
ymean = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=mean)

x = ymean[,1]
y1 = cbind(y02[,2], ymean[,2], y98[,2])
y2 = cbind(ymin[,2], ymean[,2], ymax[,2])

plotAreaCI(x,y2, ylim=c(0,2), xlab="time", ylab="variable")
plotAreaCI(x,y1, ylim=c(0,2), poly.col="blue", add=TRUE)

或

plotAreaCI(x,y2, ylim=c(0,2), xlab="time", ylab="variable", nice.x = TRUE)
plotAreaCI(x,y1, ylim=c(0,2), mean.lwd=2, poly.col="blue", add=TRUE)

函数plotAreaCI定义为：

plotAreaCI = function(x, y, add=FALSE, nice.x = FALSE,
                          xlim=NULL, ylim=NULL,
                          mean.col="black", mean.lwd=1.5,
                          poly.col="gray", poly.lty=3,
                          xlab=NULL, ylab=NULL, main="",
                          ...) {
      isFactorX = isClass("factor", x)
      if(isFactorX) {
        x.label = x
        x = as.numeric(x)
      }
      if(is.null(xlim)) xlim=range(x, na.rm=TRUE)
      if(is.null(ylim)) ylim=range(y, na.rm=TRUE)
      x.pol = c(x, rev(x), x[1])
      y.pol = c(y[,1], rev(y[,3]), y[,1][3])
      if(!add) {
        plot.new()
        plot.window(xlim=xlim, ylim=ylim, ...)
        if(!nice.x & isFactorX) {
          axis(1, at=x, labels=x.label)
        } else {
          xticks = axTicks(1)
          if(isFactorX) {
            xticks = xticks[xticks>=1]
            axis(1, at=xticks, labels=x.label[xticks])
          } else {
            axis(1)
          }
        }
            axis(2, las=1)
        box()
        title(xlab=xlab, ylab=ylab, main=main)
      }
      polygon(x.pol, y.pol, col=poly.col, lty=poly.lty)
      lines(x, y[,2], col=mean.col, lwd=mean.lwd)
      return(invisible())
    }

【讨论】：

谢谢，效果很好。我猜 plot(type="n") 主要需要传递变量参数（用于指定 x/ylab 等？）。
好的，根据这个线程 plot() 中的 type="n" 不起作用，因为它正在处理一个因素。我有点需要一种解决方法。除此之外，在你的帮助下，我现在有一个 2 种色调的走廊，最大和平均线，这很好。 r.789695.n4.nabble.com/…
我没有海'x'是一个因素，我正在编辑代码并添加一些临时解决方案。
我认为“x”轴并不是一个真正的因素，这就是为什么这种图表是合适的，但找到了一个将因素转换为数字并返回的解决方案。还有情节（type =“n”），我不再使用它了，谢谢。

【解决方案2】：

很高兴在这里见到 Debian 老前辈 :) 你的回答已经很不错了。由于我碰巧在时间序列方面做了很多工作，所以我想我会使用出色的 zoo 和 xts 包添加一个变体。后者建立在前者之上，除其他外，我们可以在此处使用 period.apply() 函数以及 endpoints() 函数来获取两小时的汇总数据。

所以在顶部我会使用

library(zoo)                                # for zoo objects
library(xts)                                # for period.apply

gcdata <- read.table("http://bernd.eckenfels.net/view/gc1001.ygc.csv",
                     header=TRUE, sep=",", dec=".")
timestamps <- gcdata$Timestamp + 
              as.POSIXct(strptime("2012-01-01 00:00:00", 
                         format="%Y-%m-%d %H:%M:%S"))
gcdatazoo <- zoo(gcdata[-1], order.by=timestamps)    # as zoo object

创建一个zoo 对象。你的功能仍然存在：

plotAreaCorridor <- function(x, y, col.poly1="lightgray", col.poly2="gray",...) {
    x.pol <- c(x, rev(x), x[1])
    y.pol <- c(y[,1], rev(y[,5]),y[,1][1])
    plot(x, y[,6]+1, type="n", ...) 
    polygon(x.pol, y.pol, col=col.poly1, lty=0)

    x.pol <- c(x, rev(x), x[1])
    y.pol <- c(y[,2], rev(y[,4]), y[,1][1])
    polygon(x.pol, y.pol, col=col.poly2, lty=0)

    lines(x, y[,3], col="blue") # median
    lines(x, y[,6], col="red")  # max

    invisible(NULL)
}

然后我们可以稍微简化一下：

agg <- period.apply(gcdatazoo[,"Pause.s."],               # to which data
                    INDEX=endpoints(gcdatazoo, "hours", k=2), # every 2 hours
                    FUN=function(x) quantile(x,               # what fun.
                                             probs=c(5,20,50,80,95,100)/100)) 

#v99 = q99(gcdata$Pause.s.)        # what is q99 ?
v99 <- mean(agg[,5])                  # mean of 95-th percentile?
plotAreaCorridor(index(agg),          # use time index as x axis
                 coredata(agg),       # and matrix part of zoo object as data
                 ylim=c(0,max(agg[,5])*1.5),
                 ylab="Quantiles of GC events",
                 main="NewPar Collection Activity")
abline(h=median(gcdatazoo[,"Pause.s."]), col="lightblue")
abline(h=v99, col="grey")
labeltxt <- paste("99%=",round(v99,digits=3),"s n=", nrow(gcdatazoo),sep="")
text(x=index(agg)[20], y=1.5*v99, labeltxt, col="grey", pos=3)  # or legend()

这给了

轴现在是自动的，仅在跨度小于周时显示工作日；这可以根据需要覆盖。

【讨论】：

酷，danke Dirk（q99 是第 99 个百分位数）

【解决方案3】：

这是我用来绘制实验室分析物（在本例中为收缩压）随时间变化的代码：

 SBP.qtr.mat <- aggregate(set1HLI$SBP, 
                          list(  year(set1HLI$Drawdt)+0.25* quarter(set1HLI$Drawdt)), 
                           quantile, prob=c(0.1,0.25,0.5,0.75, 0.9,0.95, 0.975), na.rm=TRUE)
 matplot(SBP.qtr.mat[,1], SBP.qtr.mat$x, type="pl")

不应该太难适应你的问题......或者你可以发布一个可重现的例子来使用。这给出了单个 data.frame 中的第 10、25、50、75、90、95 和 97.5 个百分位数，matplot 处理此类对象的绘图。

灰色区域？，...通常的方法是绘制一个多边形，从下界向外，在最右边“转弯”，然后在高边返回，然后在左侧连接回来. polygon 参数设置为 x, y。您可以将 col 参数设置为“灰色”。

制作“2 小时”序列，您可以将数据框合并到其中或与cut.POSIXt" as a breaks argument , there is the option of using multiples of time units withseq.POSIXt` 一起使用：

> seq(ISOdate(1910,1,1), ISOdate(1999,1,1), "10 years")
[1] "1910-01-01 12:00:00 GMT" "1920-01-01 12:00:00 GMT" "1930-01-01 12:00:00 GMT" "1940-01-01 12:00:00 GMT"
[5] "1950-01-01 12:00:00 GMT" "1960-01-01 12:00:00 GMT" "1970-01-01 12:00:00 GMT" "1980-01-01 12:00:00 GMT"
[9] "1990-01-01 12:00:00 GMT"

我没有看到它记录在案，但您可以将间隔的倍数与 cut.POSIXt 一起使用：

> str( cut( seq(ISOdate(1910,1,1), ISOdate(1999,1,1), "years"), "10 years") )
 Factor w/ 9 levels "1910-01-01","1920-01-01",..: 1 1 1 1 1 1 1 1 1 1 ...
> str( cut( seq(ISOdate(1910,1,1), ISOdate(1999,1,1), "years"), "5 years") )
 Factor w/ 18 levels "1910-01-01","1915-01-01",..: 1 1 1 1 1 2 2 2 2 2 ...

【讨论】：

这个例子是可以重现的，不是吗？
我猜是的。起初我没有意识到这一点。您可以按时使用我的 cmets 来修复 Ricardo 的代码。

【解决方案4】：

我目前还没有到达以下脚本（仍然需要查看 DWin 的更高级答案）。它现在看起来有点像我在寻找，但代码仍然很丑陋（例如我不知道如何对齐标签以及如何获得正确的 xlab 标签）：

plotAreaCorridor = function(x, y, col.poly1="lightgray", col.poly2="gray",...) {
   x.pol = c(x, rev(x), x[1])
   y.pol = c(y[,1], rev(y[,5]),y[,1][1])
   plot(x, y[,6]+1, type="n", ...) # ugly since type="n" does not work for factor
   polygon(x.pol, y.pol, col=col.poly1, lty=0)

   x.pol = c(x, rev(x), x[1])
   y.pol = c(y[,2], rev(y[,4]), y[,1][1])
   polygon(x.pol, y.pol, col=col.poly2, lty=0)

   lines(x, y[,3], col="blue") # median
   lines(x, y[,6], col="red")  # max

   return(invisible())
}
pause = gcdata$Pause.s.
hours = droplevels(cut(gcdata$date, breaks="hours")) # can I have 2 hours?
agg = aggregate(pause ~ hours, FUN=quantile, probs=c(5,20,50,80,95,100)/100)
x = agg$hours
ys = agg$pause
q99 <- function(x, ...) {  x <- quantile(x,probs=c(0.99)) }  
v99 = q99(gcdata$Pause.s.)
vmed = median(gcdata$Pause.s.)
plotAreaCorridor(x, ys,ylim=c(0,v99*1.5))
abline(h=vmed, col="lightblue")
abline(h=v99, col="grey")
label=paste("99%=",round(v99,digits=3),"s n=", length(gcdata$date),sep="")
text(x=30, y=v99, label, col="grey", pos=3)
title("NewPar Collection Activity")

【讨论】：