permutation test

Posted by ygc on October 26, 2012 Leave a comment (0) Go to comments

在关于统计-生物学家需要知道的五件事中,说到了非参数统计的重要性,确实生物实验很多数据是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

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)
}