自助法和置换检验

Posted by Sein on September 11, 2018

自助法(bootstrap)在统计学的发展中出现的较晚,Bradley Efron于1979年发表。而且在很长的时间都受到质疑。一方面自助法在非参数估计方法中表现十分出色,甚至在大多数情况下,精度接近甚至优于参数估计方法。另一方面自助法的原理和实施与其他基于秩的非参方法完全不同。

在当下,数据科学已经将自助法推向前所未有的高度,更激进的观点认为其将取代统计学中基于假设和数理推导的各类估计和检验方法。

在介绍自助法之前,需要简单的说明参数方法,以便对比:

a. 区间估计和假设检验的理论基础都是中心极限定理

b. 中心极限定理假设抽样分布服从正态分布,实际这个假设并不精确。当然这并非是一个漏洞,因为即使如此,参数方法仍然能得到一个不错的结果。

c. 除此之外,假设检验在大多数情况下还需要引入t分布,一方面减弱适用性,另一方面通过样本参数估计抽样分布参数也会降低精度。

d. 参数方法最被诟病的是假设检验结果受到两类错误的限制,阈值在实际使用中很难确定。

现在我们来看自助法的实现:假设样本规模为n ,任意指定某一样本统计量。

1. 从数据中有放回的随机抽取规模为n的样本

2. 计算并记录该样本的统计量

3. 大量重复步骤1~2,例如R

4. 通过获取的R个结果,我们能获得统计量的分布情况

自助法实际上通过对样本的重抽样,来构建统计量的相对频率分布,并且认为此分布是总体对应的统计量抽样分布的近似。而参数方法则是假设统计量抽样分布形状。在实际应用中,认为从样本本身获取的估计要比假设更可靠。并且样本规模越小越明显。

实际上自助法不补偿小规模样本,不创建新数据,也不填补已有数据集中的缺口。只会展示在从原始样本中作抽样时,大量额外样本所具有的行为。

值得一提得是,自助法也可以在多变量数据中发挥作用,直接使用“数据行”作为抽样单元。进而可在自助数据上运行模型,估计模型参数的稳定性(或变异性),或是改进模型的预测能力。

基于自助法估计置信区间简单直接,只需要在自助法的基础上加入下列步骤:

5. 对于x%置信度,对R个样本统计量结果的相对频率分布两端切尾x%的一半

6. 切尾点就是x%自助法置信区间的区间端点

自助法作为非参数方法,和参数方法相比,有不需要对总体分布做出假设的优势。另外自助法的样本统计量不受假设分布的影响,不仅能使用均值、方差、比例,还可以是中位数等任意统计量,甚至可以扩展到某统计量和某一固定值的关系。也就是说基于自助法能够实现单样本的假设检验(同样也适用于配对样本的假设检验)。步骤如下:

1. 提出零假设,确认统计量和某个观测值

2. 根据指定统计量实现自助法流程,获得R个统计量结果

3. 计算比观测值更极端的统计量的数量(注意单尾/双尾的区别),再除以R即可获得p

自助法实现:

在R语言中实现自助法,需要实现两个步骤:首先创建统计量的计算函数,然后用boot库中的boot()函数实现有放回的随机抽样

> library(boot)

> stat_fun <- function(x, ide) median(x[idx])

> boot_obj <- boot(data_name, R = 10000, statistic = stat_fun)

也可以通过sample函数来构造整个过程

> sample(data,100,replace=T)

对照A/B实验中经常出现的两个独立样本的假设检验,自助法发展出了对应的非参数方法,也就是置换检验

就如同参数估计是参数假设检验的基础,置换检验也是基于自助法的非参数方法,在实际的应用中,置换检验涉及两组或多组样本。零假设是针对某一样本统计量,无论是那一组,统计量的值都没有差别。“置换”的方法就是混洗:

1. 有A、B、(C、D......)多个样本,计算统计量的观测值并记录

2. 将所有样本随机混洗,组成一个新的数据集,无放回的随机抽取一个规模和A相等的新样本

3. 从剩下的数据中无放回的随机抽取一个规模与B相等的新样本

4. 如果有C组、D组甚至更多,执行同样的操作

5. 通过新获取的多组样本以相同的方式计算统计量并记录,此过程称作置换迭代

6. 大量重复上述步骤R次,获得R个统计量结果

7. 计算比观测值更极端的统计量的个数(注意单尾/双尾的区别),再除以R即可获得p

除了上述随机混洗的方法之外,还有两种变体。穷尽置换检验和自助置换检验,穷尽置换检验会尝试所有可能的分组,适合规模较小的样本;自助置换检验在步骤2、3中都进行有放回的随机抽样。但是两种变体由于计算太过复杂不被实践关注。

相对于参数假设检验,置换检验最大的优点在于其广泛的适用性。统计量可以任意选取,样本可以是数值,也可以是二元的,两个样本规模可以相差很大,不需要假设总体符合正态分布,也不对两个样本对应的总体方差作任何要求。

置换检验实现:

在R语言中实现置换检验可以直接搭建流程。

> perm_fun <- function(x, n1, n2) {

+ n <- n1 + n2

+ idx_b <- sample(1: n, n1)

+ idx_a <- setdiff(1: n, idx_b)

+ mean_diff <- mean(x[idx_b]) - mean(x[idx_a])

+ return(mean_diff)

+ }

>

> v <- mean_b - mean_a

> perm_diffs <- rep(0, 1000)

> for(i in 1: 1000)

> perm_diffs[i] = perm_fun(data_name, 21, 15)

> hist(perm_diffs)

> abline(v)

> p <- mean(perm_diffs > v)

先创建估计量函数,再用for循环进行R次(1000次),并绘制置换分布的直方图,计算p值。

既然置换检验能处理多组样本,也就很容易实现方差分析,具体步骤如下:

1. 已知有多组数据,计算并记录各组均值间的方差,作为观测值

2. 将所有数据混洗,并随机抽取相同组数相同规模的新样本

3. 计算各组组均值并记录

4. 计算各种均值间的方差并记录

5. 大量重复步骤2~4(例如1000次)

6. 重抽样方差超过观测值方差的比率就是p

除了组间因子,还可以增加其他因子,构成双向方差分析,实现步骤类似

同样的,基于置换检验的方差分析不对各组的分布、方差做出要求。但是各组之间的独立性是必须的

基于置换检验的方差分析实现:

在R语言中lmPerm库提供了aovp函数能够实现此过程,输出结果包括p值,迭代次数,和ANOVA

> library(lmPerm)

> summary(aovp(Time ~ Page, data=four_sessions))

方差分析用于数值型数据,对于分类型数据的多组样本,置换检验也能实现卡方检验。具体步骤如下:

现有一个2×3列联表,第1行记作1,共计n个;第二行记作0,共计m个,n+m=3000。记录n/3为1的预期计数,n为1的观测值的频数

1. 将n个1和m个0合并成一个数据集,即数据集的规模为3000

2. 将数据混洗,随机抽取3组规模为1000的新样本,并计算每组样本中1的个数

3. 计算各组中混洗计数和预期计数的平方差,并求和

4. 大量重复步骤2~3(例如1000次)

5. 重抽样偏差的平方和大于观测值的频数就是p

基于置换检验的卡方检验实现:

使用R语言的chisq.test函数能够计算重抽样的卡方统计量

> chisq.test(data, simulate.p.value=TRUE)

值得注意的是卡方检验不能应用在频数小于5的情况,如果不能合并,就需要使用费舍尔精确检验,使用R语言的fisher.test函数

> fisher.test(data)

虽然方差分析、卡方检验、费舍尔精确检验能够同时检验多个样本,避免多重检验问题,但是只能验证观测差异是否产生于随机性

有一点需要注意的是,无论是自助法还是置换检验,都没有突破α值的问题。它们提供了一种对抽样分布的估计,从而提高了检验效能,但是无法完全确定α的大小,使得在x%的置信度下保证检验的正确,所以对p值的使用仍然需要谨慎对待。更有效的方法是使用验证集,或者是交叉验证

自助法并非没有使用限制,除去对计算能力有较高的要求外,使用自助法的时候样本规模不能太小(当然同样的情况参数方法也无法使用),另一个要求是样本要有足够的代表性。通过随机抽样获得的样本是无偏的,在样本量足够大的情况下,可以认为是满足的。尤其是当前我们面对的数据要远超过统计学对大样本的定义。另外在估计精度和检验效能上,如果样本足够大,参数方法是优于自助法的,但是差别并不影响一般使用