R語言中plyr包
前言
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函數，即表示輸入一個數據框，輸出也是一個數據框。
plyr包是Hadley Wickham大神為解決split – apply – combine問題而寫的一個包，其動機在與提供超越for循環和內置的apply函數族的一個一攬子解決方案。使用plyr包可以針對不同的數據類型，在一個函數內同時完成split – apply – combine三個步驟。
plyr 的功能已經遠遠超出數據整容的范圍，Hadley在plyr中應用了split-apply-combine的數據處理哲學，即：先將數據分離，然后應用某些處理函數，最后將結果重新組合成所需的形式返回。某些人士喜歡用“揉”來表述這樣的數據處理；“揉”，把數據當面團搗來搗去，很哲，磚家們的磚頭落下來，拍死人絕不償命。
先別哲了，來點實際的：plyr的函數命名方式比較規律，很容易記憶和使用。比如 a開頭的函數aaply, adply 和 alply 將數組（array）分別轉成數組、數據框和列表；daply, ddply 和 dlply 將數據框分別轉成數組、數據框和列表；而laply, ldaply, llply將列表（list）分別轉成數組、數據框和列表。
目錄
1. 數據轉換: split – apply – combine 模式
2. 知識點_1
3. 知識點_2
4. 知識點_3
5. 知識點_4
6. 參考資料
數據轉換: split – apply – combine 模式
按：這一篇是讀Hadley Wickham的文章The Split-Apply-Combine Strategy for Data Analysis的筆記。
在數據分析中，有許多問題可以由類似的類型和方法步驟解決，可稱之為模式，設計模式或者分析模式。下面要討論的是數據轉換的一個常用模式：split – apply – combine。其解決之道，在R語言中，有3種方式：
（1） for 顯式循環,但是這種方式的缺點也很明顯，代碼長，易出錯，也難以並行化；
（2） 拜R語言的向量計算特點所賜，在R當中，大多數問題不需要用顯示循環方式，而代之以base包中的apply函數族及其它的一些函數，直接對向量，數組，列表和數據框實現循環的操作。
（3） Hadley Wickham大神覺得apply族還是不夠簡潔，所以開發了pylr包，以更少的代碼來解決split – apply – combine問題。
-
- split – apply – combine 模式
大型數據集通常是高度結構化的，結構使得我們可以按不同的方式分組，有時候我們需要關注單個組的數據片斷，有時需要聚合不同組內的信息，並相互比較。
因此對數據的轉換，可以采用split – apply – combine模式來進行處理：
split：把要處理的數據分割成小片斷；
apply：對每個小片斷獨立進行操作；
combine：把片斷重新組合。
利用這種模式解決問題在很多數據分析或編程問題中都會出現:
Python當中，是map( )，filter( )，reduce( )；
Google 有MapReduce；
R 當中是split( )，*apply( )，aggregate( )…，以及plyr包。
-
- 分划：split函數
在R當中，split這個步驟是由split( )，subset( )等等函數完成的。
下面主要介紹split這個函數。
函數split（）可以按照分組因子，把向量，矩陣和數據框進行適當的分組。它的返回值是一個列表，代表分組變量每個水平的觀測。這個列表可以使用sapply（）,lappy（）進行處理（apply – combine步驟），得到問題的最終結果。
split( )的基本用法是：group <- split(X,f)
其中X 是待分組的向量，矩陣或數據框。f是分組因子。
例1：對向量分組
> library(MASS)
#使用Cars93數據集，利用其中的Origin變量（兩個水平），對Price變量分組
> g<-split(Cars93$Price,Cars93$Origin)
#分組結果是個列表：
$USA
[1] 15.7 20.8 23.7 26.3 34.7 40.1 13.4 11.4 15.1 15.9 16.3 16.6 18.8 38.0
[15] 18.4 15.8 29.5 9.2 11.3 13.3 19.0 15.6 25.8 12.2 19.3 7.4 10.1 11.3
[29] 15.9 14.0 19.9 20.2 20.9 34.3 36.1 14.1 14.9 13.5 16.3 19.5 20.7 14.4
[43] 9.0 11.1 17.7 18.5 24.4 11.1
$`non-USA`
[1] 15.9 33.9 29.1 37.7 30.0 8.4 12.5 19.8 12.1 17.5 8.0 10.0 10.0 13.9
[15] 47.9 28.0 35.2 8.3 11.6 16.5 19.1 32.5 31.9 61.9 10.3 26.1 11.8 15.7
[29] 19.1 21.5 28.7 8.4 10.9 19.5 8.6 9.8 18.4 18.2 22.7 9.1 19.7 20.0
[43] 23.3 22.7 26.7
#計算組的長度和組內均值
> sapply(g,length)
USA non-USA
48 45
> sapply(g,mean)
USA non-USA
18.57292 20.50889
#用lapply也可以，返回值是列表
> lapply(g,mean)
$USA
[1] 18.57292
$`non-USA`
[1] 20.50889
例2：對矩陣分組（按列）
> m<-cbind(x=1:10,y=11:20)
> split(m,col(m))
$`1`
[1] 1 2 3 4 5 6 7 8 9 10
$`2`
[1] 11 12 13 14 15 16 17 18 19 20
例3：對數據框進行分組
#還是利用Cars93,它就是個數據框
g<-split(Cars93,Cars93$Origin)
#分組結果
> summary(g)
Length Class Mode
USA 27 data.frame list
non-USA 27 data.frame list
split還有一個逆函數，unsplit,可以讓分組完好如初。在base包里和split功能接近的函數有cut(對屬性數據分划)，strsplit(對字符串分划)以及subset（對向量，矩陣或數據框按給定條件取子集）等。
-
- base包中的apply()函數族
接下來的apply – combine步驟主要由apply函數族完成。對於R的初學者來說，apply函數族以其強大的向量運算功能，把for循環語句送到了被遺忘的角落，這種震撼效果，無異於天賜神器。但是，這些冠以*apply的函數，也經常讓人使用混亂，找不着方向。
那么apply包括哪些函數呢？在 R console中鍵入??apply,可以看到在base包中包含如下函數：（apply在別的包里還有起到各自功能的相關函數，我們只看涉及數據轉換的這些個）
apply ：Apply Functions Over Array Margins
by ：Apply a Function to a Data Frame Split by Factors
eapply ：Apply a Function Over Values in an Environment
lapply ：Apply a Function over a List or Vector
mapply ：Apply a Function to Multiple List or Vector Arguments
rapply ：Recursively Apply a Function to a List
tapply ：Apply a Function Over a Ragged Array
除此之外，還有可作為lapply變形的sapply，vapply和 replicate，共計10個函數。具體的內容，見apply()函數族的詳細講解！！！
知識點_1
（1）下面我們看看如何使用ldply函數將ath1121501.db包中的KEGG列表數據轉成數據框：
> library(ath1121501.db)
> keggs <- as.list(ath1121501PATH[mappedkeys(ath1121501PATH)])
> head(ldply(keggs, paste, collapse='; '))
.id V1
1 261579_at 00190
2 261569_at 04712
3 261583_at 00010; 00020; 00290; 00620; 00650; 01100; 01110
4 261574_at 00903; 00945; 01100; 01110
5 261043_at 00051; 00520; 01100
6 261044_at 04122
（2）對於簡單的問題，plyr和apply函數的效果差不多
> m <- matrix(c(1:4,1,4,1:6),ncol=3)
> apply(m,1,mean)
[1] 1.666667 3.333333 3.000000 4.000000
> aaply(m,1,mean)
1 2 3 4
1.666667 3.333333 3.000000 4.000000
plyr包的函數較多，不再一一介紹，更多用法請參考它的在線幫助，Hadley 也寫了很詳細的tutorial：http://plyr.had.co.nz/09-user/
知識點_2
下面首先來用一個簡單的例子說明一下用法。還是用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]))
# 錯誤在於數據框的統計應該用colMeans
ddply(iris,.(Species),function(df) colMeans(df[1:4]))
# 這個例子並沒有顯示出什么優勢,但我們可以學習到plyr的基礎使用方式，即對於**ply這種格式的函數，其第一個數據參數是數據源，其他參數可以用?來查詢.
初看起來plyr包所具備的功能並不很出彩，下面我們看一個略為復雜例子。還是用iris數據，我們希望對每一種花做一個簡單回歸。
# 首先定義回歸函數
model <- function(x) {
lm(Petal.Length~Petal.Width,data=x)
}
# 如果用普通的函數則需要如下的分割、計算、整合三個步驟共四條命令
# Split函數的作用是將數據框按照指定字段分組，但不做后續計算。lapply函數可以對每組數據都執行同樣的算法。Split和lapply兩者結合可以實現本案例。
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包中還有兩個比較特別的函數，分別是rply和mply，它們分別對應的是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)
下面是其他操作函數
-
each()：each(min, max)等價於function(x) c(min = min(x), max = max(x))。
-
colwise()：colwise(median)將計算列的中位數。
-
arrange()：超級順手的函數，可以方便的給dataframe排序。
-
rename()：又是一個handy的函數，按變量名而不是變量位置重命名。
-
count()：返回unique值，等價於length(unique(**))。
-
match_df()：方便的配合count()等，選出符合條件的行，有點像merge(…,all=F)的感覺。
-
join()：對於習慣SQL的童鞋，可能比merge()用起來更順手吧（當然也更快一點），不過靈活性還是比不上merge()。
知識點_3
pylr包的使用
- （1）命名規則
| array | data frame | list | nothing | |
| array | aaply | adply | alply | a_ply |
| data frame | daply | ddply | dlply | d_ply |
| list | laply | ldply | llply | l_ply |
| n replicates | raply | rdply | rlply | r_ply |
| function arguments |
maply | mdply | mlply | m_ply |
| array | data frame | list | nothing | |
| array | apply | adply | alply | a_ply |
| data frame | daply | aggregate | by | d_ply |
| list | sapply | ldply | lapply | l_ply |
| n replicates | replicate | rdply | replicate | r_ply |
| function arguments |
mapply | mdply | mapply | m_ply |
命名規則：前三行是基本類型。
根據輸入類型和輸出類型：a=array，d=data frame，l=list，_ 表示輸出放棄。第一個字母表示輸入，第2個字母表示輸出。
后兩行是對應apply族的replicates和mapply這兩個函數，分別表示n次重復和多元函數參數的情況，第2個字母還是表示輸出類型。
從命名特點來看，我們不需要列出每個函數的情況了，只要從輸入和輸出兩方面分別討論即可。
-
（2）參數說明
a*ply(.data, .margins, .fun, ..., .progress = "none") d*ply(.data, .variables, .fun, ..., .progress = "none") l*ply(.data, .fun, ..., .progress = "none")
這些函數有兩到三個主要的參數，依賴於輸入的類型：
參數.data是我們要用來分片-計算-合並的；參數.margins或者.variables描述了分片的方式；參數.fun表示用來處理的函數，其它更多的參數（...）是傳遞給處理函數的；參數.progress用來控制顯示一個進度條。
- (3) 輸入
輸入類型有三種，每一種類型給出了如何進行分片的不同方法。
簡單來說：
a*ply( )：數組（包括矩陣和向量）按維數分為低維的片。
d*ply( )：數據框被變量組合分成子集。
l*ply( )：列表的每個元素就是一個分片。
因此，對輸入數據集的分片，不是取決於數據的結構，而是取決於所采用的方法。
1. 一個對象采用 a*ply( )分片必須對應dim( )且接受多維的索引；
2. 采用 d*ply( )分片，要利用split( )並強制轉換為列表；
3. 采用 l*ply( )，需要用 length() and [[。
所以數據框可以被傳遞給 a * ply( )，可以象2維的矩陣那樣處理，也可以傳遞給l*ply( )，被視為一個向量的列表。
三種類型各自的特點：
（a）： 輸入數組（ a*ply() ）
a*ply()的分片特點在於.margins參數，它和apply很相似。對於2維數組， .margins 可以取1，2，或者c(1:2)，對應2維數組的3種分片方式。
• .margins = 1: Slice up into rows.
• .margins = 2: Slice up into columns.
• .margins = c(1,2): Slice up into individual cells.
對於2維數組，則有3種分片方式：
對於3維數組，則有7種分片方式：
.margins對應更高維的情況，可能會面臨一種爆發式的組合。
（b）輸入數據框（ d*ply() ）
使用d*ply時，需要特別指定分組所用的變量或變量函數，它們會被首先計算，然后才是整個數據框。
有下面幾種指定方式：
• .(var1)。按照變量var1的值來對數據框分組
• 多重變量 .(a,b,c)。將按照三個變量的交互值來分組。
這種形式輸出的時候，有點復雜。如果輸出為數組，則數組會有三個維度，分別以 a,b,c 的值作為維數名。如果輸出為數據框，將會包含 a,b,c 取值的三個額外的列。如果輸出為列表，則列表元素名為按周期分割的 a,b,c 的值。
• 作為列名的字符向量：c("var1", "var2")。
• 公式~ var1 + var2。
(c) 輸入列表（ l*ply() ）
l * ply() 不需要描述如何分片的函數，因為列表本身就是按照元素的分划。使用l * ply()相當於a*ply()作用於一維數組的效果。
知識點_4
plyr包案例集
本文給出幾個關於使用plyr包來進行數據處理的例子。重點在於plyr包的應用，而不是深入地探討數據的分析。
- 1. 點球成金
此例來源於 Wickham 的論文[ Hadley Wickham ：The Split-Apply-Combine Strategy for Data Analysis Journal of Statistical Software，April 2011, Volume 40, Issue 1. ]。案例里的圖形是我作的。
布拉德.皮特在《點球成金》里用數據方法發掘棒球運動員的價值，我們也可以來試試。plyr 包的baseball數據集包括了1887-2007年間1228位美國職業棒球運動員15年以上的擊球記錄。
> library(plyr)
> data(baseball)
> dim(baseball)
[1] 21699 22
下面的分析，是要研究擊球手的表現。
關注四個變量：
id：運動員的身份；
year：紀錄的年份；
rbi：跑壘得分，運動員在賽季內的跑動數目；
ad：輪到擊球或者面對投手的次數。
先來計算“職業年份”，就是運動員開始職業生涯以來的年數。
對某一位運動員：
> baberuth <- subset(baseball, id == "ruthba01")
> # 在baberuth里加入cyear
> baberuth <- transform(baberuth, cyear = year - min(year) + 1)
對所有的運動員：
> baseball <- ddply(baseball, .(id), transform, cyear = year - min(year) + 1)
先看看Babe Ruth這位運動員的模式，畫rbi/ab（每次擊球的跑動）的時間序列圖：
> library(ggplot2)
> cyear <- baberuth$cyear
> ra <- (baberuth$rbi)/(baberuth$ab)
> a <- data.frame(cyear,ra)
> p <- ggplot(a,aes(cyear,ra))
> p + geom_line()
> # 或者（畫圖代碼）
> # qplot(cyear, rbi / ab, data = baberuth, geom = "line")
下面使用d_ply函數按下面的方式對各位運動員們畫圖並保存為pdf格式：
（1）按照運動員的平均rbi / ab值對運動員排序，這樣確保如果某個圖畫不出的時候可以畫出其它的。
（2）只取ab大於25的運動員
> baseball <- subset(baseball, ab >= 25)
> xlim <- range(baseball$cyear, na.rm=TRUE)
> ylim <- range(baseball$rbi / baseball$ab, na.rm=TRUE)
> plotpattern <- function(df) {
+ qplot(cyear, rbi / ab, data = df, geom = "line", xlim = xlim, ylim = ylim)
+ }
> pdf("paths.pdf", width = 8, height = 4)
> d_ply(baseball, .(reorder(id, rbi / ab)), failwith(NA, plotpattern), .print = TRUE)
> dev.off()
這些圖里能看出大概的線性趨勢。來對每個運動員擬合一下線性模型。先看Babe Ruth的：
> model <- function(df) {
+ lm(rbi / ab ~ cyear, data = df)
+ }
> model(baberuth)
Call:
lm(formula = rbi/ab ~ cyear, data = df)
Coefficients:
(Intercept) cyear
0.203200 0.003413
然后用dlply對所有的運動員：
> bmodels <- dlply(baseball, .(id), model)
這樣共有1145個模型的列表，下面提取這些模型的系數（slope and intercept）和R^2，
> rsq <- function(x) summary(x)$r.squared
> bcoefs <- ldply(bmodels, function(x) c(coef(x), rsquare = rsq(x)))
> names(bcoefs)[2:3] <- c("intercept", "slope")
> head(bcoefs)
>
id intercept slope rsquare
1 aaronha01 0.18329371 0.0001478121 0.000862425
2 abernte02 0.00000000 NA 0.000000000
3 adairje01 0.08670449 -0.0007118756 0.010230121
4 adamsba01 0.05905337 0.0012002168 0.030184694
5 adamsbo03 0.08867684 -0.0019238835 0.108372596
6 adcocjo01 0.14564821 0.0027382939 0.229040266
把所有模型的R^2繪直方圖：（代碼有問題？？？沒畫出來？？？待修改）
> qplot(rsquare, data=bcoefs, geom="histogram", binwidth=0.01)
可以看出這些模型整體來說擬合的效果是很差的。
下面把那些擬合效果好的模型找出來：先利用系數進行“合並”，再挑出R^2為1的模型：
> baseballcoef <- merge(baseball, bcoefs, by = "id")
> subset(baseballcoef, rsquare > 0.999)$id
>
[1] "bannifl01" "bannifl01" "bedrost01" "bedrost01"
[5] "burbada01" "burbada01" "carrocl02" "carrocl02"
[9] "cookde01" "cookde01" "davisma01" "davisma01"
[13] "jacksgr01" "jacksgr01" "lindbpa01" "lindbpa01"
[17] "oliveda02" "oliveda02" "penaal01" "penaal01"
[21] "powerte01" "powerte01" "splitpa01" "splitpa01"
[25] "violafr01" "violafr01" "wakefti01" "wakefti01"
[29] "weathda01" "weathda01" "woodwi01" "woodwi01"
所有這些擬合好的模型都只有兩個數據點。
換個角度看一下這些模型。
繪制斜率和截距的散點圖：（代碼有問題？？？沒畫出來？？？待修改）
> p <- ggplot(bcoefs, aes(slope, intercept))
> p <- p + geom_point(aes( size = rsquare), alpha = 0.6)
放大局部：（代碼有問題？？？沒畫出來？？？待修改）
> s <- subset(baseballcoef,slope>-0.01&slope<0.01)
> p <- ggplot(s, aes(slope, intercept))
> p <- p + geom_point(aes( size = rsquare), alpha = 0.6)
可以看出，斜率和截距之間是負相關的，且擬合壞的模型其估計值都接近於0。
- 2. 起個好名字
plyr包還有很多對於數據集操作的函數。本文涉及summarise，arrange，mutate
這個案例來自Wickham的課程主頁（http://stat405.had.co.nz/）。
數據集Baby name包括了美國從1880年到2008年間排名前1000位的嬰兒姓名。數據集在課程主頁下載。
library(plyr)
library(ggplot2)
bnames <- read.csv("bnames2.csv.bz2")
births <- read.csv("births.csv")
# 先取一個名字看一看，
steve<-subset(bnames,name == "Steve")
qplot(year, prop, data = steve, geom = "line")
蘋果教主生於1955年，看來是個當時的流行名字。
# 在steve里求概括統計量
summarise(steve,least = year[prop == min(prop)],most = year[prop == max(prop)])
least most
1 2008 1959
# 按比例降序重排steve
arrange(steve, desc(prop))
# 在steve中增加一列per1000，取1000 * prop的整數
mutate(steve, per1000 = round(1000 * prop))
下面來探索一下名字的趨勢，
這涉及到使用外部數據集，下面合並bnames和 births
bnames2 <- join(bnames, births, by = c("year", "sex"))
bnames2 <- mutate(bnames2, n = round(prop * births))
看看出生趨勢
qplot(year, births, data = births, geom = "line",color = sex)
# 那些年，有多少個steve出生呢
steve1<- subset(bnames2, name == "Steve")
summarise(steve1, sum(n))
sum(n)
1 237613
# 對每個名字都算一下：（這就是split-apply-combine步驟了）
counts <- ddply(bnames2, "name", summarise, n=sum(n))
把名字按數量重排一下：
best<-arrange(counts, desc(n))
大眾款的名字是：
head(best)
name n
1 James 5043259
2 John 5036828
3 Robert 4771447
4 Michael 4226596
5 Mary 4111514
6 William 3966170
稀有版：
tail(best)
name n
6777 Delina 4
6778 Delle 4
6779 Dema 4
6780 Dollye 4
6781 Eithel 4
6782 Elzada 4