有没有一种有效的方法来计算预聚合数据 (R) 的百分位数？

Question

首先：这是我在这里的第一个问题，我对 R 也比较陌生。所以，如果这是一个愚蠢的问题或错误的提问方式，我很抱歉。

我有一个这样的数据框：

df <- data.frame(Website = c("A", "A", "A", "B", "B", "B"),
             seconds = c(1,12,40,3,5,14),
             visitors = c(200000,100000,12000,250000,180000,90000))


> df
  Website seconds visitors
       A       1   200000
       A      12   100000
       A      40    12000
       B       3   250000
       B       5   180000
       B      14    90000

如何解读数据：网站 A 有 200000 名访问者在网站上仅停留 1 秒，100000 名访问者仅停留 12 秒，以此类推。 实际上，数据有大约一百个不同的网站，每个网站的秒数从 0 到大约 900 秒不等（分别有大量访问者）。

现在，我想计算访问持续时间（每个网站）的百分位数或至少四分位数。

我已经在这里找到并尝试了这个解决方案：https://stackoverflow.com/a/53882909 但是，这种解决方案效率非常低，因为它会生成包含数百万行的数据帧（并且处理时间很长）。

我现在的问题是：有没有更快（更有效的方法）从这些预先汇总的数据中计算百分位数？

Answer 1

我相信这会更快。首先创建一个函数来计算您指定的分位数。然后将数据拆分成列表，使用sapply：

quant <- function(x, p=c(.25, .50, .75)) {
        v <- c(0, cumsum(x$visitors)/sum(x$visitors))
        s <- c(0, x$seconds)
        approx(v, s, p)$y
}
df.split <- split(df, df$Website)
p <- c(.1, .2, .3, .4, .5, .6, .7, .8, .9)
stats <- t(sapply(df.split, quant, p=p))
colnames(stats) <- as.character(p)
round(stats, 1)
#   0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
# A 0.2 0.3 0.5 0.6 0.8 0.9 3.0 6.5 9.9
# B 0.6 1.2 1.9 2.5 3.1 3.7 4.3 4.8 8.8

为了更好地了解这里发生的情况，我们绘制了一个显示网站 A 数据的图表：

test1 <- df[1:3, ]
test1$cumvis <- cumsum(test1$visitors)
barplot(test1$seconds, test1$visitors, space=0, xlim=c(0, 325000))
axis(1, seq(0, 300000, 50000), c("0", "50K", "100K", "150K", "200K",
     "250K", "300K"), xpd=NA)
axis(3, seq(0, sum(test1$visitors), by=31200), seq(0, 1, by=.1), lty=1)
lines(c(0, test1$cumvis), c(0, test1$seconds), col="red", lwd=2)
lines(c(0, test1$cumvis-.5*test1$visitors, tail(test1$cumvis, 1)),
     c(0, test1$seconds, tail(test1$seconds,  1)), col="blue", lwd=2)

该图将数据显示为灰色矩形。底部 x 轴显示累计访问次数，顶部 x 轴显示累计比例。我们可以将矩形视为分布，或者我们可以假设矩形是近似底层分布的样本。我建议的解决方案采用红线并使用approx 函数在数据点之间使用线性插值来估计沿该曲线的秒数。

同样的方法可以用于曲线的不同定义，其中数据点放置在每个矩形的中间，蓝色曲线。我也会为这种方法提供代码。也可以在不复制原始数据的情况下估计分位数。

首先一个函数来估计沿蓝线的分位数：

quant2 <- function(x, p=c(.25, .50, .75)) {
        v <- c(0, cumsum(x$visitors)-(.5*x$visitors)/sum(x$visitors), 1)
        s <- c(0, x$seconds, tail(x$seconds, 1))
        approx(v, s, p)$y
}
p <- c(.1, .2, .3, .4, .5, .6, .7, .8, .9)
stats <- t(sapply(df.split, quant2, p=p))
colnames(stats) <- as.character(p)
round(stats, 1)
#   0.1 0.2  0.3  0.4 0.5  0.6  0.7  0.8  0.9
# A 4.0 8.0 12.0 16.0  20 24.0 28.0 32.0 36.0
# B 1.4 2.8  4.2  5.6   7  8.4  9.8 11.2 12.6

估计值更高，因为蓝线高于红线。

最后，我们可以简单地使用矩形而无需任何插值。基本上，我们在数据点的边界设置中断，并使用这些中断来确定哪些比例属于哪些观察组（秒）。

quant3 <- function(x, p=c(.25, .50, .75)){
    v <- c(0, cumsum(x$visitors)/sum(x$visitors))
    limits <- cut(p, breaks=v, include.lowest=TRUE, labels=x$seconds)
    limits <- as.numeric(as.character(limits))
}
p <- 0:10/10
stats <- t(sapply(df.split, quant3, p=p))
colnames(stats) <- as.character(p)
stats
#   0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9  1
# A 1   1   1   1   1   1   1  12  12  12 40
# B 3   3   3   3   3   5   5   5   5  14 14

所以对于网站 A，1 秒是分位数 0 - .6 的值。