第7章 基本统计分析
基本内容
- 描述性统计分析
- 频数表和列联表
- 相关系数和协方差
- t检验
- 非参数统计
7。1 描述性统计分析
我们将关注 连续变量的中心趋势、变化性、和分布形状的方法。
我们关注的重点是每加仑油的行驶英里数(mpg)，马力(hp)，车重(wt)
> myvars<-c("mpg","hp","wt")
> head(mtcars[myvars])
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
7。1。1 方法云集
使用summary()函数来获取描述性统计量。
代码7-1 通过summary()计算描述性统计量
> myvars<-c("mpg","hp","wt")
> summary(mtcars[myvars])
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
summary()函数提供了最小值和最大值、四分位数和数值型变量的均值，以及因子向量和逻辑向量的频数统计。
可以使用apply() 或 sapply()函数计算所选择的任意描述性统计变量。
对sapply()函数，使用格式为：
sapply(x,FUN,options)
其中，x是数据框（或矩阵），FUN为一个任意函数。若指定了options，将被传递给FUN，可以插入的典型函数有，mean()、sd()、var()、min()、max()、median()、length()、range()、quantile()
函数fivenum()可以返回图基五数总括，最小值，下四分位数，中位数，上四分位数，最大值。
代码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
如果不只希望忽略缺失值，那么就应当使用
sapply(mtcars[myvars],mystats,na.omit=TRUE)
7.1.2更多方法
关于计算描述性统计方法，其中还包括Hmisc、pastecs、psych
Hmisc包中的describe()函数可以返回变量和观测值的数量、缺失值和唯一值的数目、平均值、分位数，以及五个最大的值和五个最小的值。
代码7-3 通过Hmisc包中的describe()函数计算描述性统计变量
> library(Hmisc)
> myvars<-c("mpg","hp","wt")
> describe(mtcars[myvars])
mtcars[myvars]
3 Variables 32 Observations
-----------------------------------------------------------------------------mpg
n missing distinct Info Mean Gmd .05 .10 .25 .50
32 0 25 0.999 20.09 6.796 12.00 14.34 15.43 19.20
.75 .90 .95
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 .25 .50
32 0 22 0.997 146.7 77.04 63.65 66.00 96.50 123.00
.75 .90 .95
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 .25 .50
32 0 29 0.999 3.217 1.089 1.736 1.956 2.581 3.325
.75 .90 .95
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=FALSE(不是默认的)，则返回正态分布统计量，包括偏度和峰度（以及他们的统计显著程度）和Shapiro-Wilk正态检验结果。这里使用p值来计算平均数的置信区间（默认置信度为0.95）
代码7-4 通过pastecs包中的stat.desc()函数计算描述性统计量
> library(pastecs)
> options(digits = 3)
> myvars<-c("mpg","hp","wt")
> stat.desc(mtcars[myvars])
mpg hp wt
nbr.val 32.00 32.000 32.000
nbr.null 0.00 0.000 0.000
nbr.na 0.00 0.000 0.000
min 10.40 52.000 1.513
max 33.90 335.000 5.424
range 23.50 283.000 3.911
sum 642.90 4694.000 102.952
median 19.20 123.000 3.325
mean 20.09 146.688 3.217
SE.mean 1.07 12.120 0.173
CI.mean.0.95 2.17 24.720 0.353
var 36.32 4700.867 0.957
std.dev 6.03 68.563 0.978
coef.var 0.30 0.467 0.304
我们还可以通过psych包中的describe()函数，来计算非缺失值的数量、平均数、标准值、中位数、截尾数、绝对中位差、最小值、最大值、值域、偏度和平均值的标准误。
代码7-5
> library(psych)
> library(psych)
> myvars<-c("mpg","hp","wt")
> describe(mtcars[myvars])
vars n mean sd median trimmed mad min max range skew kurtosis se
mpg 1 32 20.09 6.03 19.20 19.70 5.41 10.40 33.90 23.50 0.61 -0.37 1.07
hp 2 32 146.69 68.56 123.00 141.19 77.10 52.00 335.00 283.00 0.73 -0.14 12.12
psych包和Himsc包均有describe()函数，R如何知道使用哪个？那就是最后载入的程序优先。
7.1.3 分组计算描述性统计量
在比较多组的个体或观测时，关注焦点在各组的描述性统计信息。我们使用aggregate()函数来分组描述获取描述性统计量。
代码7-6 使用aggregate()分组获得描述性统计量
> myvars<-c("mpg","hp","wt")
> aggregate(mtcars[myvars],by=list(am=mtcars$am),mean)
am mpg hp wt
1 0 17.1 160 3.77
2 1 24.4 127 2.41
> aggregate(mtcars[myvars],by=list(am=mtcars$am),sd)
am mpg hp wt
1 0 3.83 53.9 0.777
2 1 6.17 84.1 0.617
wt 3 32 3.22 0.98 3.33 3.15 0.77 1.51 5.42 3.91 0.42 -0.02 0.17
aggregate()允许在每次调用中使用平均数、标准差这样的单返回值函数。
它无法一次返回若干个统计量，要完成这项任务，可以使用by()函数，格式：
by(data,INDICES,FUN)
其中，data是一个数据框或矩阵，NDICES是一个因子或因子组成的列表，定义了分组，FU那是任意函数。
代码7-7 使用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.000 19.0000 19.000
mean 17.147 160.2632 3.769
stdev 3.834 53.9082 0.777
skew 0.014 -0.0142 0.976
kurtosis -0.803 -1.2097 0.142
-----------------------------------------------------------------------------
mtcars$am: 1
mpg hp wt
n 13.0000 13.000 13.000
mean 24.3923 126.846 2.411
stdev 6.1665 84.062 0.617
skew 0.0526 1.360 0.210
kurtosis -1.4554 0.563 -1.174
dstats()调用了mystats()函数，将其应用于数据框的每一栏中，再通过by()函数则可以得到am中的每一水平的概括统计量。
7.1.4 分组计算的扩展
doBy包和psych包提供了描述性统计量的函数。doBy包中的函数summaryBy()使用格式：
summaryBy(formula,data=dataframe,FUN=function)
formula接受以下格式：
var1+var2+...varN ~ groupvar1+groupvar2+...+groupvarN
在~左侧的变量是需要分析的数值型变量，而右侧的变量是类别型的分组变量。function可为任何内建或自编的R函数。
代码7-8 使用doBy包中的summaryBy()分组计算概述统计量
> library(doBy)
> summaryBy(mpg+hp+wt~am,data=mtcars,FUN=mystats)
am mpg.n mpg.mean mpg.stdev mpg.skew mpg.kurtosis hp.n hp.mean hp.stdev hp.skew hp.kurtosis wt.n
1 0 19 17.1 3.83 0.0140 -0.803 19 160 53.9 -0.0142 -1.210 19
2 1 13 24.4 6.17 0.0526 -1.455 13 127 84.1 1.3599 0.563 13
wt.mean wt.stdev wt.skew wt.kurtosis
1 3.77 0.777 0.976 0.142
2 2.41 0.617 0.210 -1.174
psych包中的describeBy()函数可计算和describe()相同的描述性统计量，只按照一个或多个分组的变量分层。
代码7-9 使用psych()包中的describeBy()分组计算概述统计量
> library(doBy)
> summaryBy(mpg+hp+wt~am,data=mtcars,FUN=mystats)
am mpg.n mpg.mean mpg.stdev mpg.skew mpg.kurtosis hp.n hp.mean hp.stdev hp.skew hp.kurtosis wt.n
1 0 19 17.1 3.83 0.0140 -0.803 19 160 53.9 -0.0142 -1.210 19
2 1 13 24.4 6.17 0.0526 -1.455 13 127 84.1 1.3599 0.563 13
wt.mean wt.stdev wt.skew wt.kurtosis
1 3.77 0.777 0.976 0.142
2 2.41 0.617 0.210 -1.174
7-9
> library(psych)
> myvars<-c("mpg","hp","wt")
> describeBy(mtcars[myvars],list(am=mtcars$am))
$`0`
vars n mean sd median trimmed mad min max range skew kurtosis se
mpg 1 19 17.15 3.83 17.30 17.12 3.11 10.40 24.40 14.00 0.01 -0.80 0.88
hp 2 19 160.26 53.91 175.00 161.06 77.10 62.00 245.00 183.00 -0.01 -1.21 12.37
wt 3 19 3.77 0.78 3.52 3.75 0.45 2.46 5.42 2.96 0.98 0.14 0.18
$`1`
vars n mean sd median trimmed mad min max range skew kurtosis se
mpg 1 13 24.39 6.17 22.80 24.38 6.67 15.00 33.90 18.90 0.05 -1.46 1.71
hp 2 13 126.85 84.06 109.00 114.73 63.75 52.00 335.00 283.00 1.36 0.56 23.31
wt 3 13 2.41 0.62 2.32 2.39 0.68 1.51 3.57 2.06 0.21 -1.17 0.17
attr(,"call")
by.data.frame(data = x, INDICES = group, FUN = describe, type = type)
分析：describeBy()函数不允许指定任意函数，所以它的普适性较低。
7.2 频数表和列联表
下面的数据来源于vcd包中的Arthritis数据集。数据展示了一项关节炎的新疗法的临床效果。
前几个观测如下：
> library(vcd)
载入需要的程辑包：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
7.2.1 生成频数表
R中提供了创建频数表和列联表的若干方法。重要函数如下：
table(var1,var2,... varN) 使用N个类别型变量（因子）创建一个N维列联表
xtabs(formula,data) 根据一个公式和一个矩阵或数据框创建一个N维向量
prop.table(table,margins) 依margins定义边际列表将表中的条目表示为分数形式
margin.table(table,margins) 依margins定义边际列表将表中的条目的和
addmargins(table,margins) 将概述边为margins（默认是求和结果）放入表中
ftable(table) 创建了一个紧凑的“平铺”式列联表
1. 一维列联表
> mytable<-with(Arthritis,table(Improved))
> mytable
Improved
None Some Marked
42 14 28
可以使用prop.table()将这些频数转化为比例值
> prop.table(mytable)
Improved
None Some Marked
0.500 0.167 0.333
使用prop.table()*100 转化为百分比
> prop.table(mytable)*100
Improved
None Some Marked
50.0 16.7 33.3
2.二维列联表
对于二维列联表，table()函数的使用格式为：
mytable<-table(A,B)
其中，A是变量，B是列变量，xtabs()函数还可以使用公式风格的输入创建列联表，格式
mytable<-xtabs(~ A+B,data=mydata)
其中，mydata是一个矩阵或数据框。
对于Arthritis数据，有：
> 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()语句中的第一个变量
接受安慰剂的个体中有显著改善有16%
接受治疗的个体中病情有显著改善的有51%
列和列比例可以这样计算：
> 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
下标2指代table()语句中的第二个变量
各单元格所占比例可按如下语句：
> prop.table(mytable)
Improved
Treatment None Some Marked
Placebo 0.34523810 0.08333333 0.08333333
Treated 0.15476190 0.08333333 0.25000000
可以使用addmargin()函数为这些表格添加边际和。
> addmargins(mytable)
Improved
Treatment None Some Marked Sum
Placebo 29 7 7 43
Treated 13 7 21 41
Sum 42 14 28 84
> addmargins(prop.table(mytable))
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()时，默认行为是为表中所有的变量创建边际和。
仅添加了各行的和，类似地：
> addmargins(prop.table(mytable,1),2)
Improved
Treatment None Some Marked Sum
Placebo 0.6744186 0.1627907 0.1627907 1.0000000
Treated 0.3170732 0.1707317 0.5121951 1.0000000
添加了各列的和。
> addmargins(prop.table(mytable,2),1)
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
使用gmodels包中的CrossTable()函数是创建二维联表的一种方法。
代码7-10 使用CrossTable生成二维列联表
> library(gmodels)
> CrossTable(Arthritis$Treatment,Arthritis$Improved)
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 | |
--------------------|-----------|-----------|-----------|-----------|
3.多维列联表
table()和 xtabs()都可以基于三个或更多的类别型变量生成多维列联表
margins.table() prop.table() addmargins() 函数可以推广到高于二维的情况
代码7-11 三维列联表
#生成了三维分组各单元格频数
> 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
#治疗情况（Treatment）*改善情况（Improved）的边际频数，由不同性别的单元格加和而成。
> margin.table(mytable,c(1,3))
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
#治疗情况（Treatment）*性别（Sex）的各类改善情况比例
> ftable(prop.table(mytable,c(1,2)))
Improved None Some Marked
Treatment Sex
Placebo Female 0.5938 0.2188 0.1875
Male 0.9091 0.0000 0.0909
Treated Female 0.2222 0.1852 0.5926
Male 0.5000 0.1429 0.3571
> ftable(addmargins(prop.table(mytable,c(1,2)),3))
Improved None Some Marked Sum
Treatment Sex
Placebo Female 0.5938 0.2188 0.1875 1.0000
Male 0.9091 0.0000 0.0909 1.0000
Treated Female 0.2222 0.1852 0.5926 1.0000
Male 0.5000 0.1429 0.3571 1.0000
列联表可以告诉你组成表格的各种变量组合的频数或比例。
7.2.2 独立性检验
1. 卡方独立性检验
可以使用chisq.test() 函数对二维表的行变量和列变量进行卡方独立性检验。
代码7-12 卡方独立性检验
> library(vcd)
> 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(~Improved+Sex,data = Arthritis)
> chisq.test(mytable)
Pearson's Chi-squared test
data: mytable
X-squared = 4.8407, df = 2, p-value = 0.08889
Warning message:
In chisq.test(mytable) : Chi-squared近似算法有可能不准
可以用fisher.test()函数进行Fisher精确检验。
Fisher精确检验的原假设是：边界固定的列联表中行与列是相互独立的。
其调用格式为fisher.test(mytable)，其中mytable是一个列联表。例如：
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()函数可用来进行Cochran-Mantel-Haenszel 检验。
检验假设：两个名义变量在第三个变量的每一层中都是条件独立的。下列代码可以检验治疗情况和改善情况在性别的每一水平下是否独立。
> 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.2.3 相关性的度量
vcd包中的assosctats()函数可以用来计算二维列联表的phi系数，列联系数和Cramer's V系数。
代码7-13 二维列联表的相关性度量
> library(vcd)
> mytable<-xtabs(~Treatment+Improved,data=Arthritis)
> assocstats(mytable)
X^2 df P(> X^2)
Likelihood Ratio 13.530 2 0.0011536
Pearson 13.055 2 0.0014626
Phi-Coefficient : NA
Contingency Coeff.: 0.367
Cramer's V : 0.394
7.2.4 结果的可视化
使用条形图来进行一维频数的可视化
7.3相关系数
相关系数可以用来描述定量变量之间的关系。相关系数符号（正负号）表明关系的方向（正相关或负相关），其值的大小表示关系的强弱程度（完全不相关是为0，完全相关时为1）
7.3.1相关的类型
1.Pearson、Spearman、Kendall相关
Pearson积差相关系数衡量了两个定量变量之间的线性相关程度。
Spearman等级相关系数则衡量分级定序变量之间的相关程度。
Kendall相关系数是一种非参数的等级相关度量。
cor()函数可以计算这三种相关系数，而cov()函数可以用来计算协方差。两个函数的参数有很多，其中与相关系数的计算有关的参数可以简化为：
cor(x,use= ,method= )
其中，x是矩阵或数据框
use指定缺失数据的处理方式。
method指定相关系数的类型。可选类型为Pearson、Spearman、Kendall相关
默认参数为，use="everything" , method="pearson"
代码7-14 协方差和相关系数
#计算方差和协方差
> states<-state.x77[,1:6]
> options(digits=3)
> cov(states)
Population Income Illiteracy Life Exp Murder HS Grad
Population 19931684 571230 292.868 -407.842 5663.52 -3551.51
Income 571230 377573 -163.702 280.663 -521.89 3076.77
Illiteracy 293 -164 0.372 -0.482 1.58 -3.24
Life Exp -408 281 -0.482 1.802 -3.87 6.31
Murder 5664 -522 1.582 -3.869 13.63 -14.55
HS Grad -3552 3077 -3.235 6.313 -14.55 65.24
#计算Pearson积差相关系数
> cor(states)
Population Income Illiteracy Life Exp Murder HS Grad
Population 1.0000 0.208 0.108 -0.0681 0.344 -0.0985
Income 0.2082 1.000 -0.437 0.3403 -0.230 0.6199
Illiteracy 0.1076 -0.437 1.000 -0.5885 0.703 -0.6572
Life Exp -0.0681 0.340 -0.588 1.0000 -0.781 0.5822
Murder 0.3436 -0.230 0.703 -0.7808 1.000 -0.4880
HS Grad -0.0985 0.620 -0.657 0.5822 -0.488 1.0000
#计算Spearman等级相关系数
> cor(states,method = "spearman")
Population Income Illiteracy Life Exp Murder HS Grad
Population 1.000 0.125 0.313 -0.104 0.346 -0.383
Income 0.125 1.000 -0.315 0.324 -0.217 0.510
Illiteracy 0.313 -0.315 1.000 -0.555 0.672 -0.655
Life Exp -0.104 0.324 -0.555 1.000 -0.780 0.524
Murder 0.346 -0.217 0.672 -0.780 1.000 -0.437
HS Grad -0.383 0.510 -0.655 0.524 -0.437 1.000
上述代码是在默认的情况的得到的结果是一个方阵（即所有变量之间两两相关）。
我们看到收入和高中毕业率之间存在很强的正相关；而文盲和预期寿命之间存在很强的负相关。
我们计算非方形的相关矩阵。例如：
> x<-states[,c("Population","Income","Illiteracy","HS Grad")]
> y<-states[,c("Life Exp","Murder")]
> cor(x,y)
Life Exp Murder
Population -0.0681 0.344
Income 0.3403 -0.230
Illiteracy -0.5885 0.703
HS Grad 0.5822 -0.488
2. 偏相关
偏相关是指在控制一个或多个定量或者变量的时，另外两个定量变量之间的相互关系。
你可以使用ggm包中的por()函数计算偏相关系数。格式：
pcor(u,s)
其中，u为数值向量。
s为变量的协方差阵。
> library(ggm)
> colnames(states)
[1] "Population" "Income" "Illiteracy" "Life Exp" "Murder" "HS Grad"
> pcor(c(1,5,2,3,6),cov(states))
[1] 0.346
上面例子中，控制了收入，文盲率和高中毕业率的影响时，人口和谋杀率之间的相关系数为0.346
7.3.2 相关性的显著性检验
在计算好相关系数后，对它进行统计性显著性检验。常用的原假设为变量之间不相关。可以使用cor.test() 函数对单个的Pearson、Spearman、Kendall相关系数进行检验。格式：
cor.tese(x,y,alternative = ,method= )
其中，x和y为要检验相关性的变量（取值为“two.side”、“less”、“greater”）；
alternative则用来指定进行双侧检验或单侧检验；
method指定要计算的相关类型（Pearson、Spearman、Kendall）
默认情况是，假设alternative="two.side"（总体相关系数不等于0）
代码7-15 检验某种相关系数的显著性
> cor.test(states[,3],states[,5])
Pearson's product-moment correlation
data: states[, 3] and states[, 5]
t = 7, df = 50, p-value = 1e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.528 0.821
sample estimates:
cor
0.703
上面的代码 检验了预期寿命和谋杀率的Peasson相关系数为0的原假设。假设总体的相关度为0，则预计在1000万次中只会有少于一次的机会见到0.703这样大的样本相关度。
corr.test() 函数可以为Pearson、Spearman、Kendall 相关呢计算相关矩阵和显著性水平。
代码7-16 通过corr.test计算相关矩阵并进行显著性检验
> library(psych)
> corr.test(x=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”（表示对缺失值执行行成对删除）
参数use=的取值可为“comple”（表示对缺失值执行删除）
可以看出人口数量和高中毕业率的相关系数为（-0.1），并不显著的不为0
其他显著性检验
psych包中的pcor.test()函数可以用来检验在控制一个或多个额外变量之间的条件独立性。格式：
pcor.test(r,q,n)
其中，r是由pcor()函数计算得到的偏相关系数，q为控制的变量数（以数值表示位置），n为样本大小。
psych包中的r.test()函数提供了多种实用的显著性检验方法。此函数可用来检验：
- 某种相关系数的显著性；
- 两个独立相关系数的差异是否显著
- 两个基于一个共享变量得到的非独立相关系数的差异是否显著
- 两个基于完全不同的变量得到的非独立相关系数的差异是否显著
7.3.3 相关系数的可视化以相关系数表示的二元关系可以通过散点图和散点图矩阵进行可视化。
7.4 t检验
7.4.1 独立样本的t检验
一个针对两组的独立样本t检验可以用于检验的两个总体的均值相等的假设。这里假设的两组数据是独立的，并且是从正态总体中抽得。格式：
t.test(y ~ x ,data)
其中，y是一个数值型变量，x是一个二分变量
调用格式：
t.test(y1,y2)
其中，y1和y2是数值型向量（即各组的结构变量）
data的取值为一个包含了这些变量的矩阵或数据框
下列代码，使用了一个假设方差不等的双侧检验。比较了南方和非南方各州的监禁概率：
> library(MASS)
> t.test(Prob~So,data = UScrime)
Welch Two Sample t-test
data: Prob by So
t = -4, df = 20, p-value = 7e-04
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.0385 -0.0119
sample estimates:
mean in group 0 mean in group 1
0.0385 0.0637
7.4.2 非独立样本的t检验
在两组的观测之间相关时，获得一个非独立组织设计。前—后测试设计或重复测量设计同样也会产生非独立的组。
非独立样本的t检验假定组间的差异呈正态分布。格式：
t.test(y1,y2,paired=TRUE)
其中，y1和y2为两个非独立组的数值向量。结果：
> library(MASS)
> sapply(UScrime[c("U1","U2")],function(x)(c(mean=mean(x),sd=sd(x))))
U1 U2
mean 95.5 33.98
sd 18.0 8.45
> with(UScrime,t.test(U1,U2,paired = TRUE))
Paired t-test
data: U1 and U2
t = 30, df = 50, p-value <2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
57.7 65.3
sample estimates:
mean of the differences
61.5
差异的均值（61.5）足够大，保证拒绝年长和年轻男性的平均失业率相同的假设。
7.4.3 多于两组的情况
使用方差分析（ANOVA）
7.5 组间差异的非参数检验
如果数据无法满足t检验和ANOVA参数假设，可以使用非参数的方法。
7.5.1 两组的比较
若两组数据独立，使用Wilcoxon秩和检验来评估观测是否要从相同的概率分布中抽得。格式：
wilcox.test(y~x,data)
其中，y是数值型变量，x是一个二分变量。格式：
wilcox.test(y1,y2)
其中，y1，y2为各组的结果变量。
使用Mann-Whitney U检验回答上一节中关于监禁的问题：
> with(UScrime,by(Prob,So,median))
So: 0
[1] 0.0382
----------------------------------------------------------------------------
So: 1
[1] 0.0556
> wilcox.test(Prob~So,data = UScrime)
Wilcoxon rank sum test
data: Prob by So
W = 80, p-value = 8e-05
alternative hypothesis: true location shift is not equal to 0
Wilcox符号秩检验是非独立样本t检验的一种非参数替代方法。适用于两组成对数据和无法保证正态的性假设的情境。
> sapply(UScrime[c("U1","U2")],median)
U1 U2
92 34
> with(UScrime,wilcox.test(U1,U2,paired = TRUE))
Wilcoxon signed rank test with continuity correction
data: U1 and U2
V = 1000, p-value = 2e-09
alternative hypothesis: true location shift is not equal to 0
含参的t检验和与其他作用相同的非参数检验得到了相同的结论。当t检验的假设合理时，参数检验的功效更强。
7.5.2 多于两组的比较
若无法满足ANOVA设计的假设，可以使用非参数方法来评估组间的差异。
如果各组独立，则Kruskal-Wallis检验将是一种实用的方法
格式： kruskal-wallis.test(y~A,data)
其中，y是一个数值型结果变量，A是一个拥有两个或多个的分组变量
如果各组不独立，那么Friedman检验会更合适：
friedman.test(y~A|B,data)
其中，y是数值型结果变量，A是一个分组变量，B是一个用以认定匹配测试的区组变量。
> states<-data.frame(state.region,state.x77)
> kruskal.test(Illiteracy~state.region,data = states)
Kruskal-Wallis rank sum test
data: Illiteracy by state.region
Kruskal-Wallis chi-squared = 20, df = 3, p-value = 5e-05
显著性检验的结果意味着美国的四个地区的文盲率各不相同（p<0.001）