PH525x series - Inference for Hi

2019-12-02 本文已影响0人 3between7

不算练习章节，PH525x 系列教程中共有4章是介绍高维数据的，本篇文章是前三章的笔记，最后一章单独做笔记：

第一章没什么干货，就是介绍了下一般的高通量数据的保存格式，很简单，略。。。所以从第二章开始做笔记喽~~

p值属于随机变量

在显著性检验中，p值不仅是随机变量，而且当零假设为真时， $p$ 值还满足均匀分布，举例说明：

set.seed(1)
population = unlist( read.csv("femaleControlsPopulation.csv") )
N <- 12
B <- 10000
pvals <- replicate(B,{
  control = sample(population,N)
  treatment = sample(population,N)
  t.test(treatment,control)$p.val 
  })
hist(pvals)

hop.png

rowttest函数

因为高通量数据比较大，使用t.test()函数运算速度会比较慢，原文推荐了rowttest()函数，使用方法如下：

library(devtools)
install_github("genomicsclass/GSE5859Subset")

library(GSE5859Subset)
data(GSE5859Subset)

install_bioc("genefilter")
library(genefilter)
g <- sampleInfo$group
results <- rowttests(geneExpression,factor(g))
head(results)

##            statistic           dm    p.value
##1007_s_at -2.11988375 -0.186257601 0.04553344
##1053_at    2.26498355  0.226080913 0.03370683
##117_at     1.54736766  0.096183548 0.13604026
##121_at    -0.54071279 -0.034530859 0.59413846
##1255_g_at -0.03995292 -0.001775578 0.96849102
##1294_at    1.80428848  0.154955035 0.08489586

错误率（Error Rates）

在进行多重比较时，由于要在同一个数据集上检验多个结论，会导致犯第一类错误的概率增加，先来介绍下下面这个表格：

	Called significant	Not called significant	Total
Null True	$V$	$m_0 - V$	$m_0$
Alternative True	$S$	$m_1 - S$	$m_1$
Total	$R$	$m - R$	$m$

为了便于理解，我们以分析差异表达基因的例子来解释上述这个表格中：

$m$ ：检验总次数，也就是基因的个数；
$m_0$ ：零假设为真的基因的个数，也就是组间无差异的基因的个数；
$m_1$ ：零假设为假的基因的个数，也就是组间确实有差异的那些基因的个数；
$R$ ：推论为差异显著的基因总个数，所以 $m - R$ 就是推论为不显著的基因的总个数；
$V$ ：犯 $I$ 类错误或假阳性的基因的个数，也就是本身是无差异的，却得出差异的推论的基因的个数；
$S$ ：真阳性的基因的个数，也就是本身确实有差异，而且检验后也得出差异显著地推论的基因的个数；
$m_1 - R$ ：犯二类错误或者说假阴性的基因个数，也就是本身是有差异的，但是却没有得出差异显著地推论的基因的个数；
$m_0 - V$ ：真阴性的基因个数，也就是本身无差异，也没有得出错误推论的基因个数。

如果仅检验一个基因的话，所谓的 $p$ 值其实就是， $m = m_0 = 1$ 时， $V = 1$ 的概率，而统计效力（ $power$ ）就是当 $m = m_1 = 1$ 时， $S = 1$ 的概率。

Bonferroni校正

族系误差率（Family Wide Error Rate,FWER）

在多重假设检验中，至少犯一次 $I$ 类错误的概率就是FWER，当 $m = 1$ 时，这个概率其实就是 $p$ 值。而事实上，FWER与1非常接近：

$\begin{align} Pr(at\:least\:one\:rejection) & = 1- Pr(no\:rejection) \\ & = 1- ∏_{i=1}^{10000}Pr(p_i > 0.01) \\ & = 1 - 0.99^{10000} ≈ 1 \end{align}$

显然，FWER的概率太高，远远大于了 $0.05$ 。

校正

而Bonferroni校正就是通过控制FWER来降低多重假设检验的错误率，具体办法，就是将显著性阈值设为 $k$ ，其值为 $\alpha/m$ ，这样的话就可以控制多重假设检验整体犯I类错误的概率低于预先设定的显著性水平α（证明过程略）。但是要注意，这个方法有些过于严格保守，用FDR可以替代这种方法。

Benjamini-Hochberg校正

错误发现率（FDR）

FDR，即错误发现率，其定义如下：
$FDR = \bar Q$ 而： $Q ≡ V/R$
注意，当 $R = 0$ 与 $V = 0$ 时， $Q = 0$ ；当 $R = 0$ 时， $V = 0$ 。

校正

BH校正的大致思想其实就是将FDR的值控制在低于阈值 $\alpha$ ，具体过程：

首先，对所有的p值从小到大排序，并记作 $p_{(1)},p_{(2)},\cdots,p_{(m)}$ ，其对应的空假设为 $H_{(1)},H_{(2)},\cdots,H_{(m)}$ 。若想控制FDR不超过 $α$ ，需要找到最大的正整数 $k$ ，使得： $p(k)≤ \frac{k∗α}m$
然后，拒绝 $1≤i≤k$ 时的所有空假设 $H_{(1)},H_{(2)},\cdots,H_{(i)},\cdots,H_{(k)}$ 。这样就能从统计学上保证FDR不超过 $α$ ，保证多重假设检验整体犯I类错误的概率低于预先设定的显著性水平 $α$ （证明过程略）。另外，BH是有效的条件是要求 $m$ 个检验是相互独立的。

将公式进行变化后得到的 $q$ 值：
$q = \frac{p_{(k)}*m}k$

便是校正后的 $p$ 值。举例说明如何求出 $k$ 值：

set.seed(1)
population = unlist( read.csv("femaleControlsPopulation.csv") )

alpha <- 0.05
N <- 12
m <- 10000
p0 <- 0.90 ##10% of diets work, 90% don't
m0 <- m*p0
m1 <- m-m0
nullHypothesis <- c( rep(TRUE,m0), rep(FALSE,m1))
delta <- 3

pvals <- sapply(1:m, function(i){
  control <- sample(population,N)
  treatment <- sample(population,N)
  if(!nullHypothesis[i]) treatment <- treatment + delta
  t.test(treatment,control)$p.value
})

i = seq(along=pvals)
k <- max( which( sort(pvals) < i/m*alpha) )
cutoff <- sort(pvals)[k]
cat("k =",k,"p-value cutoff=",cutoff)

## k = 11 p-value cutoff= 3.763357e-05

p.adjust()函数

在R中，可以使用函数p.adjust()实现校正，其返回值是矫正后的 $p$ 值，它的p.adjust.methods中有很多可选的矫正方法：

p.adjust.methods

## [1] "holm"       "hochberg"   "hommel"     "bonferroni" "BH"        
## [6] "BY"         "fdr"        "none"

注意，这里面的"fdr"就是"BH"。