7.1 描述性统计分析
> vars <- c("mpg", "hp", "wt")
> head(mtcars[vars])
mpg hp wt
Mazda RX4 21.0 110 2.620
Mazda RX4 Wag 21.0 110 2.875
Datsun 710 22.8 93 2.320
Hornet 4 Drive 21.4 110 3.215
Hornet Sportabout 18.7 175 3.440
Valiant 18.1 105 3.460
以上述数据集为例，对于基础安装R的用户，可以使用summary()函数来获取描述性统计量。这个函数提供了最小值、最大值、四分位数和数值型变量的均值，以及因子向量和逻辑型向量的频数统计。
> summary(mtcars[vars])
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
分组计算描述性统计量
你可以使用aggregate()函数来分组获取描述性统计量。
> aggregate(mtcars[vars],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
> aggregate(mtcars[vars],by=list(am=mtcars$am),sd)
am mpg hp wt
1 0 3.833966 53.90820 0.7774001
2 1 6.166504 84.06232 0.6169816
注意list的使用。如果使用的是list(mtcars$am)，则am列将被标注为Group.1而不是am。你使用这个赋值指定了一个更有帮助的列标签。如果有多个分组变量，你可以使用一下语句：
by=list(name1=groupvar1, name2=groupvar2,...)
而aggregate函数仅允许每次调用中使用平均数、标准差这样的单返回值函数。如果想要一次返回若干个统计量，可以使用by()函数。
7.2 频数表和列联表
下面的示例中，假设A、B和C代表类别型变量。
> library(vcd)
> library(grid)
> head(Arthritis)
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
1. 生成频数表
1.1 一维列联表
可以使用table()函数生成简单的频数统计表。
> mytable <- with(Arthritis, table(Improved))
> mytable
Improved
None Some Marked
42 14 28
可以使用prop.table()将这些频数转换为比例值，或使用prop.table()*100来转换为百分比：
> prop.table(mytable)
Improved
None Some Marked
0.5000000 0.1666667 0.3333333
1.2 二维列联表
对于二维表，table的使用格式为：
mytable <- table(A, B)
其中，A是行变量，B是列变量。
除此之外，xtabs()函数还可使用公式风格的输入创建列联表。
mytable <- xtabs(~ A + B, data=mydata)
其中，mydata是一个矩阵或者数据框，要进行交叉分类的变量应出现在公式的右侧：
> mytable <- xtabs(~ Treatment+Improved, data=Arthritis)
> mytable
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
可以使用margin.table()和prop.table()函数分别生成边际频数和比例。
> margin.table(mytable,1)
Treatment
Placebo Treated
43 41
> prop.table(mytable,1)
Improved
Treatment None Some Marked
Placebo 0.6744186 0.1627907 0.1627907
Treated 0.3170732 0.1707317 0.5121951
下标1指代table语句的第一个变量，列和列比例可以这样计算：
> margin.table(mytable,2)
Improved
None Some Marked
42 14 28
> prop.table(mytable,2)
Improved
Treatment None Some Marked
Placebo 0.6904762 0.5000000 0.2500000
Treated 0.3095238 0.5000000 0.7500000
各单元格所占比例可用如下语句获取：
> prop.table(mytable)
Improved
Treatment None Some Marked
Placebo 0.34523810 0.08333333 0.08333333
Treated 0.15476190 0.08333333 0.25000000
可以使用addmargins()函数为这些表格添加边际和。
3. 多维列联表
除了上面二维列联表里介绍的方法之外，还有ftable()函数可以以一种紧凑而吸引人的方式输出多维列联表。
> mytable <- xtabs(~Treatment+Sex+Improved, data=Arthritis)
> mytable
, , 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
> ftable(mytable)
Improved None Some Marked
Treatment Sex
Placebo Female 19 7 6
Male 10 0 1
Treated Female 6 5 16
Male 7 2 5
> margin.table(mytable,1)
Treatment
Placebo Treated
43 41
> margin.table(mytable,2)
Sex
Female Male
59 25
> margin.table(mytable,3)
Improved
None Some Marked
42 14 28
> margin.table(mytable,c(1,3))
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
> ftable(prop.table(mytable,c(1,2)))
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
> ftable(addmargins(prop.table(mytable,c(1,2)),3))
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
>
2. 独立性检验
1. 卡方独立性检验
可以使用chisq.test()函数对二维列表的行变量和列变量进行卡方独立性检验。
> mytable <- xtabs(~Treatment+Improved, data=Arthritis)
> chisq.test(mytable)
Pearson's Chi-squared test
data: mytable
X-squared = 13.055, df = 2, p-value = 0.001463
> mytable <- xtabs(~Treatment+Sex, data=Arthritis)
> chisq.test(mytable)
Pearson's Chi-squared test with Yates' continuity correction
data: mytable
X-squared = 0.38378, df = 1, p-value = 0.5356
Treatment和Improved有相关性，和Sex没有相关性。
2. Fisher精确检验
使用fisher.test()函数进行Fisher精确检验，Fisher精确检验的原假设是：边界固定的列联表中行和列是相互独立的。
> mytable <- xtabs(~Treatment+Improved, data=Arthritis)
> fisher.test(mytable)
Fisher's Exact Test for Count Data
data: mytable
p-value = 0.001393
alternative hypothesis: two.sided
3. Cochran-Mantel-Haenszel检验
mantelhaen.test()函数可用来进行该检验，其原假设是：两个名义变量在第三个变量的每一层中都是条件独立的。
> mytable <- xtabs(~Treatment+Improved+Sex, data=Arthritis)
> mantelhaen.test(mytable)
Cochran-Mantel-Haenszel test
data: mytable
Cochran-Mantel-Haenszel M^2 = 14.632, df = 2, p-value = 0.0006647
7.3 相关
1. 相关的类型
Pearson，Spearman和Kendall相关
Pearson积差关系数衡量了两个定量变量之间的线性相关程度。
Spearman等级相关系数则衡量分级定序变量之间的相关程度。
Kendall's Tau相关系数也是一种非参数的等级相关度量。
cor()函数可以计算这三种相关系数，而cov()函数可用来计算协方差。
cor(x, use= ,method= )
x指一个矩阵或者数据框；
use是指定数据缺失的处理方式；
method是指定相关系数的类型，默认为Pearson类型。
> states <- state.x77[,1:6]
> cov(states)
Population Income Illiteracy Life Exp Murder HS Grad
Population 19931683.7588 571229.7796 292.8679592 -407.8424612 5663.523714 -3551.509551
Income 571229.7796 377573.3061 -163.7020408 280.6631837 -521.894286 3076.768980
Illiteracy 292.8680 -163.7020 0.3715306 -0.4815122 1.581776 -3.235469
Life Exp -407.8425 280.6632 -0.4815122 1.8020204 -3.869480 6.312685
Murder 5663.5237 -521.8943 1.5817755 -3.8694804 13.627465 -14.549616
HS Grad -3551.5096 3076.7690 -3.2354694 6.3126849 -14.549616 65.237894
> cor(states)
Population Income Illiteracy Life Exp Murder HS Grad
Population 1.00000000 0.2082276 0.1076224 -0.06805195 0.3436428 -0.09848975
Income 0.20822756 1.0000000 -0.4370752 0.34025534 -0.2300776 0.61993232
Illiteracy 0.10762237 -0.4370752 1.0000000 -0.58847793 0.7029752 -0.65718861
Life Exp -0.06805195 0.3402553 -0.5884779 1.00000000 -0.7808458 0.58221620
Murder 0.34364275 -0.2300776 0.7029752 -0.78084575 1.0000000 -0.48797102
HS Grad -0.09848975 0.6199323 -0.6571886 0.58221620 -0.4879710 1.00000000
> cor(states,method="spearman")
Population Income Illiteracy Life Exp Murder HS Grad
Population 1.0000000 0.1246098 0.3130496 -0.1040171 0.3457401 -0.3833649
Income 0.1246098 1.0000000 -0.3145948 0.3241050 -0.2174623 0.5104809
Illiteracy 0.3130496 -0.3145948 1.0000000 -0.5553735 0.6723592 -0.6545396
Life Exp -0.1040171 0.3241050 -0.5553735 1.0000000 -0.7802406 0.5239410
Murder 0.3457401 -0.2174623 0.6723592 -0.7802406 1.0000000 -0.4367330
HS Grad -0.3833649 0.5104809 -0.6545396 0.5239410 -0.4367330 1.0000000
在默认情况下，得到的结果是一个方阵（所有变量之间两两计算相关）。你同样可以计算非方形的相关矩阵。
> x <- states[,c("Population", "Income", "Illiteracy", "HS Grad")]
> y <- states[,c("Life Exp", "Murder")]
> cor(x, y)
Life Exp Murder
Population -0.06805195 0.3436428
Income 0.34025534 -0.2300776
Illiteracy -0.58847793 0.7029752
HS Grad 0.58221620 -0.4879710
偏相关
偏相关是控制一个或多个定量变量时，另外两个定量变量之间的相关系数。可以使用ggm包中的pcor()函数计算偏相关系数。
> pcor(c(1,5,2,3,6), cov(states))
[1] 0.3462724
7.3 相关性的显著性检验
常用的原假设为变量间不相关（即总体的相关系数为0）。你可以使用cor.test()函数对单个的Pearson、Spearman和Kendall相关系数进行检验。简化后的格式为：
cor.test(x,y, alternative= , method= ,)
x和y为要检验相关性的变量；
alternative用来指定进行双侧检验或单侧检验；
method用来指定要计算的相关类型。
> cor.test(states[,3],states[,5])
Pearson's product-moment correlation
data: states[, 3] and states[, 5]
t = 6.8479, df = 48, p-value = 1.258e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.5279280 0.8207295
sample estimates:
cor
0.7029752
但是，cor.test()一次只能检验一种相关关系。psych包中提供的corr.test()函数可以一次做更多事情。
> library(psych)
> corr.test(states, use="complete")
Call:corr.test(x = states, use = "complete")
Correlation matrix
Population Income Illiteracy Life Exp Murder HS Grad
Population 1.00 0.21 0.11 -0.07 0.34 -0.10
Income 0.21 1.00 -0.44 0.34 -0.23 0.62
Illiteracy 0.11 -0.44 1.00 -0.59 0.70 -0.66
Life Exp -0.07 0.34 -0.59 1.00 -0.78 0.58
Murder 0.34 -0.23 0.70 -0.78 1.00 -0.49
HS Grad -0.10 0.62 -0.66 0.58 -0.49 1.00
Sample Size
[1] 50
Probability values (Entries above the diagonal are adjusted for multiple tests.)
Population Income Illiteracy Life Exp Murder HS Grad
Population 0.00 0.59 1.00 1.0 0.10 1
Income 0.15 0.00 0.01 0.1 0.54 0
Illiteracy 0.46 0.00 0.00 0.0 0.00 0
Life Exp 0.64 0.02 0.00 0.0 0.00 0
Murder 0.01 0.11 0.00 0.0 0.00 0
HS Grad 0.50 0.00 0.00 0.0 0.00 0
To see confidence intervals of the correlations, print with the short=FALSE option
参数use=的取值可为"pairwise"或"complete"（分别表示对缺失值执行成对删除或行删除）。
参数method=的取值是三种方法，默认为pearson。
7.4 t检验
1. 独立样本的t检验
一个针对两组的独立样本t检验可以用于检验两个总体的均值相等的假设。这里假设两组数据是独立的，并且是从正态总体中抽的的。
t.test(y ~ x, data)
> library(MASS)
> t.test(Prob ~ So, data=UScrime)
Welch Two Sample t-test
data: Prob by So
t = -3.8954, df = 24.925, p-value = 0.0006506
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.03852569 -0.01187439
sample estimates:
mean in group 0 mean in group 1
0.03851265 0.06371269
2. 非独立样本的t检验
在两组观测之间相关时，你获得的是一个非独立组设计。非独立样本的t检验假定组间的差异呈正态分布。
7.5 组间差异的非参数检验