DALS017-高维数据推断2-统计原理

title: DALS017-高维数据推断2-统计原理
date: 2019-08-17 12:0:00
type: "tags"
tags:

高维数据
categories:
Data Analysis for the life sciences

前言

这一部分是《Data Analysis for the life sciences》的第6章线性模型的第2小节，这一部分的主要内容涉及高维数据统计的一些原理，相应的R markdown文档可以参考作者的Github。

高维数据的p值

在前面我们了解了，当我们在分析高维数据时，p值就不再是一个很好的统计指标了。这是因为，我们同一次检测了许多特征值(feature)。这种检测手段被称为多重比较( multiple comparison)，或多重检测(multiple testing)，或多重性问题(multiplicity)。在此情况下，p值就不再适用。另外，当我们同时检测多个假设问题时，仅仅基于一个小p值的阈值，例如0.01，这就很容易导假阳性。针对这种情况，我们需要定义一个新的术语来研究高通量数据。

为了处理多重比较的问题，我们广泛使用的方法就是定义一个程序(procedure，也可以说是一种算法，也可以翻译为校正等等，总之，表达的是一个意思)，然后用它来估计或控制(control)计算过程中的错误率(rate error)。我们这里所说的控制(control)的意思是说，我们会采用这个程序来保证错误率(error rate)低于某个提前设定的值。通过参数或截止值(cutoff)来进行设定的这个程序通常比较灵活，它会让我们能够控制特异性(specificity)和灵敏度(sensitivity)，这种程序的一个典型功能如下所示：

计算每个基因的p值；
计算出p值小于 $\alpha$ 的所有显著性基因。

这里需要注意的是，当我们改变 $\alpha$ 值时，会调整相应的特异性(specificity)和灵敏度(sensitivity)。

接着我们来定义错误率(error rate)，它会让我们对统计过程进行估计和控制。

错误率

在这一部分中，我们会了解到I类错误与II类错误，这两类错误分别代表假阳性(false positives)与假阴性(false negatives)。通常，特异性(specificity)与I类错误有关，灵敏性(sensitivity)与II类错误有关。

在一次高通量实验里，我们会犯第I类错误和第II类错误。我们参考了Benjamini-Hochberg的论文，做了以下表格，总结了这些错误，如下所示：

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

为了说明这个表格中的内容，我们来打个比方，假设我们检测了10000个基因，这就意味着我们需要做检测的次数为 $m=10000$ 。

那些零假设是真的基因数目（多数是我们不感兴趣的基因，零假设为真，也就是说这些基因在对照和实验组中都没有显著性差异）为 $m_{0}$ ，那些零假设为假的基因数目为 $m_{1}$ ，零假设为假，也就是说，替代假设(alternative hypothesis)为真（不一定严谨，反正就是说零假设为假）。通常来说，我们感兴趣的是尽可能地检测到那些替代假设为真的基因（真阳性），避免检测到那些零假设为真的基因（假阳性）。对于多数高通量实验来说，我们会假设 $m_{0}$ 远大于 $m_{1}$ （这句话我的理解就是，在一次高通量实验中，没差异基因的数目 $m_{0}$ 要大于有差异基因的数目 $m_{1}$ ）。

例如我们检测了10000个基因，对其中约有100个基因感兴趣。这也就是说， $m_1 \leq 100$ 并且 $m_0 \geq 19,900$ 。

在这一章中，我们指的特征值(feature)就是我们的检测值。在遗传学中，这些特征值可以是基因(genes)，转录本(transcripts)，结合位点(binding sites)，CpG岛和SNPs。

在上面的那个表格中， $R$ 表示的是经过检测后，有显著性差异的特征值的数目总和，而 $m-R$ 则表示不显著的基因数目。表格中剩下的部分表示的是一些实际上未知的重要的量。

$V$ 表示I类错误或假阳性。更具体地来说就是， $V$ 表示了那些零假设为真的特征值的数目。
$S$ 表示的是真阳性的数目。具体地来说就是， $S$ 表示替代假设为真的特征值的数目。

$m_{1}-S$ 表示了II类错误或假阴性， $m_{0}-V$ 表示真阴性。

这里需要牢记的是，当我们只做一次检测时，p就仅仅是当 $V=1$ ， $m=m_{0}=1$ 时的概率。功效(power)就是当 $S=1$ ， $m=m_{1}=1$ 时的概率。在这种非常简单的案例里，我们并不制作上面类似的表格，但是我们会说明一下，如何定义表格中的术语，从而帮助我们处理高维数据。

数据案例

现在看一个案例。在这个案例中我们会使用小鼠数据进行Monte Carlo模拟，从而模拟一种情况，在这种情况里，我们会检测10000种对小鼠体重无影响的减肥饲料(fad diets)。这就是说，在零假设下，这些饲料对小鼠体重没影响为真，也就是说 $m=m_{0}=10000$ ，并且 $m_{1}=0$ ，现在我们先进行一个样本数目为12的计算，并且我们认为p值小于 $\alpha=0.05$ 显著，如下所示：

set.seed(1)
population = unlist( read.csv("femaleControlsPopulation.csv") )
alpha <- 0.05
N <- 12
m <- 10000
pvals <- replicate(m,{
  control = sample(population,N)
  treatment = sample(population,N)
  t.test(treatment,control)$p.value
})
sum(pvals < 0.05) ##This is R

结果如下所示：

> sum(pvals < 0.05) ##This is R
[1] 462

从结果我们可以看出，这个假阳性（462个）还是比较高的，这是要在多数分析中是要避免的。

下面我们来看一个更加复杂的代码，这段代码会进行人为设定10%的饮食有效，平均效应(average effect size)为 $\Delta=3$ 盎司（约85克）。仔细研究这段代码可以帮助我们理解上面的那个表格，现在我们先来定义这样的数据，其中10%有效：

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

现在我们做10000次模拟统计，每次都采用t检验，并记录下来：

set.seed(1)
calls <- sapply(1:m, function(i){
  control <- sample(population,N)
  treatment <- sample(population,N)
  if(!nullHypothesis[i]) treatment <- treatment + delta
  ifelse( t.test(treatment,control)$p.value < alpha,
          "Called Significant",
          "Not Called Significant")
})

由于我们知道哪些是数据是有差异的（毕竟是自己人为生成的，保存在了nullHypothesis中)，我们现在计算一下上面的表格，代码如下所示：

null_hypothesis <- factor( nullHypothesis, levels=c("TRUE","FALSE"))
table(null_hypothesis,calls)

结果如下所示：

> null_hypothesis <- factor( nullHypothesis, levels=c("TRUE","FALSE"))
> table(null_hypothesis,calls)
               calls
null_hypothesis Called Significant Not Called Significant
          TRUE                 421                   8579
          FALSE                520                    480

结果中的第1列就是 $V$ 与 $S$ ，需要注意的是， $V$ 与 $S$ 是随机变量，如果我们再次运行这段代码以，这些值就会改变，如下所示：

B <- 10 ##number of simulations
VandS <- replicate(B,{
  calls <- sapply(1:m, function(i){
    control <- sample(population,N)
    treatment <- sample(population,N)
    if(!nullHypothesis[i]) treatment <- treatment + delta
    t.test(treatment,control)$p.val < alpha
  })
  cat("V =",sum(nullHypothesis & calls), "S =",sum(!nullHypothesis & calls),"\n")
  c(sum(nullHypothesis & calls),sum(!nullHypothesis & calls))
})

结果如下所示：

V = 410 S = 564 
V = 400 S = 552 
V = 366 S = 546 
V = 382 S = 553 
V = 372 S = 505 
V = 382 S = 530 
V = 381 S = 539 
V = 396 S = 554 
V = 380 S = 550 
V = 405 S = 569

针对这种情况，我们就定义了错误率(error rate)。例如，我们可以估计 $V$ 大于0的概率。它可以解释为，在10000次检验中，我们出现I类错误的概率。在上述模拟数据中， $V$ 在每次模拟中都大于1，因此我们怀疑这个概率实际上就是1（这里的1就是“V大于0”这个事件的概率，换句话讲，就是，在这个检验中，必然会出现假阳性）。

当 $m=1$ 时，这个概率就等于p值。当我们进行多重检验模拟时，我们称之为多重比较谬误(Family Wide Error Rate,FWER)，它与我们广泛使用的一个检验方法有关，即Bonferroni校正( Bonferroni Correction)。

Bonferroni校正

现在我们了解一下FWER是如何运用的，在实际的分析中，我们会选择一个程序(procedure)来保证FWER小于某个提前设定好的阈值，例如0.05，通常情况下，就使用 $\alpha$ 来表示。

现在我们来描述一个这样的程序：拒绝所有p值小于0.01的假设。

为了说明这个目的，我们假设所有的检验都是独立的（在10000个饲料检验的实验里，这个假设是成立的；但是在一些遗传学实验里，这个假设有可能不成立，因为某些基因会共同发挥作用）。每次检验得到的p值我们用 $p_1,\dots,p_{10000}$ 表示。这些独立的随机变量如下所示：
$\begin{align*} \mbox{Pr}(\mbox{at least one rejection}) &= 1 -\mbox{Pr}(\mbox{no rejections}) \\ &= 1 - \prod_{i=1}^{10000} \mbox{Pr}(p_i>0.01) \\ &= 1-0.99^{10000} \approx 1 \end{align*}$
如果我们要模拟这个过程，代码如下所示：

B<-10000
minpval <- replicate(B, min(runif(10000,0,1))<0.01)
mean(minpval>=1)

结果如下所示：

> mean(minpval>=1)
[1] 1

因此，我们计算出来的FWER是1，这并非我们想要的结果。如果我们想要它低于 $\alpha=0.05$ ，那么我们的统计就是失败的。

我们如何做才能使得一个错误出现的概率低于 $\alpha$ 呢，我们可以使用上面的公式推导一下，通过选择更严格的截止值(cutoff），之前的是0.01，从而将我们的错误降低至至少5%的水平，如下所示：
$\mbox{Pr}(\mbox{at least one rejection}) = 1-(1-k)^{10000}$
解这个方程，我们就得到了 $1-(1-k)^{10000}=0.01 \implies k = 1-0.99^{1/10000} \approx 1e-6$

这就是我们给出的一个程序案例。这实际上就是Sikdak过程。尤其是，当我们定义一个说明，例如拒绝所有p值小于 0.000005的零假设。然后，我们知道了p值是一个随机变量，我们会使用统计理论来计算，如果我们遵循这个程序，平均会犯多少错误。更确切地讲就是，我们可以计算出这些错误的边界，也就是说，这些错误小于某些预定值的比例。

Sidak校正的一个问题是，这个校正假设所有的检验都是独立的。因此当这个假设时成立的，它只能控制FWER。百Bonferroini校正则更为通用，因为即使每个检测不独立，它也能控制FWER。。与Sidak校正一样，我们首先来看一下：
$FWER = \mbox{Pr}(V>0) \leq \mbox{Pr}(V>0 \mid \mbox{all nulls are true})$

现在使用前面表格的那种表示方法为：
$\mbox{Pr}(V>0) \leq \mbox{Pr}(V>0 \mid m_1=0)$
Bonferoni校正会设定 $k=\alpha/m$ ，因此可以写为如下形式：
$\begin{align*} \mbox{Pr}(V>0 \,\mid \, m_1=0) &= \mbox{Pr}\left( \min_i \{p_i\} \leq \frac{\alpha}{m} \mid m_1=0 \right)\\ &\leq \sum_{i=1}^m \mbox{Pr}\left(p_i \leq \frac{\alpha}{m} \right)\\ &= m \frac{\alpha}{m}=\alpha \end{align*}$
将FWER控制在0.05水平是一种非常保守的方法，现在我们使用前面计算的p值来计算一下：

set.seed(1)
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
})

只需要关于p值是否小于0.05/10000即可，如下所示：

sum(pvals < 0.05/10000)

结果为：

> sum(pvals < 0.05/10000)
[1] 2

当使用了Bonferroni校正后，即使我们进行了10000次饮食实验，却只有发现2个假阳性的结果，Bonferroni是一种非常严格的校正。

错误发现率(FDR)

在许多情况下，要求FWER是0.05没多大意义，因为它太严格了。例如，我们常见的做法就是先进行初步的小型实验来确定小数候选基因。这种做法被称为之发现驱动的实验或项目。我们也许会寻找一个未知的致病基因，而不仅仅是对候选基因进行采用更多的样本进行后续研究。如果我们开发一个程序，例如一个10个基因列表，从中发现1到2个为重要的基因，这个实验就算非常成功。小样本实验，实现 $FWER\leq0.05$ 的唯一途径就是使用一个空的基因列表。在前面我们已经看到，虽然只有1000种包含有效，但是最终我们只得到2个饮食有效这样一个结果，如果把样本数目降低到6，这个结果有可能就是0，如下所示：

set.seed(1)
pvals <- sapply(1:m, function(i){
  control <- sample(population,6)
  treatment <- sample(population,6)
  if(!nullHypothesis[i]) treatment <- treatment + delta
  t.test(treatment,control)$p.value
})
sum(pvals < 0.05/10000)

结果如下所示：

> sum(pvals < 0.05/10000)
[1] 0

由于我们要求 $FWER\leq 0.05$ ，因此我们实际上就是0功效（灵敏度）。在许多方法中，这种情况称为特异性(specificity)过强(over-kill)。替代FWER的另外一种方法就是FDR，即错误发现率(false discover rate)。FDR背后的思想就是集中关注Q值，即 $Q \equiv V/R$ ，当 $R=0,V=0$ 时， $Q=0$ 。其中，当 $R=0$ （没有显著性）就表明， $V=0$ （没有假阳性）。因此 $Q$ 是一个随机变量，它的范围是0到1，通过计算Q的平均值，我们可以定义一个比值(rate)，它所表示的是意思是，在显著的基因里，假阳性的基因占的比例。为了更好地理解这个概含，我们计算Q的程序要求调用所有的p值都小于0.05。

向量化代码

在R中可以使用sapply()系列函数来加快运行速度，前面已经看到了这个函数的使用方法，现在使用传统的代码来看一下Q值的计算方法，如下所示：

library(genefilter) ##rowttests is here
set.seed(1)
##Define groups to be used with rowttests
g <- factor( c(rep(0,N),rep(1,N)) )
B <- 1000 ##number of simulations
Qs <- replicate(B,{
  ##matrix with control data (rows are tests, columns are mice)
  controls <- matrix(sample(population, N*m, replace=TRUE),nrow=m)
  ##matrix with control data (rows are tests, columns are mice)
  treatments <- matrix(sample(population, N*m, replace=TRUE),nrow=m)
  ##add effect to 10% of them
  treatments[which(!nullHypothesis),]<-treatments[which(!nullHypothesis),]+delta
  ##combine to form one matrix
  dat <- cbind(controls,treatments)
  calls <- rowttests(dat,g)$p.value < alpha
  R=sum(calls)
  Q=ifelse(R>0,sum(nullHypothesis & calls)/R,0)
  return(Q)
})

结果如下所示：

> head(Qs)
[1] 0.4513274 0.4063786 0.4568627 0.4490414 0.4468314 0.4315569
> tail(Qs)
[1] 0.4390756 0.4718657 0.4395373 0.4425711 0.4536391 0.4284284

控制FDR

上述代码使用Monte Carlo模拟计算了1000次10000个样本的实验，并保存了Q值，现在看一下Q值的直方图，如下所示：

library(rafalib)
mypar(1,1)
hist(Qs) ##Q is a random variable, this is its distribution

image

FDR就是Q值的平均值，如下所示：

FDR=mean(Qs)
print(FDR)

如下所示：

> print(FDR)
[1] 0.4463354

这里的计算结果表明，FDR值比较高。这是因为对于90%的统计检验而言，零假设是真。这就暗示了，由于p值是cutoff值为0.05，100个检验中，大概有4到5个是假阳性。再加上我们没有考虑到替代所设为真时的所有情况，因此FDR值就比较高。那么我们如何控制它呢，如果我们想要更低的FDR值，比如5%怎么办？

为了用图形说明FDR为什么这么高，我们来绘制p值的直方图。我们使用一个较大的m值从直方图中获取更多的数据。再绘制一条水平线，表示当NULL为真时，从 $m_{0}$ 个案例（指的基因或其它检测值， $m_{0}$ 就是没有统计学差异的基因）中到的阳性结果的均匀分布，如下所示：

set.seed(1)
controls <- matrix(sample(population, N*m, replace=TRUE),nrow=m)
treatments <- matrix(sample(population, N*m, replace=TRUE),nrow=m)
treatments[which(!nullHypothesis),]<-treatments[which(!nullHypothesis),]+delta
dat <- cbind(controls,treatments)
pvals <- rowttests(dat,g)$p.value
h <- hist(pvals,breaks=seq(0,1,0.05))
polygon(c(0,0.05,0.05,0),c(0,0,h$counts[1],h$counts[1]),col="grey")
abline(h=m0/20)

image

第1个柱子（灰色）表示提p值小于0.05的基因的数目，从水平线可知，大概有一半p值小于0.05的基因是假阳性，这与前面提到的FDR=0.5是一致的。我们来看一下当柱子的宽度为0.01时FDR的值会更低，但同时我们的显著性差异数目也会降低。

h <- hist(pvals,breaks=seq(0,1,0.01))
polygon(c(0,0.01,0.01,0),c(0,0,h$counts[1],h$counts[1]),col="grey")
abline(h=m0/100)

image

当柱子的宽度变为0.01时，我们可以看到一个更小的p值cutoff，并且检测到到特征值的数目也会下降（降低了灵敏性sensitivity)，但是，FDR也会降低（提高了特异性specificity）。我们如何设定这个cutoff呢？其中一个方法就是设定一个FDR水平 $\alpha$ ，然后设置控制错误率的程序即可： $FDR \leq \alpha$ 。

Benjamini-Hochberg

对于任意给定的 $\alpha$ ，Benjamini-Hochberg方法都非常适用，这种方法可以让使用者地每个检验的p值进行校正，也能很好地定义一个程序。

Benjamini-Hochberg方法的原理是，先按照升序对p值进行排列，即 $p_{(1)},\dots,p_{(m)}$ ，然后定义一个 $k$ 用来表示最大的秩 $i$ ，它所对应的p值为：
$p_{(i)} \leq \frac{i}{m}\alpha$
这个程序会拒绝那些p值小于或等于 $p_{(k)}$ 的检验。现在看一个案例，如何选择 $k$ 来计算p值，如下所示：

alpha <- 0.05
i = seq(along=pvals)
mypar(1,2)
plot(i,sort(pvals))
abline(0,i/m*alpha)
##close-up
plot(i[1:15],sort(pvals)[1:15],main="Close-up")
abline(0,i/m*alpha)

image

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

结果如下所示：

> cat("k =",k,"p-value cutoff=",cutoff)
k = 11 p-value cutoff= 3.763357e-05

我们可以从数学上证明到这个程序可以将FDR控制在5%以下，具体的算法可以参考Benjamini-Hochberg在1995年的论文。这种新的程序计算的结果就是将原来得到的2个有统计学差异的数值提高到了11个。如果我们将FDR的值设为50%，那么这个数字会提高到1063。FWER无法提供这种灵活性，因为只要检测值的数量增加，都会造成FWER的值为1。

这里我们需要注意的是，在R中，已经有了计算FDR的函数，前面的那些复杂代码只是为了说明这种算法，在R中我们可以使用下面的代码进行计算，如下所示：

fdr <- p.adjust(pvals, method="fdr")
mypar(1,1)
plot(pvals,fdr,log="xy")
abline(h=alpha,v=cutoff) ##cutoff was computed above

image

我们也可以使用Monte Carlo模拟来确认一下FDR的值实际上小于0.05，如下所示：

alpha <- 0.05
B <- 1000 ##number of simulations. We should increase for more precision
res <- replicate(B,{
  controls <- matrix(sample(population, N*m, replace=TRUE),nrow=m)
  treatments <- matrix(sample(population, N*m, replace=TRUE),nrow=m)
  treatments[which(!nullHypothesis),]<-treatments[which(!nullHypothesis),]+delta
  dat <- cbind(controls,treatments)
  pvals <- rowttests(dat,g)$p.value
  ##then the FDR
  calls <- p.adjust(pvals,method="fdr") < alpha
  R=sum(calls)
  Q=ifelse(R>0,sum(nullHypothesis & calls)/R,0)
  return(c(R,Q))
})
Qs <- res[2,]
mypar(1,1)
hist(Qs) ##Q is a random variable, this is its distribution

image

上述是Q值的直方图（假阳性除以显著性特征值的数目），现在看一下FDR值，如下所示：

FDR=mean(Qs)
print(FDR)

结果如下所示：

> FDR=mean(Qs)
> print(FDR)
[1] 0.03813818

这个FDR的值小于0.05，这个结果是符合预期的，因为我们需要保守一点，从而确保任何值的 $m_{0}$ 的FDR都小于0.05，例如当每个假设检验都是零的极端情况，例如 $m=m_{0}$ 。如果你重新模拟上述情况，你会发现FDR会增加。

我们需要注意下面的值：

> Rs <- res[1,]
> mean(Rs==0)*100
[1] 0.7

在模拟中，这个比例是0.7%，我们没有调用任何有显著性差异的基因。

在R中，p.adjust()函数可以选择一些算法来控制FDR，如下所示：

> p.adjust.methods
[1] "holm"       "hochberg"   "hommel"     "bonferroni" "BH"         "BY"         "fdr"        "none"

上面的这些方法不仅仅是估计错误率的不同方法，并且它们的估计值也不同，也就是说FWER与FDR不同。

直接FDR与q值

这里我们先回顾一下由John D.Storey在J. R. Statist. Soc. B(2002)中提到的结果。Storey与Benjamini-Hochberg方法的不同之处在于，前者不再提前设定 $\alpha$ 水平。因为在一些高通量实验中，我们感兴趣的仅仅是找到一个基因列表用于验证这些基因，我们会事先考虑将那些p值小于0.01的基因都进行验证。我们人随后会考虑估计的错误率。通常使用这些方法，我们会确保我们的 $R>0$ 。在上述定义的FDR里，当 $R=V=0$ 时，我们指定 $Q=0$ ，因此我们可以计算FDR如下所示：
$\mbox{FDR} = E\left( \frac{V}{R} \mid R>0\right) \mbox{Pr}(R>0)$
在Storey提出的方法里，我们需要构建一个非空列表，也就是说 $R>0$ ，那么我们计算阳性FDR(positive FDR)的公式如下所示：
$\mbox{pFDR} = E\left( \frac{V}{R} \mid R>0\right)$
第二点不同之处在于，Benjamini和Hochberg的程度是在最差的情况下进行控制，最差的情况是指零假设都为真（ $m=m_{0}$ 的情况），Storey的方法则是让我们从数据中估计 $m_{0}$ 。因为在高通量实验中，我们已经获得了如此多的数据，使Storey的算法成为了可能。这种算法的大致思路就是，挑选一个相对高值的p值cut-off，将它称为 $\lambda$ ，并且假设那些p值大于 $\Lambda$ 的检验在多数情况下其零假设是成立的，因此我们可以计算出估计值 $\pi_0 = m_0/m$ 为：
$\hat{\pi}_0 = \frac{\#\left\{p_i > \lambda \right\} }{ (1-\lambda) m }$
还有其它比这种算法更复杂的算法，但是它们的思路都是一样的，例如当我们设定 $\lambda=0.1$ 时，我们计算一下p值，如下所示：

hist(pvals,breaks=seq(0,1,0.05),freq=FALSE)
lambda = 0.1
pi0=sum(pvals> lambda) /((1-lambda)*m)
abline(h= pi0)
print(pi0) ##this is close to the trye pi0=0.9

image

> print(pi0) ##this is close to the trye pi0=0.9
[1] 0.9311111

根据这个估计，我们可以改变一下Benjamini-Hochberg程序，例如选择一个 $k$ 作为最大值（ $k$ 这里是最大的 $i$ ），因此就有如下公式：
$\hat{\pi}_0 p_{(i)} \leq \frac{i}{m}\alpha$
我们还有一种方法可以计算校正后的p值，即对每个检验计算q值(q-value)，如果一个特征计算的p值为 $p$ ，那么q值就是一系列含有最小p值为 $p$ 的特征值的估计pFDR。

除此之外，我们还有一种方法可以计算校正后的p值，即对每个检验计算q值(q-value)，如果一个特征最终计算的p值为 $p$ ，那么q值就是一系列含有尽可能最小与 $p$ 相等的p值的特征值的估计pFDR（原文是：However, instead of doing this, we compute a q-value for each test. If a feature resulted in a p-value of $p$ , the q-value is the estimated pFDR for a list of all the features with a p-value at least as small as $p$ .）。

上面这段我也不理解，后来查中文资料，根据Benjamini-Hochberg的算法，q值的定义如下所示：

对于m个独立的假设检验，它们对就的p值为： $p_{i},i=1,2,\cdots,m$

按照升序的方法对这些p值进行排序，得到：

$p_{(1)} \leq p_{(1)} \cdots \leq p_{(m)}$

对于给定的统计学检验水平 $\alpha in(0,1]$ ，找到最大的 $k$ ，使得：

$p_{(k)} \leq \frac{a*k}{m}$

对于排序靠前的 $k$ 个假设检验，认为它们是真阳性，看下面的案例：

现在我们做了6个基因的检验，它们的p值如下所示：

Gene	p-value
G1	p1=0.053
G2	p2=0.001
G3	p3=0.045
G4	p4=0.03
G5	p5=0.02
G6	p6=0.01

现在取检测水平 $\alpha=0.05$ ，先把p值按从小到大的升序排列：

Gene	p-value
G2	p2=0.001
G6	p6=0.01
G5	p5=0.02
G4	p4=0.03
G3	p3=0.045
G1	p1=0.053

取检测水平 $\alpha=0.05$ ，现在我们找到一个值 $k$ ，这个 $k$ 满足以下公式：
$p_{(k)} \leq \frac{a*k}{m}$
当 $k=1$ 时， $p_{(1)}=0.001 < 0.05*1/6=0.008333333$

当 $k=2$ 时， $p_{(2)}=0.01<0.05*2/6=0.01666667$

当 $k=3$ 时， $p_{(3)}=0.02<0.05*3/6=0.025$

当 $k=4$ 时， $p_{(4)}=0.03<0.05*4/6=0.03333333$

当 $k=5$ 时， $p_{(5)}=0.045>0.05*5/6=0.04166667$

当 $k=6$ 时， $p_{(6)}=0.053>0.05*6/6=0.05$

从上面的计算过程可以知道，这个 $k$ 等于4，也就是说，在FDR<0.05的情况下，G2，G6，G5，G4有差异。

到这里，只是说明了，G2，G6，G5，G4是有差异的，但是q值还没有算出来，继续计算：

前面我们已经把原始的p值按照从小到大的顺序排列好了，也就是[1] 0.001 0.010 0.020 0.030 0.045 0.053，其中最大的p值是0.053，它校正后还是这个值，也就是说这个值是校正后的最大p值，次大的p值是0.045，这个值需要校正，它排序是第5，那么校正的公式就是所有p值的数目（一共是6个p值）除以秩（这里是5），再乘以p值大小，也就是0.045*6/5=0.054，但是，这个值已经大于原来最大的p值（0.053）了，因此这个把它校正后的值也作为0.053，再看倒数第3个值，即0.03的校正p值，即0.03*6/4=0.045，它小于0.053，因此它校正后可以是0.045，现在依次校正剩下的值：

0.02*6/3=0.04，0.01*6/2=0.03，0.001*6/1=0.006，也就是说校正后的p值（就是q值），按从小到大的顺序排列就是：0.006，0.03，0.04，0.045，0.053，0.053，从结果中我们可以发现，前4个是有差异的，它们都小于0.05，也就是说G2，G6，G5，G4有差异。

其实公式就是：
$q^{(i)}=p_{(k)}^{(i)} * \frac{\text {Total Gene } N u m b e r}{\operatorname{rank}\left(p^{(i)}\right)}=p_{(k)}^{(i)} * \frac{m}{k}$
以上是BH法进行q值计算的过程，R中可以使用p.adjust()函数计算q值，所示：

p1 <- 0.053
p2 <- 0.001
p3 <- 0.045
p4 <- 0.03
p5 <- 0.02
p6 <- 0.01

p0 <- c(p1,p2,p3,p4,p5,p6)
p.adjust(p0,method = "BH")
p.adjust(sort(p0),method = "BH")
p.adjust(sort(p0),method = "fdr")
sort(p.adjust(p0,method = "BH"))

结果如下所示：

> p.adjust(p0,method = "BH")
[1] 0.053 0.006 0.053 0.045 0.040 0.030
> p.adjust(sort(p0),method = "BH")
[1] 0.006 0.030 0.040 0.045 0.053 0.053
> p.adjust(sort(p0),method = "fdr")
[1] 0.006 0.030 0.040 0.045 0.053 0.053
> sort(p.adjust(p0,method = "BH"))
[1] 0.006 0.030 0.040 0.045 0.053 0.053

从上面的计算结果看来：

method="BH"等于method="fdr"，结果没有区别；
在使用p.adjust()函数计算q值时，不用对原始数据进行排序，如果想要结果按从小到大的序列排列，那么就排序原始值，或计算后的值均可。

回到原文。

在R中，qvalue包中的qvalue()函数可以用来计算q值，如下所示：

library(qvalue)
res <- qvalue(pvals)
qvals <- res$qvalues
plot(pvals,qvals)

image

估计一下 $\hat{\pi_{0}}$ ，如下所示：

> res$pi0
[1] 0.8813727

练习

P274

基础探索数据分析

与低通量数据相比，高通量数据的一个被低估的优点就是它很容易展现数据。例如当我们有了海量的高通量数据后，我们很容易发现那些在少量数据时并不明显的问题。研究这些数据的一个强有力工具就是探索性数据分析(EDA,exploratory data analysis)。这一部分我们就来了解这方面的内容，可以参考作者的Github上的原文。

火山图

我们使用前面数据做的t检验结果来看一下火山图，如下所示：

library(genefilter)
library(GSE5859Subset)
data(GSE5859Subset)
g <- factor(sampleInfo$group)
results <- rowttests(geneExpression,g)
pvals <- results$p.value

我们还可以从一个数据集中生成一些p值，这些数据集中的零假设为真，如下所示：

m <- nrow(geneExpression)
n <- ncol(geneExpression)
randomData <- matrix(rnorm(n*m),m,n)
nullpvals <- rowttests(randomData,g)$p.value

我们前面提到过，当我们报告效应大小(effect size)时，如果仅仅计算p值则会造成一些错误。我们可以通过画一个火山图来展示高通量数据的统计结果。火山图背后的思想是，它能展示出所有的特征值（这里指的是基因表达谱）。火山图的y轴是-log(base 10) p-value，x轴是效应大小(effect size)。当我们经过-log(base 10)转换后，那些有着极显著差异的特征值就会出现在火山图的上方。这里使用log转换是因为，我们可用log转换可以更好地区分一些非常小的数据，例如区分0.01和 $10^6$ ，此时我们来绘制一个火山图，如下所示：

plot(results$dm,-log10(results$p.value),
xlab="Effect size",ylab="- log (base 10) p-values")

image

p值直方图

另外一种看整体数据的思路就是绘制p值的直方图。当我们的数据完全无效时，那么p值的直方图是服从均匀分布的，而我们的数据有效时，我们可以发现较小的p值频率较高，如下所示：

library(rafalib)
mypar(1,2)
hist(nullpvals,ylim=c(0,1400))
hist(pvals,ylim=c(0,1400))

image

左侧的p值是无效数据产生的p值，它服从均匀分布，右侧的p值则是则是我们基因表达谱的数据。

当我们预期大多数假设都是无效时，我们就不会看到服从均匀分布的p值，它也许是一些异常属性的指标，例如相关的样本。如果我们对结果时行置换检验，并计算出p值后，如果样本是独立的，那么我们应该会看到均匀分布，但是，我们的数据却看不到这个结果：

permg <- sample(g)
permresults <- rowttests(geneExpression,permg)
hist(permresults$p.value)

image

在后面部分中，我们会了解到这个数据集中的每个检验并不是独立的，因此这里用于检验的假设（我们使用的是t检验，而t检验的前提就是样本独立）是不正确的。

箱线图与直方图

高通量数据实验中，我们会检测每个实验单元的数千个特征值，我们前面也提到了，使用箱线图可以辅助我们来判断这些数据的质量。例如，如果一个样本的分布完全不同于剩余的样本，那么我们就可以认为，这个样本存在着一定问题。虽然一个完全的完全的变化可能是由于真正的生物学差异引起的，但是多数情况下，这就是技术因素造成的。这里我们使用Bioconductor中的基因表达数据，为了模拟出一组异常的数据，我们会对其中的一个样本进行log2转换，如下所示：

#BiocManager::install("Biobase")
#devtools::install_github("genomicsclass/GSE5859")
library(Biobase)
library(genefilter)
load("GSE5859.rda")
data(GSE5859)
ge <- exprs(e) ##ge for gene expression
ge[,49] <- ge[,49]/log2(exp(1)) ##immitate error

library(rafalib)
mypar(1,1)
boxplot(ge,range=0,names=1:ncol(e),col=ifelse(1:ncol(ge)==49,1,2))

运行过重中发现，GSE5859这个包无法安装，因此可以去官网下载原始文件，直接加载到RStudio中，另外在加载data(GSE5859)时会出错，不用管，运行就行，图形如下所示：

image

从上图我们可以看到，样本数据大，很难看清楚中间的箱子形状，但是我们很容易发现有一个样本异常，除此之外，我们可以用另外一种方式来展示一下数据，这个数据不画箱子，如下所示：

qs <- t(apply(ge,2,quantile,prob=c(0.05,0.25,0.5,0.75,0.95)))
matplot(qs,type="l",lty=1)

image

这种图形可以称为kaboxplot，因为这是由Karl Broman首次使用的，它绘制的是0.05，0.25，0.5，0.75和0.95分位数的图形。

我们也可以绘制直方图。因为我们的数据很多，因此可以使用很窄的间隔(bin)与柱子，然后将这些柱子进行平滑处理，最终绘制成平滑直方图(smooth histogram)。我们重新再校正这些平滑曲线的高度，那么在 $x_{0}$ 处的高度与基本单元构成的面积内，它们的点的数目就都比较接近，如下所示：

元区域内的点的数目就与这个基本单元的面积接近积，如下所示：

mypar(1,1)
shist(ge,unit=0.5)

image

MA图

散点图(scatterplot)与相关性(correlation)不是检测重复性问题的最佳工具。检测重复性更好的工作就是检测两次之间的差值，如果重复性好，那么这些差值应该是一样的。因此，一种更好的绘图工具就是将散点图旋转一下，这个散点图的y轴上差值，x轴是平均值。这种图最初被叫做Bland-Altman图，但在遗传学中它经常被称为MA图。MA的名称来源于图形的内容，M表示减(minus)，表示两值之差（有的时候是log值之差），A表示平均(Average)，现在绘制一下MA图

x <- ge[,1]
y <- ge[,2]
mypar(1,2)
plot(x,y)
plot((x+y)/2,x-y)

image

左图是常规的相关图，右图是MA图，需要注意的是，当我们把左图进行旋转后，这些数据的差异的SD就非常直观了：

sd(y-x)
[1] 0.2025465

左图的散点图显示这两个样本的相关性很强，但是显示的信息有限。

参考资料

如何理解与计算FDR？（第二版）

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 206,214评论 6赞 481
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 88,307评论 2赞 382
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 152,543评论 0赞 341
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 55,221评论 1赞 279
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 64,224评论 5赞 371
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 49,007评论 1赞 284
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 38,313评论 3赞 399
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,956评论 0赞 259
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 43,441评论 1赞 300
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,925评论 2赞 323
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 38,018评论 1赞 333
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,685评论 4赞 322
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 39,234评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 30,240评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,464评论 1赞 261
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 45,467评论 2赞 352
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,762评论 2赞 345

DALS017-高维数据推断2-统计原理

前言

高维数据的p值

错误率

数据案例

Bonferroni校正

错误发现率(FDR)

向量化代码

控制FDR

Benjamini-Hochberg

直接FDR与q值

练习

基础探索数据分析

火山图

p值直方图

箱线图与直方图

MA图

参考资料

推荐阅读更多精彩内容