构建均值/方差汇总表的快速/优雅方式答案

【问题标题】：quick/elegant way to construct mean/variance summary table构建均值/方差汇总表的快速/优雅方式
【发布时间】：2011-09-16 18:58:03
【问题描述】：

我可以完成这项任务，但我觉得必须有一种“最好的”（最流畅、最紧凑、最清晰的代码、最快？）的方式来完成它，但到目前为止还没有想出来......

对于一组指定的分类因素，我想按组构建一个均值和方差表。

生成数据：

set.seed(1001)
d <- expand.grid(f1=LETTERS[1:3],f2=letters[1:3],
                 f3=factor(as.character(as.roman(1:3))),rep=1:4)
d$y <- runif(nrow(d))
d$z <- rnorm(nrow(d))

想要的输出：

  f1 f2  f3    y.mean      y.var
1  A  a   I 0.6502307 0.09537958
2  A  a  II 0.4876630 0.11079670
3  A  a III 0.3102926 0.20280568
4  A  b   I 0.3914084 0.05869310
5  A  b  II 0.5257355 0.21863126
6  A  b III 0.3356860 0.07943314
... etc. ...

使用aggregate/merge：

library(reshape)
m1 <- aggregate(y~f1*f2*f3,data=d,FUN=mean)
m2 <- aggregate(y~f1*f2*f3,data=d,FUN=var)
mvtab <- merge(rename(m1,c(y="y.mean")),
      rename(m2,c(y="y.var")))

使用ddply/summarise（可能是最好的，但无法使其工作）：

mvtab2 <- ddply(subset(d,select=-c(z,rep)),
                .(f1,f2,f3),
                summarise,numcolwise(mean),numcolwise(var))

结果

Error in output[[var]][rng] <- df[[var]] : 
  incompatible types (from closure to logical) in subassignment type fix

使用melt/cast（也许是最好的？）

mvtab3 <- cast(melt(subset(d,select=-c(z,rep)),
          id.vars=1:3),
     ...~.,fun.aggregate=c(mean,var))
## now have to drop "variable"
mvtab3 <- subset(mvtab3,select=-variable)
## also should rename response variables

不会 (?) 在 reshape2 中工作。向某人解释...~. 可能很棘手！

【问题讨论】：

标签： r aggregate plyr reshape2

【解决方案1】：

这是使用data.table的解决方案

library(data.table)
d2 = data.table(d)
ans = d2[,list(avg_y = mean(y), var_y = var(y)), 'f1, f2, f3']

【讨论】：

【解决方案2】：

我有点困惑。这不起作用：

mvtab2 <- ddply(d,.(f1,f2,f3),
            summarise,y.mean = mean(y),y.var = var(y))

这给了我这样的东西：

   f1 f2  f3    y.mean       y.var
1   A  a   I 0.6502307 0.095379578
2   A  a  II 0.4876630 0.110796695
3   A  a III 0.3102926 0.202805677
4   A  b   I 0.3914084 0.058693103
5   A  b  II 0.5257355 0.218631264

格式正确，但看起来值与您指定的不同。

编辑

以下是使用 numcolwise 制作您的版本的方法：

mvtab2 <- ddply(subset(d,select=-c(z,rep)),.(f1,f2,f3),summarise,
                y.mean = numcolwise(mean)(piece),
                y.var = numcolwise(var)(piece))

您忘记将实际数据传递给numcolwise。然后是ddply 的小技巧，每个部分在内部都被称为piece。（不应依赖 Hadley 在 cmets 中指出的内容，因为它可能会在 plyr 的未来版本中发生变化。）

【讨论】：

价值观问题是因为我并没有真正做到set.seed 我声称的地方——后来添加了它（哎呀，现在已修复，所以你关于“不同”的评论不再有意义——对不起）。这看起来很有效。我需要弄清楚numcolwise 是怎么回事——我认为它会自动处理数据框中只剩下一个数字变量的事实，但也许我对此感到困惑。
很难接受，但这是第一个答案，我喜欢优雅（无需定义匿名函数），即使速度较慢。
你不应该依赖任何内部变量的使用——当我最终弄清楚如何修复 plyr 中的作用域时，它会中断
@hadley Bummer！我实际上认为这个小技巧在很多情况下都很方便。
如果您发现自己这样做了，这可能表明您不应该使用摘要。不过，我确实需要更多地考虑代词。

【解决方案3】：

（我投票给 Joshua 的。）这是一个 Hmisc::summary.formula 解决方案。这对我来说的好处是它与 Hmisc::latex 输出“通道”很好地集成在一起。

summary(y ~ interaction(f3,f2,f1), data=d, method="response", 
                    fun=function(y) c(mean.y=mean(y) ,var.y=var(y) ))
#-----output----------
y    N=108

+-----------------------+-------+---+---------+-----------+
|                       |       |N  |mean.y   |var.y      |
+-----------------------+-------+---+---------+-----------+
|interaction(f3, f2, f1)|I.a.A  |  4|0.6502307|0.095379578|
|                       |II.a.A |  4|0.4876630|0.110796695|

截断输出以显示乳胶 -> PDF -> png 输出：

【讨论】：

【解决方案4】：

@joran 对ddply 的回答非常准确。以下是我将如何使用aggregate。请注意，我避免使用公式界面（它比较慢）。

aggregate(d$y, d[,c("f1","f2","f3")], FUN=function(x) c(mean=mean(x),var=var(x)))

【讨论】：

【解决方案5】：

我有点沉迷于速度比较，尽管在这种情况下它们对我来说基本上是无关紧要的......

joran_ddply <- function(d) ddply(d,.(f1,f2,f3),
                                 summarise,y.mean = mean(y),y.var = var(y))
joshulrich_aggregate <- function(d) {
  aggregate(d$y, d[,c("f1","f2","f3")],
            FUN=function(x) c(mean=mean(x),var=var(x)))
}

formula_aggregate <- function(d) {
  aggregate(y~f1*f2*f3,data=d,
            FUN=function(x) c(mean=mean(x),var=var(x)))
}
library(data.table)
d2 <- data.table(d)
ramnath_datatable <- function(d) {
  d[,list(avg_y = mean(y), var_y = var(y)), 'f1, f2, f3']
}


library(Hmisc)
dwin_hmisc <- function(d) {summary(y ~ interaction(f3,f2,f1), 
                   data=d, method="response", 
                   fun=function(y) c(mean.y=mean(y) ,var.y=var(y) ))
                         }


library(rbenchmark)
benchmark(joran_ddply(d),
          joshulrich_aggregate(d),
          ramnath_datatable(d2),
          formula_aggregate(d),
          dwin_hmisc(d))

aggregate 是最快的（甚至比data.table 还要快，这让我很惊讶，尽管要聚合的表格更大，情况可能会有所不同），即使使用公式界面也是如此......）

                     test replications elapsed relative user.self sys.self
5           dwin_hmisc(d)          100   1.235 2.125645     1.168    0.044
4    formula_aggregate(d)          100   0.703 1.209983     0.656    0.036
1          joran_ddply(d)          100   3.345 5.757315     3.152    0.144
2 joshulrich_aggregate(d)          100   0.581 1.000000     0.596    0.000
3   ramnath_datatable(d2)          100   0.750 1.290878     0.708    0.000

（现在我只需要 Dirk 站出来发布一个Rcpp 解决方案，它比其他任何东西都快 1000 倍......）

【讨论】：

我检查了有 2700 行的较大表，发现 data.table 比基于 aggregate 的解决方案高出 1.5 倍。
@Ramnath 如果表更大，您的解决方案的改进是 21 倍快约书亚的聚合。
这并不奇怪，因为data.table 已经优化了多少。完全归功于 Matt Dowle 和 Arun :)

【解决方案6】：

我发现doBy package 有一些非常方便的功能来处理这样的事情。例如，函数?summaryBy 非常方便。考虑：

> summaryBy(y~f1+f2+f3, data=d, FUN=c(mean, var))
   f1 f2  f3    y.mean       y.var
1   A  a   I 0.6502307 0.095379578
2   A  a  II 0.4876630 0.110796695
3   A  a III 0.3102926 0.202805677
4   A  b   I 0.3914084 0.058693103
5   A  b  II 0.5257355 0.218631264
6   A  b III 0.3356860 0.079433136
7   A  c   I 0.3367841 0.079487973
8   A  c  II 0.6273320 0.041373836
9   A  c III 0.4532720 0.022779672
10  B  a   I 0.6688221 0.044184575
11  B  a  II 0.5514724 0.020359289
12  B  a III 0.6389354 0.104056229
13  B  b   I 0.5052346 0.138379070
14  B  b  II 0.3933283 0.050261804
15  B  b III 0.5953874 0.161943989
16  B  c   I 0.3490460 0.079286849
17  B  c  II 0.5534569 0.207381592
18  B  c III 0.4652424 0.187463143
19  C  a   I 0.3340988 0.004994589
20  C  a  II 0.3970315 0.126967554
21  C  a III 0.3580250 0.066769484
22  C  b   I 0.7676858 0.124945402
23  C  b  II 0.3613772 0.182689385
24  C  b III 0.4175562 0.095933470
25  C  c   I 0.3592491 0.039832864
26  C  c  II 0.7882591 0.084271963
27  C  c III 0.3936949 0.085758343

所以函数调用很简单，易于使用，而且我会说，优雅。

现在，如果您主要关心的是速度，那么它似乎是合理的——至少对于较小的任务（请注意，无论出于何种原因，我都无法让 ramnath_datatable 函数工作）：

                     test replications elapsed relative user.self 
4           dwin_hmisc(d)          100    0.50    2.778      0.50 
3    formula_aggregate(d)          100    0.23    1.278      0.24 
5       gung_summaryBy(d)          100    0.34    1.889      0.35 
1          joran_ddply(d)          100    1.34    7.444      1.32 
2 joshulrich_aggregate(d)          100    0.18    1.000      0.19

【讨论】：

【解决方案7】：

我遇到了这个问题，发现基准测试是用小表完成的，所以很难判断哪种方法在 100 行时更好。

我还稍微修改了数据以使其“未排序”，这将是一种更常见的情况，例如数据在数据库中。我已经添加了更多 data.table 试验，以查看是否预先设置键更快。似乎在这里，预先设置密钥并不能提高性能，因此 ramnath 解决方案似乎是最快的。

set.seed(1001)
d <- data.frame(f1 = sample(LETTERS[1:3], 30e5, replace = T), f2 = sample(letters[1:3], 30e5, replace = T),
                f3 = sample(factor(as.character(as.roman(1:3))), 30e5, replace = T), rep = sample(1:4, replace = T))

d$y <- runif(nrow(d))
d$z <- rnorm(nrow(d))

str(d)

require(Hmisc)
require(plyr)
require(data.table)
d2 = data.table(d)
d3 = data.table(d)

# Set key of d3 to compare how fast it is if the DT is already keyded
setkey(d3,f1,f2,f3)

joran_ddply <- function(d) ddply(d,.(f1,f2,f3),
                                 summarise,y.mean = mean(y),y.var = var(y))

formula_aggregate <- function(d) {
  aggregate(y~f1*f2*f3,data=d,
            FUN=function(x) c(mean=mean(x),var=var(x)))
}

ramnath_datatable <- function(d) {
  d[,list(avg_y = mean(y), var_y = var(y)), 'f1,f2,f3']
}

key_agg_datatable <- function(d) {
  setkey(d2,f1,f2,f3)
  d[,list(avg_y = mean(y), var_y = var(y)), 'f1,f2,f3']
}

one_key_datatable <- function(d) {
  setkey(d2,f1)
  d[,list(avg_y = mean(y), var_y = var(y)), 'f1,f2,f3']
}

including_3key_datatable <- function(d) {
  d[,list(avg_y = mean(y), var_y = var(y)), 'f1,f2,f3']
}

dwin_hmisc <- function(d) {summary(y ~ interaction(f3,f2,f1), 
                                   data=d, method="response", 
                                   fun=function(y) c(mean.y=mean(y) ,var.y=var(y) ))
}

require(rbenchmark)
benchmark(joran_ddply(d),
          joshulrich_aggregate(d),
          ramnath_datatable(d2),
          including_3key_datatable(d3),
          one_key_datatable(d2),
          key_agg_datatable(d2),
          formula_aggregate(d),
          dwin_hmisc(d)
          )

#                         test replications elapsed relative user.self sys.self
#                dwin_hmisc(d)          100 1757.28  252.121   1590.89   165.65
#         formula_aggregate(d)          100  433.56   62.204    390.83    42.50
# including_3key_datatable(d3)          100    7.00    1.004      6.02     0.98
#               joran_ddply(d)          100  173.39   24.877    119.35    53.95
#      joshulrich_aggregate(d)          100  328.51   47.132    307.14    21.22
#        key_agg_datatable(d2)          100   24.62    3.532     19.13     5.50
#        one_key_datatable(d2)          100   29.66    4.255     22.28     7.34
#        ramnath_datatable(d2)          100    6.97    1.000      5.96     1.01

【讨论】：

【解决方案8】：

这是使用 Hadley Wickham 的新 dplyr 库的解决方案。

library(dplyr)
d %>% group_by(f1, f2, f3) %>% 
summarise(y.mean = mean(y), z.mean = mean(z))

【讨论】：