统计学 | bootstrap 和 permutation test 的区别与联系？

如果我有两组数字，想做统计检验其差异。
但是样本太小，不确定原始分布，只能做非参数检验。

0. 准备数据集

a1, a2 之间不显著，t test p=0.11
a1，a2B 之间显著，t test p=0.023

a0=1:20; a0
a1=a0[1:7]; a1
a2=a0[c(1:7) + 2.999999999]; print(a2); t.test(a1, a2) #0.11
a2B=a0[c(1:7) + 3]; print(a2B); t.test(a1, a2B) #0.023

1. bootstrap 法：错误用法

每组抽样，算均值。
抽样n次。
计算两组之间一系列均值的统计学差异…

结果发现，随着抽样次数n的增高，t.test 差异越来越显著。
这种方法是错的！

因为没有拿到稳定的p值。甚至相当于变相增加了样本量，而实际上样本量是固定的。

# 各抽样3次
my.boot = function(a1, a2, n=3){#n=3arr1=c()arr2=c()len=length(a1)/2for(i in 1:n){arr1=c(arr1, mean(sample(a1, len, replace = F)) )arr2=c(arr2, mean(sample(a2, len, replace = F)) )}hist(arr1, n=20, col="#FF000055", main=n)hist(arr2, n=20, col="#00FFFF55", add=T)pval=t.test(arr1, arr2)$p.valuec(n, pval)
}

(1) 不显著组：也显著了

my.boot(a1, a2, 3)
my.boot(a1, a2, 30)
my.boot(a1, a2, 300)
my.boot(a1, a2, 3000)
my.boot(a1, a2, 30000)

在这里插入图片描述

(2) 显著的组

set.seed(42)
my.boot(a1, a2B, 3)
my.boot(a1, a2B, 30)
my.boot(a1, a2B, 300)
my.boot(a1, a2B, 3000)
my.boot(a1, a2B, 30000)

在这里插入图片描述

2. bootstrap 法：对统计学的模糊启发

bootstrap法是用来求mean±ci的。

如果非要求p值，应该每次抽样都计算组件p值，
如果不显著，理论上每次抽样计算的p值是均匀分布，且p=1特别高频。
如果显著，每次抽样计算的p值会向着0倾斜，且p=1频率较低。

my.boot2 = function(a1, a2, n=3){#n=3arrP=c()len=length(a1)/2for(i in 1:n){#arr1=c(arr1, mean(sample(a1, len, replace = F)) )#arr2=c(arr2, mean(sample(a2, len, replace = F)) )arr1=sample(a1, len, replace = F)arr2=sample(a2, len, replace = F)arrP=c(arrP, t.test(arr1, arr2)$p.value)}hist(arrP, n=100, main=n)
}

(1) 不显著组：应该是均匀分布

set.seed(42)
par(mfrow=c(1,4))
my.boot2(a1, a2, 300)
my.boot2(a1, a2, 3000)
my.boot2(a1, a2, 30000)
my.boot2(a1, a2, 60000)

在这里插入图片描述

p=1有一个类似旗杆的bar: 记忆：像是杨戬和孙悟空比试时，孙悟空变化出来的房子和露馅的旗杆。
原理：p<P的概率等于p，则该分布是均匀分布: P(t<T) = P(p<P) = P

我没怎么看懂，谁懂数学，可以证明一下，留言给个链接。。。。

(2) 显著组：p值的分布向着0倾斜

set.seed(42)
par(mfrow=c(1,4))
my.boot2(a1, a2B, 300)
my.boot2(a1, a2B, 3000)
my.boot2(a1, a2B, 30000)
my.boot2(a1, a2B, 60000)

在这里插入图片描述

3. permutation test

所有数据打乱
前一部分作为a组，后一部分作为b组。
求每组的差
重复以上过程n次，做分布
判断真实差值在以上分布中的显著性

随着n的提升，p值慢慢稳定。
这才是统计学该有的样子。

##
my.permut=function(a1, a2, n){dif.real = mean(a1) - mean(a2)dat0=c( a1, a2 )dif=c()for(i in 1:n){# 打乱dat=sample(dat0)a1B=dat[1:length(a1)]a2B=dat[(length(a1)+1):length(dat)]dif=c(dif, mean(a1B) - mean(a2B) )}pval=pnorm(dif.real, mean= mean(dif), sd=sd(dif) )hist(dif, n=100, main=sprintf("n=%s\np=%s", n, pval ))abline(v=dif.real, col="red", lty=2)
}

(1) 不显著组：还是不显著

set.seed(42)
par(mfrow=c(2,3))
my.permut(a1, a2, 3)
my.permut(a1, a2, 30)
my.permut(a1, a2, 300)
my.permut(a1, a2, 3000)
my.permut(a1, a2, 30000)
my.permut(a1, a2, 60000)

在这里插入图片描述

(2) 显著组：依旧显著

set.seed(42)
par(mfrow=c(2,3))
my.permut(a1, a2B, 3)
my.permut(a1, a2B, 30)
my.permut(a1, a2B, 300)
my.permut(a1, a2B, 3000)
my.permut(a1, a2B, 30000)
my.permut(a1, a2B, 60000)

在这里插入图片描述