对于基础安装,你可以使用summary()函数来获取描述性统计量。代码清单7-1展示了一个 示例。
## mpg hp wt
## Min. :10.40 Min. : 52.0 Min. :1.513
## 1st Qu.:15.43 1st Qu.: 96.5 1st Qu.:2.581
## Median :19.20 Median :123.0 Median :3.325
## Mean :20.09 Mean :146.7 Mean :3.217
## 3rd Qu.:22.80 3rd Qu.:180.0 3rd Qu.:3.610
## Max. :33.90 Max. :335.0 Max. :5.424
使用apply()函数或sapply()函数计算所选择的任意描述性统计量。对于sapply()函数,其使用格式为: sapply(x, FUN, options) 其中的x是你的数据框(或矩阵),FUN为一个任意的函数。如果指定了options,它们将被传递 给FUN。你可以在这里插入的典型函数有mean()、sd()、var()、min()、max()、median()、length()、range()和quantile()。函数fivenum()可返回图基五数总括(Tukey’s five-number summary,即最小值、下四分位数、中位数、上四分位数和最大值)。令人惊讶的是,基础安装并没有提供偏度和峰度的计算函数,不过你可以自行添加。代码清 单7-2中的示例计算了若干描述性统计量,其中包括偏度和峰度。常用于定制化的描述性统计分析
#通过sapply()计算描述性统计量
mystats <- function(x, na.omit=FALSE){
if (na.omit)
x <- x[!is.na(x)]
m <- mean(x)
n <- length(x)
s <- sd(x)
skew <- sum((x-m)^3/s^3)/n
kurt <- sum((x-m)^4/s^4)/n - 3
return(c(n=n, mean=m, stdev=s, skew=skew, kurtosis=kurt))
}
myvars <- c("mpg", "hp", "wt")
sapply(mtcars[myvars], mystats)
## mpg hp wt
## n 32.000000 32.0000000 32.00000000
## mean 20.090625 146.6875000 3.21725000
## stdev 6.026948 68.5628685 0.97845744
## skew 0.610655 0.7260237 0.42314646
## kurtosis -0.372766 -0.1355511 -0.02271075
Hmisc包中的describe()函数可返回变量和观测的数量、缺失值和唯一值的数目、平均值、 分位数,以及五个最大的值和五个最小的值。
## Warning: package 'Hmisc' was built under R version 3.5.3
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.5.3
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
## mtcars[myvars]
##
## 3 Variables 32 Observations
## ---------------------------------------------------------------------------
## mpg
## n missing distinct Info Mean Gmd .05 .10
## 32 0 25 0.999 20.09 6.796 12.00 14.34
## .25 .50 .75 .90 .95
## 15.43 19.20 22.80 30.09 31.30
##
## lowest : 10.4 13.3 14.3 14.7 15.0, highest: 26.0 27.3 30.4 32.4 33.9
## ---------------------------------------------------------------------------
## hp
## n missing distinct Info Mean Gmd .05 .10
## 32 0 22 0.997 146.7 77.04 63.65 66.00
## .25 .50 .75 .90 .95
## 96.50 123.00 180.00 243.50 253.55
##
## lowest : 52 62 65 66 91, highest: 215 230 245 264 335
## ---------------------------------------------------------------------------
## wt
## n missing distinct Info Mean Gmd .05 .10
## 32 0 29 0.999 3.217 1.089 1.736 1.956
## .25 .50 .75 .90 .95
## 2.581 3.325 3.610 4.048 5.293
##
## lowest : 1.513 1.615 1.835 1.935 2.140, highest: 3.845 4.070 5.250 5.345 5.424
## ---------------------------------------------------------------------------
pastecs包中有一个名为stat.desc()的函数,它可以计算种类繁多的描述性统计量。使用 格式为: stat.desc(x, basic=TRUE, desc=TRUE, norm=FALSE, p=0.95) 其中的x是一个数据框或时间序列。若basic=TRUE(默认值),则计算其中所有值、空值、缺失 值的数量,以及最小值、最大值、值域,还有总和。若desc=TRUE(同样也是默认值),则计算 中位数、平均数、平均数的标准误、平均数置信度为95%的置信区间、方差、标准差以及变异系 数。最后,若norm=TRUE(不是默认的),则返回正态分布统计量,包括偏度和峰度(以及它们 的统计显著程度)和Shapiro-Wilk正态检验结果。这里使用了p值来计算平均数的置信区间(默认 置信度为0.95)。
## Warning: package 'pastecs' was built under R version 3.5.3
## mpg hp wt
## nbr.val 32.0000000 32.0000000 32.0000000
## nbr.null 0.0000000 0.0000000 0.0000000
## nbr.na 0.0000000 0.0000000 0.0000000
## min 10.4000000 52.0000000 1.5130000
## max 33.9000000 335.0000000 5.4240000
## range 23.5000000 283.0000000 3.9110000
## sum 642.9000000 4694.0000000 102.9520000
## median 19.2000000 123.0000000 3.3250000
## mean 20.0906250 146.6875000 3.2172500
## SE.mean 1.0654240 12.1203173 0.1729685
## CI.mean.0.95 2.1729465 24.7195501 0.3527715
## var 36.3241028 4700.8669355 0.9573790
## std.dev 6.0269481 68.5628685 0.9784574
## coef.var 0.2999881 0.4674077 0.3041285
psych包也拥有一个名为describe()的函数,它可以计算非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均值的标准误。代码清单7-5中有一个示例。
## Warning: package 'psych' was built under R version 3.5.3
##
## Attaching package: 'psych'
## The following object is masked from 'package:Hmisc':
##
## describe
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
## vars n mean sd median trimmed mad min max range skew
## mpg 1 32 20.09 6.03 19.20 19.70 5.41 10.40 33.90 23.50 0.61
## hp 2 32 146.69 68.56 123.00 141.19 77.10 52.00 335.00 283.00 0.73
## wt 3 32 3.22 0.98 3.33 3.15 0.77 1.51 5.42 3.91 0.42
## kurtosis se
## mpg -0.37 1.07
## hp -0.14 12.12
## wt -0.02 0.17
注意
在前面的示例中,psych包和Hmisc包均提供了名为describe()的函数。R如何知道该 使用哪个呢?简言之,如代码清单7-5所示,最后载入的程序包优先。在这里,psych在 Hmisc之后被载入,然后显示了一条信息,提示Hmisc包中的describe()函数被psych 包中的同名函数所屏蔽(masked)。键入describe()后,R在搜索这个函数时将首先找 到psych包中的函数并执行它。如果你想改而使用Hmisc包中的版本,可以键入 Hmisc::describe(mt)。这个函数仍然在那里。你只是需要给予R更多信息以找到它。
##使用aggregate()分组获取描述性统计量
在比较多组个体或观测时,关注的焦点经常是各组的描述性统计信息,而不是样本整体的描 述性统计信息。同样地,在R中完成这个任务有若干种方法。我们将以获取变速箱类型各水平的 描述性统计量开始。 在第5章中,我们讨论了整合数据的方法。你可以使用aggregate()函数(5.6.2节)来分组 获取描述性统计量,如代码清单7-6所示。
#代码清单7-6 使用aggregate()分组获取描述性统计量
myvars <- c("mpg", "hp", "wt")
aggregate(mtcars[myvars], by=list(am=mtcars$am), mean)
## am mpg hp wt
## 1 0 17.14737 160.2632 3.768895
## 2 1 24.39231 126.8462 2.411000
## am mpg hp wt
## 1 0 3.833966 53.90820 0.7774001
## 2 1 6.166504 84.06232 0.6169816
注意
list(am=mtcars(am)的使用。如果使用的是list(mtcars)am),则am列将被标注为 Group.1而不是am。你使用这个赋值指定了一个更有帮助的列标签。如果有多个分组变量,可以 使用by=list(name1=groupvar1, name2=groupvar2, … , nameN=groupvarN)这样的 语句。
遗憾的是,aggregate()仅允许在每次调用中使用平均数、标准差这样的单返回值函数。 它无法一次返回若干个统计量。要完成这项任务,可以使用by()函数。格式为: by(data, INDICES, FUN) 其中data是一个数据框或矩阵,INDICES是一个因子或因子组成的列表,定义了分组,FUN是任 意函数。
#使用by()分组计算描述性统计量
dstats <- function(x)sapply(x, mystats)
myvars <- c("mpg", "hp", "wt")
by(mtcars[myvars], mtcars$am, dstats)
## mtcars$am: 0
## mpg hp wt
## n 19.00000000 19.00000000 19.0000000
## mean 17.14736842 160.26315789 3.7688947
## stdev 3.83396639 53.90819573 0.7774001
## skew 0.01395038 -0.01422519 0.9759294
## kurtosis -0.80317826 -1.20969733 0.1415676
## --------------------------------------------------------
## mtcars$am: 1
## mpg hp wt
## n 13.00000000 13.0000000 13.0000000
## mean 24.39230769 126.8461538 2.4110000
## stdev 6.16650381 84.0623243 0.6169816
## skew 0.05256118 1.3598859 0.2103128
## kurtosis -1.45535200 0.5634635 -1.1737358
这里,dstats()调用了代码清单7-2中的mystats()函数,将其应用于数据框的每一栏中。 再通过by()函数则可得到am中每一水平的概括统计量。
doBy包中summaryBy()函数的使用格式为: summaryBy(formula, data=dataframe, FUN=function)
, 其中的formula接受以下的格式:
var1 + var2 + var3 + ... + varN ~ groupvar1 + groupvar2 + ... + groupvarN
, 在~左侧的变量是需要分析的数值型变量,而右侧的变量是类别型的分组变量。function可为 任何内建或用户自编的R函数。使用7.2.1节中创建的mystats()函数的一个示例如代码清单7-8 所示。
## Warning: package 'doBy' was built under R version 3.5.3
## am mpg.n mpg.mean mpg.stdev mpg.skew mpg.kurtosis hp.n hp.mean
## 1 0 19 17.14737 3.833966 0.01395038 -0.8031783 19 160.2632
## 2 1 13 24.39231 6.166504 0.05256118 -1.4553520 13 126.8462
## hp.stdev hp.skew hp.kurtosis wt.n wt.mean wt.stdev wt.skew
## 1 53.90820 -0.01422519 -1.2096973 19 3.768895 0.7774001 0.9759294
## 2 84.06232 1.35988586 0.5634635 13 2.411000 0.6169816 0.2103128
## wt.kurtosis
## 1 0.1415676
## 2 -1.1737358
psych包中的describeBy()函数可计算和describe()相同的描述性统计量,只是按照一 个或多个分组变量分层。
##
## Descriptive statistics by group
## am: 0
## vars n mean sd median trimmed mad min max range skew
## mpg 1 19 17.15 3.83 17.30 17.12 3.11 10.40 24.40 14.00 0.01
## hp 2 19 160.26 53.91 175.00 161.06 77.10 62.00 245.00 183.00 -0.01
## wt 3 19 3.77 0.78 3.52 3.75 0.45 2.46 5.42 2.96 0.98
## kurtosis se
## mpg -0.80 0.88
## hp -1.21 12.37
## wt 0.14 0.18
## --------------------------------------------------------
## am: 1
## vars n mean sd median trimmed mad min max range skew
## mpg 1 13 24.39 6.17 22.80 24.38 6.67 15.00 33.90 18.90 0.05
## hp 2 13 126.85 84.06 109.00 114.73 63.75 52.00 335.00 283.00 1.36
## wt 3 13 2.41 0.62 2.32 2.39 0.68 1.51 3.57 2.06 0.21
## kurtosis se
## mpg -1.46 1.71
## hp 0.56 23.31
## wt -1.17 0.17