按类别逐列计算加权中位数（matrixStats）

Question

假设我有两个数据框“values”和“weights”，我想按类别（A，B，C）列计算按权重加权的年度中位数（year1，year2）：

values <- data.frame(TICKER=c("A","A","B","B","B","C","C","C","C"), year1=c(1,2,3,4,5,6,7,8,9), year2=c(9,8,7,6,5,4,3,2,1))
weights <- data.frame(TICKER=c("A","A","B","B","B","C","C","C","C"), year1=c(0.3,0.7,0.25,0.25,0.5,0.1,0.1,0.6,0.2), year2=c(0.6,0.4,0.3,0.5,0.2,0.4,0.2,0.1,0.3))

为此，我想使用ddply和weightedMedian函数（matrixStats包）。

output <- ddply(values, .(TICKER), colwise(weightedMedian(values, weights), na.rm=TRUE))

但是，我收到了错误信息：

"(list) object cannot be coerced to type 'double'"

有人知道如何调整代码以获得可行的解决方案吗？

我尝试将数据框转换为矩阵（通过as.matrix），因为weightedMedian需要矩阵作为输入。然而，这并没有帮助。 到目前为止，我找到的唯一解决方案是使用子集循环（但这非常缓慢且不太优雅）。

output <- matrix(data=0, nrow=3, ncol=2)
for (i in 2:ncol(values)){
 for (j in 1:length(unique(values$TICKER))){
  values.j <- subset(values, values$TICKER == as.character(unique(values$TICKER)[j]))
  weights.j <-  subset(weights, weights$TICKER == as.character(unique(values$TICKER)[j]))
  output[j,(i-1)] <- weightedMedian(values.j[,i], weights.j[,i], na.rm=TRUE)
}}

任何帮助都将不胜感激。非常感谢。

- Marcel

---

嗨 - 不，真实数据包含约70个周期（列）和大约15,000个类别。 - Marcel

Answer 1

除了OP提到的weightedMedian函数外，Hmisc包还提供了一个更通用的wtd.quantile函数。

我将两个数据框拆分成列表，并使用嵌套的sapply将这些函数应用于两个年份变量。通过比较下面的结果，可以看出weightedMedian产生了所需的结果。

为了准备数据，将值和权重分别拆分成它们的TICKER列表。

# split values and weights into lists by category
valuesList <- split(values, values$TICKER)
weightsList <- split(weights, values$TICKER)

如果我在上述代码中使用 OP 问题中的 weightedMedian，我会得到以下结果：

library(matrixStats)
sapply(names(valuesList),
  function(i) sapply(names(valuesList[[i]])[-1],
                function(j) weightedMedian(valuesList[[i]][[j]],
                                           w=weightsList[[i]][[j]])))

        A        B C
year1 1.7 4.333333 8
year2 8.6 6.125000 3

---

另一个包，Hmisc，有一个加权分位数函数wtd.quantile。

# load Hmisc package
library(Hmisc)

sapply(names(valuesList),
   function(i) sapply(names(valuesList[[i]])[-1],
                   function(j) {
                     wtd.quantile(valuesList[[i]][[j]],
                                  weights=weightsList[[i]][[j]], probs=0.5)}))

这将返回

myMedians 
          A B C
year1.50% 2 5 9
year2.50% 9 7 4

经过检查，matrixStats 的结果似乎更加合理。例如，TICKER == C，year == 2 不应该是4。

- lmo

Answer 2

如果您想留在plyr/weightedStats上下文中，我建议首先合并两个data.frames，然后结合您的值变量的已知列索引使用ldply和ddply。

df <- data.frame(values,wt=weights)

output <- lapply(names(values)[-1], 
   function(i) ddply(df,.(TICKER), 
     function(x) setNames(weightedMedian(x=x[,i],w=x[,match(i,names(x))+ncol(x)/2]),i)))

然后你可以采取以下方法将结果合并到单个数据框中。

do.call('join',output)

它会给你

  TICKER    year1 year2
1      A 1.700000 8.600
2      B 4.333333 6.125
3      C 8.000000 3.000

顺便说一下，您出现错误的原因是您只对values数据框进行了“切片”，而在weightedMedian中调用了整个weights数据框。

- jkt