《R语言实战》自学笔记69-重抽样和自助法

Python019

《R语言实战》自学笔记69-重抽样和自助法,第1张

数据准备

许多实际情况中统计假设(假定观测数据抽样自正态分布或者其他性质较好的理论分布)并不一定满足,比如数据抽样于未知或混合分布、样本量过小、存在离群点、基于理论分布设计合适的统计检验过于复杂且数学上难以处理等情况,这时基于随机化和重抽样的统计方法就可派上用场。

置换检验的定义

置换检验(Permutation test),也称随机化检验或重随机化检验,是Fisher于20世纪30年代提出的一种基于大量计算(computationally intensive),利用样本数据的全(或随机)排列,进行统计推断的方法,因其对总体分布自由,应用较为广泛,特别适用于总体分布未知的小样本资料,以及某些难以用常规方法分析资料的假设检验问题。

置换检验的原理

1、提出原假设,比如XX处理后结果没有变化

2、计算统计量,如两组的均值之差,记作t0

3、将所有样本放在一起,然后随机排序进行分组,再计算其统计量t1

4、重复第3步骤,直至所有排序可能性都齐全(比如有A组有n样本,B组有m样本,则总重复次数相当于从n+m中随机抽取n个的次数),得到一系列的统计量(t1-tn)

5、最后将这些统计量按照从小到大排序,构成抽样分布,再看t0是否落在分布的置信区间内(如95%置信区间),这时候可计算一个P值(如果抽样总体1000次统计量中大于t0的有10个,则估计的P值为10/1000=0.01),落在置信区间外则拒绝原假设

6、如果第3步骤是将所有可能性都计算了的话,则是精确检验;如果只取了计算了部分组合,则是近似结果,这时一般用蒙特卡罗模拟(Monte Carlo simulation)的方法进行置换检验

7、置换检验和参数检验都计算了统计量,但是前者是跟置换观测数据后获得的经验分布进行比较,后者则是跟理论分布进行比较。

请牢记:置换检验都是使用伪随机数来从所有可能的排列组合中进行抽样(当做近似检验时)。因此,每次检验的结果都有所不同。

coin包提供了一个进行置换检验的一般性框架。通过该包,你可以回答如下问题。

响应值与组的分配独立吗?

两个数值变量独立吗?

两个类别型变量独立吗?

表12-2列出来的每个函数都是如下形式:

function_name(formula, data, distribution=)

其中:

 formula描述的是要检验变量间的关系。示例可参见表12-2;

 data是一个数据框;

 distribution指定经验分布在零假设条件下的形式,可能值有exact,asymptotic和

approximate。

若distribution = "exact",那么在零假设条件下,分布的计算是精确的(即依据所有可能的排列组合)。当然,也可以根据它的渐进分布(distribution = "asymptotic")或蒙特卡洛重抽样(distribution = "approxiamate(B = #)")来做近似计算,其中#指所需重复的次数。

distribution = "exact"当前仅可用于两样本问题。

传统t检验表明存在显著性差异(p <0.05),而精确检验却表明差异并不显著(p >0.072)。

第7章我用自己的数据进行了t检验,对比一下传统t检验和置换检验,结果如下:

两种检验方式下结果都是显著的

Wilcoxon-Mann-Whitney U检验

coin包规定所有的类别型变量都必须以因子形式编码。

wilcox.test()默认计算的也是精确分布。

K样本检验的置换检验

通过chisq_test()或cmh_test()函数,我们可用置换检验判断两类别型变量的独立性。 当数据可根据第三个类别型变量进行分层时,需要使用后一个函数。若变量都是有序型,可使用lbl_test()函数来检验是否存在线性趋势。

卡方独立性检验

卡方独立性检验的置换检验

你可能会有疑问,为什么需要把变量Improved从一个有序因子变成一个分类因子?(好问题!)这是因为,如果你用有序因子,coin()将会生成一个线性与线性趋势检验,而不是卡方检验。

结果解读:两种检验下p值都是小于0.05,说明Treatment和Improved之间相互不独立

自己数据的演示

结果解读:p值均为1,表明nitrogen和variety相互独立。

spearman_test()函数提供了两数值变量的独立性置换检验。

当处于不同组的观测已经被分配得当,或者使用了重复测量时,样本相关检验便可派上用场。

对于两配对组的置换检验,可使用wilcoxsign_test()函数;多于两组时,使用friedman_test()函数。

自己数据演示

lmPerm包可做线性模型的置换检验。比如lmp()和aovp()函数即lm()和aov()函数的修改版,能够进行置换检验,而非正态理论检验。

lmp()和aovp()函数的参数与lm()和aov()函数类似,只额外添加了perm =参数。

perm =选项的可选值有"Exact"、"Prob"或"SPR"。Exact根据所有可能的排列组合生成精确检验。Prob从所有可能的排列中不断抽样,直至估计的标准差在估计的p值0.1之下,判停准则由可选的Ca参数控制。SPR使用贯序概率比检验来判断何时停止抽样。注意,若观测数大于10,perm = "Exact"将自动默认转为perm = "Prob",因为精确检验只适用于小样本问题。

简单线性回归的置换检验

R语言实战的例子:

多项式回归的置换检验

R语言实战的例子:

自己数据集的例子:

R语言实战的例子:

自己数据集的例子:

当两种方法所得结果不一致时,你需要更加谨慎地审视数据,这很可能是因为违反了正态性假设或者存在离群点。

R语言实战的例子:

自己数据集的例子:

R语言实战的例子:

自己数据集的例子:

R语言实战的例子:

自己数据集的例子:

值得注意的是,当将aovp()应用到方差分析设计中时,它默认使用唯一平方和法(SAS也称为类型III平方和)。每种效应都会依据其他效应做相应调整。R中默认的参数化方差分析设计使用的是序贯平方和(SAS是类型I平方和)。每种效应依据模型中先出现的效应做相应调整。对于平衡设计,两种方法结果相同,但是对于每个单元格观测数不同的不平衡设计,两种方法结果则不同。不平衡性越大,结果分歧越大。若在aovp()函数中设定seqs = TRUE,可以生成你想要的序贯平方和。

你可能已经注意到,基于正态理论的检验与上面置换检验的结果非常接近。在这些问题中数据表现非常好,两种方法结果的一致性也验证了正态理论方法适用于上述示例。

当然,置换检验真正发挥功用的地方是处理非正态数据(如分布偏倚很大)、存在离群点、样本很小或无法做参数检验等情况。不过,如果初始样本对感兴趣的总体情况代表性很差,即使是置换检验也无法提高推断效果。

置换检验主要用于生成检验零假设的p值,它有助于回答“效应是否存在”这样的问题。不过,置换方法对于获取置信区间和估计测量精度是比较困难的。幸运的是,这正是自助法大显神通的地方。

所谓自助法,即从初始样本重复随机替换抽样,生成一个或一系列待检验统计量的经验分布。 无需假设一个特定的理论分布,便可生成统计量的置信区间,并能检验统计假设。

倘若你假设均值的样本分布不是正态分布,该怎么办呢?可使用自助法。

(1)从样本中随机选择10个观测,抽样后再放回。有些观测可能会被选择多次,有些可能一直都不会被选中。

(2)计算并记录样本均值。

(3)重复1和2一千次。

(4)将1000个样本均值从小到大排序。

(5)找出样本均值2.5%和97.5%的分位点。此时即初始位置和最末位置的第25个数,它们就限定了95%的置信区间。

boot包扩展了自助法和重抽样的相关用途。你可以对一个统计量(如中位数)或一个统计量向量(如一列回归系数)使用自助法。

一般来说,自助法有三个主要步骤。

(1)写一个能返回待研究统计量值的函数。如果只有单个统计量(如中位数),函数应该返回一个数值;如果有一列统计量(如一列回归系数),函数应该返回一个向量。

(2)为生成R中自助法所需的有效统计量重复数,使用boot()函数对上面所写的函数进行处理。

(3)使用boot.ci()函数获取第(2)步生成的统计量的置信区间。

主要的自助法函数是boot(),它的格式为:

bootobject <- boot(data=, statistic=, R=, ...)

参数见下表:

boot()函数调用统计量函数R次,每次都从整数1:nrow(data)中生成一列有放回的随机指标,这些指标被统计量函数用来选择样本。统计量将根据所选样本进行计算,结果存储在bootobject中。

你可以用bootobject t0和bootobject t来获取这些元素。

一旦生成了自助样本,可通过print()和plot()来检查结果。如果结果看起来还算合理, 使用boot.ci()函数获取统计量的置信区间。格式如下:

boot.ci(bootobject, conf=, type= )

type参数设定了获取置信区间的方法。perc方法(分位数)展示的是样本均值,bca将根据偏差对区间做简单调整。

回归的R平方值

1000次自助抽样

输出结果

结果可视化

95%的置信区间获取

回归系数向量函数

自助抽样1000次

获得车重和发动机排量95%的置信区间

置换检验和自助法并不是万能的,它们无法将烂数据转化为好数据。当初始样本对于总体情况的代表性不佳,或者样本量过小而无法准确地反映总体情况,这些方法也是爱莫能助。

参考资料:

你有分类数据然后想要检验是否这些数据值的频数分布是否与预期不符,或者是否组间的频数分布有(显著)差异。

频数检验通常解决两类问题:

通常用于解决这样问题的统计检验方法,分为 精确检验 近似检验 两种。

注意 :精确二项检验仅能用于有两个水平的单变量。Fisher精确检验仅能用于二维列联表(比如,当存在一个独立变量和一个非独立变量时它可以使用;但不能用于两个独立变量和一个非独立变量的情况)。

想要检验配对或被试内效应,我们可以使用McNemar检验。使用该检验必须满足存在两个水平的独立变量和两个水平的非独立变量。

想要检验有重复测量的两个变量独立性,我们可以使用Cochran-Mantel-Haenszel 检验。

假设你有下面的数据,其中每一行代表一个记录:

相比于以 记录 的数据框存储,你的数据可能是 计数 的数据框,或者是一个列联表。本文提到的分析必须使用列联表,你可以参见 this page 获取更多解决方案信息。

想要检验假设:结果列result(忽略条件condition)中的两个值在总体中几乎相等(50%-50%)。

想要检验有不同期望频率的样本(比如下面一个0.75,一个0.25):

如果你想要从检验结果中提取信息,可以将结果保存进一个变量,然后用 str() 函数查看变量信息,接着把你想要的部分取出来。例如:

精确二项检验仅能用于存在两个值的单变量数据。

如果你想要从检验结果中提取信息,可以将结果保存进一个变量,然后用 str() 函数查看变量信息,接着把你想要的部分取出来。例如:

想要检验控制和处理组结果的频数差异,使用2维列联表。

对 2x2 列表,默认使用 Yates’s continuity correction 。这个检验对小样本进行更加保守地估计,设置选项 correct=FALSE 使用无校正的Pearson卡方检验。

对于小样本而言Fisher精确检验更为适合。小样本的2x2列表非常典型,样本更多、更复杂的列表计算强度非常大。当然,用R进行比较复杂的计算也是没有太大问题的。

Cochran-Mantel-Haenszel 检验 (或称为 Mantel-Haenszel 检验))用于检验重复测量两离散变量的独立性。通常使用 2x2xK列表表示,K是测量条件的次数。比如你想要指导是否一个处理(C vs. D)是否影响了恢复的概率(yes or no)。假设该处理一天监控测量三次——早上、中午和晚上,而你想要你的检验能够控制它。那么你可以使用CMH检验对2x2x3列联表进行操作,第三个变量是你想要控制的变量。

R中的CMH检验可以处理比2x2xK维度更高的数据,例如你处理3x3xK列联表。

在接下来的例子里有三个变量:Location,Allele和Habitat。问题是——当控制location变量时,Allel(94或非94)和Habitat(marine或estuarine)两个变量是否独立。

注意上面的数据是 计数 的数据框,而不是像之前的例子是 记录 的数据框。这里我们使用 xtabs() 函数将它转换为列联表。

根据检验结果,当控制Location变量时Allele与Habitat变量存在相关(p=.025)。

注意列联表的前两个维度处理是一致的,所以前后顺序变化都不会影响结果。而最后一个变量变化会导致结果的不同,下面是一个实例。

McNemar检验概念上是频数数据的一个被试内检验。例如,假设你想要检验是否一个处理增加了一个人对某个问题反应“yes”的概率,而且你只有每个人处理前和处理后的数据。标准的卡方检验将不合适,因为它假设了组别是独立的。取而代之,我们可以使用McNemar检验。该检验仅适用于当存在一个独立变量的两次测量时。用于McNemar的列联表与用于卡方检验的非常相似,但结构上是不同的。

假设你有下面的数据。每个对象有处理前和后的反应。

如果你的数据不是宽格式,必须要进行转换(参见 this page 获取更多信息)。

接下来从数据框的pre和post列生成列联表:

执行检验:

对于小样本,它会使用连续校正。我们可以使用 精确 校正的McNemar检验替换这种校正方式,前者更加的精确,可通过 exact2x2 包获取。

原文链接: http://www.cookbook-r.com/Statistical_analysis/Frequency_tests/#cochran-mantel-haenszel-test