专栏首页小白鱼的生统笔记置换检验概述及其在R中的实现

置换检验概述及其在R中的实现

置换检验及其在R中的实现

某些差异分析方法建立在观测数据服从较好的特定理论分布的基础上,如常见的参数检验T检验、方差分析等,假定观测数据抽样自正态分布;再如分析差异表达基因时最常用的负二项分布模型,假定整体基因表达具有负二项分布特征。在这些情况下,根据数据的理论分布,即可比较容易建立假设检验和总体参数的置信区间估计。

但在许多情况中,统计检验并不一定满足,比如数据抽样于未知或混合分布、样本量过小、存在离群点、基于理论分布设计合适的统计检验过于复杂且数学上难以处理等情况,这时基于随机化的统计方法就可以派上用场。

置换检验概述

置换检验(Permutation Test),也称随机化检验或重随机化检验,提供了稳健的检验框架。

置换检验的基本思想

为理解置换检验的逻辑,可首先考虑如下问题。

在两种处理条件的试验中,10个个体被随机分配到其中一种条件(A或B)中,相应的结果变量也被记录,如下所示。

如果使用T检验分析,可假设数据抽样自正态分布,然后使用双尾T检验来验证结果。此时,零假设为A处理的总体均值与B处理的总体均值相等,根据数据计算t统计量,将其与理论分布进行比较,如果观测到的t统计值落在理论分布值的95%置信区间外,将会拒绝零假设,断言在0.05显著性水平下两组的总体均值不相等。

置换检验的思路与之不同,为了检验两种处理方式的差异,其工作方式如下:

(1)首先计算观测数据的t统计量,称为t0;

(2)随机置换数据,即将10个观测变量取出,再随机分配5个变量给组A,另外5个给组B,并计算置换后数据的t统计量;如此随机置换N次,获得N个置换后数据的t统计量,分别称为tn(n为第N次置换次数);

(3)如此可获得tn的经验分布,如果t0落在tn经验分布中间95%区域的外侧,或者说t0大于(或小于)95%的tn,则在0.05的显著性水平下,拒绝两个处理组的总体均值相等的零假设。

以该图为例,Observed t-ratio就是t0,横轴的t-ratio就代表了tn的频数分布(经验分布),考虑到t0出现在经验分布的右侧,并由于tn代表一种随机化的统计量,因此通常认为若t0>95%的tn,则表明t0是显著的。

并且此时t0<tn的概率就是p值,即

其中nx是置换数据后所获得统计量的零分布期望值高于观测值(上图红色区域)的置换次数,N是置换总次数。该公式还解释了为什么通常在置换检验中,使用的置换次数总设置以“9”结尾(例如99、199、499,而不使用100、200、500):由观测数据所得的初始统计量也被认为是零分布的一部分,因此“+1”。

二者相比,尽管传统的T检验和置换检验过程都计算了相同的t统计量,但置换方法并不是将统计量与理论分布进行比较,而是将其与置换观测数据后获得的经验分布进行比较,根据统计量的极端性判断是否有足够的理由拒绝零假设。如果数据偏离正态分布,存在离群值,或者感觉对于标准的参数方法来说数据集太小,那么置换检验则是一个不错的选择。

关于随机化排列

由于置换检验通过随机置换数据的方式实现评估观测数据结构,因此可知每次置换后数据集的变量分布都是不一致的,即每次所得的随机统计量都有所不同,所以对于同一数据,每次通过置换检验所得的p值会略有差异。再次以该图为例,也就是每次所得的红色区域范围会有差别,但是当置换次数足够大时,p值的差异将很小。

此外,对于执行置换检验的软件,如R,也可以通过设置随机数种子,实现结果重现的目的。

当置换次数达到一定数值时,数据中所有可能的随机排列情况均已被考虑完全,即置换次数达到饱和。此时经验分布依据的是数据所有可能的排列组合,这种情况下的置换检验也常被称为“精确检验”。

随着样本量的增加,获取所有可能排列的时间开销会非常大。这种情况下,可以使用蒙特卡洛模拟,从所有可能的排列中进行抽样,获得一个近似值。

置换检验的广泛应用

不难看出,除了上述演示的代替T检验,置换检验的逻辑还可以延伸至很多统计检验或线性模型等方法上来,为它们提供非参数的检验框架,用以评估结果是否是合理的。

依靠基础的抽样分布理论知识,置换检验提供了另外一种强大的可选检验思路,由于其不受数据分布特征的限制,因此它的应用更为灵活。

尽管如此,置换检验主要发挥作用的地方仍是处理数据偏倚程度较大(如偏离非正态性)、存在离群点、样本很小或无法做参数检验等情况。并且如果初始样本对总体情况代表性很差,即使置换检验也无法提高推断效果。此外,置换检验主要用于生成检验零假设的p值,有助于回答“效应是否存在”,但对于获取置信区间和估计测量精度则比较困难。

置换检验的思想在数十年前就已被提出,但由于置换数据的计人工计算比较费力,限制的它的使用。直到计算机普及后,该方法才得到了真正的应用价值。

在前述已写过的推文中,很多分析中都借鉴了置换检验的思想。

例如在基于距离的差异检验方法中,提到的几种方法均首先计算数据的一种初始统计量,然后再随机置换数据计算置换后数据的统计量获得经验分布,最后比较初始统计量在经验分布中的区间位置,评估检验的显著性。

再例如,非约束排序中被动添加解释变量一文中,将外部解释变量以多元回归的方式被动拟合到非约束排序轴中,可获得拟合的R2,其代表了外部解释变量与非约束轴的关联程度。之后再通过置换检验的方式,获得随机置换后数据的拟合R2,比较观测R2与随机R2的关系,若观测R2明显大于95%以上的随机R2,则代表这种关联是可靠的。与此类似,约束排序中确立约束轴的有效性时,也广泛使用了置换检验。

再例如协惯量分析,统计量RV代表了两个数据矩阵之间的pearson相关性,对于两矩阵结构相似程度的有效性的评估,置换检验就提供了一种便捷的手段。通过比较观测数据集的RV和置换数据后所得的RV,若观测RV明显大于95%以上的随机RV,则代表两组数据集的潜在相关是合理的。

以及更多的方法,不再多说。

R语言中的置换检验

接下来展示置换检验在R语言中的实现方法示例。

使用置换检验替代传统分析方法的R包

R目前有一些非常全面的包可用来做置换检验。

coin包

例如,coin包对于独立性问题提供了一个非常全面的置换检验框架,相对于多种传统的统计检验方法,提供的可选置换检验替代函数。

lmPerm包

lmPerm包则专门用来做方差分析和回归分析的置换检验。它提供了线性模型的置换方法,包括回归(lmp()函数,使用方法类似线性回归函数lm())和方差分析(aovp()函数,使用方法类似方差分析函数aov()),而非正态理论的检验。

其它常见R包

例如perm包,与coin包中的部分功能类似;corrperm包,提供了有重复测量的相关性的置换检验;logregperm包提供了Logistic回归的置换检验;glmperm,涵盖了广义线性模型的置换检验;以及其它一系列的更多R包等。

自写置换检验过程实现的一个示例

其实在理解了置换检验的原理后(不难理解,对不对),也可以根据这种思想,自写代码去解决一些问题,对它进行广泛的应用。

例如这里以计算某数据中,变量间的Pearson相关系数为例。R函数cor()只能计算相关系数,但无法给出相关系数是否显著,我们借助置换检验的思想,首先计算观测数据的Pearson相关系数(cor0),然后对数据集随机置换并计算置换后数据的Pearson相关系数(corN),最后根据|corN|>|cor0|的概率获得p值。

同时也通过psych包中带显著性检验的相关系数计算方法,比较手写置换检验结果和psych包中函数的计算结果是否一致。

#数据获取自 http://blog.sciencenet.cn/blog-3406804-1173120.html
dat <- read.table('data.txt', sep = '\t', row.names = 1, header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
 
#只包含细菌类群丰度数据,计算细菌类群之间的相关性
dat_phylum <- dat[8:17]
var_name <- names(dat_phylum)
 
##自写过程实现
#计算观测数据的相关系数
dat_phylum <- scale(dat_phylum)
cor0 <- cor(dat_phylum, method = 'pearson')
 
p_num <- cor0
p_num[abs(p_num)>0] <- 1
 
#随机置换数据 999 次,计算每次获得的相关系数,并统计 |corN|>|cor0| 的频数
set.seed(123)
 
for (i in 1:999) {
        random <- matrix(sample(dat_phylum), ncol = 10)
        colnames(random) <- var_name
        corN <- cor(random, method = 'pearson')
        
        corN[abs(corN) >= abs(cor0)] <- 1
        corN[abs(corN) < abs(cor0)] <- 0
        p_num <- p_num + corN
}
 
#p 值矩阵,即 |corN|>|cor0| 的概率
p <- p_num/1000
p
 
##使用 R 内置函数的相关性系数计算并获得检验
library(psych)
 
cor_test <- corr.test(dat_phylum, method = 'pearson')
 
cor_test$r
cor_test$p
 
##相关图比较
#左图为手写的置换检验结果,右图为 psych 包获得的结果
#仅显著的相关系数标以背景色
library(corrplot)
 
layout(matrix(c(1,2), 1, 2, byrow = TRUE))
corrplot(cor0, method = 'square', type = 'lower', p.mat = p, sig.level = 0.05, insig = 'blank', addCoef.col = 'black', diag = FALSE, number.cex = 0.8, tl.cex = 0.8)
corrplot(cor_test$r, method = 'square', type = 'lower', p.mat = cor_test$p, sig.level = 0.05, insig = 'blank', addCoef.col = 'black', diag = FALSE, number.cex = 0.8, tl.cex = 0.8)

手写置换检验结果(左图)和psych包中函数的计算结果(右图)是一致的。

参考资料

Robert I. Kabacoff. R语言实战(第二版)(王小宁 刘撷芯 黄俊文 等 译). 人民邮电出版社, 2016.

本文分享自微信公众号 - 小白鱼的生统笔记(gh_5f751e893315),作者:鲤小白

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-01-18

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • RDA、CCA的R2校正及约束轴的显著性检验概述

    RDA排序涵括多元回归分析,关联程度的大小通过回归决定系数R2来反映,R2量化了线性模型中由环境变量所解释的物种组成的变差。RDA模型的解释变差反映了响应变量变...

    用户7585161
  • R语言执行两组间差异分析Wilcox秩和检验

    在进行两组数据间的差异分析时,我们通常会想到使用t检验。但若数据不满足执行t检验的参数假设(例如数据分布不符合正态性,变量在本质上就严重偏倚或呈现有序关系),无...

    用户7585161
  • R包ggraph绘制弧线图(arc diagram)

    在定义上,弧线图(arc diagram)是一种特殊的网络图。和常规的网络图相似,使用节点(node)表示网络中的元素,存在关系的元素之间通过边(edge)连接...

    用户7585161
  • Springsecurity之SecurityContextHolderStrategy

    注:下面分析的版本是spring-security-4.2.x,源码的github地址是: https://github.com/spring-projects...

    克虏伯
  • 如何在小程序中获取用户信息

    在以前的文章中,我们介绍了小程序的登录鉴权功能,方便开发者去获取用户的appid和session_key以便确认用户的身份。但是,仅仅通过appid和sessi...

    it大叔
  • Python 利用socket 实现 s

        主机A---->主机B(80)--->主机C(22), A通过B的80访问主机C131 

    py3study
  • 给大数据入门小伙伴的几个小挑战No.28

    我是小蕉。 子曰:视其所以,观其所由,察其所安,人焉廋哉?人焉廋哉? 子曰:不患无位,患所以立;不患莫己知,求为可知也。 ---- 今天突然神来之笔,有小伙伴...

    大蕉
  • CIO学习:深入了解腾讯大数据平台

    目前腾讯数据平台部的技术团队规模和结构是怎样的? 目前我们数据平台部共有200多人。整个数据平台是按照基础平台、核心应用、产品包装和质量监控的思...

    小莹莹
  • 苹果iphone安装企业根证书(Lync登录)

    苹果手机安装Skype for Business(lync),如果没有购买公网SSL证书,是不允许登录的。解决方法就是手动导入企业内部CA颁发的根证书,如果直接...

    杨强生
  • Pandas进阶修炼120题|第五期

    以上就是Pandas进阶修炼120题第五期全部内容,也是该系列最后一期的内容,如果对本期内容有任何疑问或者更好的方法欢迎给我留言。我会结合所有读者给出的新方法对...

    刘早起

扫码关注云+社区

领取腾讯云代金券