这里是佳奥!继图形学习后,我们开始统计分析的部分。
在数据被组织成合适的形式后,我们也开始使用图形探索数据,而下一步通常就是使用数值描述每个变量的分布,接下来则是两两探索所选择变量之间的关系。其目的是回答如下问题:
1、各车型的油耗如何?特别是,在对车型的调查中,每加仑汽油行驶英里数的分布是什么样的?(均值、标准差、中位数、值域等。)2、在进行新药实验后,用药组和安慰剂组的治疗结果(无改善、一定程度的改善、显著的改善)相比如何?实验参与者的性别是否对结果有影响?3、收入和预期寿命的相关性如何?它是否明显不为零?4、美国的某些地区是否更有可能因为你犯罪而将你监禁?不同地区的差别是否在统计上显著?
本篇我们将评述用于生成基本的描述性统计量和推断统计量的R函数。
首先,我们将着眼于定量变量的位置和尺度的衡量方式。
然后我们将学习生成类别型变量的频数表和列联表的方法(以及连带的卡方检验)。
接下来,我们将考察连续型和有序型变量相关系数的多种形式。
最后,我们将转而通过参数检验(t检验)和非参数检验(Mann—Whitney U检验、Kruskal—Wallis检验)方法研究组间差异。
那么我们开始吧!
1 描述性统计分析
本节中,我们将关注分析连续型变量的中心趋势、变化性和分布形状的方法。
我们将使用第1篇中Motor Trend 杂志的车辆路试(mtcars)数据集。
我们的关注焦点是每加仑汽油行驶英里数(mpg)、马力(hp)和车重(wt)。
我们将首先查看所有32种车型的描述性统计量,然后按照变速箱类型(am)和汽缸数(cyl)考察描述性统计量。
变速箱类型是一个以0表示自动挡、1表示手动挡来编码的二分变量,而汽缸数可为4、5或6。
1.1 方法云集
让我们从基础安装中包含的函数开始,然后查看那些用户贡献包中的扩展函数。
对于基础安装,你可以使用summary()函数来获取描述性统计量:
summary( )函数提供了最小值、最大值、四分位数和数值型变量的均值,以及因子向量和逻辑型向量的频数统计。
可以使用第5篇中的apply( )函数或sapply( )函数计算所选择的任意描述性统计量。
对于sapply()函数,其使用格式为:
其中的x是数据框(或矩阵),FUN为一个任意的函数。
如果指定了options,它们将被传递给FUN。可以在这里插入的典型函数有mean、sd、var、min、max、median、length、range和quantile。
函数fivenum( )可返回图基五数总括(Tukey’s five-number summary,即最小值、下四分位数、中位数、上四分位数和最大值)。
但是基础安装并没有提供偏度和峰度的计算函数,需要自己安装。
下列代码的示例计算了若干描述性统计量,其中包括偏度和峰度:
对于样本中的车型,每加仑汽油行驶英里数的平均值(mean)为20.1,标准差为6.0。分布呈现右偏(偏度+0.61),并且较正态分布稍平(峰度-0.37)。
请注意,如果只希望单纯地忽略缺失值,那么应当使用sapply(mtcars[vars], mystats,na.omit=TRUE)。
na.omit( ):移除所有含有缺失值的观测(行删除)。
扩展
若干用户贡献包都提供了计算描述性统计量的函数,其中包括Hmisc、pastecs和psych。
Hmisc包中的describe()函数可返回变量和观测的数量、缺失值和唯一值的数目、平均值、分位数,以及五个最大的值和五个最小的值。
通过Hmisc包中的describe()函数计算描述性统计量:
pastecs包中有一个名为stat.desc( )的函数,它可以计算种类繁多的描述性统计量。使用格式为:
其中的x是一个数据框或时间序列。
若basic=TRUE(默认值),则计算其中所有值、空值、缺失值的数量,以及最小值、最大值、值域,还有总和。
若desc=TRUE(同样也是默认值),则计算中位数、平均数、平均数的标准误、平均数置信度为95%的置信区间、方差、标准差以及变异系数。
最后,若norm=TRUE(不是默认的),则返回正态分布统计量,包括偏度和峰度(以及它们的统计显著程度)和Shapiro–Wilk正态检验结果。
这里使用了p值来计算平均数的置信区间(默认置信度为0.95)。
通过pastecs包中的stat.desc()函数计算描述性统计量:
psych包也拥有一个名为describe( )的函数,它可以计算非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均值的标准误。
通过psych包中的describe( )计算描述性统计量:
注意
在前面的示例中,psych包和Hmisc包均提供了名为describe( )的函数。
R如何知道该使用哪个呢?简言之,最后载入的程序包优先。
在这里,psych在Hmisc之后被载入,然后显示了一条信息,提示Hmisc包中的describe( )函数被psych包中的同名函数所屏蔽(masked)。
键入describe( )后,R在搜索这个函数时将首先找到psych包中的函数并执行它。
如果你想改而使用Hmisc包中的版本,可以键入Hmisc::describe(mt)。这个函数仍然在那里。你只是需要给予R更多信息以找到它。
1.2 分组计算描述性统计量
使用aggregate( )分组获取描述性统计量:
注意list(am=mtcars$am)的使用。
如果使用的是list(mtcars$am),则am列将被标注为Group.1而不是am。使用这个赋值指定了一个更有帮助的列标签。
如果有多个分组变量,可以使用by=list(name1=groupvar1, name2=groupvar2, ... , groupvarN) 这样的语句。
不过,aggregate( )仅允许在每次调用中使用平均数、标准差这样的单返回值函数。
它无法一次返回若干个统计量。
要完成这项任务,可以使用by( )函数。格式为:
其中data是一个数据框或矩阵,INDICES是一个因子或因子组成的列表,定义了分组,FUN是任意函数。
使用by( )分组计算描述性统计量:
(开始报错,,,,,)啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊
问了舍友以后,得到的解答是,mean( )函数无法计算数据框,而我输入的mtcars为数据框格式,修改运用apply( )函数以后可以运行了。没事用用?mean。
扩展
doBy包和psych包也提供了分组计算描述性统计量的函数。
doBy包中summaryBy()函数的使用格式为:
其中的formula接受以下的格式:
在~左侧的变量是需要分析的数值型变量,而右侧的变量是类别型的分组变量。
function可为任何内建或用户自编的R函数。
使用doBy包中的summaryBy()分组计算概述统计量:
psych包中的describe.by( )函数可计算和describe相同的描述性统计量,只是按照一个或多个分组变量分层。
使用psych包中的describe.by( )分组计算概述统计量:
与前面的示例不同,describe.by( )函数不允许指定任意函数,所以它的普适性较低。
若存在一个以上的分组变量,可以使用list(groupvar1, groupvar2, ... , groupvarN)来表示它们。但这仅在分组变量交叉后不出现空白单元时有效。
最后,可以使用reshape包灵活地按组导出描述性统计量。
首先,使用:
融合数据框。其中的dataframe包含着数据,y是一个向量,指明了要进行概述的数值型变量(默认使用所有变量),而g是由一个或多个分组变量组成的向量。
然后使用:
重铸数据。分组变量以+号分隔,这里的variable只取其字面含义(即仅表示重铸后数据框中的变量variable),而FUN是一个任意函数。
在本节的最后一个例子中,我们将运用数据重塑的方法来取得由变速箱类型与汽缸数形成的每个亚组的描述性统计量。
我们要获取的描述性统计量是样本大小、平均数和标准差。
通过reshape包分组计算概述统计量:
书本认为这种方式最为简洁动人。
当然,数据分析人员对于展示哪些描述性统计量以及结果采用什么格式都有着自己的偏好,这也许就是有如此多不同方法的原因。
1.3 结果的可视化
分布特征的数值刻画的确很重要,但是这并不能代替视觉呈现。
对于定量变量,我们有直方图、密度图、箱线图和点图。
它们都可以让我们洞悉那些依赖于观察一小部分描述性统计量时忽略的细节。
2 频数表和列联表
在本节中,我们将着眼于类别型变量的频数表和列联表,以及相应的独立性检验、相关性的度量、图形化展示结果的方法。
我们除了使用基础安装中的函数,还将连带使用vcd包和gmodels包中的函数。
下面的示例中,假设A、B和C代表类别型变量。
本节中的数据来自vcd包中的Arthritis数据集,表示了一项风湿性关节炎新疗法的双盲临床实验的结果。前几个观测是这样的:
治疗情况(安慰剂治疗、用药治疗)、性别(男性、女性)和改善情况(无改善、一定程度的改善、显著改善)均为类别型因子。
接下来,我们将使用此数据创建频数表和列联表(交叉的分类)。
2.1 生成频数表
用于创建和处理列联表的函数:
接下来,我们将逐个使用以上函数来探索类别型变量。
我们首先考察简单的频率表,接下来是二维列联表,最后是多维列联表。
第一步是使用table( )或xtabs( )函数创建一个表,然后使用其他函数处理它。
1.一维列联表
可以使用table()函数生成简单的频数统计表。
可以用prop.table( )将这些频数转化为比例值:
或使用prop.table( )*100转化为百分比:
这里可以看到,有50%的研究参与者获得了一定程度或者显著的改善(16.7 + 33.3)。
2.二维列联表
对于二维列联表,table()函数的使用格式为:
其中的A是行变量,B是列变量。除此之外,xtabs()函数还可使用公式风格的输入创建列联表,格式为:
其中的mydata是一个矩阵或数据框。
总的来说,要进行交叉分类的变量应出现在公式的右侧(即~符号的右方),以+作为分隔符。
若某个变量写在公式的左侧,则其为一个频数向量(在数据已经被表格化时很有用)。
对于Arthritis数据,有:
可以使用margin.table( )和prop.table( )函数分别生成边际频数和比例。行和与行比例可以这样计算:
下标1指代table( )语句中的第一个变量。
观察表格可以发现,与接受安慰剂的个体中有显著改善的16%相比,接受治疗的个体中的51%的个体病情有了显著的改善。
列和与列比例可以这样计算:
这里的下标2指代table( )语句中的第二个变量。
各单元格所占比例可用如下语句获取:
可以使用addmargins( )函数为这些表格添加边际和。
例如,以下代码添加了各行的和与各列的和:
在使用addmargins( )时,默认行为是为表中所有的变量创建边际和。作为对照:
仅添加了各行的和。类似地:
添加了各列的和。在表中可以看到,有显著改善患者中的25%是接受安慰剂治疗的。
注意
table( )函数默认忽略缺失值(NA)。要在频数统计中将NA视为一个有效的类别,请设定参数useNA="ifany"。
使用gmodels包中的CrossTable( )函数是创建二维列联表的第三种方法。
CrossTable()函数仿照SAS中PROC FREQ或SPSS中CROSSTABS的形式生成二维列联表。
使用CrossTable生成二维列联表:
CrossTable( )函数有很多选项,可以做许多事情:计算(行、列、单元格)的百分比;指定小数位数;进行卡方、Fisher和McNemar独立性检验;计算期望和(皮尔逊、标准化、调整的标准化)残差;将缺失值作为一种有效值;进行行和列标题的标注;生成SAS或SPSS风格的输出。
3.多维列联表
table( )和xtabs( )都可以基于三个或更多的类别型变量生成多维列联表。
margin.table()、prop.table()和addmargins()函数可以自然地推广到高于二维的情况。
另外,ftable()函数可以以一种紧凑而吸引人的方式输出多维列联表。
三维列联表:
第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。例如:
2.2 独立性检验
R提供了多种检验类别型变量独立性的方法。本节中描述的三种检验分别为卡方独立性检验、Fisher精确检验和Cochran-Mantel–Haenszel检验。
1.卡方独立性检验
使用chisq.test( )函数对二维表的行变量和列变量进行卡方独立性检验:
在结果1中,患者接受的治疗和改善的水平看上去存在着某种关系(p < 0.01)。
而患者性别和改善情况之间却不存在关系(p > 0.05) 。这里的p值表示从总体中抽取的样本行变量与列变量是相互独立的概率。由于1的概率值很小,所以我们拒绝了治疗类型和治疗结果相互独立的原假设。
由于2的概率不够小,故没有足够的理由说明治疗结果和性别之间是不独立的。
产生警告信息的原因是,表中的6个单元格之一(男性 - 一定程度上的改善)有一个小于5的值,这可能会使卡方近似无效。
2.Fisher精确检验
可以使用fisher.test ()函数进行Fisher精确检验。
Fisher精确检验的原假设是:边界固定的列联表中行和列是相互独立的。
其调用格式为fisher.test(mytable),其中的mytable是一个二维列联表。
与许多统计软件不同的是,这里的fisher.test( )函数可以在任意行列数大于等于2的二维列联表上使用,但不能用于2×2的列联表。
3.Cochran-Mantel-Haenszel检验
mantelhaen.test( )函数可用来进行Cochran—Mantel—Haenszel卡方检验,其原假设是,两个名义变量在第三个变量的每一层中都是条件独立的。
下列代码可以检验治疗情况和改善情况在性别的每一水平下是否独立。此检验假设不存在三阶交互作用(治疗情况×改善情况×性别)。
结果表明,患者接受的治疗与得到的改善在性别的每一水平下并不独立(即,分性别来看,用药治疗的患者较接受安慰剂的患者有了更多的改善)。
2.3 相关性的度量
vcd包中的assocstats( )函数可以用来计算二维列联表的phi系数、列联系数和Cramer’s V系数。
二维列联表的相关性度量:
总体来说,较大的值意味着较强的相关性。vcd包也提供了一个kappa()函数,可以计算混淆矩阵的Cohen’s kappa值以及加权的kappa值。
2.4 结果的可视化
通常,我们会使用条形图进行一维频数的可视化。
vcd包中拥有优秀的、用于可视化多维数据集中类别型变量间关系的函数,可以绘制马赛克图和关联图。
最后,ca包中的对应分析函数允许使用多种几何表示可视地探索列联表中行和列之间的关系。
2.5 将表转换为扁平格式
在已经拥有一个列联表却需要原始的数据时怎么办?举例来说,假设有以下列联表:
但是需要的格式:
R中的许多统计函数接受的是后一种格式而不是前一种格式。
通过table2flat将表转换为扁平格式:
此函数可以接受一个R中的表格(行列数任意)并返回一个扁平格式的数据框。你也可以使用这个函数来输入已发表的研究中的表格。举例来说,假设你在某期刊上看到了表7-2,并想以扁平格式将其保存到R中。
使用table2flat( )函数转换已发表的数据:
3 相关
相关系数可以用来描述定量变量之间的关系。相关系数的符号(+)表明关系的方向(正相关或负相关),其值的大小表示关系的强弱程度(完全不相关时为0,完全相关时为1)。
本节中,我们将关注多种相关系数和相关性的显著性检验。
我们将使用R基础安装中的state.x77数据集,它提供了美国50个州在1977年的人口、收入、文盲率、预期寿命、谋杀率和高中毕业率数据。
3.1 相关的类型
R可以计算多种相关系数,包括Pearson相关系数、Spearman相关系数、Kendall相关系数、偏相关系数、多分格(polychoric)相关系数和多系列(polyserial)相关系数。
1.Pearson、Spearman和Kendall相关
Pearson积差相关系数衡量了两个定量变量之间的线性相关程度。
Spearman等级相关系数则衡量分级定序变量之间的相关程度。
Kendall’s Tau相关系数也是一种非参数的等级相关度量。
cor( )函数可以计算这三种相关系数,而cov( )函数可用来计算协方差。
两个函数的参数有很多,其中与相关系数的计算有关的参数可以简化为:
cor和cov的参数:
默认参数为use="everything"和method="pearson"
协方差和相关系数:
首个语句计算了方差和协方差,
第二个语句则计算了Pearson积差相关系数,
而第三个语句计算了Spearman等级相关系数。
举例来说,我们可以看到收入和高中毕业率之间存在很强的正相关,而文盲率和预期寿命之间存在很强的负相关。
请注意,在默认情况下得到的结果是一个方阵(所有变量之间两两计算相关)。
我们同样可以计算非方形的相关矩阵。观察以下示例:
当对某一组变量与另外一组变量之间的关系感兴趣时,cor( )函数的这种用法是非常实用的。
注意,上述结果并未指明相关系数是否显著不为0(即,根据样本数据是否有足够的证据得出总体相关系数不为0的结论)。
由于这个原因,需要对相关系数进行显著性检验。
2.偏相关
偏相关是指在控制一个或多个定量变量时,另外两个定量变量之间的相互关系。
可以使用ggm包中的pcor( )函数计算偏相关系数。函数调用格式为:
其中的u是一个数值向量,前两个数值表示要计算相关系数的变量下标,其余的数值为条件变量(即要排除影响的变量)的下标。
S为变量的协方差阵。
这个示例有助于阐明用法:
本例中,在控制了收入、文盲率和高中毕业率的影响时,人口和谋杀率之间的相关系数为0.346。
3.其他类型的相关
polycor包中的hetcor( )函数可以计算一种混合的相关矩阵,其中包括数值型变量的Pearson积差相关系数、数值型变量和有序变量之间的多系列相关系数、有序变量之间的多分格相关系数以及二分变量之间的四分相关系数。
多系列、多分格和四分相关系数都假设有序变量或二分变量由潜在的正态分布导出。
3.2 相关性的显著性检验
在计算好相关系数以后,如何对它们进行统计显著性检验呢?
常用的原假设为变量间不相关(即总体的相关系数为0)。
可以使用cor.test()函数对单个的Pearson、Spearman和Kendall相关系数进行检验。
简化后的使用格式为:
其中的x和y为要检验相关性的变量,
alternative则用来指定进行双侧检验或单侧检验(取值为"two.side"、"less"或"greater"),
而method用以指定要计算的相关类型("pearson"、"kendall"或"spearman")。
当研究的假设为总体的相关系数小于0时,请使用alternative="less"。
在研究的假设为总体的相关系数大于0时,应使用alternative="greater"。
在默认情况下,假设为alternative="two.side"(总体相关系数不等于0)
检验某种相关系数的显著性:
这段代码检验了预期寿命和谋杀率的Pearson相关系数为0的原假设。假设总体的相关度为0,则预计在一千万次中只会有少于一次的机会见到0.703这样大的样本相关度(即p = 1.258e-08)。由于这种情况几乎不可能发生,所以你可以拒绝原假设,从而支持了要研究的猜想,即预期寿命和谋杀率之间的总体相关度不为0。
遗憾的是,cor.test每次只能检验一种相关关系。但幸运的是,psych包中提供的corr.test()函数可以一次做更多事情。corr.test()函数可以为Pearson、Spearman或Kendall相关计算相关矩阵和显著性水平
通过corr.test计算相关矩阵并进行显著性检验:
参数use=的取值可为"pairwise"或"complete"(分别表示对缺失值执行成对删除或行删除)。参数method=的取值可为"pearson"(默认值)、"spearman"或"kendall"。这里可以看到,人口数量和高中毕业率的相关系数(-0.10)并不显著地不为0(p = 0.5)。
其他显著性检验
在4.1节中,我们关注了偏相关系数。在多元正态性的假设下,ggm包中的pcor.test( )函数可以用来检验在控制一个或多个额外变量时两个变量之间的条件独立性。使用格式为:
其中的r是由pcor( )函数计算得到的偏相关系数,
q为要控制的变量数(以数值表示位置),
n为样本大小。
与此同时,psych包中的r.test( )函数提供了多种实用的显著性检验方法。此函数可用来检验:某种相关系数的显著性;两个独立相关系数的差异是否显著;两个基于一个共享变量得到的非独立相关系数的差异是否显著;两个基于完全不同的变量得到的非独立相关系数的差异是否显著。
3.3 相关关系的可视化
以相关系数表示的二元关系可以通过散点图和散点图矩阵进行可视化。
4 t检验
在研究中最常见的行为就是对两个组进行比较。
接受某种新药治疗的患者是否较使用某种现有药物的患者表现出了更大程度的改善?
某种制造工艺是否较另外一种工艺制造出的不合格品更少?
两种教学方法中哪一种更有效?
如果结果变量是类别型的,那么可以直接使用3节中阐述的方法。
这里我们将关注结果变量为连续型的组间比较,并假设其呈正态分布。
为了阐明方法,我们将使用MASS包中的UScrime数据集。
它包含了1960年美国47个州的刑罚制度对犯罪率影响的信息。我们感兴趣的结果变量为Prob(监禁的概率)、U1(14~24岁年龄段城市男性失业率)和U2(35~39岁年龄段城市男性失业率)。类别型变量So(指示该州是否位于南方的指示变量)将作为分组变量使用。
4.1 独立样本的t检验
如果你在美国的南方犯罪,是否更有可能被判监禁?
我们比较的对象是南方和非南方各州,因变量为监禁的概率。
一个针对两组的独立样本t检验可以用于检验两个总体的均值相等的假设。
这里假设两组数据是独立的,并且是从正态总体中抽得。检验的调用格式为:
其中的y是一个数值型变量,x是一个二分变量。调用格式或为:
其中的y1和y2为数值型向量(即各组的结果变量)。
可选参数data的取值为一个包含了这些变量的矩阵或数据框。
与其他多数统计软件不同的是,这里的t检验默认假定方差不相等,并使用Welsh的修正自由度。
可以添加一个参数var.equal=TRUE以假定方差相等,并使用合并方差估计。
默认的备择假设是双侧的(即均值不相等,但大小的方向不确定)。
可以添加一个参数alternative="less"或alternative="greater"来进行有方向的检验。
在下列代码中,我们使用了一个假设方差不等的双侧检验,比较了南方(group 1)和非南方(group 0)各州的监禁概率:
我们可以拒绝南方各州和非南方各州拥有相同监禁概率的假设p-value = 0.0006506(p < .001)。
注意
由于结果变量是一个比例值,你可以在执行t检验之前尝试对其进行正态化变换。
在本例中,所有对结果变量合适的变换(Y/1-Y、log(Y/1-Y)、arcsin(Y)、arcsin(sqrt(Y))都会将检验引向相同的结论。
4.2 非独立样本的t检验
再举个例子:较年轻(14~24岁)男性的失业率是否比年长(35~39岁)男性的失业率更高?
在这种情况下,这两组数据并不独立。
不能说亚拉巴马州的年轻男性和年长男性的失业率之间没有关系。
在两组的观测之间相关时,获得的是一个非独立组设计(dependent groups design)。
前—后测设计(pre-post design)或重复测量设计(repeated measures design)同样也会产生非独立的组。
非独立样本的t检验假定组间的差异呈正态分布。对于本例,检验的调用格式为:
其中的y1和y2为两个非独立组的数值向量。结果如下:
差异的均值(61.5)足够大,可以拒绝年长和年轻男性的平均失业率相同的假设。
年轻男性的失业率更高。事实上,若总体均值相等,获取一个差异如此大的样本的概率小于0.000 000 000 000 000 22( 2.2e-16)。
4.3 多于两组的情况
如果能够假设数据是从正态总体中独立抽样而得的,那么可以使用方差分析(ANOVA)。
5 组间差异的非参数检验
如果数据无法满足t检验或ANOVA的参数假设,可以转而使用非参数方法。
5.1 两组的比较
若两组数据独立,可以使用Wilcoxon秩和检验(更广为人知的名字是Mann–Whitney U检验)来评估观测是否是从相同的概率分布中抽得的(即,在一个总体中获得更高得分的概率是否比另一个总体要大)。
调用格式为:
其中的y是数值型变量,而x是一个二分变量。调用格式或为:
其中的y1和y2为各组的结果变量。
可选参数data的取值为一个包含了这些变量的矩阵或数据框。
默认进行一个双侧检验。可以添加参数exact来进行精确检验,指定alternative="less"或alternative="greater"进行有方向的检验。
如果使用Mann–Whitney U检验回答上一节中关于监禁率的问题,将得到这些结果:
可以再次拒绝南方各州和非南方各州监禁率相同的假设(p < 0.001)。
Wilcoxon符号秩检验是非独立样本t检验的一种非参数替代方法。
它适用于两组成对数据和无法保证正态性假设的情境。
调用格式与Mann–Whitney U检验完全相同,不过还可以添加参数paired=TRUE。
让我们用它解答上一节中的失业率问题:
再次得到了与配对t检验相同的结论。
在本例中,含参的t检验和与其作用相同的非参数检验得到了相同的结论。
当t检验的假设合理时,参数检验的功效更强(更容易发现存在的差异)。
而非参数检验在假设非常不合理时(如对于等级有序数据)更适用。
5.2 多于两组的比较
在要比较的组数多于两个时,必须转而寻求其他方法。
考虑4节中的state.x77数据集。它包含了美国各州的人口、收入、文盲率、预期寿命、谋杀率和高中毕业率数据。
如果想比较美国四个地区(东北部、南部、中北部和西部)的文盲率,应该怎么做呢?
这称为单向设计(one-waydesign),我们可以使用参数或非参数的方法来解决这个问题。
如果无法满足ANOVA设计的假设,那么可以使用非参数方法来评估组间的差异。
如果各组独立,则Kruskal—Wallis检验将是一种实用的方法。
如果各组不独立(如重复测量设计或随机区组设计),那么Friedman检验会更合适。
Kruskal–Wallis检验的调用格式为:
其中的y是一个数值型结果变量,A是一个拥有两个或更多水平的分组变量(grouping variable)。
(若有两个水平,则它与Mann–Whitney U检验等价。)
而Friedman检验的调用格式为:
其中的y是数值型结果变量,A是一个分组变量,而B是一个用以认定匹配观测的区组变量(blockingvariable)。
在以上两例中,data皆为可选参数,它指定了包含这些变量的矩阵或数据框。
让我们利用Kruskal–Wallis检验回答文盲率的问题。
首先,将地区的名称添加到数据集中。这些信息包含在随R基础安装分发的state.region数据集中:
现在就可以进行检验了:
显著性检验的结果意味着美国四个地区的文盲率各不相同(p <0.001)。
虽然可以拒绝不存在差异的原假设,但这个检验并没有告诉你哪些地区显著地与其他地区不同。
要回答这个问题,可以使用Mann–Whitney U检验每次比较两组数据。
一种更为优雅的方法是在控制犯第一类错误的概率(发现一个事实上并不存在的差异的概率)的前提下,执行可以同步进行的多组比较,这样可以直接完成所有组之间的成对比较。
npmc包提供了所需要的非参数多组比较程序。
此包中的npmc()函数接受的输入为一个两列的数据框,其中一列名为var(因变量),另一列名为class(分组变量)。
非参数多组比较:
调用了npmc的语句生成了六对统计比较结果(东北部对南部、东北部对中北部、东北部对西部、南部对中北部、南部对西部,以及中北部对西部) 。
可以从双侧的p值(p.value.2s)看出南部与其他三个地区显著不同,而其他三个地区之间并没有什么不同。
在处可以看到南部的文盲率中间值更高。
注意,npmc在计算积分时使用了随机数,所以每次计算的结果会有轻微的不同。
6 组间差异的可视化
在4节和5节中,我们关注了进行组间比较的统计方法。
使用视觉直观地检查组间差异,同样是全面的数据分析策略中的一个重要组成部分。
它允许我们评估差异的量级、甄别出任何会影响结果的分布特征(如偏倚、双峰或离群点)并衡量检验假定的合理程度。
R中提供了许多比较组间数据的图形方法,其中包括第六篇中讲解的箱线图(简单箱线图、含凹槽的箱线图、小提琴图)、节中叠加的核密度图,以及在第9篇中讨论的评估检验假定的图形方法。
7 小结
在本章中,我们评述了R中用于生成统计概要和进行假设检验的函数。
我们关注了样本统计量和频数表、独立性检验和类别型变量的相关性度量、定量变量的相关系数(和连带的显著性检验)以及两组或更多组定量结果变量的比较。
下一篇中,我们将探索一元回归和多元回归,讨论的焦点在于如何理解一个预测变量(一元回归)或多个预测变量(多元回归)与某个被预测变量或效标变量(criterion variable)之间的关系。
终于,我们将迈向R语言的中级方法,那么我们,下篇再见!
网友评论