6.10 数据聚合
分组求和 https://stackoverflow.com/questions/1660124
主要是分组统计
apropos("apply")## [1] "apply" "dendrapply" "eapply" "frollapply" "kernapply"
## [6] "lapply" "mapply" "rapply" "sapply" "tapply"
## [11] "vapply"# 分组求和 colSums colMeans max
unique(iris$Species)## [1] setosa versicolor virginica
## Levels: setosa versicolor virginica# 分类求和
# colSums(iris[iris$Species == "setosa", -5])
# colSums(iris[iris$Species == "virginica", -5])
colSums(iris[iris$Species == "versicolor", -5])## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 296.8 138.5 213.0 66.3# apply(iris[iris$Species == "setosa", -5], 2, sum)
# apply(iris[iris$Species == "setosa", -5], 2, mean)
# apply(iris[iris$Species == "setosa", -5], 2, min)
# apply(iris[iris$Species == "setosa", -5], 2, max)
apply(iris[iris$Species == "setosa", -5], 2, quantile)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0% 4.3 2.300 1.000 0.1
## 25% 4.8 3.200 1.400 0.2
## 50% 5.0 3.400 1.500 0.2
## 75% 5.2 3.675 1.575 0.3
## 100% 5.8 4.400 1.900 0.6aggregate: Compute Summary Statistics of Data Subsets
# 按分类变量 Species 分组求和
# aggregate(subset(iris, select = -Species), by = list(iris[, "Species"]), FUN = sum)
aggregate(iris[, -5], list(iris[, 5]), sum)## Group.1 Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 setosa 250.3 171.4 73.1 12.3
## 2 versicolor 296.8 138.5 213.0 66.3
## 3 virginica 329.4 148.7 277.6 101.3# 先确定位置，假设有很多分类变量
ind <- which("Species" == colnames(iris))
# 分组统计
aggregate(iris[, -ind], list(iris[, ind]), sum)## Group.1 Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 setosa 250.3 171.4 73.1 12.3
## 2 versicolor 296.8 138.5 213.0 66.3
## 3 virginica 329.4 148.7 277.6 101.3按照 Species 划分的类别，分组计算，使用公式表示形式，右边一定是分类变量，否则会报错误或者警告，输出奇怪的结果，请读者尝试运行aggregate(Species ~ Sepal.Length, data = iris, mean)。公式法表示分组计算，~ 左手边可以做加 + 减 - 乘 * 除 / 取余 %% 等数学运算。下面以数据集 iris 为例，只对 Sepal.Length 按 Species 分组计算
aggregate(Sepal.Length ~ Species, data = iris, mean)## Species Sepal.Length
## 1 setosa 5.006
## 2 versicolor 5.936
## 3 virginica 6.588与上述分组统计结果一样的命令，在大数据集上， 与 aggregate 相比，tapply 要快很多，by 是 tapply 的包裹，处理速度差不多。读者可以构造伪随机数据集验证。
# tapply(iris$Sepal.Length, list(iris$Species), mean)
with(iris, tapply(Sepal.Length, Species, mean))## setosa versicolor virginica
## 5.006 5.936 6.588by(iris$Sepal.Length, iris$Species, mean)## iris$Species: setosa
## [1] 5.006
## ------------------------------------------------------------
## iris$Species: versicolor
## [1] 5.936
## ------------------------------------------------------------
## iris$Species: virginica
## [1] 6.588对所有变量按 Species 分组计算
aggregate(. ~ Species, data = iris, mean)## Species Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 setosa 5.006 3.428 1.462 0.246
## 2 versicolor 5.936 2.770 4.260 1.326
## 3 virginica 6.588 2.974 5.552 2.026对变量 Sepal.Length 和 Sepal.Width 求和后，按 Species 分组计算
aggregate(Sepal.Length + Sepal.Width ~ Species, data = iris, mean)## Species Sepal.Length + Sepal.Width
## 1 setosa 8.434
## 2 versicolor 8.706
## 3 virginica 9.562对多个分类变量做分组计算，在数据集 ChickWeight 中 Chick和Diet都是数字编码的分类变量，其中 Chick 是有序的因子变量，Diet 是无序的因子变量，而 Time 是数值型的变量，表示小鸡出生的天数。
# 查看数据
str(ChickWeight)## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame': 578 obs. of 4 variables:
## $ weight: num 42 51 59 64 76 93 106 125 149 171 ...
## $ Time : num 0 2 4 6 8 10 12 14 16 18 ...
## $ Chick : Ord.factor w/ 50 levels "18"<"16"<"15"<..: 15 15 15 15 15 15 15 15 15 15 ...
## $ Diet : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "formula")=Class 'formula' language weight ~ Time | Chick
## .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## - attr(*, "outer")=Class 'formula' language ~Diet
## .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## - attr(*, "labels")=List of 2
## ..$ x: chr "Time"
## ..$ y: chr "Body weight"
## - attr(*, "units")=List of 2
## ..$ x: chr "(days)"
## ..$ y: chr "(gm)"查看数据集ChickWeight的前几行
head(ChickWeight)## weight Time Chick Diet
## 1 42 0 1 1
## 2 51 2 1 1
## 3 59 4 1 1
## 4 64 6 1 1
## 5 76 8 1 1
....str(ChickWeight)## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame': 578 obs. of 4 variables:
## $ weight: num 42 51 59 64 76 93 106 125 149 171 ...
## $ Time : num 0 2 4 6 8 10 12 14 16 18 ...
## $ Chick : Ord.factor w/ 50 levels "18"<"16"<"15"<..: 15 15 15 15 15 15 15 15 15 15 ...
## $ Diet : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "formula")=Class 'formula' language weight ~ Time | Chick
....对于数据集ChickWeight中的有序变量Chick，aggregate 会按照既定顺序返回分组计算的结果
aggregate(weight ~ Chick, data = ChickWeight, mean)## Chick weight
## 1 18 37.00000
## 2 16 49.71429
## 3 15 60.12500
## 4 13 67.83333
## 5 9 81.16667
....aggregate(weight ~ Diet, data = ChickWeight, mean)## Diet weight
## 1 1 102.6455
## 2 2 122.6167
## 3 3 142.9500
## 4 4 135.2627分类变量没有用数字编码，以 CO2 数据集为例，该数据集描述草植对二氧化碳的吸收情况，Plant 是具有12个水平的有序的因子变量，Type表示植物的源头分别是魁北克(Quebec)和密西西比(Mississippi)，Treatment表示冷却(chilled)和不冷却(nonchilled)两种处理方式，conc表示周围环境中二氧化碳的浓度，uptake表示植物吸收二氧化碳的速率。
# 查看数据集
head(CO2)## Plant Type Treatment conc uptake
## 1 Qn1 Quebec nonchilled 95 16.0
## 2 Qn1 Quebec nonchilled 175 30.4
## 3 Qn1 Quebec nonchilled 250 34.8
## 4 Qn1 Quebec nonchilled 350 37.2
## 5 Qn1 Quebec nonchilled 500 35.3
## 6 Qn1 Quebec nonchilled 675 39.2str(CO2)## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame': 84 obs. of 5 variables:
## $ Plant : Ord.factor w/ 12 levels "Qn1"<"Qn2"<"Qn3"<..: 1 1 1 1 1 1 1 2 2 2 ...
## $ Type : Factor w/ 2 levels "Quebec","Mississippi": 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment: Factor w/ 2 levels "nonchilled","chilled": 1 1 1 1 1 1 1 1 1 1 ...
## $ conc : num 95 175 250 350 500 675 1000 95 175 250 ...
## $ uptake : num 16 30.4 34.8 37.2 35.3 39.2 39.7 13.6 27.3 37.1 ...
## - attr(*, "formula")=Class 'formula' language uptake ~ conc | Plant
## .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## - attr(*, "outer")=Class 'formula' language ~Treatment * Type
## .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv>
## - attr(*, "labels")=List of 2
## ..$ x: chr "Ambient carbon dioxide concentration"
## ..$ y: chr "CO2 uptake rate"
## - attr(*, "units")=List of 2
## ..$ x: chr "(uL/L)"
## ..$ y: chr "(umol/m^2 s)"对单个变量分组统计
aggregate(uptake ~ Plant, data = CO2, mean)## Plant uptake
## 1 Qn1 33.22857
## 2 Qn2 35.15714
## 3 Qn3 37.61429
## 4 Qc1 29.97143
## 5 Qc3 32.58571
## 6 Qc2 32.70000
## 7 Mn3 24.11429
## 8 Mn2 27.34286
## 9 Mn1 26.40000
## 10 Mc2 12.14286
## 11 Mc3 17.30000
## 12 Mc1 18.00000aggregate(uptake ~ Type, data = CO2, mean)## Type uptake
## 1 Quebec 33.54286
## 2 Mississippi 20.88333aggregate(uptake ~ Treatment, data = CO2, mean)## Treatment uptake
## 1 nonchilled 30.64286
## 2 chilled 23.78333对多个变量分组统计，查看二氧化碳吸收速率uptake随类型Type和处理方式Treatment
aggregate(uptake ~ Type + Treatment, data = CO2, mean)## Type Treatment uptake
## 1 Quebec nonchilled 35.33333
## 2 Mississippi nonchilled 25.95238
## 3 Quebec chilled 31.75238
## 4 Mississippi chilled 15.81429tapply(CO2$uptake, list(CO2$Type, CO2$Treatment), mean)## nonchilled chilled
## Quebec 35.33333 31.75238
## Mississippi 25.95238 15.81429by(CO2$uptake, list(CO2$Type, CO2$Treatment), mean)## : Quebec
## : nonchilled
## [1] 35.33333
## ------------------------------------------------------------
## : Mississippi
## : nonchilled
## [1] 25.95238
## ------------------------------------------------------------
## : Quebec
## : chilled
## [1] 31.75238
## ------------------------------------------------------------
## : Mississippi
## : chilled
## [1] 15.81429在这个例子中 tapply 和 by 的输出结果的表示形式不一样，aggregate 返回一个 data.frame 数据框，tapply 返回一个表格 table，by 返回特殊的数据类型 by。
Function by is an object-oriented wrapper for tapply applied to data frames.
# 分组求和
# by(iris[, 1], INDICES = list(iris$Species), FUN = sum)
# by(iris[, 2], INDICES = list(iris$Species), FUN = sum)
by(iris[, 3], INDICES = list(iris$Species), FUN = sum)## : setosa
## [1] 73.1
## ------------------------------------------------------------
## : versicolor
## [1] 213
## ------------------------------------------------------------
## : virginica
## [1] 277.6by(iris[1:3], INDICES = list(iris$Species), FUN = sum)## : setosa
## [1] 494.8
## ------------------------------------------------------------
## : versicolor
## [1] 648.3
## ------------------------------------------------------------
## : virginica
## [1] 755.7by(iris[1:3], INDICES = list(iris$Species), FUN = summary)## : setosa
## Sepal.Length Sepal.Width Petal.Length
## Min. :4.300 Min. :2.300 Min. :1.000
## 1st Qu.:4.800 1st Qu.:3.200 1st Qu.:1.400
## Median :5.000 Median :3.400 Median :1.500
## Mean :5.006 Mean :3.428 Mean :1.462
## 3rd Qu.:5.200 3rd Qu.:3.675 3rd Qu.:1.575
## Max. :5.800 Max. :4.400 Max. :1.900
## ------------------------------------------------------------
## : versicolor
## Sepal.Length Sepal.Width Petal.Length
## Min. :4.900 Min. :2.000 Min. :3.00
## 1st Qu.:5.600 1st Qu.:2.525 1st Qu.:4.00
## Median :5.900 Median :2.800 Median :4.35
## Mean :5.936 Mean :2.770 Mean :4.26
## 3rd Qu.:6.300 3rd Qu.:3.000 3rd Qu.:4.60
## Max. :7.000 Max. :3.400 Max. :5.10
## ------------------------------------------------------------
## : virginica
## Sepal.Length Sepal.Width Petal.Length
## Min. :4.900 Min. :2.200 Min. :4.500
## 1st Qu.:6.225 1st Qu.:2.800 1st Qu.:5.100
## Median :6.500 Median :3.000 Median :5.550
## Mean :6.588 Mean :2.974 Mean :5.552
## 3rd Qu.:6.900 3rd Qu.:3.175 3rd Qu.:5.875
## Max. :7.900 Max. :3.800 Max. :6.900by(iris, INDICES = list(iris$Species), FUN = summary)## : setosa
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.300 Min. :1.000 Min. :0.100
## 1st Qu.:4.800 1st Qu.:3.200 1st Qu.:1.400 1st Qu.:0.200
## Median :5.000 Median :3.400 Median :1.500 Median :0.200
## Mean :5.006 Mean :3.428 Mean :1.462 Mean :0.246
## 3rd Qu.:5.200 3rd Qu.:3.675 3rd Qu.:1.575 3rd Qu.:0.300
## Max. :5.800 Max. :4.400 Max. :1.900 Max. :0.600
## Species
## setosa :50
## versicolor: 0
## virginica : 0
##
##
##
## ------------------------------------------------------------
## : versicolor
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## Min. :4.900 Min. :2.000 Min. :3.00 Min. :1.000 setosa : 0
## 1st Qu.:5.600 1st Qu.:2.525 1st Qu.:4.00 1st Qu.:1.200 versicolor:50
## Median :5.900 Median :2.800 Median :4.35 Median :1.300 virginica : 0
## Mean :5.936 Mean :2.770 Mean :4.26 Mean :1.326
## 3rd Qu.:6.300 3rd Qu.:3.000 3rd Qu.:4.60 3rd Qu.:1.500
## Max. :7.000 Max. :3.400 Max. :5.10 Max. :1.800
## ------------------------------------------------------------
## : virginica
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.900 Min. :2.200 Min. :4.500 Min. :1.400
## 1st Qu.:6.225 1st Qu.:2.800 1st Qu.:5.100 1st Qu.:1.800
## Median :6.500 Median :3.000 Median :5.550 Median :2.000
## Mean :6.588 Mean :2.974 Mean :5.552 Mean :2.026
## 3rd Qu.:6.900 3rd Qu.:3.175 3rd Qu.:5.875 3rd Qu.:2.300
## Max. :7.900 Max. :3.800 Max. :6.900 Max. :2.500
## Species
## setosa : 0
## versicolor: 0
## virginica :50
##
##
## Group Averages Over Level Combinations of Factors 分组平均
str(warpbreaks)## 'data.frame': 54 obs. of 3 variables:
## $ breaks : num 26 30 54 25 70 52 51 26 67 18 ...
## $ wool : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
## $ tension: Factor w/ 3 levels "L","M","H": 1 1 1 1 1 1 1 1 1 2 ...head(warpbreaks)## breaks wool tension
## 1 26 A L
## 2 30 A L
## 3 54 A L
## 4 25 A L
## 5 70 A L
## 6 52 A Lave(warpbreaks$breaks, warpbreaks$wool)## [1] 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704
## [9] 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704
## [17] 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704 31.03704
## [25] 31.03704 31.03704 31.03704 25.25926 25.25926 25.25926 25.25926 25.25926
## [33] 25.25926 25.25926 25.25926 25.25926 25.25926 25.25926 25.25926 25.25926
## [41] 25.25926 25.25926 25.25926 25.25926 25.25926 25.25926 25.25926 25.25926
## [49] 25.25926 25.25926 25.25926 25.25926 25.25926 25.25926with(warpbreaks, ave(breaks, tension, FUN = function(x) mean(x, trim = 0.1)))## [1] 35.6875 35.6875 35.6875 35.6875 35.6875 35.6875 35.6875 35.6875 35.6875
## [10] 26.3125 26.3125 26.3125 26.3125 26.3125 26.3125 26.3125 26.3125 26.3125
## [19] 21.0625 21.0625 21.0625 21.0625 21.0625 21.0625 21.0625 21.0625 21.0625
## [28] 35.6875 35.6875 35.6875 35.6875 35.6875 35.6875 35.6875 35.6875 35.6875
## [37] 26.3125 26.3125 26.3125 26.3125 26.3125 26.3125 26.3125 26.3125 26.3125
## [46] 21.0625 21.0625 21.0625 21.0625 21.0625 21.0625 21.0625 21.0625 21.0625# 分组求和
with(warpbreaks, ave(breaks, tension, FUN = function(x) sum(x)))## [1] 655 655 655 655 655 655 655 655 655 475 475 475 475 475 475 475 475 475 390
## [20] 390 390 390 390 390 390 390 390 655 655 655 655 655 655 655 655 655 475 475
## [39] 475 475 475 475 475 475 475 390 390 390 390 390 390 390 390 390# 分组求和
with(iris, ave(Sepal.Length, Species, FUN = function(x) sum(x)))## [1] 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3
## [13] 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3
## [25] 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3
## [37] 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3 250.3
## [49] 250.3 250.3 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8
## [61] 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8
## [73] 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8
## [85] 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8 296.8
## [97] 296.8 296.8 296.8 296.8 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4
## [109] 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4
## [121] 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4
## [133] 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4 329.4
## [145] 329.4 329.4 329.4 329.4 329.4 329.4