From 《Data Manipulation with R》- table系列函数
对于构建简单表格或者列联表，table函数可以实现这样的功能。对于更复杂的需求，可以联合apply、sweep、mapply、sapply、lapply等函数对数据框进行操作。 1.1 table table函数中的参数可以是向量vector，或者是列表list、数据框data frame。该函数返回的值是与一个数组，数组的维度与向量中不同元素个数一致，数组维度名即是用于作表的元素名。 pets = c("dog", "cat", "duck", "chicken", "duck", "cat", "dog")
(tt = table(pets))
## pets
## cat chicken dog duck
## 2 1 2 2
tt["duck"]
## duck
## 2
table函数结果可以转换为数据框 as.data.frame(tt)
## pets Freq
## 1 cat 2
## 2 chicken 1
## 3 dog 2
## 4 duck 2
table函数中的参数可以是多个向量，下例将引用state.x77和state.region两个数据 hiinc = state.x77[, "Income"] > median(state.x77[, "Income"])
stateinc = table(state.region, hiinc)
stateinc
## hiinc
## state.region FALSE TRUE
## Northeast 4 5
## South 12 4
## North Central 5 7
## West 4 9
将stateinc转换为数据框 as.data.frame(stateinc)
## state.region hiinc Freq
## 1 Northeast FALSE 4
## 2 South FALSE 12
## 3 North Central FALSE 5
## 4 West FALSE 4
## 5 Northeast TRUE 5
## 6 South TRUE 4
## 7 North Central TRUE 7
## 8 West TRUE 9
如果table函数中的参数引用的是数据框格式，那么table将视数据框每一列为单独的变量，例如： x = data.frame(a = c(1, 2, 2, 1, 2, 2, 1), b = c(1, 2, 2, 1, 1, 2, 1), c = c(1,
1, 2, 1, 2, 2, 1))
x
## a b c
## 1 1 1 1
## 2 2 2 1
## 3 2 2 2
## 4 1 1 1
## 5 2 1 2
## 6 2 2 2
## 7 1 1 1
as.data.frame(table(x))
## a b c Freq
## 1 1 1 1 3
## 2 2 1 1 0
## 3 1 2 1 0
## 4 2 2 1 1
## 5 1 1 2 0
## 6 2 1 2 1
## 7 1 2 2 0
## 8 2 2 2 2
1.2 addmargins 有时候需要显示列联表每一行或每一列的汇总情况，addmargins函数可以实现该功能；如果要指定某一个维度(行或列)增加汇总，可以使用margin参数，margin=1表示增加行汇总，margin=2表示增加列汇总。使用sum函数时会默认创建汇总，如果是其他计算函数，可以增加FUN参数。例如： # 首先,创建一个列联表
tt = table(infert$education, infert$parity)
tt
##
## 1 2 3 4 5 6
## 0-5yrs 3 0 0 3 0 6
## 6-11yrs 42 42 21 12 3 0
## 12+ yrs 54 39 15 3 3 2
# 其次，增加一行汇总
tt1 = addmargins(tt, 1)
tt1
##
## 1 2 3 4 5 6
## 0-5yrs 3 0 0 3 0 6
## 6-11yrs 42 42 21 12 3 0
## 12+ yrs 54 39 15 3 3 2
## Sum 99 81 36 18 6 8
# 增加行列汇总
tt12 = addmargins(tt, c(1, 2))
tt12
##
## 1 2 3 4 5 6 Sum
## 0-5yrs 3 0 0 3 0 6 12
## 6-11yrs 42 42 21 12 3 0 120
## 12+ yrs 54 39 15 3 3 2 116
## Sum 99 81 36 18 6 8 248
dimnames(tt12)
## [[1]]
## [1] "0-5yrs" "6-11yrs" "12+ yrs" "Sum"
##
## [[2]]
## [1] "1" "2" "3" "4" "5" "6" "Sum"
1.3 prop.table 如果列联表中需要利用比率替代计数，可以使用sweep函数，利用计数除以行列的汇总值。prop.table函数可以便捷的得到这样的需求。prop.table函数中margin参数忽略，则所有比率加总和为1；margin=1则每一行比率加总为1；margin=2则每一列比率加总为1。例： prop.table(tt, 2)
##
## 1 2 3 4 5 6
## 0-5yrs 0.0303 0.0000 0.0000 0.1667 0.0000 0.7500
## 6-11yrs 0.4242 0.5185 0.5833 0.6667 0.5000 0.0000
## 12+ yrs 0.5455 0.4815 0.4167 0.1667 0.5000 0.2500
1.4 ftable 对于多于两个维度的表格可以采用flattened模式的表进行性展示，以UCBAdmissions数据为例,该数据有三个维度，table格式是由一系列二维表格呈现的，ftable函数可以更紧凑地展现该数据： ftable(UCBAdmissions)
## Dept A B C D E F
## Admit Gender
## Admitted Male 512 353 120 138 53 22
## Female 89 17 202 131 94 24
## Rejected Male 313 207 205 279 138 351
## Female 19 8 391 244 299 317
1.5 xtable xtable函数可以得到类似于table函数的结果，但是它的参数是公式形式的 xtabs(~state.region + hiinc)
## hiinc
## state.region FALSE TRUE
## Northeast 4 5
## South 12 4
## North Central 5 7
## West 4 9
xtable中公式左边是对公式右边变量的计数，例： x = data.frame(a = c(1, 2, 2, 1, 2, 2, 1), b = c(1, 2, 2, 1, 1, 2, 1), c = c(1,
1, 2, 1, 2, 2, 1))
dfx = as.data.frame(table(x))
xtabs(Freq ~ a + b + c, data = dfx)
## , , c = 1
##
## b
## a 1 2
## 1 3 0
## 2 0 1
##
## , , c = 2
##
## b
## a 1 2
## 1 0 0
## 2 1 2
2.数据聚合(Aggregation)的映射路线 Aggregation主要涉及三个方面： a.数据划分的组别是如何定义的? b.数据的本质是什么? c.理想的最终结果是什么? 根据数据属性、格式，数据分组映射可以有以下几中方式： a.根据list的元素进行分组 b.根据矩阵的行或列进行分组 c.根据分组变量进行分组 3.作用于向量或列表的函数(lapply、sapply、replicate) 尽管R中大部分函数都是直接作用域向量(vector)中每个元素，列表(list)略有不同。很多R函数返回结果的类型是list格式的，R对列表(list)中每个层级(level)的处理又与向量(vector)相同。 3.1 sapply lapply和sapply函数是处理向量和列表的常用函数，前者返回一个列表格式的结果，后者通常返回向量或者矩阵。 3.1.1 sapply可以结合strsplit函数计算字符串中单词的个数： text = c("R is a free environment for statistical analysis", "It compiles and runs on a variety of platforms",
"Visit the R home page for more information")
result = strsplit(text, " ")
result
## [[1]]
## [1] "R" "is" "a" "free" "environment"
## [6] "for" "statistical" "analysis"
##
## [[2]]
## [1] "It" "compiles" "and" "runs" "on" "a"
## [7] "variety" "of" "platforms"
##
## [[3]]
## [1] "Visit" "the" "R" "home" "page"
## [6] "for" "more" "information"
由于字符串向量中每个字符串中都有多个单词，strsplit函数得到的是一个列表格式的结果，我们使用length函数只会得到列表中levels的个数： length(result)
## [1] 3
想要查询每个level的单词个数，我们可以使用lapply或者sapply函数，它们会对列表的每个level进行运算： # nwords = lapply(result,length)#返回list格式
nwords = sapply(result, length) #返回vector格式
nwords
## [1] 8 9 8
3.1.2 sapply对数据框(data frames)同样适用，数据框的列保持原来类型，以ChickWeight数据为例： # class函数分析数据框性质
class(ChickWeight)
## [1] "nfnGroupedData" "nfGroupedData" "groupedData" "data.frame"
# sapply函数分析数据中每个变量的属性
sapply(ChickWeight, class)
## $weight
## [1] "numeric"
##
## $Time
## [1] "numeric"
##
## $Chick
## [1] "ordered" "factor"
##
## $Diet
## [1] "factor"
需要注意的是，因为变量Chick是有序因子类型(ordered & factor)，sapply函数返回的是list格式。如果sapply返回vector或array格式将丢失部分信息。 sapply可以用来剥离数据框中的变量，如创建一个只包含数值型变量的数据框： df[, sapply(df,class)== 'numeric'] 3.1.3 sapply或lapply还可以用于一些重复性工作，作为除loops的另一种选择。 使用这些函数的时候，需要注意返回结果的类型，尽量避免用向量或矩阵存储循环的结果。比如，我们要生成随机矩阵，并计算矩阵中任意变量间最大的相关系数。首先，我们创建一个随机数矩阵，并计算变量间最大的相关系数： maxcor = function(i, n = 10, m = 5) {
mat = matrix(rnorm(n * m), n, m)
corr = cor(mat)
diag(corr) = NA
max(corr, na.rm = TRUE)
}
sapply函数总会向其应用函数参数中传递参数，所以，maxcor中增加了一个哑参数(i),它可以用于自定义maxcor中参数m、n；又因为相关系数矩阵中对角线上总是1，这里预先将对角线元素设为空(NA)。假设生成1000个100×5的矩阵，并计算最大相关系数的平均值: maxcors = sapply(1:1000, maxcor, n = 100)
mean(maxcors)
## [1] 0.153
3.2 replicate replicate函数可以用于一些简单的模拟，该函数包含一个重复次数的参数和一个用于计算模拟数据统计量的表达式(非函数)。例如： t.test(rnorm(10), rnorm(10))$statistic #t检验表达式
## t
## -0.01107
tsim = replicate(10000, t.test(rnorm(10), rnorm(10))$statistic)
quantile(tsim, c(0.5, 0.75, 0.9, 0.95, 0.99))
## 50% 75% 90% 95% 99%
## 0.002184 0.704160 1.385977 1.745731 2.515936
4.作用矩阵和数组的常见函数(apply) R中的apply函数可以对数组或矩阵的各个维度上的数据进行操作。apply包含三个参数：待操作的数组array、指定待操作维度的索引和待运用的函数。类似sapply函数，apply还可以添加额外的一些参数。例如对于矩阵，参数1默认表示对行操作，2默认对列操作。 4.1 apply & scale apply通常与其他函数联合使用。例如，scale与apply联合使用对矩阵的每一列进行描述统计。 sstate = scale(state.x77, center = apply(state.x77, 2, median), scale = apply(state.x77,
2, mad))
类似sapply,apply通常返回向量和矩阵形式的结果。 例，假设我们要构建一个包含非缺失值数目、均值和标准差的矩阵， # 首先需要构建一个作用于列的且满足统计需求的函数：
sumfun = function(x) c(n = sum(!is.na(x)), mean = mean(x), sd = sd(x))
# 进而将该函数应用于整个数据框(矩阵、数组)：
x = apply(state.x77, 2, sumfun)
t(x)
## n mean sd
## Population 50 4246.420 4.464e+03
## Income 50 4435.800 6.145e+02
## Illiteracy 50 1.170 6.095e-01
## Life Exp 50 70.879 1.342e+00
## Murder 50 7.378 3.692e+00
## HS Grad 50 53.108 8.077e+00
## Frost 50 104.460 5.198e+01
## Area 50 70735.880 8.533e+04
例，非重叠分组情况下对向量进行计算，此时可以将该向量视作矩阵，运营apply函数在分组变量上运算。例如，求一向量中每三个相邻数的和： x = 1:12
apply(matrix(x, ncol = 3, byrow = TRUE), 1, sum)
## [1] 6 15 24 33
4.2 rowSums,colSums,rowMeans等 apply函数在计算常见统计量时非常有用，另外还有一些函数专门用于求和、均值等需求的函数：rowSums,colSums,rowMeams等，这些函数通常作用于矩阵(数据框等可以转换为矩阵的格式)，如： mns = colMeans(USJudgeRatings)
mns
## CONT INTG DMNR DILG CFMG DECI PREP FAMI ORAL WRIT PHYS RTEN
## 7.437 8.021 7.516 7.693 7.479 7.565 7.467 7.488 7.293 7.384 7.935 7.602
jscore = rowSums(USJudgeRatings >= 8)
head(jscore)
## AARONSON,L.H. ALEXANDER,J.M. ARMENTANO,A.J. BERDON,R.I. BRACKEN,J.J.
## 1 8 1 11 0
## BURNS,E.B.
## 10
4.3 sweep，mapply sweep函数可以对矩阵每一行或列进行不同的运算，该函数可以添加R中所有内嵌运算符，例如将矩阵所有元素除以每一列的最大值： maxes = apply(state.x77, 2, max)
swept = sweep(state.x77, 2, maxes, "/")
head(swept)
## Population Income Illiteracy Life Exp Murder HS Grad Frost
## Alabama 0.17053 0.5739 0.7500 0.9382 1.0000 0.6137 0.10638
## Alaska 0.01722 1.0000 0.5357 0.9417 0.7483 0.9911 0.80851
## Arizona 0.10435 0.7173 0.6429 0.9586 0.5166 0.8633 0.07979
## Arkansas 0.09954 0.5349 0.6786 0.9601 0.6689 0.5929 0.34574
## California 1.00000 0.8098 0.3929 0.9743 0.6821 0.9302 0.10638
## Colorado 0.11987 0.7734 0.2500 0.9791 0.4503 0.9495 0.88298
## Area
## Alabama 0.08952
## Alaska 1.00000
## Arizona 0.20023
## Arkansas 0.09171
## California 0.27605
## Colorado 0.18319
再如，计算各个变量大于中位数的所有数的均值 meds = apply(state.x77, 2, median) #apply计算中位数
meanmed = function(var, med) mean(var[var > med]) #计算大于中位数的数的均值函数
meanmed(state.x77[, 1], meds[1])
## [1] 7136
meanmed(state.x77[, 1], meds[2])
## [1] 9116
上面的函数只能对每个变量单独的计算，利用sweep函数可以对所有列计算： sweep(state.x77, 2, meds, meanmed)
## [1] 15570
由于数据state.x77每个变量都不一样，想要一次性得到所有变量大于其中位数的数值均值，可以使用mapply： mapply(meanmed, as.data.frame(state.x77), meds)
## Population Income Illiteracy Life Exp Murder HS Grad
## 7136.16 4917.92 1.66 71.95 10.54 59.52
## Frost Area
## 146.84 112213.40
5.作用于组变量的函数(aggregate,tapply) 5.1 aggregate aggregate函数可以用于数据框、矩阵中一列或多列计算。该函数的参数之一是带计算的矩阵列(或数据的变量),第二个参数是用于作分组变量的list，第三个参数是用于计算的函数。以鸢尾花iris数据为例，要计算三种鸢尾花(species作为分组变量)其他变量的均值： aggregate(iris[-5], iris[5], 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
由于aggregate函数的第二个参数必须是list格式的，所以，当待处理数据是数据框格式的，可以使用“[]”下标构建分组变量，中扩号得到的是list格式。当有多个分组变量时，给定的分组变量名可以自动转为返回结果的变量名，相反如果没有给定，那么返回结果将使用形如“Group.1”的变量名。 cweights = aggregate(ChickWeight$weight, ChickWeight[c("Time", "Diet")], mean)
head(cweights)
## Time Diet x
## 1 0 1 41.40
## 2 2 1 47.25
## 3 4 1 56.47
## 4 6 1 66.79
## 5 8 1 79.68
## 6 10 1 93.05
分组变量list还可以利用list函数构建： length(list(Time = ChickWeight$Time, Diet = ChickWeight$Diet))
## [1] 2
5.2 tapply 基于一个或多个分组变量对另一个向量进行运算也可以使用tapply函数 maxweight = tapply(PlantGrowth$weight, PlantGrowth$group, max)
maxweight
## ctrl trt1 trt2
## 6.11 6.03 6.31
可以利用as.table和as.data.frame将向量转换为数据框 as.data.frame(as.table(maxweight))
## Var1 Freq
## 1 ctrl 6.11
## 2 trt1 6.03
## 3 trt2 6.31
将上例中Freq替换为其他名称，可以直接使用as.data.frame.table函数 as.data.frame.table(as.table(maxweight), responseName = "MaxWeight")
## Var1 MaxWeight
## 1 ctrl 6.11
## 2 trt1 6.03
## 3 trt2 6.31
与aggregate不同，tapply除了返数值，还可以返回数值间隔 ranges = tapply(PlantGrowth$weight, PlantGrowth$group, range)
ranges
## $ctrl
## [1] 4.17 6.11
##
## $trt1
## [1] 3.59 6.03
##
## $trt2
## [1] 4.92 6.31
将上例中list转换为数据框格式可以通过先转换为矩阵格式，进而使用data.frame函数： data.frame(group = dimnames(ranges)[[1]], matrix(unlist(ranges), ncol = 2, byrow = TRUE))
## group X1 X2
## 1 ctrl 4.17 6.11
## 2 trt1 3.59 6.03
## 3 trt2 4.92 6.31
相比cbind，data.frame函数可以避免数值型被强制转换为字符型 类似aggregate函数，tapply函数也可以接多个分组变量，以CO2数据为例： meds = tapply(CO2$uptake, CO2[c("Type", "Treatment")], median)
meds
## Treatment
## Type nonchilled chilled
## Quebec 39.2 35.0
## Mississippi 28.1 17.9
在没有median参数时，得到的结果是向量的下标 inds = tapply(CO2$uptake, CO2[c("Type", "Treatment")])
inds
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [36] 3 3 3 3 3 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4
## [71] 4 4 4 4 4 4 4 4 4 4 4 4 4 4
adj.uptake = CO2$uptake - meds[inds]
还可以使用ave函数合并上面的两步 adj.uptake = CO2$uptake - ave(CO2$uptake, CO2[c("Type", "Treatment")], FUN = median)
例，考虑找出相关系数矩阵中最大的特征值： # 首先定义一个运用于单一数据框的且满足需求的函数
maxeig = function(df) eigen(cor(df))$val[1]
# 将数据框中的数据进行分割成一系列的子数据框(转换成list)
frames = split(iris[-5], iris[5])
# 运用sapply函数，调用maxeig函数
sapply(frames, maxeig)
## setosa versicolor virginica
## 2.059 2.926 2.455
# 将上面三个命令用一个表达式表示
sapply(split(iris[-5], iris[5]), function(df) eigen(cor(df))$val[1])
## setosa versicolor virginica
## 2.059 2.926 2.455
# 还可以间接地将行向量传递给tapply函数，完成需求
tapply(1:nrow(iris), iris["Species"], function(ind, data) eigen(cor(data[ind,
-5]))$val[1], data = iris)
## Species
## setosa versicolor virginica
## 2.059 2.926 2.455
# 最后by函数可以用来对分组变量进行分解：
max.e = by(iris, iris$Species, function(df) eigen(cor(df[-5]))$val[1])
as.data.frame(as.table(max.e))
## iris.Species Freq
## 1 setosa 2.059
## 2 versicolor 2.926
## 3 virginica 2.455
为说明by函数，再举一例：在CO2数据中，根据Type和Treatmen分组的变量uptake的计数、均值、标准差 sumfun = function(x) data.frame(n = length(x$uptake), mean = mean(x$uptake),
sd = sd(x$uptake))
bb = by(CO2, CO2[c("Type", "Treatment")], sumfun)
bb
## Type: Quebec
## Treatment: nonchilled
## n mean sd
## 1 21 35.33 9.596
## --------------------------------------------------------
## Type: Mississippi
## Treatment: nonchilled
## n mean sd
## 1 21 25.95 7.402
## --------------------------------------------------------
## Type: Quebec
## Treatment: chilled
## n mean sd
## 1 21 31.75 9.645
## --------------------------------------------------------
## Type: Mississippi
## Treatment: chilled
## n mean sd
## 1 21 15.81 4.059
do.call函数可以用于调用合并函数进行合并： do.call(rbind, bb)
## n mean sd
## 1 21 35.33 9.596
## 2 21 25.95 7.402
## 3 21 31.75 9.645
## 4 21 15.81 4.059
# expand.grid对do.call进行补救，添加分组变量
cbind(expand.grid(dimnames(bb)), do.call(rbind, bb))
## Type Treatment n mean sd
## 1 Quebec nonchilled 21 35.33 9.596
## 2 Mississippi nonchilled 21 25.95 7.402
## 3 Quebec chilled 21 31.75 9.645
## 4 Mississippi chilled 21 15.81 4.059
6.reshape包 reshape包提供了一些用于数据聚合(Aggregation)的方法，该包是基于可扩展的公式方法对数据进行重构，其背后的原理是利用melt函数将数据"融化"，然后再进行重铸。对数据框、列表或数组进行融化，首先需要将变量分为id变量和度量或分析变量。默认情况下，melt会将factor型变量或整数变量视为id变量，其他的变量视为分析变量。另外，也可以通过id.var和measure.var参数对两种变量进行设置。一旦一个数据被"融化"，它就可以被重构成多种形式，以state.x77为例： # 案例数据是由state.x77，state.region合并而成
states = data.frame(state.x77, state = row.names(state.x77), region = state.region,
row.names = 1:50)
由于state和region都是factor形式的，melt函数会默认其为id变量,融化后的数据只会保留id变量，其他变量合并为variable和value变量。也可以添加variable_name参数进行融化 library(reshape)
## Loading required package: plyr
##
## Attaching package: 'reshape'
##
## 下列对象被屏蔽了from 'package:plyr':
##
## rename, round_any
mstates = melt(states)
## Using state, region as id variables
利用cast函数进行数据结构的重构，函数中，公式的左边的变量将会在结果表中呈列向，公式右边的变量将会是结果表的呈行向。“.”代表对数据进行全面的统计描述，“…”表示除公式中变量以外的全部变量。所以，我们可以使用cast(mstates,…~variable)恢复数据结构。 利用reshape包进行数据聚合(Aggeregation)，需要调用计算均值等统计量的函数。举例说明，计算依region分组的其他变量的均值： cast(mstates, region ~ variable, mean)
## region Population Income Illiteracy Life.Exp Murder HS.Grad
## 1 Northeast 5495 4570 1.000 71.26 4.722 53.97
## 2 South 4208 4012 1.738 69.71 10.581 44.34
## 3 North Central 4803 4611 0.700 71.77 5.275 54.52
## 4 West 2915 4703 1.023 71.23 7.215 62.00
## Frost Area
## 1 132.78 18141
## 2 64.62 54605
## 3 138.83 62652
## 4 102.15 134463
如果向让每个变量在单独的一行显示，可以： cast(mstates, variable ~ region, mean)
## variable Northeast South North Central West
## 1 Population 5495.111 4208.125 4803.000 2.915e+03
## 2 Income 4570.222 4011.938 4611.083 4.703e+03
## 3 Illiteracy 1.000 1.738 0.700 1.023e+00
## 4 Life.Exp 71.264 69.706 71.767 7.123e+01
## 5 Murder 4.722 10.581 5.275 7.215e+00
## 6 HS.Grad 53.967 44.344 54.517 6.200e+01
## 7 Frost 132.778 64.625 138.833 1.022e+02
## 8 Area 18141.000 54605.125 62652.000 1.345e+05
我们也可以使用subset参数对满足一定条件的变量进行计算均值： cast(mstates, region ~ variable, mean, subset = variable %in% c("Population",
"Life.Exp"))
## region Population Life.Exp
## 1 Northeast 5495 71.26
## 2 South 4208 69.71
## 3 North Central 4803 71.77
## 4 West 2915 71.23
不像aggregate函数,cast接受返回向量的函数，例如要计算states数据中Population和Life.Exp变量的均值、中位数和标准差： cast(mstates, . ~ variable, c(mean, median, sd), subset = variable %in% c("Population",
"Life.Exp"))
## value Population_mean Population_median Population_sd Life.Exp_mean
## 1 (all) 4246 2838 4464 70.88
## Life.Exp_median Life.Exp_sd
## 1 70.67 1.342
cast(mstates, variable ~ ., c(mean, median, sd), subset = variable %in% c("Population",
"Life.Exp"))
## variable mean median sd
## 1 Population 4246.42 2838.50 4464.491
## 2 Life.Exp 70.88 70.67 1.342
# 添加分组变量
cast(mstates, region ~ variable, c(mean, median, sd), subset = variable %in%
c("Population", "Life.Exp"))
## region Population_mean Population_median Population_sd
## 1 Northeast 5495 3100 6080
## 2 South 4208 3710 2780
## 3 North Central 4803 4255 3703
## 4 West 2915 1144 5579
## Life.Exp_mean Life.Exp_median Life.Exp_sd
## 1 71.26 71.23 0.7439
## 2 69.71 70.07 1.0222
## 3 71.77 72.28 1.0367
## 4 71.23 71.71 1.3520
# 如果简单的颠倒region和variable的顺序，将出现一个变量对应地区和均值等统计量的情况，这将不利于进一步的操作，添加“|”可以避免这样的情况出现：
cast(mstates, variable ~ . | region, c(mean, median, sd), subset = variable %in%
c("Population", "Life.Exp"))
## $Northeast
## variable mean median sd
## 1 Population 5495.11 3100.00 6079.5651
## 2 Life.Exp 71.26 71.23 0.7439
##
## $South
## variable mean median sd
## 1 Population 4208.12 3710.50 2779.508
## 2 Life.Exp 69.71 70.07 1.022
##
## $`North Central`
## variable mean median sd
## 1 Population 4803.00 4255.00 3702.828
## 2 Life.Exp 71.77 72.28 1.037
##
## $West
## variable mean median sd
## 1 Population 2915.31 1144.00 5578.607
## 2 Life.Exp 71.23 71.71 1.352
cast(mstates, variable ~ region, c(mean, median, sd), subset = variable %in%
c("Population", "Life.Exp"))
## variable Northeast_mean Northeast_median Northeast_sd South_mean
## 1 Population 5495.11 3100.00 6079.5651 4208.12
## 2 Life.Exp 71.26 71.23 0.7439 69.71
## South_median South_sd North Central_mean North Central_median
## 1 3710.50 2779.508 4803.00 4255.00
## 2 70.07 1.022 71.77 72.28
## North Central_sd West_mean West_median West_sd
## 1 3702.828 2915.31 1144.00 5578.607
## 2 1.037 71.23 71.71 1.352
考虑多个id变量的情况， 以ChickWeight数据为例，Time、Chick、Diet作为id变量，weight作为分析变量。由于Time是数值型的，所以需要通过measure.var参数明确指定weigth为分析变量： mChick = melt(ChickWeight, measure.var = "weight")
# 创建Diet和Time维度下Weight的中位数的数据
head(cast(mChick, Diet + Time ~ variable, median))
## Diet Time weight
## 1 1 0 41
## 2 1 2 49
## 3 1 4 56
## 4 1 6 67
## 5 1 8 79
## 6 1 10 93
# 创建每个Times时点下的Diet维度下Weight的均值
cast(mChick, Diet ~ Time + variable, mean)
## Diet 0_weight 2_weight 4_weight 6_weight 8_weight 10_weight 12_weight
## 1 1 41.4 47.25 56.47 66.79 79.68 93.05 108.5
## 2 2 40.7 49.40 59.80 75.40 91.70 108.50 131.3
## 3 3 40.8 50.40 62.20 77.90 98.40 117.10 144.4
## 4 4 41.0 51.80 64.50 83.90 105.60 126.00 151.4
## 14_weight 16_weight 18_weight 20_weight 21_weight
## 1 123.4 144.6 158.9 170.4 177.8
## 2 141.9 164.7 187.7 205.6 214.7
## 3 164.5 197.4 233.1 258.9 270.3
## 4 161.8 182.0 202.9 233.9 238.6
# 创建list格式，Diet值作为list的levels
cast(mChick, Time ~ variable | Diet, mean)
## $`1`
## Time weight
## 1 0 41.40
## 2 2 47.25
## 3 4 56.47
## 4 6 66.79
## 5 8 79.68
## 6 10 93.05
## 7 12 108.53
## 8 14 123.39
## 9 16 144.65
## 10 18 158.94
## 11 20 170.41
## 12 21 177.75
##
## $`2`
## Time weight
## 1 0 40.7
## 2 2 49.4
## 3 4 59.8
## 4 6 75.4
## 5 8 91.7
## 6 10 108.5
## 7 12 131.3
## 8 14 141.9
## 9 16 164.7
## 10 18 187.7
## 11 20 205.6
## 12 21 214.7
##
## $`3`
## Time weight
## 1 0 40.8
## 2 2 50.4
## 3 4 62.2
## 4 6 77.9
## 5 8 98.4
## 6 10 117.1
## 7 12 144.4
## 8 14 164.5
## 9 16 197.4
## 10 18 233.1
## 11 20 258.9
## 12 21 270.3
##
## $`4`
## Time weight
## 1 0 41.0
## 2 2 51.8
## 3 4 64.5
## 4 6 83.9
## 5 8 105.6
## 6 10 126.0
## 7 12 151.4
## 8 14 161.8
## 9 16 182.0
## 10 18 202.9
## 11 20 233.9
## 12 21 238.6
如果id变量之间不能随意的组合，这种情况下，可以添加add.mssing=TRUE参数以保留所有可能的id变量间的组合。例如，我们首先删除ChickWeight数据中Diet和Time 变量间的一个组合 xChickWeight = subset(ChickWeight, !(Diet == 1 & Time == 4))
mxChick = melt(xChickWeight, measure.var = "weight")
head(cast(mxChick, Diet + Time ~ variable, median))
## Diet Time weight
## 1 1 0 41
## 2 1 2 49
## 3 1 6 67
## 4 1 8 79
## 5 1 10 93
## 6 1 12 106
添加add.missing=TRUE,缺失的组合将被创建，以NA表示： head(cast(mxChick, Diet + Time ~ variable, median, add.missing = TRUE))
## Diet Time weight
## 1 1 0 41
## 2 1 2 49
## 3 1 4 NA
## 4 1 6 67
## 5 1 8 79
## 6 1 10 93
前面的例子都是先对数据仅融化(melt)再进行重铸(cast)，recast函数可以合并两步： head(recast(xChickWeight, measure.var = "weight", Diet + Time ~ variable, median,
add.missing = TRUE))
## Diet Time weight
## 1 1 0 41
## 2 1 2 49
## 3 1 4 NA
## 4 1 6 67
## 5 1 8 79
## 6 1 10 93
7.R中的循环 下面以例子的形式展现各种loops应用 7.1计算矩阵每一列的均值，system.time比较几种方法的效率 dat = matrix(rnorm(1e+06), 10000, 100)
system.time(mns <- rowMeans(dat)) #rowMeans
## user system elapsed
## 0 0 0
system.time(mns <- apply(dat, 2, mean)) #apply
## user system elapsed
## 0.09 0.00 0.10
system.time({
m <- ncol(dat)
for (i in 1:m) mns[i] <- mean(dat[, i])
})
## user system elapsed
## 0.04 0.00 0.03
上面的例子都是利用向量化的思想：每一列的均值都是调用单独的mean命令进行计算的。避免对矩阵中的元素进行循环。下面的例子是对矩阵列元素进行循环加总再除以列长度计算矩阵每列均值的方法，计算效率远低于向量化运算： slowmean = function(dat) {
n = dim(dat)[1]
m = dim(dat)[2]
mns = numeric(m)
for (i in 1:n) {
sum = 0
for (j in 1:m) sum = sum + dat[j]
mns[i] = sum/n
}
return(mns)
}
system.time(mns <- slowmean(dat))
## user system elapsed
## 1.96 0.01 1.99
这个例子也给我们一个启示：每个问题都会有多种解决方式。此例还可以直接进行矩阵运算： system.time({
m = dim(dat)[1]
mns = rep(1, m) %*% dat/m
})
## user system elapsed
## 0.01 0.00 0.02
R中盲目的使用循环只会造成低效，但是如果将循环和向量化有效的结合起来，就可以发挥循环的功效。 下面的例子将从R如何存储矩阵的角度出发，间接解释R循环 # 考虑创建一个矩阵的简单任务,矩阵的每一行代表从1到100的数字，可以这样实现这一点：
system.time(m <- matrix(1:100, 10000, 100, byrow = TRUE))
## user system elapsed
## 0.03 0.00 0.03
# 利用循环
buildrow = function() {
res = NULL
for (i in 1:10000) res = rbind(res, 1:100)
res
}
system.time(buildrow())
## user system elapsed
## 96.75 5.21 102.52
两方面原因造成循环低效，1.每一个新行添加到矩阵中，res的值就变化一次，这导致每次循环R都需要重新分配内存；2.R是以列形式存储矩阵的，每增加一行意味着矩阵的每一列都变动。根据这两点，我们可以对列循环： buildcol = function() {
res = NULL
for (i in 1:10000) res = cbind(res, 1:100)
t(res)
}
system.time(buildcol())
## user system elapsed
## 58.98 1.58 61.28
随着这种方式较行循环效率有所提升，但仍然效率低下。第一中方式的效率之高是因为使用矩阵函数matrix时，生成模拟数据之前矩阵的大小就已经确定了。也可以利用这个原理构建循环： buildrow1 = function() {
res = matrix(0, 10000, 100) #确定矩阵规模
for (i in 1:10000) res[i, ] = 1:100
res
}
system.time(buildrow1())
## user system elapsed
## 0.14 0.00 0.14
即使我们不知道一个矩阵有多少行，R也会快速分配行数，然后截断矩阵得到需要的行数。例如： # 在向矩阵添加行数，构建一个空矩阵
somerow1 = function() {
res = NULL
for (i in 1:10000) if (runif(1) < 0.5)
res = rbind(res, 1:100)
res
}
system.time(somerow1())
## user system elapsed
## 19.74 0.40 20.83
# 分配一个足够大的矩阵来包含行，然后截断矩阵
somerow2 = function() {
res = matrix(0, 10000, 100)
k = 0
for (i in 1:10000) if (runif(1) < 0.5) {
k = k + 1
res[k, ] = 1:100
}
res[1:k, ]
}
system.time(somerow2())
## user system elapsed
## 0.14 0.00 0.14
看见，首先提供一个足够大的矩阵再进行内存分配的方法较循环的运行rbind函数高效很多。 如果不能分配一个明确大小的矩阵，可以利用与矩阵存储形式不同的list，向list中添加元素需要的内存更少。这种策略是构建一个行列表,最终将转换成为矩阵的行,利用do.call函数将所有行传递给rbind函数： somerow3 = function() {
res = list()
for (i in 1:10000) if (runif(1) < 0.5)
res = c(res, list(1:100))
do.call(rbind, res)
}
system.time(somerow3())
## user system elapsed
## 8.77 0.00 8.93
|
