来自: 糗世界
本文心得自：The Split-Apply-Combine Strategy for Data Analysis, Hadley Wickham, Journal of Statistical Software, April 2011, V.40.
引子：
我们常常会遇到这样的问题，数据量很大，并不需要依顺序来依次处理。合理分块处理，并最终整合起来是一个不错的选择。这也就是所谓的Split-Apply-Combine Strategy策略。这在速度上会有比做一个loop有优势，因为它可以并行处理数据。
什么时候我们需要使用到化整为零的策略呢？有以下三种情况：
- 数据需要分组处理
- 数据需要按照每行或者每列来处理
- 数据需要分级处理，和分组很类似，但是分级时需要考虑分级之间的关系。
化整为零策略有点类似于由Google推广的map-reduce策略。当然map-reduce策略的基础是网格，而这里的Split-Apply-Combine的基础完全可以是单机，甚至不支持并行处理的单机都可以。
然而，化整为零并不是一个很直观的编程过程。最直观的过程是使用Loop循环。这里使用一个例子来讲解一下如何实现化整为零策略。在plyr包中有数据ozone,它是一个三维矩阵(24X24X72)，其中最后一维72是指的6年12个月每个月的结果。也就是ozone是一个包括了连续72个月24Ｘ24的三维矩阵数据。三维分别是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包给我们提供了非常简洁的书写方式。我们来看看其主要函数的定义方式。
| Inputoutput | 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为单位
需要注意的是，这里的每一个参数都使用了.点号起始。在作者看来，这样可以避免误操作。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
