PH525x series - Inference for Hi

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

• Introduction to high-throughput data
• Inference for high-throughput data
• Multiple testing
• EDA for high-throughput data

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

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
Alternative True
Total

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

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

如果仅检验一个基因的话，所谓的值其实就是，时，的概率，而统计效力（）就是当时，的概率。

Bonferroni校正

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

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

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

• 校正

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

Benjamini-Hochberg校正

• 错误发现率（FDR）

FDR，即错误发现率，其定义如下：
而：
注意，当与时，；当时，。

• 校正

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

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

将公式进行变化后得到的值：

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

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.adjust.methods中有很多可选的矫正方法：

p.adjust.methods

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

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