数据科学中的R和Python: 用plyr包扩展apply族函数的功能

apply族函数是R语言中很有特色的一类函数，包括了apply、sapply、lapply、tapply、aggregate等等。在这篇博文里对它们进行了简略的说明。这一类函数本质上是将数据进行分割、计算和整合。它们在数据分析的各个阶段都有很好的用处。例如在数据准备阶段，我们可以按某个标准将数据分组，然后获得各组的统计描述。或是在建模阶段，为不同组的数据建立模型并比较建模结果。apply族函数与Google提出的mapreduce策略有着一致的思路。因为mapreduce的思路也是将数据进行分割、计算和整合。只不过它是将分割后的数据分发给多个处理核心进行运算。如果你熟悉了apply族函数，那么将数据转为并行运算是轻而易举的事情。plyr包则可看作是apply族函数的扩展，使之更容易运用，功能更为强大。

plyr包的主函数是**ply形式的，其中首字母可以是(d、l、a)，第二个字母可以是 (d、l、a、_)，不同的字母表示不同的数据格式，d表示数据框格式，l表示列表，a表示数组，_则表示没有输出。第一个字母表示输入的待处理的数据格式，第二个字母表示输出的数据格式。例如ddply函数，即表示输入一个数据框，输出也是一个数据框。

下面首先来用一个简单的例子说明一下用法。还是用iris数据集，其中包括了一个分类变量和四个数值变量。我们希望数据按不同类别，分别计算数值变量的均值。下面我们分别用三种方法来得到同样的结果。

library(plyr)
library(reshape2)
# 用aggregate函数进行数据汇总
aggregate(iris[1:4],list(iris$Species),mean)
# 用reshape2包进行数据汇总
data.melt <- melt(iris,id=c('Species'))
dcast(data.melt,Species~variable,mean)
# 用ddply函数进行数据汇总
ddply(iris,.(Species),function(df) mean(df[1:4]))

初看起来plyr包所具备的功能并不很出彩，下面我们看一个略为复杂例子。还是用iris数据，我们希望对每一种花做一个简单回归。

# 首先定义回归函数
model <- function(x) {
    lm(Petal.Length~Petal.Width,data=x)
}
# 如果用普通的函数则需要如下的分割、计算、整合三个步骤共四条命令
pieces <- split(iris,list(iris$Species))
models <- lapply(pieces,model)
result <- lapply(models,coef)
do.call('rbind',result)
# 用plyr包只用下面两个函数，每个函数都内置了分割、计算、整合的功能。
result1 <- dlply(iris,.(Species),model)
result2 <- ldply(result1,function(x) coef(x))

plyr包中还有两个比较特别的函数，分别是r*ply和m*ply，它们分别对应的是replicate和mapply函数。

replicate(20,mean(runif(100)))
rdply(20, mean(runif(100)))
mapply(rnorm,mean=1:5,sd=1:5, n=2)
mdply(data.frame(mean = 1:5, sd = 1:5), rnorm, n = 2)

最后我们来看一个mdply函数的应用，我们希望用神经网络包来为不同的花进行分类，使用BP神经网络需要的一个参数就是隐藏层神经元的个数。我们来尝试用1到10这十个参数运行模型十次，并观察十个建模结果的预测准确率。但我们并不需要手动运行十次。而是使用mdply函数来完成这个任务。

library(nnet)
# 确定建模函数
nnet.m <- function(...) {
  nnet(Species~.,data=iris,trace=F,...)
}
# 确定输入参数
opts <- data.frame(size=1:10,maxiter=50)
# 建立预测准确率的函数
accuracy <- function(mod,true) {
  pred <- factor(predict(mod,type='class'),levels=levels(true))
  tb <- table(pred,true)
  sum(diag(tb))/sum(tb)
}
# 用mlply函数建立包括10个元素的列表，每个元素包括了一个建模结果
models <- mlply(opts,nnet.m)
# 再用ldply函数读取列表，计算后得到最终结果
ldply(models,'accuracy',true=iris$Species)