R数据分析当中的化整为零(Split-Apply-Combine)策略

本文心得自:The Split-Apply-Combine Strategy for Data Analysis, Hadley Wickham, Journal of Statistical Software, April 2011, V.40.

引子:
我们常常会遇到这样的问题,数据量很大,并不需要依顺序来依次处理。合理分块处理,并最终整合起来是一个不错的选择。这也就是所谓的Split-Apply-Combine Strategy策略。这在速度上会有比做一个loop有优势,因为它可以并行处理数据。

什么时候我们需要使用到化整为零的策略呢?有以下三种情况:

  1. 数据需要分组处理
  2. 数据需要按照每行或者每列来处理
  3. 数据需要分级处理,和分组很类似,但是分级时需要考虑分级之间的关系。

化整为零策略有点类似于由Google推广的map-reduce策略。当然map-reduce策略的基础是网格,而这里的Split-Apply-Combine的基础完全可以是单机,甚至不支持并行处理的单机都可以。

然而,化整为零并不是一个很直观的编程过程。最直观的过程是使用Loop循环。这里使用一个例子来讲解一下如何实现化整为零策略。在plyr包中有数据ozone,它是一个三维矩阵(24X24X72),其中最后一维72是指的6年12个月每个月的结果。也就是ozone是一个包括了连续72个月24X24的三维矩阵数据。三维分别是lat,long,time。我们需要由对时间robust linear model之后的残基residuals。

> library(plyr) # need for dataset ozone
> library(MASS) # need for function rlm
> month <- ordered(rep(1:12, length=72)) #set time sequence
> #try one set
> one <- ozone[1,1,]
> model <- rlm(one ~ month - 1); model
Call:
rlm(formula = one ~ month - 1)
Converged in 9 iterations
 
Coefficients:
  month1   month2   month3   month4   month5   month6   month7   month8   month9  month10  month11  month12 
264.3964 259.2036 255.0000 252.0052 258.5089 265.3387 274.0000 276.6724 277.0000 285.0000 283.6036 273.1964 
 
Degrees of freedom: 72 total; 60 residual
Scale estimate: 4.45 
> deseas <- resid(model)


现在我们对每一组数据都做相同的处理。首先使用for loop,这样比较直观

> deseasf <- function(value) rlm(value ~ month -1) #the function
> models <- as.list(rep(NA, 24*24)) #prepare the variable
> dim(models) <- c(24, 24)
> deseas <- array(NA, c(24,24,72)) #prepare the variable
> dimnames(deseas) <- dimnames(ozone)
> for (i in seq_len(24)) { #for loop for first dimension
+ 	for(j in seq_len(24)) { #for loop for second dimension
+ 		mod <- deseasf(ozone[i, j, ]) #apply function
+ 		models[[i, j]] <- mod #save data
+ 		deseas[i, j, ] <- resid(mod) #get residure
+ 	}
+ }

接下来我们使用化整为零的策略。因为数据可以分成24X24块来处理,每一块都是单独运算,可以并行处理。而使用for loop,只能一块接一块的处理,在速度上可能没有并行处理来得快。而在R当中,有一系列相关的函数,apply, lapply, sapply, tapply, mapply, sweep。我们先得了解这些函数,然后再来应用它们。最简单的是apply。

其形式是apply(array, margin, function, …)。首先,apply的对象是矩阵array或者matrix。它的第二个参数是指的维度,如果你的array是一个二维矩阵,需要按横排的方式计算每一排的平均值,那么你的第二个参数就应该是1。如果需要按纵列的方式计算每一列的平均值,那么第二个参数就应该是2。当然还可以使用c(1,2)这样的方式来设置第二个参数,就是并行计算每个值。第三个参数是需要应用的函数。之后的…是需要传入函数的其它参数。而apply的返回值就是由function来确定的,它可能是vector, array or list。下面举个例子比较容易理解。

> x<-cbind(x1=3,x2=c(4:1,2:5))
> dimnames(x)[[1]]<-letters[1:8]
> x
  x1 x2
a  3  4
b  3  3
c  3  2
d  3  1
e  3  2
f  3  3
g  3  4
h  3  5
> apply(x,2,mean,trim =.2) #在这里,trim =.2就是mean(x, trim = 0, na.rm = FALSE, ...)函数当中的一个参数。
x1 x2 
 3  3 
> apply(x,1,mean,trim =.2)
  a   b   c   d   e   f   g   h 
3.5 3.0 2.5 2.0 2.5 3.0 3.5 4.0 
> col.sums <- apply(x, 2, sum)
> row.sums <- apply(x, 1, sum)
> rbind(cbind(x, Rtot = row.sums), Ctot = c(col.sums, sum(col.sums)))
     x1 x2 Rtot
a     3  4    7
b     3  3    6
c     3  2    5
d     3  1    4
e     3  2    5
f     3  3    6
g     3  4    7
h     3  5    8
Ctot 24 24   48
> sum.plus.y <- function(x,y){
+ 	sum(x) + y
+ }
> apply(x, 1, sum.plus.y, 3) #使用自定义函数
 a  b  c  d  e  f  g  h 
10  9  8  7  8  9 10 11

理解了apply,就可比较容易地理解lapply, sapply, vapply了。这三者针对的对象是list或者Vector。其形式为

lapply(X, FUN, ...)
sapply(X, FUN, ..., simplify = TRUE, USE.NAMES = TRUE)
vapply(X, FUN, FUN.VALUE, ..., USE.NAMES = TRUE)

我们看到,它没有了apply当中所需要的第二个参数margin,其原因就是输入对象不是array或者matrix,而是list或者Vector。先举个最简单的应用例子。

> x <- list(a = 1:10, beta = exp(-3:3), logic = c(TRUE,FALSE,FALSE,TRUE))
> x
$a
 [1]  1  2  3  4  5  6  7  8  9 10
 
$beta
[1]  0.04978707  0.13533528  0.36787944  1.00000000  2.71828183  7.38905610 20.08553692
 
$logic
[1]  TRUE FALSE FALSE  TRUE
 
> lapply(x,mean)
$a
[1] 5.5
 
$beta
[1] 4.535125
 
$logic
[1] 0.5

也就是说,其x的维度应该是一维的,当然,其下面的元素可以是多维的。那么如果x是一个矩阵呢?我们先来看例子。

> x<-cbind(x1=3,x2=c(4:1,2:5))
> dimnames(x)[[1]]<-letters[1:8]
> x
  x1 x2
a  3  4
b  3  3
c  3  2
d  3  1
e  3  2
f  3  3
g  3  4
h  3  5
> x<-as.data.frame(x)
> as.list(x)
$x1
[1] 3 3 3 3 3 3 3 3
 
$x2
[1] 4 3 2 1 2 3 4 5
 
> lapply(x,function(.ele) mean(.ele))
$x1
[1] 3
 
$x2
[1] 3
 
> sapply(x,mean)
x1 x2 
 3  3 
> vapply(x,mean,1)
x1 x2 
 3  3

从它们的说明文件我们知道,无论你传入的x是什么,它首先做的一步说是使用as.list来将其转换成一个一维的list。所以,一个data.frame传入lapply之后,它的colnames将会转换成list的names,而rownames可能会丢失。比较可知,lapply和sapply的差别在于,lapply的返回值是一个list,而sapply的返回值是一个矩阵。sapply的返回值其实就是在lapply的基础上再使用了simplify2array(x, higher=TRUE)函数,使用其结果变成一个array。为了更清楚地了解sapply和vapply,我们看下面的例子

> i39 <- sapply(3:9, seq)
> i39
[[1]]
[1] 1 2 3
 
[[2]]
[1] 1 2 3 4
 
[[3]]
[1] 1 2 3 4 5
 
[[4]]
[1] 1 2 3 4 5 6
 
[[5]]
[1] 1 2 3 4 5 6 7
 
[[6]]
[1] 1 2 3 4 5 6 7 8
 
[[7]]
[1] 1 2 3 4 5 6 7 8 9
 
> sapply(i39, fivenum)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]  1.0  1.0    1  1.0  1.0  1.0    1
[2,]  1.5  1.5    2  2.0  2.5  2.5    3
[3,]  2.0  2.5    3  3.5  4.0  4.5    5
[4,]  2.5  3.5    4  5.0  5.5  6.5    7
[5,]  3.0  4.0    5  6.0  7.0  8.0    9
> vapply(i39, fivenum,
+        c(Min. = 0, "1st Qu." = 0, Median = 0, "3rd Qu." = 0, Max. = 0))
        [,1] [,2] [,3] [,4] [,5] [,6] [,7]
Min.     1.0  1.0    1  1.0  1.0  1.0    1
1st Qu.  1.5  1.5    2  2.0  2.5  2.5    3
Median   2.0  2.5    3  3.5  4.0  4.5    5
3rd Qu.  2.5  3.5    4  5.0  5.5  6.5    7
Max.     3.0  4.0    5  6.0  7.0  8.0    9

其中,fivenum函数会返回一组数的最小值,四分位低值(lower-hinge),中值,四分位高值(upper-hinge),以及最大值。从上面的比较中,我们很清楚的看到,sapply返回值的排列形式,以list的names为colnames。可以想象,它使用的是按列填充matrix的方式输出的。而vapply是在sapply的基础上,为rownames做出了定义。

除了上面介绍的,还有tapply,mapply,sweep等。它们的定义如下。如果需要了解和掌握它们,需要熟悉上面介绍的apply, lapply, sapply以及vapply。还需要了解split。所以这里就不多加解释了(因为篇幅会很长)。

tapply(X, INDEX, FUN = NULL, ..., simplify = TRUE)
mapply(FUN, ..., MoreArgs = NULL, SIMPLIFY = TRUE,
       USE.NAMES = TRUE)
sweep(x, MARGIN, STATS, FUN="-", check.margin=TRUE, ...)

有了上面的了解,我们可以看一下最前面的for loop可以怎样改写为apply形式了。

> models <- apply(ozone, 1:2, deseasf) #这里相当于for loop当中的for(i in seq_len(24)){for(j in seq_len(24)){mod<-deseasf(ozone[i,j,]); models[[i,j]]<-mod;}}, 但是运算却是并行处理的。
> resids_list <- lapply(models, resid)
> resids <- unlist(resids_list)
> dim(resids) <- c(72, 24, 24)
> deseas <- aperm(resids, c(2, 3, 1))
> dimnames(deseas) <- dimnames(ozone)

如果,我们知道的并不是ozone这样一个24X24X72的已知确切维度的三维数组,假如我们只有一个名为ozonedf的matrix,它的5列分别为lat, long, time, month,和value。我们如果需要做上述的分析应该怎么办呢?在思路上,我们的想法可能会是先从ozonedf出发生成一个类似ozone这样子的数据,然后再使用apply,lapply这样的函数来完成就可以。第一步生成ozone这样子的数据,就是化整为零策略(Split-Apply-Combine)的第一步了。

R> deseasf_df <- function(df) {
+ rlm(value ~ month - 1, data = df)
+ }
R> pieces <- split(ozonedf, list(ozonedf$lat, ozonedf$long))
R> models <- lapply(pieces, deseasf_df)
R> results <- mapply(function(model, df) {
+ cbind(df[rep(1, 72), c("lat", "long")], resid(model))
+ }, models, pieces)
R> deseasdf <- do.call("rbind", results)

我们再来举一个airquality的例子,我们先来构造一个比较复杂的数据。

> aq<-airquality
> aq$source<-"qiuworld"
> head(aq)
  Ozone Solar.R Wind Temp Month Day   source
1    41     190  7.4   67     5   1 qiuworld
2    36     118  8.0   72     5   2 qiuworld
3    12     149 12.6   74     5   3 qiuworld
4    18     313 11.5   62     5   4 qiuworld
5    NA      NA 14.3   56     5   5 qiuworld
6    28      NA 14.9   66     5   6 qiuworld
> aq1<-aq
> aq1$source<-c("gaziou")
> head(aq1)
  Ozone Solar.R Wind Temp Month Day source
1    41     190  7.4   67     5   1 gaziou
2    36     118  8.0   72     5   2 gaziou
3    12     149 12.6   74     5   3 gaziou
4    18     313 11.5   62     5   4 gaziou
5    NA      NA 14.3   56     5   5 gaziou
6    28      NA 14.9   66     5   6 gaziou
> set.seed(123)
> x1<-runif(nrow(aq1),-0.5,0.5)
> head(x1)
[1] -0.21242248  0.28830514 -0.09102308  0.38301740  0.44046728 -0.45444350
> aq1$Temp<-aq1$Temp+x1
> head(aq1$Temp)
[1] 66.78758 72.28831 73.90898 62.38302 56.44047 65.54556
> aq<-rbind(aq,aq1)
> aq<-aq[order(aq$Wind),]
> head(aq)
    Ozone Solar.R Wind     Temp Month Day   source
53     NA      59  1.7 76.00000     6  22 qiuworld
206    NA      59  1.7 76.29892     6  22   gaziou
121   118     225  2.3 94.00000     8  29 qiuworld
274   118     225  2.3 94.14789     8  29   gaziou
126    73     183  2.8 93.00000     9   3 qiuworld
279    73     183  2.8 93.48422     9   3   gaziou

我们的任务就是比较不同source来源的每个Month的平均Temp。思路应该是先把数据按照source和Month分成小块,计算出来其Temp的平均值,然后输出。这就是一个完整而简单的Split-Apply-Combine的过程了。

> pieces <- split(aq, list(aq$source,aq$Month))
> pieces[1:2]
$gaziou.5
    Ozone Solar.R Wind     Temp Month Day source
183   115     223  5.7 78.64711     5  30 gaziou
164     7      NA  6.9 74.45683     5  11 gaziou
154    41     190  7.4 66.78758     5   1 gaziou
184    37     279  7.4 76.46302     5  31 gaziou
155    36     118  8.0 72.28831     5   2 gaziou
180    NA      NA  8.0 57.04407     5  27 gaziou
160    23     299  8.6 65.02811     5   7 gaziou
163    NA     194  8.6 68.95661     5  10 gaziou
166    11     290  9.2 66.17757     5  13 gaziou
165    16     256  9.7 68.95333     5  12 gaziou
173    11      44  9.7 62.45450     5  20 gaziou
174     1       8  9.7 59.38954     5  21 gaziou
176     4      25  9.7 61.14051     5  23 gaziou
167    14     274 10.9 68.07263     5  14 gaziou
157    18     313 11.5 62.38302     5   4 gaziou
169    14     334 11.5 64.39982     5  16 gaziou
172    30     322 11.5 67.82792     5  19 gaziou
170    34     307 12.0 65.74609     5  17 gaziou
177    32      92 12.0 61.49427     5  24 gaziou
181    23      13 12.0 67.09414     5  28 gaziou
156    12     149 12.6 73.90898     5   3 gaziou
168    18      65 13.2 57.60292     5  15 gaziou
161    19      99 13.8 59.39242     5   8 gaziou
158    NA      NA 14.3 56.44047     5   5 gaziou
159    28      NA 14.9 65.54556     5   6 gaziou
179    NA     266 14.9 58.20853     5  26 gaziou
182    45     252 14.9 80.78916     5  29 gaziou
175    11     320 16.6 73.19280     5  22 gaziou
178    NA      66 16.6 57.15571     5  25 gaziou
171     6      78 18.4 56.54206     5  18 gaziou
162     8      19 20.1 61.05144     5   9 gaziou
 
$qiuworld.5
   Ozone Solar.R Wind Temp Month Day   source
30   115     223  5.7   79     5  30 qiuworld
11     7      NA  6.9   74     5  11 qiuworld
1     41     190  7.4   67     5   1 qiuworld
31    37     279  7.4   76     5  31 qiuworld
2     36     118  8.0   72     5   2 qiuworld
27    NA      NA  8.0   57     5  27 qiuworld
7     23     299  8.6   65     5   7 qiuworld
10    NA     194  8.6   69     5  10 qiuworld
13    11     290  9.2   66     5  13 qiuworld
12    16     256  9.7   69     5  12 qiuworld
20    11      44  9.7   62     5  20 qiuworld
21     1       8  9.7   59     5  21 qiuworld
23     4      25  9.7   61     5  23 qiuworld
14    14     274 10.9   68     5  14 qiuworld
4     18     313 11.5   62     5   4 qiuworld
16    14     334 11.5   64     5  16 qiuworld
19    30     322 11.5   68     5  19 qiuworld
17    34     307 12.0   66     5  17 qiuworld
24    32      92 12.0   61     5  24 qiuworld
28    23      13 12.0   67     5  28 qiuworld
3     12     149 12.6   74     5   3 qiuworld
15    18      65 13.2   58     5  15 qiuworld
8     19      99 13.8   59     5   8 qiuworld
5     NA      NA 14.3   56     5   5 qiuworld
6     28      NA 14.9   66     5   6 qiuworld
26    NA     266 14.9   58     5  26 qiuworld
29    45     252 14.9   81     5  29 qiuworld
22    11     320 16.6   73     5  22 qiuworld
25    NA      66 16.6   57     5  25 qiuworld
18     6      78 18.4   57     5  18 qiuworld
9      8      19 20.1   61     5   9 qiuworld
 
> avgTemp<-lapply(pieces,function(.ele) mean(.ele$Temp))
> head(avgTemp)
$gaziou.5
[1] 65.63339
 
$qiuworld.5
[1] 65.54839
 
$gaziou.6
[1] 79.02871
 
$qiuworld.6
[1] 79.1
 
$gaziou.7
[1] 83.90313
 
$qiuworld.7
[1] 83.90323
 
> avgTemp<-do.call(rbind,avgTemp)
> head(avgTemp)
               [,1]
gaziou.5   65.63339
qiuworld.5 65.54839
gaziou.6   79.02871
qiuworld.6 79.10000
gaziou.7   83.90313
qiuworld.7 83.90323

由上面的过程我们可以看来,其实化整为零策略在实现起来,就是分三步走,使用split将数据化分成小块,使用lapply函数对小块进行计算,最后使用do.call使用函数将其整理成我们需要的形式。这个过程,其实使用plyr包来实现,就更为便洁了。同样是上面的操作,使用plyr的话,只需要一行即可。

> avgTemp1<-ddply(aq,.(source,Month),function(.ele) mean(.ele$Temp))
> avgTemp1
     source Month       V1
1    gaziou     5 65.63339
2    gaziou     6 79.02871
3    gaziou     7 83.90313
4    gaziou     8 83.98611
5    gaziou     9 76.89319
6  qiuworld     5 65.54839
7  qiuworld     6 79.10000
8  qiuworld     7 83.90323
9  qiuworld     8 83.96774
10 qiuworld     9 76.90000

plyr包给我们提供了非常简洁的书写方式。我们来看看其主要函数的定义方式。

Input\output Array Data frame List Discarded
Array aaply adply alply a_ply
Data frame daply ddply dlply d_ply
List laply ldply llply l_ply
a*ply(.data, .margins, .fun, ..., .progress = "none")
d*ply(.data, .variables, .fun, ..., .progress = "none")
l*ply(.data, .fun, ..., .progress = "none")

我们可以看出,函数名称的第一个字母代表输入的形式,它们分别是a->Array, d->Data frame, l->List。而第二字母代表输出的形式,它们的定义同前。对于输入为array和data frame的,函数的第二个参数为data的margins或者variables。对于margins,可以是

  • .margins = 1 #以行为单位
  • .margins = 2 #以列为单位
  • .margins = c(1,2) #以individual cell为单位

需要注意的是,这里的每一个参数都使用了.点号起始。在作者看来,这样可以避免误操作。

2 thoughts on “R数据分析当中的化整为零(Split-Apply-Combine)策略

发表评论

电子邮件地址不会被公开。 必填项已用*标注