R语言执行两组间差异分析T检验

vlambda
2020-05-17

R语言执行两组间差异分析T检验

两组间差异分析之T检验在R中实现

做统计学分析的其中一个重要目的就是寻找组间差异，在研究中，我们最常关注的问题莫过于处理组与对照组是否存在了显著不同。就两组间比较而言，t检验是常见的分析方法之一。

本文简介怎样在R中进行t检验，以实现两组间差异分析。

本文使用的作图数据的网盘链接（提取码 6268）：

https://pan.baidu.com/s/108oBkDTkuJ1HGBPX0h5Osg

文件“alpha.txt”为某16S细菌群落测序所获得的部分alpha多样性指数数据，包含3列信息：sample，样本名称；observed_species和shannon分别为两种类型的alpha多样性指数。文件“group.txt”为各样本分组信息，第一列（sample）为各样本名称；第二列（group）为各样本的分组信息。

接下来，我们期望在R中运行t检验，以查看不同分组间（两两分组之间）的各alpha多样性指数是否存在显著不同。

数据预处理及正态性假设检验

数据预处理

首先将上述两个数据表读入 R 中，并合并在一起。

library(reshape2)
 
#读入文件，合并分组信息，数据重排
alpha <- read.table('alpha.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
group  <- read.table('group.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
alpha <- melt(merge(alpha, group, by = 'sample'), id = c('sample', 'group'))
 
#我们期望查看 group1 和 group2 的 observed_species 指数是否存在显著差异
#选择要比较的分组
richness_12 <- subset(alpha, variable == 'observed_species' & group %in% c('1', '2'))
richness_12$group <- factor(richness_12$group)
head(richness_12, 10)

选取的数据框“richness_12”内容如下所示。第一列（sample），两组数据中所含样本名称；第二列（group），两组分组名称，且分组列已转化为因子类型；第三列（variable），alpha多样性指数observed_species；第四列（value），observed_species指数的数值。

正态性假设检验

t检验的一个重要前提就是数据必须符合正态分布模型。因此在执行t检验之前必须验证数据分布的正态性。若数据不符合正态分布，则t检验将无法适用于该数据（此时可以考虑转化数据，或者使用非参数的检验方法）。验证数据是否符合正态分布的方法很多，以下展示两种常见方法。

正态QQ图

以下使用car包中的qqPlot()，绘制QQ图查看数值分布。结果中，横坐标是标准的正态分布值，纵坐标是我们数据的值。如果两者基本相等，或者说所有的点都离直线很近，落在置信区间内（图中虚线部分，默认展示95%置信区间），即表明正态性假设符合得很好。

##正态 qq 图验证数据正态性
library(car)
 
#QQ-plot
qqPlot(lm(value~group, data = richness_12), simulate = TRUE, main = 'QQ Plot', labels = FALSE)

R语言执行两组间差异分析T检验

Shapiro-Wilk检验

Shapiro-Wilk类似于线性回归的方法，是检验其于回归曲线的残差，据此验证数据分布的正态性。R中提供了可用于执行Shapiro-Wilk检验的函数shapiro.test()，原假设（或称零假设）为数据集符合正态分布，若结果中p值大于0.05，则接受原假设，数据分布符合正态性。

##Shapiro-Wilk 检验，当且仅当两者 p 值均大于 0.05 时表明数据符合正态分布
shapiro <- tapply(richness_12$value, richness_12$group, shapiro.test)
shapiro
shapiro$'1'$p.value
shapiro$'2'$p.value

R语言执行两组间差异分析T检验

T检验

综上，我们的数据分布通过了正态假设检验，即可执行t检验。可分为独立样本的t检验与非独立样本的t检验。

独立样本的t检验

如果样本间是相互独立的，可选用独立样本t检验。R语言t检验函数t.test()中默认两组间相互独立（默认参数paired = FALSE），执行独立样本的t检验。同时需注意，在R中的t检验默认假定方差不相等（默认参数var.equal = FALSE），并使用Welsh的修正自由度；可以通过参数“var.equal= TRUE”假定方差相等，并使用合并方差估计。此外，t.test()默认的备择假设是双侧的（默认参数alternative = 'two.sided'），即执行双侧检验；可分别使用“alternative= 'less'”或“alternative = 'greater'”执行单侧t检验。

我们执行了一个假设方差不相等的双侧检验，如下示例。

##独立样本的 t 检验
t_test <- t.test(value~group, richness_12, paired = FALSE, alternative = 'two.sided')
t_test
t_test$p.value

由于p值远小于0.05，即拒绝了原假设（原假设两组间没有差异），group1和group2的observed_species指数间存在显著不同。

非独立样本的t检验

如果样本间并非相互独立的，可选用非独立样本t检验。例如，非独立组设计（dependent groups design），前-后测设计（per-post design），或重复测量设计（repeated measures design）等。尽管此时你选用独立t检验的方法也可以。

此时在t.test()中设定参数“paired = TRUE”即可执行非独立样本的t检验，如下示例，同样为假设方差不相等的双侧检验。

##非独立样本的 t 检验
t_test <- t.test(value~group, richness_12, paired = TRUE, alternative = 'two.sided')
t_test
t_test$p.value

由于p值远小于0.05，即拒绝了原假设（原假设两组间没有差异），group1和group2的observed_species指数间存在显著不同。

可视化展示

考虑作图将两组差异进行可视化展示。

#boxplot() 箱线图示例
boxplot(value~group, data = richness_12, col = c('blue', 'orange'), ylab = 'Observed_species', xlab = 'Group', main = 't-test: p-value < 0.001')

R语言执行两组间差异分析T检验

#ggplot2 柱形图示例
#分别计算各组中的均值以及标准差，展示为均值 ± 标准差的柱形图样式
 
library(doBy)  #使用其中的 summaryBy() 以方便按分组计算
library(ggplot2)      #ggplot2 作图
 
dat <- summaryBy(value~group, richness_12, FUN = c(mean, sd))
 
ggplot(dat, aes(group, value.mean, fill = group)) +
geom_col(width = 0.4, show.legend = FALSE) +
geom_errorbar(aes(ymin = value.mean - value.sd, ymax = value.mean + value.sd), width = 0.15, size = 0.5) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), plot.title = element_text(hjust = 0.5)) +
labs(x = 'Group', y = 'Observed_species', title = 't-test: p-value < 0.001')

R语言执行两组间差异分析T检验