permutation test

关于统计-生物学家需要知道的五件事中,说到了非参数统计的重要性,确实生物实验很多数据是messy和non-normal的,比如扫描胶片,扫描仪的设置还有图片处理软件的设置都会影响数据的分布,很多时候根本就搞不清是否是正态,生物实验的样本量通常又比较低,也没办法分析数据的分布.

保险起见,就用非参数统计,正如五件事里所说.

在这种情况下,permutation test也不失为一个好方法.

> set.seed(2012)
> ctr <- rnorm(15, 5, 2)
> trm <- rnorm(15, 6, 2)
> ctr
 [1]  3.4441635  3.8442482  6.3265121  5.1760447  7.5141573  3.7404510
 [7]  4.2080860  5.7459793  5.8142129  2.5497654  5.6246073  7.4974758
[13]  3.3182997  5.7448370 -0.4339599
> trm
 [1] 12.277560  4.365237  1.645831 10.296847  7.784031  6.853054  6.888757
 [8]  7.262824  5.963059  1.868349  6.472083 10.051960  7.878363  5.578179
[15]  7.977236


假设我们有ctr和trm两组数据,这两组数据抽样于正态分布,进行t检验,算出p值。

> t.test(ctr, trm)

	Welch Two Sample t-test

data:  ctr and trm 
t = -2.4147, df = 25.367, p-value = 0.02328
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -4.0810132 -0.3254517 
sample estimates:
mean of x mean of y 
 4.674325  6.877558 

如果用permutation test,随机变换这两组数据的标签,观察它们的均值差,我们将得到均值差值的分布,将实际的观察值对比于随机得到的差值分布,算出随机得到的差值大于实际观察值的概率。

?View Code RSPLUS
1
2
3
4
5
6
7
8
9
10
11
12
13
14
perm.test <- function(d1, d2, nIter=10000) {
	m <- mean(d2-d1)
 	pooledData <- c(d1, d2)
	n <- length(d1)
	meanDiff <- numeric(nIter)
	for (i in 1:nIter) {
		idx <- sample(1:length(pooledData), n, replace=FALSE)
		d1 <- pooledData[idx]
		d2 <- pooledData[-idx]
		meanDiff[i] <- mean(d2) - mean(d1)
	}
	p <- mean(abs(meanDiff) >= abs(m))
	return(p)
}
> perm.test(ctr, trm)
[1] 0.0218

得到的p值和t检验的p值还是很接近的。使用这个方法不需要正态和方差齐性的假设,数据的分布是由我们随机变换之后得到的。关键是它可以做各种统计量的检验,而不单单是均值。

p5rn7vb

Related Posts

Leave a Comment


NOTE - You can use these HTML tags and attributes:
<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>