p值FDR多重校正与可视化

“

小编寄语：这是来自小狐狸同学的投稿2，多才多艺心美人美的女孩。

”

1.假阳性

如果我们只执行一次检验，那么使用p值来解释没有任何问题，但是一旦测试次数增加(如，多重比较)，p<0.05作为显著性标准的适用性就会降低，因为多重比较使假阳性（false positive）的概率会大量增加。

假阳性的简化理解：一个变量进行统计检验，假阳性(第I类错误)的概率为α=5%，不犯假阳性错误的概率为95%；对5个变量一起进行检验，不犯假阳性错误的概率为0.95^5=0.77，至少出现一次第I类错误的概率就是1-0.95^5=0.23。“至少有1个错误”的概率称为总体错误率（Family-Wise Error Rate，FWER），FWER = 1- (1- α)^m，如果样本量很大很大则总体错误率无限接近100%。

由此可见，在进行多次假设检验的时候，为了避免增大犯I型错误，需要对p值进行校正。

2.多重校正P值

假如我们有多组成对的数值型数据（连续型变量）（比如有两组样本的多个变量）需要进行两组间的比较，现在拟通过t检验鉴定两组数据的平均水平是否存在显著差异。为了提高假设检验的准确性，可以借助基于p value排名的多重比较方法来减低假阳性发生的频率。最常见的p value矫正方法为FDR矫正法，获得的p value也称为“FDR adjust (false discovery rate)”或“q value” John Storey（2002）。多重校正p值方法有"holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none"，其中BH(benjaminihochberg)较容易得出阳性p值。一般情况下，我们可以认为FDR = adjusted p value= q value。

3.R实战

借助R语言（Rstudio）这款完全免费的统计工具，我们只需几行简单的代码便可实现上述目标，及数据的可视化。

#加载本次分析所用到的R包
pacman::p_load(stats,dplyr,tidyr,rstatix,ggpubr) 
#依次t检验 Wilcoxon秩和检验所在包/运用管道函数/长宽数据转换函数/统计检验/数据可视化

#---- 情景1：两独立样本的多重组间比较----
#对参与药物实验的动物施以VC或OJ两种药物，剂量分别为0.5，1.0，2。
#一段时间后，检测实验动物的牙齿长度是否存在显著差异，
#以确定药物种类和剂量对动物牙齿生长的影响。
#数据准备
rm(list = ls()) #清空环境中的所有数据
#如果是读取已有的excel数据，则运行以下代码，本文使用R包自带数据集做演示
#df= openxlsx::read.xlsx('data.xlsx')
data("ToothGrowth")
df <- ToothGrowth
View(df) #查看表格，注意V大写

test <- df %>%#test表示生成结果的数据框名称，df表示待分析的数据集
 group_by(dose) %>%#根据不同剂量（分类变量）拆分实验结果，也是分面函数的依据
 t_test(len ~ supp,#变量（牙齿长度，连续变量）~分组（施以的药物种类，分类变量）
        #做t检验，如需秩和检验则改为wilcox_test
        paired = FALSE, #不同给药组动物的数据做比较，选择非配对t检验
        alternative = "two.sided") %>%#输出双尾p value
 adjust_pvalue(method = "fdr") %>% #fdr矫正避免假阳性，也可填BH或其他
 #多重校正p值有"holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none"
 add_significance("p.adj", #根据校正后的p value换算显著性符号
                  cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                  symbols = c("***", "**", "*", "ns")) #根据需要设定显著性符号
                  #的cutoff
View(test) #查看表格，注意V大写

#使用分面函数进行可视化 
ggboxplot(df, #数据框名称
         x = "supp", #横坐标 给药种类
         y = "len", #纵坐标 牙齿长度
         color = "supp", #不同的给药种类（分类变量）使用不同的颜色
         palette = "jco", #在默认的配色板中，选择“jco”方案
         facet.by = "dose",#根据不同剂量（分类变量）拆分实验结果
         ylim = c(0, 40)#纵坐标的下限和上限
) +
 stat_pvalue_manual(test, #使用test数据框标注显著性
                    label = "p.adj", #将adjusted p value标注于标记线上
                    y.position = 35)
                    #由下图可知，两种药物的剂量均为0.5或1时，受试动物的牙齿存在显著差异；而当两种药物的剂量均为2时，受试动物的牙齿长度的均值基本一致。

#---- 情景2：配对/配伍样本多重比较----
#10名受试者参与了特定饮食试验。在试验开始，中间，结束3个时间点分别测定自尊心评分检测受试者在3个时间点的自尊心评分两两间是否存在显著差异，从而判断饮食试验对提升心理素质是否有效。
#数据准备
#rm(list = ls()) #清空环境中的所有数据
#如果是读取已有的excel数据，则运行以下代码，本文使用R包自带数据集做演示
#df= openxlsx::read.xlsx('data.xlsx')
pacman::p_load(stats,dplyr,tidyr,rstatix,ggpubr,datarium) 
data("selfesteem", package = "datarium") #从datarium包提取selfsteem数据
selfesteem
View(selfesteem) #查看表格，注意V大写

# selfesteem是一个“宽数据”（一个变量为一列），我们要把它转换为“长数据”
df <- selfesteem %>% #结果数据框df=初始数据框selfsteem
 gather(key = "time", value = "score", t1, t2, t3) %>% 
 #合并后，测试时间转换后的新列名为“time”，time一列中的结果与原先的列名
 #t1，t2，t3相匹配；t1，t2，t3列数据转换后的新列名为“value”
 convert_as_factor(id, time) 
#对转换后数据框的id和time列进行因子化，即转化为分类变量
View(df) #查看表格，注意V大写

test <- df %>%
 pairwise_t_test(score ~ time, #变量（受试者自尊心评分，连续变量）~分组
                 #（测试时间，分类变量），
                 #如需秩和检验则改为pairwise_wilcox_test
                 paired = TRUE, #由于比较的是不同时间点同样本的数据，
                 #所以需要配对
                 p.adjust.method = "fdr", #fdr矫正避免假阳性，也可填BH或其他
                 alternative = "two.sided") %>%#输出双尾p value
 add_xy_position(x = "time") %>% #自动确定ggplot系列函数显著性标签线位置，
 #如不需要可视化可将此行注释
 add_significance("p.adj", #根据校正后的p value换算显著性符号
                  cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                  symbols = c("***", "**", "*", "ns")) 
 #根据需要设定显著性符号的cutoff
View(test) #查看表格，注意V大写


#做一个简单的抖动散点+箱线图数据可视化
ggboxplot(df, x = "time", y = "score",
         color = "time",
         palette = "jco",
         add = "point") +
 stat_pvalue_manual(test)
#从下图中我们不难看出，受试者的t2时期的评分要显著高于t1，t3-t1，t3-t2的比较同理
#提示这种饮食疗法可能改善了受试者的心理健康状况。

发散思维

1.t检验和Wilcoxon秩和检验

t检验，又称学生t（Student t）检验，是由英国统计学家戈斯特（William Sealy Gosset, 1876-1937）所提出，Student为其笔名。t检验是一种检验总体均值的统计方法，当数据中仅含2组样本，且单组样本数较大时（通常样本个数≥30的样本可视为样本数较大），可用独立或配对样本t检验来检验两组之间的均值是否存在显著差异。

当数据中仅含2组样本但单组样本数较小时（通常样本个数<30的样本可视为样本数较小），或单组数据不服从正态分布时，建议改以Wilcoxon秩和检验（Wilcoxon rank-sum test）来检验总体中位数的差异。

值得注意的是，Wilcoxon检验假定样本X1, … , X n来自分布连续对称的总体X，在此假定下总体X的中位数等于均值。因此检验中位数的性质等同于检验均值。此外，两独立样本Wilcoxon检验和Mann-Whitney U 检验的p value完全一致，区别仅在于统计量的不同。因此两者在只考虑显著性的情况下，是等价的。

2.宽数据（wide data）与长数据（long data）

在数据分析中，也会遇到一个常见的问题，就是长数据和宽数据之间的相互转化。宽数据（wide data）的不同变量会各占一列。而在长数据（long data）中，各个分组的不同变量会合并到一列，不过还有一列专门的说明每一行的数值所对应的变量名称，对应的是宽数据各变量的列名。 

3. 箱线图 

进行均值或中位数的比较，我们常常会使用箱线图对数据进行可视化。箱线图的各部分含义请参下图。

由于本篇教程的侧重于对p value进行矫正，及数据的可视化，关于是使用t检验还是Wilcoxon检验来分析数据，请同学们结合自己的数据情况和上述理论进行判断。后续教程会对两种方法的选择做更加详尽的介绍。

结束语

对于t检验或Wilcoxon检验结果的可视化，大家可以进一步完善ggboxplot函数的代码，以实现与自己数据情况更为匹配的效果。例如标注线两侧tip的高度，也是通过参数设定的。在后续的教程中，我们也会对差异分析数据的可视化做更加详尽地讲解。