PH525x series - Permutation Test

假设我们无法应用任何一种标准的数学统计方法时，可以使用置换检验（permutation test）去进行统计推断。

举例说明：

#下载数据
library(downloader) ##use install.packages to install
url <- "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extdata/femaleMiceWeights.csv"
filename <- "femaleMiceWeights.csv" 
download(url, destfile=filename)
dat <- read.csv(filename)

# head(dat)
  Diet Bodyweight
1 chow      21.51
2 chow      28.14
3 chow      24.04
4 chow      23.45
5 chow      23.68
6 chow      19.79

#提取数据
library(dplyr)
control <- filter(dat,Diet=="chow") %>% select(Bodyweight) %>% unlist
treatment <- filter(dat,Diet=="hf") %>% select(Bodyweight) %>% unlist
obsdiff <- mean(treatment)-mean(control)

#obsdiff便是两组小鼠样本体重均值差异

置换检验的核心做法是：将两组样本打乱并随机分为control与case两组，而后推断null distribution（即成立的条件下样本的分布，本例中，为两组间小鼠体重无差异）。在本例中的具体过程如下：

#每组12个样本
N <- 12
#将置换检验重复1000次
avgdiff <- replicate(1000, {
    #将两组打乱
    all <- sample(c(control,treatment))
    #重新将样本分为case与control两组
    newcontrols <- all[1:N]
    newtreatments <- all[(N+1):(2*N)]
    #计算两个新组之间的差值
    return(mean(newtreatments) - mean(newcontrols))
})
hist(avgdiff)
abline(v=obsdiff, col="red", lwd=2)

• 所绘结果如下：

avgdiff.png

avgdiff中大于obsdiff的数值占比便是null distribution的p值，(Why？起初我并未理解，后来查阅了在显著性检验中p值得定义后就明白了，即：它是在零假设下，出现检验统计量的实现值及（向备选假设方向）更极端的值得概率。)，在真正计算p值时还需要在分子、分母上加1（具体原因参考Phipson and Smyth, Permutation P-values should never be zero），如下：

#the proportion of permutations with larger difference
(sum(abs(avgdiff) > abs(obsdiff)) + 1) / (length(avgdiff) + 1)

#结果
 [1] 0.07392607

需要注意的是，并没有理论支持置换检验得来的null distribution是真正的null distribution。这是因为，如果群体间真的存在差异，一些置换检验将是不平衡的，致使最后得到的null distribution会有更大的尾，这也是为什么置换检验得到的p值会比较保守，出于此，当我们样本量很小时，不能做置换检验。

如果数据中拥有隐藏的结构，那么置换检验会将原始数据中的这种结构破坏，使得最后得出的null distribution会低估尾部分布的大小。因此，还需注意的是，进行置换检验需有如下前提：样本互相之间是独立且可交换的。

遗留问题：具体到底啥样的数据算是可交换的？

本文参考