快速/优雅的方法来构建均值/方差汇总表

Question

我可以完成这个任务,但我觉得必须有一个"最好的"(最简洁,最紧凑,最清晰的代码,最快？)的方式,到目前为止还没有弄清楚...

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

生成数据:

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....~.向某人解释可能会很棘手!

Answer 1

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

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

Answer 2

(我投票支持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输出:

Answer 3

我有点不解.这不起作用:

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在评论中指出的内容,因为它可能会在未来的版本中发生变化plyr.)

Answer 4

@joran与ddply答案一致.这是我将如何做到这一点aggregate.请注意,我避免使用公式界面(它更慢).

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

Answer 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倍的解决方案......)