5 局部錯(cuò)誤發(fā)現(xiàn)率

經(jīng)典單case假設(shè)檢驗(yàn)基于對(duì)統(tǒng)計(jì)量（p值）尾部的解釋闺兢。二戰(zhàn)后桥状，多重檢驗(yàn)繼續(xù)基于p值立肘，并擴(kuò)展到大規(guī)模假設(shè)檢驗(yàn)锹引，前面3和4章進(jìn)行了介紹。然而即使控制了錯(cuò)誤發(fā)現(xiàn)率，仍然與顯著性檢驗(yàn)和一類錯(cuò)誤想去甚遠(yuǎn)款慨。
對(duì)單例假設(shè)檢驗(yàn)來說，基于尾部區(qū)域是必須的，因?yàn)閦=1.96的概率是0。大規(guī)模檢驗(yàn)中蔗崎，允許進(jìn)行局部值推斷而不包含更極端值區(qū)域。這就是局部錯(cuò)誤發(fā)現(xiàn)率扰藕。

5.1 估計(jì)局部錯(cuò)誤發(fā)現(xiàn)率

由于每個(gè)case要么null要么non-null缓苛，可以用以下描述
$\pi_0=pr\{ null \}, f_0(z) = null density$
$\pi_1=pr\{ non-null \}, f_1(z) = non-null density$
局部錯(cuò)誤發(fā)現(xiàn)率為
$fdr(z) = Pr \{ null | z \} = \pi_0f_0(z) / f(z)$
其中 $f(z)$ 是混合概率密度函數(shù)
$f(z) = \pi_0f_0(z) + \pi_1f_1(z)$
本章假設(shè)已知 $f_0(z)$ ，而且基于上章中方法估計(jì) $\pi_0$ 邓深，則只剩 $f(z)$ 需要估計(jì)未桥，而且已知了觀測(cè)值 $\mathcal z = {z_1, z_2...z_N}$

基于泊松回歸方法估計(jì)

f(z)

，假設(shè)null下z服從標(biāo)準(zhǔn)正態(tài)分布芥备，并利用中心區(qū)域?yàn)?0%估算

\pi_0

冬耿，則可估算局部錯(cuò)誤發(fā)現(xiàn)率

\hat{fdr}(z) = \hat{\pi}_0f_0(z)/\hat{f}(z)

ps：這個(gè)估計(jì)值可能會(huì)超過1，這是因?yàn)閷?duì)

\pi_0

的估計(jì)離真實(shí)差太多萌壳，或者z并不服從

N(0,1)

亦镶。
對(duì)

\hat{f}(z)

進(jìn)行積分可以得到

\hat{Fdr}(z) = \hat{\pi}_0F_0(z)/\hat{F}(z)

如果我們?nèi)?img class="math-inline" src="https://math.jianshu.com/math?formula=fdr(z)%5Cleq%200.2" alt="fdr(z)\leq 0.2" mathimg="1">，則等同于 $\frac{f_1(z)}{f_0(z)} \geq 4\frac{\pi_0}{\pi_1}$ 袱瓮，如果我們假設(shè) $\pi_0\geq0.9$ 缤骨，則要求 $\frac{f_1(z)}{f_0(z)}\geq36$ 。這被稱為對(duì)抗零假設(shè)的貝葉斯因子懂讯。

5.2 $f(z)$ 的泊松回歸估計(jì)法

對(duì) $f(z)$ 的平滑估計(jì)荷憋，采用flexible exponential family models的MLE得到台颠。
例如設(shè) $f$ 屬于J-parameter famlily
$f(z) = exp \{ \sum_{j=0}^J\beta_jz^ j \}$
其中 $\beta_0$ 取決于 $\{ \beta_1,...\beta_J \}$ 褐望，以使 $f$ 的積分為1。當(dāng)J=2時(shí)串前，會(huì)使得 $f$ 為正態(tài)分布瘫里。

Lindsey’s method是一種基于離散的z值，使用標(biāo)準(zhǔn)泊松回歸方法荡碾，估算 $\beta$ 的最大似然估計(jì)的算法：

將 $z_i$ 的取值區(qū)間 $\mathcal Z$ 谨读，按照相等的范圍 $d$ 劃分為 $K$ 段：
定義 $y_k$ 為落入對(duì)應(yīng)區(qū)間中的觀測(cè)值數(shù)量

而 $x_k$ 是對(duì)應(yīng)區(qū)間的中心值，則 $y_k$ 的期望值近似為：
$\nu _k = Ndf(x_k)$
假設(shè) $y_k$ 是來自獨(dú)立的泊松分布
$y_k \sim Poi(\nu_k)$
然后擬合回歸模型
$log(\nu_k) = \sum_{j=0}^{J}\beta_j x_k^j$
以上是standard Poisson generalized linear model (GLM)坛吁。
以此得到的 $\hat{\beta} = (\hat{\beta}_1,\hat{\beta}_2,...\hat{\beta}_J)$ 是其模型的最大似然估計(jì)值劳殖。

5.3 統(tǒng)計(jì)推斷和局部錯(cuò)誤發(fā)現(xiàn)率

視角從 $Fdr$ 切換到 $fdr$ 更符合貝葉斯習(xí)慣：從貝葉斯角度來看 $fdr=Pr \{null | z \}$ 相對(duì)于觀測(cè)尾部概率 $Fdr = Pr \{null | Z \geq z \}$ 更合適。

上圖通過非參數(shù)估計(jì)得到

然而通過區(qū)間內(nèi)的個(gè)數(shù)進(jìn)行非參數(shù)估計(jì)非常不穩(wěn)定拨脉，如果用平滑版本

\hat{y}_k = N*d*\hat{f}(x_k)

進(jìn)行估計(jì)

\hat{fdr}_{78}=1.14/\hat{y}_k=0.254

更普通的結(jié)構(gòu)
5.1節(jié)中的模型可以更一般化使1...N基因?qū)?yīng)的結(jié)構(gòu)不同：
$\pi_{i0}=pr\{case\ i\ null \}, f_{i0}(z) = case\ i\ null density$
$\pi_{i1}=pr\{ case\ i\ non-null \}, f_{i1}(z) = case\ i\ non-null density$
如果定義
$\pi_{0} = \frac{1}{N}\sum_{i = 1}^N\pi_{i0}, f_{0}(z) = \frac{1}{N}\sum_{i = 1}^N\frac{\pi_{i0}}{\pi_0}f_{i0}(z)$
$\pi_{1} = \frac{1}{N}\sum_{i = 1}^N\pi_{i1}, f_{1}(z) = \frac{1}{N}\sum_{i = 1}^N\frac{\pi_{i1}}{\pi_1}f_{i1}(z)$
則我們又回到了5.1中的兩個(gè)分組的模型哆姻。
使用先驗(yàn)知識(shí)
之前我們的推斷都是基于我們不知道 $z_i$ (第i個(gè)基因)的信息，所以只能勉強(qiáng)使用兩個(gè)分組的模型玫膀。如果我們知道 $z_i$ 的先驗(yàn)信息矛缨，則
$fdr_i(z_i) = \pi_0 f_{i0}(z_i) / f_i(z_i) = Pr\{ case\ i\ null | z_i \}$
相比于 $fdr(z_i) = \pi_0 f_{0}(z_i) / f(z_i)$ 是一個(gè)更好的模型。
可交換性
比如取 $\hat{Fdr}(3.2) = 0.108$ ，我們會(huì)報(bào)道大于等于3.2的36個(gè)基因有較大可能確實(shí)與研究?jī)?nèi)容相關(guān)箕昭。但是它們其實(shí)顯著水平并不相同灵妨，對(duì)于數(shù)值更大的來說，它們的錯(cuò)誤發(fā)現(xiàn)率低于0.108落竹。
如果采用 $fdr$ 泌霍，則問題會(huì)小一些，比如我們會(huì)認(rèn)為[3.2, 3.3)間的錯(cuò)誤發(fā)現(xiàn)率為0.25筋量，而[3.3,3.4)間的錯(cuò)誤發(fā)現(xiàn)率為0.21烹吵。
Ps：當(dāng)然如果知道單個(gè)基因的先驗(yàn)知識(shí)，可交換性就沒有意義了桨武，應(yīng)采用前一部分的方法肋拔。
伸縮性
如果研究的假設(shè)增加會(huì)怎么樣？比如前面N個(gè)基因擴(kuò)大為2N個(gè)呀酸。
對(duì) $fdr$ 來說影響并不大凉蜂，基于前面的模型可知，增大為2N后只是讓均值更趨向于期望值性誉，會(huì)讓結(jié)果更精確窿吩。
然后對(duì)于傳統(tǒng)控制FWER的方法來說，會(huì)有特別大影響错览。比如對(duì)Bonferroni方法纫雁，會(huì)導(dǎo)致閾值從 $\alpha/N$ 降低到 $\alpha/(2N)$ 。
那么對(duì) $Fdr$ 呢倾哺？如果 $z_{(1)}$ 是最小的值轧邪，而且其對(duì)應(yīng)的p值為 $p_{(1)}=F_0(z_{(1)})$ ，則 $Fdr(z_{(1)})$ 等于 $N\pi_0p_{(1)}$ 羞海。如果 $p_{(1)} \leq \frac{1}{N}\frac{q}{\pi_0}$ 忌愚，就會(huì)導(dǎo)致錯(cuò)誤發(fā)現(xiàn)率小于控制目標(biāo)q。
增大檢驗(yàn)基因數(shù)却邓，另一方面會(huì)有相關(guān)性上的影響硕糊。之前的研究可能選擇的是人為認(rèn)為最相關(guān)的基因集合，如果數(shù)量擴(kuò)大一倍會(huì)導(dǎo)致集合與研究問題的相關(guān)性下降腊徙。
更多結(jié)構(gòu)的模型
如果N個(gè)基因來自M個(gè)天然的分類简十，我們可以根據(jù)每種分類運(yùn)用locfdr算法擬合，但是在小的分類中會(huì)引入評(píng)估問題撬腾。一種更好的做法是使用以下擴(kuò)展模型：
$log\{ f(z) \} = \sum_{j=0}^{J}\beta_jz^j + \gamma_{1m}z + \gamma_{2m}z^2$
其中m代表類別螟蝙，且 $\sum \gamma_{1m} = \sum \gamma_{2m} = 0$ 。它在保持了尾部特性同時(shí)时鸵，很好的兼容了不同均值和方差的類別胶逢。會(huì)在第10章討論厅瞎。
結(jié)合Fdr和fdr
其實(shí)沒必要選擇使用Fdr或fdr，它們可以合并使用初坠。它們間是可以轉(zhuǎn)換的和簸。
貝葉斯的局限
經(jīng)驗(yàn)貝葉斯推斷的是 $\hat{fdr}(z_i)$ 即 $Pr\{ case\ i \ null | z_i\}$ ，不一定等同于
$Pr\{ case\ i \ null | z_1, z_2,...,z_N\}$ 特別是z值有相關(guān)性的情況下碟刺。會(huì)在第9章討論锁保。
假陽(yáng)性和真陽(yáng)性的期望
局部錯(cuò)誤發(fā)現(xiàn)率控制下，對(duì)應(yīng)的假陽(yáng)性為“EFP”半沽，對(duì)應(yīng)的真陽(yáng)性為“ETP”爽柒。
如果我們按個(gè)體來看，如果對(duì) $z_i$ 拒絕的閾值為 $c_i$ 者填，則
$\alpha_i = \int_{c_i}^{\inf}f_{i0}(z_i)dz_i\ \ and \ \ \beta_i = \int_{c_i}^{\inf}f_{i1}(z_i)dz_i$
所以
$EFP = \sum_{i=1}^{N}w_i\pi_{i0}\alpha_i\ \ and \ \ ETP = \sum_{i=1}^{N}w_i\pi_{i1}\beta_i\$
其中 $w_i$ 是成為第i個(gè)的先驗(yàn)概率(可以取 $1/N$ )浩村。
我們期望的是：通過調(diào)整閾值 $c=(c_1,c_2,...,c_N)$ ，在給定EFP前提下占哟，最大化ETP心墅。
由于 $\frac{\mathrmftqbjzpEFP }{\mathrm3wwcukj c_i} = \sum_{i=1}^Nw_i\pi_{i0}\frac{\mathrm3gmlk5a\alpha_i }{\mathrmt0as7vu c_i} = -\sum_{i=1}^Nw_i\pi_{i0}f_{i0}(c_i)$
同樣的 $\frac{\mathrmw7w51olETP }{\mathrmcg7tbju c_i} = -\sum_{i=1}^Nw_i\pi_{i1}f_{i1}(c_i)$ ，應(yīng)用標(biāo)準(zhǔn)拉格朗日乘子法榨乎，對(duì)最佳的 $c$ 怎燥，存在常數(shù) $\lambda$ 使得
$\pi_{i1}f_{i1}(c_i) = \lambda\pi_{i0}f_{i0}(c_i)$
由于知道先驗(yàn)知識(shí)時(shí) $fdr_i(z_i) = \pi_0 f_{i0}(z_i) / f_i(z_i)$ ，則可推導(dǎo)出
$fdr_i(c_i) = 1 / (1 + \lambda)$
因此在給定EFP前提下最大化ETP：給定fdr下蜜暑，z值等于閾值時(shí)铐姚。

5.4 power診斷

之前的討論都主要專注于控制一類錯(cuò)誤(正如fdr其名字)，本節(jié)主要討論在局部錯(cuò)誤發(fā)現(xiàn)率控制下的power診斷肛捍。
定義正確發(fā)現(xiàn)率：local true discovery rate, tdr(z):
$tdr(z) = Pr \{ non-null |z \} = 1 - fdr(z) = \pi_1f_1(z) / f(z)$
則
$\hat{tdr}(z) = 1 - \hat{fdr}(z) = \hat{\pi}_1\hat{f}_1(z) / \hat{f}(z)$
其中
$\hat{\pi}_1 = 1 - \hat{\pi}_0\ and\ \hat{f}_1(z)= \frac{\hat{f}(z) -\hat{\pi}_0f_0(z)}{1 - \hat{\pi}_0}$
如果 $y_k$ 代表落在第k個(gè)區(qū)間的基因隐绵，當(dāng)然我們不能直接區(qū)分開null和non-null，但是可以進(jìn)行估算
$y_{1k} = \hat{tdr}_k*y_k$
由于 $y_k$ 來自區(qū)間統(tǒng)計(jì)篇梭，會(huì)有較大波動(dòng)(histogram noise)氢橙，一個(gè)更好的版本是結(jié)合之前的評(píng)估的概率密度函數(shù)：
$\hat{y}_{1k} = Nd \hat{tdr}_k * \hat{f}_k$
稱為：smoothed non-null counts

上圖是前面例子的對(duì)比酝枢。
此處有一個(gè)重要的區(qū)別恬偷，

fdr(z) = \pi_0f_0(z)/f(z)

中不再假設(shè)

f_0(z)

為

N(0,1)

，而是取empirical null

\hat{f}_0(z) = N(0, 1.06^2)

這會(huì)在第6章討論帘睦。
在圖中的105個(gè)smoothed non-null中袍患，只有26.8個(gè)發(fā)生在

\hat{fdr}(z) \leq 0.2

的區(qū)域中，約占26%竣付。也就是說這項(xiàng)研究的power很低诡延。

上圖中全部non-null的cdf是

Pr_1 \{ \hat{ fdr} \leq q \} = \sum_{k: \hat{fdr}_k \leq q} \hat{y}_{1k} / \sum_k \hat{y}_{1k}

圖中還有一個(gè)模擬的高power示例（虛線）。

一個(gè)簡(jiǎn)單直接的關(guān)于power的統(tǒng)計(jì)量是 $\hat{Efdr_1}$
$\hat{ Efdr_1} = \sum_k \hat{fdr}_k\hat{y}_{1k} / \sum_k \hat{y}_{1k}$

$\hat{Efdr_1}$ 值低代表power高（non-null大部分發(fā)生在fdr低的區(qū)域）古胆，反之為power低肆良。

上表展示了一個(gè)模擬筛璧，可以發(fā)現(xiàn)增加z的個(gè)數(shù)并不會(huì)明顯影響

\hat{Efdr_1}

，只是會(huì)讓bias變小惹恃，真實(shí)

Efdr_1 = 0.329

夭谤，大部分bias會(huì)讓

\hat{\pi}_0

評(píng)估偏大，從而降低

\hat{tdr}(z)

巫糙，這是因?yàn)椴捎昧?.44中的

\hat{\pi}_0

估計(jì)法導(dǎo)致的朗儒。

在這種場(chǎng)景下，研究員經(jīng)常會(huì)發(fā)現(xiàn)自己實(shí)驗(yàn)前認(rèn)為相關(guān)的基因常常不會(huì)fdr拒絕域內(nèi)参淹，這可能是因?yàn)榈蚿ower導(dǎo)致的醉锄。如前面所講的，結(jié)合先驗(yàn)知識(shí)可能會(huì)有助改善此類問題浙值。

最后值得注意的恳不，所有本節(jié)的power診斷都是基于不需要先驗(yàn)知識(shí)的前提下，這是大規(guī)模研究的優(yōu)點(diǎn)之一开呐。

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者

人面猴
序言：七十年代末妆够，一起剝皮案震驚了整個(gè)濱河市，隨后出現(xiàn)的幾起案子负蚊，更是在濱河造成了極大的恐慌神妹，老刑警劉巖，帶你破解...
沈念sama閱讀 216,324評(píng)論 6贊 498
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件家妆，死亡現(xiàn)場(chǎng)離奇詭異鸵荠，居然都是意外死亡，警方通過查閱死者的電腦和手機(jī)伤极，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 92,356評(píng)論 3贊 392
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進(jìn)店門蛹找，熙熙樓的掌柜王于貴愁眉苦臉地迎上來，“玉大人哨坪，你說我怎么就攤上這事庸疾。” “怎么了当编？”我有些...
開封第一講書人閱讀 162,328評(píng)論 0贊 353
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵届慈，是天一觀的道長(zhǎng)。經(jīng)常有香客問我忿偷，道長(zhǎng)金顿，這世上最難降的妖魔是什么？我笑而不...
開封第一講書人閱讀 58,147評(píng)論 1贊 292
?港島之戀（遺憾婚禮）
正文為了忘掉前任鲤桥，我火速辦了婚禮揍拆，結(jié)果婚禮上，老公的妹妹穿的比我還像新娘茶凳。我一直安慰自己嫂拴，他們只是感情好播揪，可當(dāng)我...
茶點(diǎn)故事閱讀 67,160評(píng)論 6贊 388
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布。她就那樣靜靜地躺著筒狠，像睡著了一般剪芍。火紅的嫁衣襯著肌膚如雪。梳的紋絲不亂的頭發(fā)上窟蓝，一...
開封第一講書人閱讀 51,115評(píng)論 1贊 296
城市分裂傳說
那天罪裹，我揣著相機(jī)與錄音，去河邊找鬼运挫。笑死状共，一個(gè)胖子當(dāng)著我的面吹牛，可吹牛的內(nèi)容都是我干的谁帕。我是一名探鬼主播峡继，決...
沈念sama閱讀 40,025評(píng)論 3贊 417
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼，長(zhǎng)吁一口氣：“原來是場(chǎng)噩夢(mèng)啊……” “哼匈挖！你這毒婦竟也來了碾牌？” 一聲冷哼從身側(cè)響起，我...
開封第一講書人閱讀 38,867評(píng)論 0贊 274
萬(wàn)榮殺人案實(shí)錄
序言：老撾萬(wàn)榮一對(duì)情侶失蹤儡循，失蹤者是張志新（化名）和其女友劉穎舶吗，沒想到半個(gè)月后，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體择膝，經(jīng)...
沈念sama閱讀 45,307評(píng)論 1贊 310
?護(hù)林員之死
正文獨(dú)居荒郊野嶺守林人離奇死亡誓琼，尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點(diǎn)故事閱讀 37,528評(píng)論 2贊 332
?白月光啟示錄
正文我和宋清朗相戀三年，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了肴捉。大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片腹侣。...
茶點(diǎn)故事閱讀 39,688評(píng)論 1贊 348
活死人
序言：一個(gè)原本活蹦亂跳的男人離奇死亡，死狀恐怖齿穗，靈堂內(nèi)的尸體忽然破棺而出傲隶，到底是詐尸還是另有隱情，我是刑警寧澤窃页，帶...
沈念sama閱讀 35,409評(píng)論 5贊 343
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布跺株，位于F島的核電站，受9級(jí)特大地震影響腮出，放射性物質(zhì)發(fā)生泄漏帖鸦。R本人自食惡果不足惜芝薇，卻給世界環(huán)境...
茶點(diǎn)故事閱讀 41,001評(píng)論 3贊 325
男人毒藥：我在死后第九天來索命
文/蒙蒙一胚嘲、第九天我趴在偏房一處隱蔽的房頂上張望。院中可真熱鬧洛二，春花似錦馋劈、人聲如沸攻锰。這莊子的主人今日做“春日...
開封第一講書人閱讀 31,657評(píng)論 0贊 22
一樁弒父案妓雾，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽(yáng)娶吞。三九已至，卻和暖如春械姻，著一層夾襖步出監(jiān)牢的瞬間妒蛇，已是汗流浹背。一陣腳步聲響...
開封第一講書人閱讀 32,811評(píng)論 1贊 268
情欲美人皮
我被黑心中介騙來泰國(guó)打工楷拳，沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留绣夺，地道東北人。一個(gè)月前我還...
沈念sama閱讀 47,685評(píng)論 2贊 368
代替公主和親
正文我出身青樓欢揖，卻偏偏與公主長(zhǎng)得像陶耍，于是被迫代替她去往敵國(guó)和親。傳聞我的和親對(duì)象是個(gè)殘疾皇子她混，可洞房花燭夜當(dāng)晚...
茶點(diǎn)故事閱讀 44,573評(píng)論 2贊 353

5 局部錯(cuò)誤發(fā)現(xiàn)率

5.1 估計(jì)局部錯(cuò)誤發(fā)現(xiàn)率

5.2 的泊松回歸估計(jì)法

5.3 統(tǒng)計(jì)推斷和局部錯(cuò)誤發(fā)現(xiàn)率

5.4 power診斷

推薦閱讀更多精彩內(nèi)容

5.2 $f(z)$ 的泊松回歸估計(jì)法