vlambda博客
学习文章列表

不符合正态性前提的数据，无法应用t检验去比较差异。我们考虑使用非参数的方法作为替代，对于两组数据的比较，可使用wilcoxon检验。类似于t检验，根据样本间是否独立，wilcoxon检验分为wilcox秩和检验以及wilcox符号秩和检验。

假设样本间并非相互独立的，可考虑wilcox符号秩和检验，它是非独立样本t检验的一种非参数替代方法。例如，非独立组设计（dependent groups design）、重复测量设计（repeated measures design）等。尽管此时你选用独立样本的wilcox秩和检验方法也是可行的，这种分析方法本身并没问题（仅仅是在统计算法上存在一些不同，相较之下可能wilcox符号秩和检验会更合适一些）。

此时在wilcox.test()中设定参数“paired = TRUE”即可执行wilcox符号秩和检验。

##wilcox 符号秩和检验，我们执行了一个双侧检验
wilcox_test <- wilcox.test(value~group, shannon_23, paired = TRUE, alternative = 'two.sided')
wilcox_test
wilcox_test$p.value

根据p值（0.039，低于0.05）可知group2和group3的shannon指数间存在显著不同。

考虑作图将两组差异进行可视化展示。例如，一个简单的箱线图示例。

#boxplot() 箱线图
boxplot(value~group, data = shannon_23, col = c('blue', 'orange'), ylab = 'Shannon', xlab = 'Group', main = 'wilcox test: p-value = 0.00295')

相较于参数分析的t检验，wilcox秩和检验不必事先验证数据分布的正态性，因此理论上来讲，只要是两组数据间的差异分析均可使用wilcox秩和检验去完成，因此其适用范围相较于t检验也更广泛。在数据量较大的情况下（可能会存在部分数据满足t检验分析的条件，而另一部分数据则不满足，无法做到全部使用t检验去实现），可以考虑使用循环逐一挑选分组后，直接执行wilcox秩和检验进行各两两分组间的差异分析。尽管这种方法比较“粗暴”，但也不失为一种备选方案。

如下将展示一个批处理示例。

网盘附件中提供了另一示例数据集“gene.txt”。表格中每一行为一种基因，每一列为一个样本，交叉区域为各基因在各样本中的相对丰度。接下来，我们期望通过wilcox秩和检验，找到在group1和group2组中，具有丰度差异的基因。

##wilcox 检验批处理示例
library(doBy) #使用其中的 summaryBy() 以方便按分组计算均值、中位数
 
#读取数据
gene <- read.table('gene.txt', sep = '\t', row.names = 1, header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
group <- read.table('group.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
result <- NULL
 
#使用循环，逐一对各基因进行两组间 wilcox 秩和检验
for (n in 1:nrow(gene)) {
gene_n <- data.frame(t(gene[n,subset(group, group %in% c(1, 2))$sample]))
gene_id <- names(gene_n)[1]
names(gene_n)[1] <- 'gene'

gene_n$sample <- rownames(gene_n)
gene_n <- merge(gene_n, group, by = 'sample', all.x = TRUE)

gene_n$group <- factor(gene_n$group)
p_value <- wilcox.test(gene~group, gene_n)$p.value
if (!is.na(p_value) & p_value < 0.05) {
stat <- summaryBy(gene~group, gene_n, FUN = c(mean, median))
result <- rbind(result, c(gene_id, as.character(stat[1,1]), stat[1,2], stat[1,3], as.character(stat[2,1]), stat[2,2], stat[2,3], p_value))
}
}
 
#输出统计结果
result <- data.frame(result)
names(result) <- c('gene_id', 'group1', 'mean1', 'median1', 'group2', 'mean2', 'median2', 'p_value')
result$p_adjust <- p.adjust(result$p_value, method = 'BH') #推荐加个 p 值校正的过程
write.table(result, 'gene.wilcox.txt', sep = '\t', row.names = FALSE, quote = FALSE)