在R中对数据子集进行计算

Question

我希望能够对数据框中PERMNO列中的每个公司编号进行计算，以下是摘要：

> summary(companydataRETS)
     PERMNO           RET           
 Min.   :10000   Min.   :-0.971698  
 1st Qu.:32716   1st Qu.:-0.011905  
 Median :61735   Median : 0.000000  
 Mean   :56788   Mean   : 0.000799  
 3rd Qu.:80280   3rd Qu.: 0.010989  
 Max.   :93436   Max.   :19.000000  

到目前为止，我的解决方案是创建一个包含所有可能公司编号的变量。

compns <- companydataRETS[!duplicated(companydataRETS[,"PERMNO"]),"PERMNO"]

然后使用并行计算的 foreach 循环调用我的函数 get.rho()，该函数执行所需的计算。

rhos <- foreach (i=1:length(compns), .combine=rbind) %dopar% 
      get.rho(subset(companydataRETS[,"RET"],companydataRETS$PERMNO == compns[i]))

我已经对我的一部分数据进行了测试，一切正常。问题是我有7200万个观察值，即使让电脑整夜运行，它仍然无法完成。

我是R的新手，所以我想我的代码结构可以改进，并且有更好的（更快，计算量更小）方法来执行相同的任务（也许使用apply或with，但我不理解）。 有什么建议吗？

- Vivi

---

必须提及 data.table 软件包。 - joran

没有一个好的答案回答你的主要问题，但是你公司名称的唯一列表可能更容易获取为unique(companydataRETS$PERMNO)。 - Matt Parker

@joran 我会看一下 data.table 包。 - Vivi

@joran，我安装了data.table并修改了代码以使用它，system.time()的结果从使用foreach的(43.925 12.413 56.337)减少到使用data.table的(0.229 0.047 0.276)，这太不可思议了！真正符合我的期望。你认为你能在我的问题中发布一个答案吗？我稍后可以用代码修改来完成它，但我希望积分归你... - Vivi

不需要了，我贡献的只是指针。你自己发布代码并进行时间比较，然后把分数给自己。我有很多分数。 :) - joran

Answer 1

根据Joran的建议，我研究了data.table库。修改代码如下：

library(data.table) 
companydataRETS <- data.table(companydataRETS)
setkey(companydataRETS,PERMNO)

rhos <- foreach (i=1:length(compns), .combine=rbind) %do% 
      get.rho(companydataRETS[J(compns[i])]$RET)

我按照原本的方式运行了代码（使用subset），还有一次使用data.table，其中变量compns只包括数据集中的30个公司，这是两个版本的system.time()输出:

使用subset:

用户.......系统.....经过
43.925 ... 12.413...... 56.337

使用data.table

用户.......系统.....经过
0.229..... 0.047....... 0.276

（出于某种原因，在原始代码中使用%do%而不是%dopar%让它运行得更快。对于subset的system.time()是在这种情况下较快的%do%。）

我让原始代码跑了一整夜，5个小时后还没完成，所以我放弃并终止了它。通过这个小修改，我在不到5分钟内得到了结果（我想大约3分钟左右）！

编辑

使用data.table甚至有更简单的方法来做到这一点，不需要使用foreach，只需将上面代码的最后一行替换为

rhos <- companydataRETS[ , get.rho(RET), by=PERMNO]

- Vivi

---

很高兴它能够正常工作。如果您可以提供更多有关您不理解的冲突的细节（例如错误消息），那么我们可以提供帮助。请先检查最新的实时NEWS（也链接在?data.table的顶部）以查看问题是否已经被报告并修复。否则，请就此提出一个新问题。这对我来说没有任何印象。 - Matt Dowle

@MatthewDowle，错误信息与我之前尝试使用库doMC并在使用source()加载一些函数后得到的错误相同。我现在尝试重现错误，但没有发生。这可能与data.table无关。我将编辑上面的答案以删除评论，如果再次发生，我会发布另一个问题。顺便说一句，我做了你上面建议的crantastic事情。 - Vivi

%do%和%dopar%的速度取决于你使用的并行后端 - 并行运行具有将数据拆分、分配到不同节点以及重新组合的所有开销。如果你没有很多用于并行处理的节点，或者你应用的函数非常简单，那么它有时会比%do%更慢。 - Matt Parker

Answer 2

有很多方法可以做到这样的事情，您的foreach解决方案是其中之一。仅从您提供的代码来看，我只能猜测最合适的解决方案...

然而，我认为您的代码中最大的减速因素实际上是您的get.rho函数，而不是循环或子集。如果您想分享该函数，我敢打赌您会得到一些惊人的答案，这些答案既可以加快速度，又可以澄清一些“R-isms”。

话虽如此，还有许多替代方法可以完成您正在进行的操作。

plyr包专门用于此类计算。它使用拆分-应用-组合策略。函数的前两个字母表示输入和输出数据类型。

由于您正在输入一个数据框并输出一个数据框，因此要使用ddply函数：

library(plyr)
ddply(companydataRETS, .(PERMNO), summarise, get.rho(RET))

如果你不在使用Windows系统，你可以轻松地使用多线程计算

library(doMC)
registerDoMC()
ddply(companydataRETS, .(PERMNO), summarise, get.rho(RET), .parallel=TRUE)

tapply 也是一个完美的选择：

tapply(companydataRETS$RET, companydataRET$PERMNO, get.rho)

data.table 包，如评论中所述，也非常擅长这个，但是我会把那段代码留给读者作为练习。

然而，正如我上面所说，如果你的 get.rho 函数很慢，无论你用什么样的子集和循环技巧，计算都会花费很长时间。

---

在帖子中编辑函数代码：

如果这是时间序列数据或可以视为此类数据，则有许多包和函数可进行此类滞后比较。我对它们不是很熟悉，但快速浏览谷歌和CRAN任务视图将为您提供出色的选项概述。

我没有进行基准测试，但我认为可以安全地假设您的代码最慢的部分是lm调用。仅对数据的sample进行操作而不是全部集合将大大加快速度。但我相信有人会有更好、更完整的解决方案。

- Justin

---

我没有使用Windows，但我正在使用doMC进行并行计算。我的get.rho()函数很糟糕，肯定可以改进。我将在主要问题中包含它。 - Vivi

针对建议在一个拥有7200万行数据集中使用plyr的做法，给出-1分。特别是因为plyr是专门为此而设计的？ - Matt Dowle

@MatthewDowle 好的。我本意是想探讨可能性，而不是给出完美答案。但是，你说得对，在plyr中处理7200万行数据是一个非常糟糕的想法。 - Justin