R语言-组间多重比较

梦里梦见梦不见的497

已于 2023-12-24 21:55:20 修改

阅读量230

点赞数

分类专栏：统计文章标签： r语言

于 2023-12-24 10:54:11 首次发布

原文链接：https://zhuanlan.zhihu.com/p/613817808

版权

统计专栏收录该内容

6 篇文章 0 订阅

订阅专栏

今天需要对ROC曲线下面积AUC进行多重比较，查阅了资料，觉得木舟笔记 - 知乎写的很好，为了以后方便查找，就复制了过来，以下笔记转自木舟笔记 - 知乎的《R实战｜多组间两两比较：从统计方法（多重比较）到结果可视化》

一个引子

Q：多组间的两两比较能不能用「成组T检验」？

A：可以但「必须校正」。

不校正会怎样？

I类错误和II类错误

建立假设是进行假设检验的第一步，通常我们会先建立一个原假设，或者也叫零假设或无效假设（null hypothesis），记为H0，例如某两个（或多个）总体参数相等，或总体参数之差为0。与原假设对立的为备择假设，也称对立假设（alternative hypothesis），记为H1，例如某两个（或多个）总体参数不相等，或总体参数之差不为0。通常备择假设包括「大于或者小于」两种情况，故一般为「双侧检验」。若凭借专业知识有充分把握认为只存在大于或小于两者中的一种可能，则可采用「单侧检验」。

P值，是用来判定假设检验结果的一个参数。如果P值很小，说明原假设H0的发生概率很小，可认为是小概率事件，当P值小到一定程度时，我们就有理由拒绝原假设H0的成立。**但需要注意的是，P值的大小并不能代表所检验的差异的大小，也就是说P值越小，并不能说明差异越大「。这一点很容易引起误解，因此我们在报告结果的时候，提倡使用“」差异有统计学显著性**”的描述，而非“有显著性差异”。

那么，P值一般要小到什么程度才能被认为是小概率事件呢？此时我们就要设立一个「检验水准」，即α，它确定了小概率事件的标准。通常设定α=0.05或0.01，但α的取值并非一成不变，可以根据研究目的的不同给予不同的设置。当P≤α时，在设定α的检验水准下，可认为原假设H0为小概率事件，因此拒绝H0，接受备择假设H1，差异有统计学显著性。

假设检验是基于「抽样样本」来进行结果推断的，而抽样样本只是总体的一小部分，从总体中抽取不同的样本，可能会得出不同的结果，因此我们通常希望抽样样本是一个能够很好地反映总体特征的具有代表性的样本。但由于抽样误差的存在，在进行假设检验根据P值做出推断时具有一定的概率性，因此所得的结论就不一定完全正确，这就是我们常见的假设检验的陷阱：「I类错误和II类错误」。

「I类错误」，也称为「假阳性错误」，就是说实际上总体并无差异，原假设H0是成立的，但是通过假设检验P≤α，在设定α的检验水准下，拒绝了H0，认为有差异，出现了假阳性的现象。前面提到的检验水准α，就是预先设定允许犯I类错误概率的最大值，此时犯I类错误的概率即为α。

「II类错误」，也称为「假阴性错误」，就是说实际上原假设H0不成立，但是通过假设检验P＞α，在设定α的检验水准下，不拒绝H0，得出了阴性的结论，此时犯II类错误的概率为β。

回到原来的多组间的两两检验这个问题来。

对于三组样本组内两两检验，我们共需要比较3次。若每次比较的检验水准 α = 0.05，那么每次拒绝H0犯I型错误为 「0.05」，正确的概率为 「1-0.05 = 0.95」，「3」次t检验拒绝H0「均正确」的概率为0.95x0.95x0.95 = 0.857，3次犯I型错误概率为1-0.857 = 「0.143 > 0.05」。同理，当有4个处理，需要做6个比较，则6个比较中至少有1个被误判的概率为 1−(1−0.05)^6=「0.2649」 。因此，对于多组间的两两比较我们应该使用多重比较。

示例数据和代码

详见：

R实战｜多组间两两比较：从统计方法（多重比较）到结果可视化mp.weixin.qq.com/s/LV2W5aen9QSRcs4wB7MgDA编辑

多重比较

「多重比较（multiple comparisons）」是指「方差分析后」对「各组样本均数之间是否有显著差异」的假设检验的统称，用于方差分析后进一步确定具体哪两个样本均数间有差异。

在重复进行检验时，为了防止检验统计量轻易落入拒绝域，我们需要使用能让检验更加严格的校正方法。校正方法可根据校正对象分为 3 类。

多重比较分类

常用方法示例

1.最小显著差异法(least significance difference,LSD)

rm(list = ls())
library(ggplot2)
library(tidyverse)
library(agricolae)
library(car)
library(reshape2)

# example data
df <- read.csv("crop.data.csv", 
                      header = TRUE, 
                      colClasses = c("factor", "factor", "factor", "numeric"))
head(df)
summary(df)

# 正态性检验
## qq图 点沿着线分布 则满足正态分布
qqPlot(lm(df$yield ~ df$block, data=df), 
       simulate=TRUE, main="Q-Q Plot", 
       lables=FALSE)

# 分组正态性检验
## Shapiro–Wilk正态性检验，p > 0.05 则正态分布
tapply(df$yield,df$block,shapiro.test)

# 方差齐性检验
## Bartlett检验，p > 0.05 则方差齐
bartlett.test(df$yield ~ df$block, data=df)

qqplot

若数据满足「正态分布且方差齐」，则继续「方差分析」。若不满足可对数据进行预处理或使用「非参数检验」。

# 方差分析
variance<-aov(yield ~ block, data=df)
summary(variance)

> summary(variance)
            Df Sum Sq Mean Sq F value  Pr(>F)   
block        3   5.61  1.8693   4.732 0.00409 **
Residuals   92  36.35  0.3951                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

p<0.05 说明组间存在差异。下面进行多重检验。

# 进行多重比较
MC <- LSD.test(variance, 'block', p.adj = 'bonferroni')
print(MC$groups)

> print(MC$groups)
     yield groups
2 177.3169      a
4 177.1760     ab
1 176.8564     ab
3 176.7126      b

「标记字母法」的结果怎么看？

「标记字母法标记法则：」凡是存在「相同字母」的组别之间都「不存在显著差异」，字母「完全不同」的才存在「显著差异」。

因此，只有组2和组3有差异，其他组间均没有差异。

可视化

# 标记字母法数据准备
mark <- data.frame(MC$groups)
mark$group = row.names(mark)
head(mark)
# 绘制
ggplot(df,aes(block, yield, color = block))+
  geom_boxplot()+
  geom_jitter()+
  scale_color_manual(values = c('#eb4b3a', "#48bad0", "#1a9781","#355783"))+
  geom_text(data=mark,
            aes(x=group,y=yield+2,label=groups),
            color="black",
            size = 5,
            fontface="bold")+#添加字母标记
  xlab("")+
  theme_classic()+
  theme(axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title.y = element_text(size = 14,face = 'bold'))
ggsave('test.pdf',width = 5,height = 5)

fig1

2.Bonferroni校正

回到开头提到的问题，多重比较下几种常用的校正P值的方法：Bonferroni、Benjamini & Hochberg（BH或fdr）法等。

在满足正态分布和方差齐后，我们可以使用校正的t检验。

library(rstatix)
# 批量t检验：
stat.test <- df %>%
  t_test(
    yield ~ block,
    p.adjust.method = "bonferroni"
  )
# p.adjust.methods
# c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY","fdr", "none")
stat.test

> stat.test
# A tibble: 6 × 10
  .y.   group1 group2    n1    n2 statistic    df     p p.adj p.adj.signif
* <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <dbl> <chr>       
1 yield 1      2         24    24    -2.51   46.0 0.016 0.095 ns          
2 yield 1      3         24    24     0.817  45.8 0.418 1     ns          
3 yield 1      4         24    24    -1.73   45.9 0.09  0.538 ns          
4 yield 2      3         24    24     3.38   45.6 0.001 0.009 **          
5 yield 2      4         24    24     0.754  46.0 0.455 1     ns          
6 yield 3      4         24    24    -2.59   45.6 0.013 0.078 ns

可视化

stat.test <- stat.test %>% 
  add_xy_position(x='block',dodge = 1)

# 绘制
ggplot(df,aes(block, yield, color = block))+
  geom_boxplot()+
  geom_jitter()+
  scale_color_manual(values = c('#eb4b3a', "#48bad0", "#1a9781","#355783"))+
  stat_pvalue_manual(
    stat.test,  label = "p = {p.adj}", tip.length = 0
  ) +
  xlab("")+
  theme_classic()+
  theme(axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title.y = element_text(size = 14,face = 'bold'))