R语言多重比较示例：Bonferroni校订法和Benjamini & Hochberg法

原文连接：http://tecdat.cn/?p=21825

=

假设检验的基本原理是小几率原理，即咱们认为小几率事件在一次试验中实际上不可能发生。

=

多重比较的问题

当同一研究问题下进行屡次假设检验时，再也不符合小几率原理所说的“一次试验”。若是在该研究问题下只要有检验是阳性的，就对该问题下阳性结论的话，对该问题的检验的犯一类错误的几率就会增大。若是同一问题下进行n次检验，每次的检验水准为α（每次假阳性几率为α），则n次检验至少出现一次假阳性的几率会比α大。假设每次检验独立的条件下该几率可增长至

​
常见的多重比较情景包括：

• 多组间比较
• 多个主要指标
• 临床试验中期中分析
• 亚组分析

=

控制多重比较谬误（Familywise error rate)：Bonferroni矫正

Bonferroni法获得的矫正P值=P×n
Bonferroni法很是简单，它的缺点在于很是保守（大概是各类方法中最保守的了），尤为当n很大时，通过Bonferroni法矫正后总的一类错误可能会远远小于既定α。

控制错误发现率：Benjamini & Hochberg法

简称BH法。首先将各P值从小到大排序，生成顺序数
排第k的矫正P值=P×n/k
另外要保证矫正后的各检验的P值大小顺序不发生变化。

怎么作检验

R内置了一些方法来调整一系列p值，以控制多重比较谬误（Familywise error rate)或控制错误发现率。

Holm、Hochberg、Hommel和Bonferroni方法控制了多重比较谬误（Familywise error rate)。这些方法试图限制错误发现的几率（I型错误，在没有实际效果时错误地拒绝无效假设），所以都是相对较保守的。

方法BH（Benjamini-Hochberg，与R中的FDR相同）和BY（Benjamini & Yekutieli）控制错误发现率，这些方法试图控制错误发现的指望比例。
 
请注意，这些方法只须要调整p值和要比较的p值的数量。这与Tukey或Dunnett等方法不一样，Tukey和Dunnett也须要基础数据的变异性。Tukey和Dunnett被认为是多重比较谬误（Familywise error rate)方法。
 
要了解这些不一样调整的保守程度，请参阅本文下面的两个图。
 
关于使用哪一种p值调整度量没有明确的建议。通常来讲，你应该选择一种你的研究领域熟悉的方法。此外，可能有一些逻辑容许你选择如何平衡犯I型错误和犯II型错误的几率。例如，在一项初步研究中，你可能但愿保留尽量多的显著值，来避免在将来的研究中排除潜在的显著因素。另外一方面，在危及生命而且治疗费用昂贵的医学研究中，得出一种治疗方法优于另外一种治疗方法的结论以前，你应该有很高的把握。

 具备25个p值的多重比较示例

====================

### --------------------------------------------------------------
### 多重比较示例
### --------------------------------------------------------------

Data = read.table(Input,header=TRUE)

按p值排序数据

Data = Data[order(Data$Raw.p),]

检查数据是否按预期的方式排序

​

执行p值调整并添加到数据框

Data$Bonferroni =
      p.adjust(Data$Raw.p,
               method = "bonferroni")

Data$BH =
      p.adjust(Data$Raw.p,
               method = "BH")

Data$Holm =
      p.adjust(Data$ Raw.p,
               method = "holm")

Data$Hochberg =
      p.adjust(Data$ Raw.p,
               method = "hochberg")

Data$Hommel =
      p.adjust(Data$ Raw.p,
               method = "hommel")

Data$BY =
      p.adjust(Data$ Raw.p,
               method = "BY")

Data

​

绘制图表

=====

plot(X, Y,
        xlab="原始的p值",
        ylab="矫正后的P值"
        lty=1,
        lwd=2

 ​

调整后的p值与原始的p值的图为一系列的25个p值。虚线表示一对一的线。

5个p值的多重比较示例

================

### --------------------------------------------------------------
### 多重比较示例，假设示例
### --------------------------------------------------------------
Data = read.table(Input,header=TRUE)

执行p值调整并添加到数据帧

Data$Bonferroni =
      p.adjust(Data$Raw.p,
               method = "bonferroni")

Data$BH =
      signif(p.adjust(Data$Raw.p,
               method = "BH"),
             4)

Data$Holm =
      p.adjust(Data$ Raw.p,
               method = "holm")

Data$Hochberg =
      p.adjust(Data$ Raw.p,
               method = "hochberg")

Data$Hommel =
      p.adjust(Data$ Raw.p,
               method = "hommel")

Data$BY =
      signif(p.adjust(Data$ Raw.p,
               method = "BY"),
             4)

Data

 ​

=

绘制(图表)

=======

plot(X, Y,
        type="l",

​

调整后的p值与原始p值在0到0.1之间的一系列5个p值的绘图。请注意，Holm和Hochberg的值与Hommel相同，所以被Hommel隐藏。虚线表示一对一的线。

---

​

最受欢迎的看法

1.Matlab马尔可夫链蒙特卡罗法（MCMC）估计随机波动率（SV，Stochastic Volatility） 模型

2.基于R语言的疾病制图中自适应核密度估计的阈值选择方法

3.WinBUGS对多元随机波动率模型：贝叶斯估计与模型比较

4.R语言回归中的hosmer-lemeshow拟合优度检验

5.matlab实现MCMC的马尔可夫切换ARMA – GARCH模型估计

6.R语言区间数据回归分析

7.R语言WALD检验 VS 似然比检验

8.python用线性回归预测股票价格

9.R语言如何在生存分析与Cox回归中计算IDI，NRI指标