R语言实现bootstrap和jackknife检验方法

2023-02-23 11:08:02Python016

R语言实现bootstrap和jackknife检验方法,第1张

写在最前面：

首先需要说一下，本文的bootstrap和jackknife都算是蒙特卡罗方法(Monte Carlo method)的一种。应用广泛的的MCMC链(马尔可夫链蒙特卡洛方法Markov chain Monte Carlo)也是蒙特卡罗与马尔可夫链的结合。简单来说，蒙特卡罗方法就是从已知样本的分布中随机抽取新的样本集进行评估，然后放回，再次抽取的方法。根据具体方法的不同，抽取样本集的手段也不同。

bootstrap抽样方法将观测到的样本视为一个有限的总体，是唯一的信息来源，从中有放回的随机抽样来评估总体特征，以及对抽样总体进行推断统计。bootstrap 也分参数bootstrap和非参数bootstrap，前者的分布已完全知道。但在生信领域一般没有这种情况。所以下面讨论的是非参数bootstrap。

直接上例子：

假设现在有bootstrap包中的law数据集如下，

现在我们要计算LSAT成绩（美国法学入学考试）和GPA之间的相关系数。但因为样本量太少了，所以我们使用bootstrap重复抽样评估其标准误。

200次循环抽样后，计算得se.R标准误为0.1474629

得到如下的图：

1e6次循环抽样后，计算得se.R标准误为0.1333802

得到如下的图：

如果用bootstrap包的bootstrap函数会快一些：

bootstrap函数的用法： bootstrap(抽取样本范围,重复次数,进行bootstrap的函数,bootstrap的数据集)

偏差定义为bootstrap结果（多个数值）与原数据统计结果（单个数值）的均值：

得到bias大约为0.001817608，比较小

换一个包，boot包

这里用了三种方法计算置信区间：basic、正态和百分数。样本相关系数分布接近正态，则正态置信区间接近百分数区间。此外还有“Better Bootstrap Confivendence Interval” 更好的bootstrap置信区间，称为BCa区间，使用偏差和偏度对百分数置信区间进行矫正。设置type="bca"即可。

简单的说，bootstrap是从原有真实样本中有放回地抽取n个。jacknife就是每次都抽取n-1个样本，也就是每次只剔除一个原样本。

同样地，如果以bootstrap包中的law数据进行演示：

Jackknife计算的bias为-0.006473623。 这里jackknife的偏差公式相比于bootstrap有一个(n-1)系数，推导就不写了。

标准误se为0.1425186，与bootstrap得出的比较接近。

当统计量不太平滑的时候，Jacknife有很大误差。比如说对中位数进行统计，其变化很大。在进行Jacknife之后最好再跑一次bootstrap，看看是否相差很大。

居然还能这么嵌套着玩，针对每次bootstrap形成的数列向量计算jackknife的标准差，这样可以看出bootstrap若干次取样之间的差异。

算出来分别为0.1344824和0.08545141。后者较小，表面bootstrap取样之间的variance较小。

简单来说就是一种数据分割检验的方法，将数据分割为K份，称为"K-fold"交叉检验，每次第i个子集作为测试集来评估模型，其余的用来构建模型。Admixture使用的就是这个原理。Jackknife也属于Cross Validation的应用之一。

现在我创建一个这样的alignment：

这棵树长这样，符合遗传距离：

进行bootstrap：

phylogeny的bootstrap是对每一个节点都进行bootstrap取样并建树，比如说在9号节点，查看其bootstrap子集建的树符合系统发育关系((human2,human4,human3)(human8,human1,human6,human7,human5))的百分比（不管内部怎么样，先看这个节点）。发现Node1支持率是100（1000次都符合）。而后移到下一个节点，并且只看节点内部的分支支持率是多少。

其实原理都比较简单，计算bootstrap也会有专门的软件。

参考资料：

1）中科大张伟平教授课件

2） https://ecomorph.wordpress.com/2014/10/09/phylogenetic-trees-in-r-4/

数据准备

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

置换检验的定义

置换检验(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包扩展了自助法和重抽样的相关用途。你可以对一个统计量(如中位数)或一个统计量向量(如一列回归系数)使用自助法。

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