对于基础安装，你可以使用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更多信息以找到它。
##2.1使用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
与前面的示例不同，describeBy()函数不允许指定任意函数，所以它的普适性较低。若存在一个以上的分组变量，你可以使用list(name1=groupvar1, name2=groupvar2, … ,nameN=groupvarN)来表示它们。但这仅在分组变量交叉后不出现空白单元时有效。
我们将着眼于类别型变量的频数表和列联表，以及相应的独立性检验、相关性的 度量、图形化展示结果的方法。我们除了使用基础安装中的函数，还将连带使用vcd包和gmodels 包中的函数。下面的示例中，假设A、B和C代表类别型变量。 数据来自vcd包中的Arthritis数据集。这份数据来自Kock & Edward （1988），表 示了一项风湿性关节炎新疗法的双盲临床实验的结果。前几个观测是这样的：
## Warning: package 'vcd' was built under R version 3.5.2
## Loading required package: grid
## ID Treatment Sex Age Improved
## 1 57 Treated Male 27 Some
## 2 46 Treated Male 29 None
## 3 77 Treated Male 30 None
## 4 17 Treated Male 32 Marked
## 5 36 Treated Male 46 Marked
## 6 23 Treated Male 58 Marked
治疗情况（安慰剂治疗、用药治疗）、性别（男性、女性）和改善情况（无改善、一定程度 的改善、显著改善）均为类别型因子①。我们将使用此数据创建频数表和列联表（交 叉的分类）。 ## 3.1生成频数表 R中提供了用于创建频数表和列联表的若干种方法。其中最重要的函数已列于表7-1中。
## Improved
## None Some Marked
## 42 14 28
## Improved
## None Some Marked
## 0.5000000 0.1666667 0.3333333
## Improved
## None Some Marked
## 50.00000 16.66667 33.33333
对于二维列联表，table()函数的使用格式为：
mytable <- table(A, B)
其中的A是行变量，B是列变量。
除此之外，xtabs()函数还可使用公式风格的输入创建列联表， 格式为：
mytable <- xtabs(~ A + B, data=mydata)
其中的mydata是一个矩阵或数据框。总的来说，要进行交叉分类的变量应出现在公式的右侧（即 ~符号的右方），以+作为分隔符。若某个变量写在公式的左侧，则其为一个频数向量（在数据已经被表格化时很有用）。
## Improved
## Treatment None Some Marked
## Placebo 29 7 7
## Treated 13 7 21
可以使用margin.table()和prop.table()函数分别生成边际频数和比例。行和与行比 例可以这样计算：
## Treatment
## Placebo Treated
## 43 41
## Improved
## Treatment None Some Marked
## Placebo 0.6744186 0.1627907 0.1627907
## Treated 0.3170732 0.1707317 0.5121951
下标1指代table()语句中的第一个变量。观察表格可以发现，与接受安慰剂的个体中有显 著改善的16%相比，接受治疗的个体中的51%的个体病情有了显著的改善。 列和与列比例可以这样计算：
## Improved
## None Some Marked
## 42 14 28
## Improved
## Treatment None Some Marked
## Placebo 0.6904762 0.5000000 0.2500000
## Treated 0.3095238 0.5000000 0.7500000
各单元格所占比例可用如下语句获取：
## Improved
## Treatment None Some Marked
## Placebo 0.34523810 0.08333333 0.08333333
## Treated 0.15476190 0.08333333 0.25000000
你可以使用addmargins()函数为这些表格添加边际和。例如，以下代码添加了各行的和与 各列的和：
## Improved
## Treatment None Some Marked Sum
## Placebo 29 7 7 43
## Treated 13 7 21 41
## Sum 42 14 28 84
## Improved
## Treatment None Some Marked Sum
## Placebo 0.34523810 0.08333333 0.08333333 0.51190476
## Treated 0.15476190 0.08333333 0.25000000 0.48809524
## Sum 0.50000000 0.16666667 0.33333333 1.00000000
在使用addmargins()时，默认行为是为表中所有的变量创建边际和。作为对照：
## Improved
## Treatment None Some Marked Sum
## Placebo 0.6744186 0.1627907 0.1627907 1.0000000
## Treated 0.3170732 0.1707317 0.5121951 1.0000000
仅添加了各行的和。类似地，
## Improved
## Treatment None Some Marked
## Placebo 0.6904762 0.5000000 0.2500000
## Treated 0.3095238 0.5000000 0.7500000
## Sum 1.0000000 1.0000000 1.0000000
注意 table()函数默认忽略缺失值（NA）。要在频数统计中将NA视为一个有效的类别，请设定 参数useNA=“ifany”。
使用CrossTable生成二维列联表 使用gmodels包中的CrossTable()函数是创建二维列联表的第三种方法。CrossTable() 函数仿照SAS中PROC FREQ或SPSS中CROSSTABS的形式生成二维列联表。
## Warning: package 'gmodels' was built under R version 3.5.3
##
##
## Cell Contents
## |-------------------------|
## | N |
## | Chi-square contribution |
## | N / Row Total |
## | N / Col Total |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 84
##
##
## | Arthritis$Improved
## Arthritis$Treatment | None | Some | Marked | Row Total |
## --------------------|-----------|-----------|-----------|-----------|
## Placebo | 29 | 7 | 7 | 43 |
## | 2.616 | 0.004 | 3.752 | |
## | 0.674 | 0.163 | 0.163 | 0.512 |
## | 0.690 | 0.500 | 0.250 | |
## | 0.345 | 0.083 | 0.083 | |
## --------------------|-----------|-----------|-----------|-----------|
## Treated | 13 | 7 | 21 | 41 |
## | 2.744 | 0.004 | 3.935 | |
## | 0.317 | 0.171 | 0.512 | 0.488 |
## | 0.310 | 0.500 | 0.750 | |
## | 0.155 | 0.083 | 0.250 | |
## --------------------|-----------|-----------|-----------|-----------|
## Column Total | 42 | 14 | 28 | 84 |
## | 0.500 | 0.167 | 0.333 | |
## --------------------|-----------|-----------|-----------|-----------|
##
##
CrossTable()函数有很多选项，可以做许多事情：计算（行、列、单元格）的百分比；指 定小数位数；进行卡方、Fisher和McNemar独立性检验；计算期望和（皮尔逊、标准化、调整的 标准化）残差；将缺失值作为一种有效值；进行行和列标题的标注；生成SAS或SPSS风格的输出。 参阅help(CrossTable)以了解详情。
如果有两个以上的类别型变量，那么你就是在处理多维列联表。table() 和 xtabs() 都可以基于三个或更多的类别型变量生成多维列联表。 margin.table()、prop.table()和addmargins()函数可以自然地推广到高于二维的情况。 另外，ftable()函数可以以一种紧凑而吸引人的方式输出多维列联表。
## , , Improved = None
##
## Sex
## Treatment Female Male
## Placebo 19 10
## Treated 6 7
##
## , , Improved = Some
##
## Sex
## Treatment Female Male
## Placebo 7 0
## Treated 5 2
##
## , , Improved = Marked
##
## Sex
## Treatment Female Male
## Placebo 6 1
## Treated 16 5
## Improved None Some Marked
## Treatment Sex
## Placebo Female 19 7 6
## Male 10 0 1
## Treated Female 6 5 16
## Male 7 2 5
## Treatment
## Placebo Treated
## 43 41
## Sex
## Female Male
## 59 25
## Improved
## None Some Marked
## 42 14 28
## Improved
## Treatment None Some Marked
## Placebo 29 7 7
## Treated 13 7 21
## Improved None Some Marked
## Treatment Sex
## Placebo Female 0.59375000 0.21875000 0.18750000
## Male 0.90909091 0.00000000 0.09090909
## Treated Female 0.22222222 0.18518519 0.59259259
## Male 0.50000000 0.14285714 0.35714286
## Improved None Some Marked Sum
## Treatment Sex
## Placebo Female 0.59375000 0.21875000 0.18750000 1.00000000
## Male 0.90909091 0.00000000 0.09090909 1.00000000
## Treated Female 0.22222222 0.18518519 0.59259259 1.00000000
## Male 0.50000000 0.14285714 0.35714286 1.00000000
第1部分代码生成了三维分组各单元格的频数。这段代码同时演示了如何使用ftable()函 数输出更为紧凑和吸引人的表格。 第2部分代码为治疗情况（Treatment）、性别（Sex）和改善情况（Improved）生成了边 际频数。由于使用公式~Treatement+Sex+Improve创建了这个表，所以Treatment需要通过下标1来引用，Sex通过下标2来引用，Improve通过下标3来引用。 第3部分代码为治疗情况（Treatment）× 改善情况（Improved）分组的边际频数，由不 同性别（Sex）的单元加和而成。每个Treatment×Sex组合中改善情况为None、Some和Marked 患者的比例由4部分给出。在这里可以看到治疗组的男性中有36%有了显著改善，女性为59%。总而 言之，比例将被添加到不在prop.table()调用中的下标上（本例中是第三个下标，或称 Improve）。在最后一个例子中可以看到这一点，你在那里为第三个下标添加了边际和。 如果想得到百分比而不是比例，可以将结果表格乘以100。例如：
## Improved None Some Marked Sum
## Treatment Sex
## Placebo Female 59.375000 21.875000 18.750000 100.000000
## Male 90.909091 0.000000 9.090909 100.000000
## Treated Female 22.222222 18.518519 59.259259 100.000000
## Male 50.000000 14.285714 35.714286 100.000000
列联表可以告诉你组成表格的各种变量组合的频数或比例，不过你可能还会对列联表中的变 量是否相关或独立感兴趣。