第7章--基本统计分析
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 组间差异的非参数检验
