午夜精品aaa,欧美中文一区,久久青草日本精品久久

title: DALS017-高維數據推斷2-統(tǒng)計原理
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小節(jié)，這一部分的主要內容涉及高維數據統(tǒng)計的一些原理，相應的R markdown文檔可以參考作者的Github。

高維數據的p值

在前面我們了解了，當我們在分析高維數據時，p值就不再是一個很好的統(tǒng)計指標了。這是因為，我們同一次檢測了許多特征值(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)，它會讓我們對統(tǒng)計過程進行估計和控制。

錯誤率

在這一部分中，我們會了解到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ǒng)計，每次都采用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個飼料檢驗的實驗里，這個假設是成立的；但是在一些遺傳學實驗里，這個假設有可能不成立，因為某些基因會共同發(fā)揮作用）。每次檢驗得到的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$ ，那么我們的統(tǒng)計就是失敗的。

我們如何做才能使得一個錯誤出現的概率低于 $\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值是一個隨機變量，我們會使用統(tǒng)計理論來計算，如果我們遵循這個程序，平均會犯多少錯誤。更確切地講就是，我們可以計算出這些錯誤的邊界，也就是說，這些錯誤小于某些預定值的比例。

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次飲食實驗，卻只有發(fā)現2個假陽性的結果，Bonferroni是一種非常嚴格的校正。

錯誤發(fā)現率(FDR)

在許多情況下，要求FWER是0.05沒多大意義，因為它太嚴格了。例如，我們常見的做法就是先進行初步的小型實驗來確定小數候選基因。這種做法被稱為之發(fā)現驅動的實驗或項目。我們也許會尋找一個未知的致病基因，而不僅僅是對候選基因進行采用更多的樣本進行后續(xù)研究。如果我們開發(fā)一個程序，例如一個10個基因列表，從中發(fā)現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，即錯誤發(fā)現率(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()系列函數來加快運行速度，前面已經看到了這個函數的使用方法，現在使用傳統(tǒng)的代碼來看一下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%的統(tǒng)計檢驗而言，零假設是真。這就暗示了，由于p值是cutoff值為0.05，100個檢驗中，大概有4到5個是假陽性。再加上我們沒有考慮到替代所設為真時的所有情況，因此FDR值就比較高。那么我們如何控制它呢，如果我們想要更低的FDR值，比如5%怎么辦？

為了用圖形說明FDR為什么這么高，我們來繪制p值的直方圖。我們使用一個較大的m值從直方圖中獲取更多的數據。再繪制一條水平線，表示當NULL為真時，從 $m_{0}$ 個案例（指的基因或其它檢測值， $m_{0}$ 就是沒有統(tǒng)計學差異的基因）中到的陽性結果的均勻分布，如下所示：

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

當柱子的寬度變?yōu)?.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)}$ 的檢驗?，F在看一個案例，如何選擇 $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個有統(tǒng)計學差異的數值提高到了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}$ 。如果你重新模擬上述情況，你會發(fā)現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)}$

對于給定的統(tǒng)計學檢驗水平 $\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值還沒有算出來，繼續(xù)計算：

前面我們已經把原始的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，從結果中我們可以發(fā)現，前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"，結果沒有區(qū)別；
在使用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

基礎探索數據分析

與低通量數據相比，高通量數據的一個被低估的優(yōu)點就是它很容易展現數據。例如當我們有了海量的高通量數據后，我們很容易發(fā)現那些在少量數據時并不明顯的問題。研究這些數據的一個強有力工具就是探索性數據分析(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值則會造成一些錯誤。我們可以通過畫一個火山圖來展示高通量數據的統(tǒng)計結果。火山圖背后的思想是，它能展示出所有的特征值（這里指的是基因表達譜）?；鹕綀D的y軸是-log(base 10) p-value，x軸是效應大小(effect size)。當我們經過-log(base 10)轉換后，那些有著極顯著差異的特征值就會出現在火山圖的上方。這里使用log轉換是因為，我們可用log轉換可以更好地區(qū)分一些非常小的數據，例如區(qū)分0.01和 $10^6$ ，此時我們來繪制一個火山圖，如下所示：

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

image

p值直方圖

另外一種看整體數據的思路就是繪制p值的直方圖。當我們的數據完全無效時，那么p值的直方圖是服從均勻分布的，而我們的數據有效時，我們可以發(fā)現較小的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中的基因表達數據，為了模擬出一組異常的數據，我們會對其中的一個樣本進行l(wèi)og2轉換，如下所示：

#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))

運行過重中發(fā)現，GSE5859這個包無法安裝，因此可以去官網下載原始文件，直接加載到RStudio中，另外在加載data(GSE5859)時會出錯，不用管，運行就行，圖形如下所示：

image

從上圖我們可以看到，樣本數據大，很難看清楚中間的箱子形狀，但是我們很容易發(fā)現有一個樣本異常，除此之外，我們可以用另外一種方式來展示一下數據，這個數據不畫箱子，如下所示：

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}$ 處的高度與基本單元構成的面積內，它們的點的數目就都比較接近，如下所示：

元區(qū)域內的點的數目就與這個基本單元的面積接近積，如下所示：

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

左圖是常規(guī)的相關圖，右圖是MA圖，需要注意的是，當我們把左圖進行旋轉后，這些數據的差異的SD就非常直觀了：

sd(y-x)
[1] 0.2025465

左圖的散點圖顯示這兩個樣本的相關性很強，但是顯示的信息有限。

參考資料

如何理解與計算FDR？（第二版）

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

DALS017-高維數據推斷2-統(tǒng)計原理

DALS017-高維數據推斷2-統(tǒng)計原理

前言

高維數據的p值

錯誤率

數據案例

Bonferroni校正

錯誤發(fā)現率(FDR)

向量化代碼

控制FDR

Benjamini-Hochberg

直接FDR與q值

練習

基礎探索數據分析

火山圖

p值直方圖

箱線圖與直方圖

MA圖

參考資料

相關閱讀更多精彩內容

友情鏈接更多精彩內容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九 欧美,1769亚洲,黄色成人av

DALS017-高維數據推斷2-統(tǒng)計原理

前言

高維數據的p值

錯誤率

數據案例

Bonferroni校正

錯誤發(fā)現率(FDR)

向量化代碼

控制FDR

Benjamini-Hochberg

直接FDR與q值

練習

基礎探索數據分析

火山圖

p值直方圖

箱線圖與直方圖

MA圖

參考資料

相關閱讀更多精彩內容

友情鏈接更多精彩內容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av