Permutation Test

常规步骤：一般多组间样本均数的比较是先做正态性检验和方差齐性，然后在各组都满足正态性和方差齐性的前提下做方差分析。

如果数据呈现正态分布，但方差不齐，则可以选择近似方法(Welch法 和 Brown-Forsythe法)。

如果不满足正态性检验，则选择非参数方法Kruskal-Wallis检验或其他方法。

(来源：陈平雁主编的《IBM SPSS 19统计软件应用教程》)。

Permutation Test。该方法可用于：

1. 完全随机设计资料的分析(两或多样本的均数比较、率比较和等级分布比较)

2. 配对设计资料的分析(配对设计的定量资料差值均值与0比较、配对四格表资料、配对设计等级资料)

3. 随机区组设计资料的分析

4. 相关分析

5. 线性回归分析等

(来源：饶克勤主编的《卫生统计方法与应用进展 第2卷》)。

 Permutation Test的基本思想  

从特定的H0假设出发，根据所研究的问题构造一个统计量，并利用手头样本，按排列组合的原理，导出检验统计量的理论抽样分布；若难以导出确切的理论分布，则采用抽样模拟的方法估计其近似分布。然后求出从该分布中获得手头样本及更极端样本的概率(P值)，并界定此概率，做出推断性结论。

注意：假设现在有18个对象平均分布在3组，则其确切的组合数是18!/(6!*6!*6!) = 17,153,136(约1715万)。如果是成千上万的数据，则计算量之大可见一斑。所以，当计算量大时，一般会采用抽样模拟的方法估计近似分布。

实施步骤(以两样本均数比较为例)：

(1) 建立假设，确定检验水准

H0: μ1=μ2; H1: μ1≠μ2; α=0.05(双侧)

(2) 构造统计量D，并计算现有样本统计量D(obs)

两样本均数的比较，则构造样本均数之差作为统计量

(3) 在H0假设成立的条件下，通过计算机大量计算得到统计量D的“经验抽样分布”。

(4) 计算概率P值

在H0假设成立的前提下，P值为“经验抽样分布”中D值大于等于(或小于等于)现有样本统计量D(obs)的概率。即

P=P(|D|≥|D(obs)|)=#( |D|≥|D(obs)|)/R

 

其中，分母R为确切排列组合数或随机重复的次数，分子为分母中|D|≥|D(obs)|的次数。

(5) 根据小概率原理做出推断性结论

  案例分析  

数据背景：美国各地区(东北部A、南部B、中北部C和西部D)的文盲率，其中每个地区包含多个州，每个州有一个文盲率。所以在比较四个地区(A/B/C/D)的文盲率时文盲率属于数值型变量。

目的：对比四组之间的文盲率是否有差异

方法：单因素方差分析、Kruskal-Wallis检验和Permutation Test

#加载R包

library(tidyverse) 

#获得数据，%>%表示管道操作符，表示参数直接传递给下一个函数，mutate用于增添变量，rename用于修改变量名称，select用于选择变量

states <- as.data.frame(cbind(state.region,state.x77[,3])) %>%
mutate(region = factor(state.region, labels= c("A", "B", "C", "D"))) %>%
rename(illiteracy = V2) %>%
select(region, illiteracy)

#查看数据结构

summary(states)

#绘制箱式图，看组间文盲率的差异

states %>% ggplot( aes (x= region, y = illiteracy,color = region)) + geom_boxplot() + geom_point()

#绘制各组的频数分布直方图，结果可见D组比较偏态，binwidth=0.3表示间距是0.3，facet_wrap是分面函数，指按照地区分别绘制，theme_bw()是主题函数，表示让背景显示白色

states %>% ggplot(aes(x = illiteracy, color= region )) +

geom_histogram(binwidth = 0.3) + facet_wrap(~ region) + theme_bw()

#调用FSA包的Summarize函数，获得每组的描述性指标

library(FSA)

Summarize(illiteracy ~ region, data =states)

#分组正态性检验，发现D组不满足正态性检验，by = list(分组变量)表示按组分别计算

a <- aggregate(states$illiteracy, by = list(states$region),shapiro.test)

#查看结果

a$x

 

#1.进行Kruskal-Wallis检验

kruskal.test(illiteracy ~ region, data =states)

 

#1.1基于Kruskal-Wallis检验的两两比较，bon表示采用bonferroni矫正

kw.pairtest <- dunnTest(illiteracy ~ region,data = states, method="bon")

#提取两两比较的结果

kw.pairtest <- kw.pairtest$res 

#2进行Permutation Test

#此包可进行Permutation Test

library(coin)

#设置种子。由于Permutation Test会模拟抽样，所以每次的结果会存在一定的差异，设置种子有利于结果的重现。
set.seed(1234)

oneway_test(illiteracy ~ region, data = states, distribution= approximate(B=9999))

#2.1基于PermutationTest的两两比较

library(rcompanion)
PT = pairwisePermutationTest(illiteracy ~ region, data =states, method="bon")

#3进行单因素方差分析

#加载R包

library(multcomp) 

#拟合模型

oneway.fit <- aov(illiteracy ~ region, data = states) 

#查看整体模型是否有意义

summary(oneway.fit)

#3.1基于单因素方差分析的两两比较

TukeyHSD(oneway.fit)

表1 三种方法的结果对比

组间	diff	PT.adj.P	KW.adj.P	Anova.adj.P
A-B	-0.74	0.016	0.110	0.002
A-C	0.30	0.046	0.600	0.464
A-D	-0.02	1.000	1.000	0.999
B-C	1.04	<0.001	<0.001	<0.001
B-D	0.71	0.027	0.0038	<0.001
C-D	-0.32	0.518	1.000	0.314

理由：在数据不满足正态分布和结合图1箱式图的情况下，单因素方差分析的结果不太可靠。秩和检验由于舍弃原始数据的具体信息而导致检验功效低。相对而言，Permutation Test 充分利用样本数据的信息，且对原始数据无正态性、方差齐性等要求，故检验效能高。