5 局部错误发现率

经典单case假设检验基于对统计量（p值）尾部的解释。二战后，多重检验继续基于p值，并扩展到大规模假设检验，前面3和4章进行了介绍。然而即使控制了错误发现率，仍然与显著性检验和一类错误想去甚远。
对单例假设检验来说，基于尾部区域是必须的，因为z=1.96的概率是0。大规模检验中，允许进行局部值推断而不包含更极端值区域。这就是局部错误发现率。

5.1 估计局部错误发现率

由于每个case要么null要么non-null，可以用以下描述
$\pi_0=pr\{ null \}, f_0(z) = null density$
$\pi_1=pr\{ non-null \}, f_1(z) = non-null density$
局部错误发现率为
$fdr(z) = Pr \{ null | z \} = \pi_0f_0(z) / f(z)$
其中 $f(z)$ 是混合概率密度函数
$f(z) = \pi_0f_0(z) + \pi_1f_1(z)$
本章假设已知 $f_0(z)$ ，而且基于上章中方法估计 $\pi_0$ ，则只剩 $f(z)$ 需要估计，而且已知了观测值 $\mathcal z = {z_1, z_2...z_N}$

基于泊松回归方法估计

f(z)

，假设null下z服从标准正态分布，并利用中心区域为50%估算

\pi_0

，则可估算局部错误发现率

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

ps：这个估计值可能会超过1，这是因为对

\pi_0

的估计离真实差太多，或者z并不服从

N(0,1)

。
对

\hat{f}(z)

进行积分可以得到

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

如果我们取 $fdr(z)\leq 0.2$ ，则等同于 $\frac{f_1(z)}{f_0(z)} \geq 4\frac{\pi_0}{\pi_1}$ ，如果我们假设 $\pi_0\geq0.9$ ，则要求 $\frac{f_1(z)}{f_0(z)}\geq36$ 。这被称为对抗零假设的贝叶斯因子。

5.2 $f(z)$ 的泊松回归估计法

对 $f(z)$ 的平滑估计，采用flexible exponential family models的MLE得到。
例如设 $f$ 属于J-parameter famlily
$f(z) = exp \{ \sum_{j=0}^J\beta_jz^ j \}$
其中 $\beta_0$ 取决于 $\{ \beta_1,...\beta_J \}$ ，以使 $f$ 的积分为1。当J=2时，会使得 $f$ 为正态分布。

Lindsey’s method是一种基于离散的z值，使用标准泊松回归方法，估算 $\beta$ 的最大似然估计的算法：

将 $z_i$ 的取值区间 $\mathcal Z$ ，按照相等的范围 $d$ 划分为 $K$ 段：
定义 $y_k$ 为落入对应区间中的观测值数量

而 $x_k$ 是对应区间的中心值，则 $y_k$ 的期望值近似为：
$\nu _k = Ndf(x_k)$
假设 $y_k$ 是来自独立的泊松分布
$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)$ 是其模型的最大似然估计值。

5.3 统计推断和局部错误发现率

视角从 $Fdr$ 切换到 $fdr$ 更符合贝叶斯习惯：从贝叶斯角度来看 $fdr=Pr \{null | z \}$ 相对于观测尾部概率 $Fdr = Pr \{null | Z \geq z \}$ 更合适。

上图通过非参数估计得到

然而通过区间内的个数进行非参数估计非常不稳定，如果用平滑版本

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

进行估计

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

更普通的结构
5.1节中的模型可以更一般化使1...N基因对应的结构不同：
$\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中的两个分组的模型。
使用先验知识
之前我们的推断都是基于我们不知道 $z_i$ (第i个基因)的信息，所以只能勉强使用两个分组的模型。如果我们知道 $z_i$ 的先验信息，则
$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)$ 是一个更好的模型。
可交换性
比如取 $\hat{Fdr}(3.2) = 0.108$ ，我们会报道大于等于3.2的36个基因有较大可能确实与研究内容相关。但是它们其实显著水平并不相同，对于数值更大的来说，它们的错误发现率低于0.108。
如果采用 $fdr$ ，则问题会小一些，比如我们会认为[3.2, 3.3)间的错误发现率为0.25，而[3.3,3.4)间的错误发现率为0.21。
Ps：当然如果知道单个基因的先验知识，可交换性就没有意义了，应采用前一部分的方法。
伸缩性
如果研究的假设增加会怎么样？比如前面N个基因扩大为2N个。
对 $fdr$ 来说影响并不大，基于前面的模型可知，增大为2N后只是让均值更趋向于期望值，会让结果更精确。
然后对于传统控制FWER的方法来说，会有特别大影响。比如对Bonferroni方法，会导致阈值从 $\alpha/N$ 降低到 $\alpha/(2N)$ 。
那么对 $Fdr$ 呢？如果 $z_{(1)}$ 是最小的值，而且其对应的p值为 $p_{(1)}=F_0(z_{(1)})$ ，则 $Fdr(z_{(1)})$ 等于 $N\pi_0p_{(1)}$ 。如果 $p_{(1)} \leq \frac{1}{N}\frac{q}{\pi_0}$ ，就会导致错误发现率小于控制目标q。
增大检验基因数，另一方面会有相关性上的影响。之前的研究可能选择的是人为认为最相关的基因集合，如果数量扩大一倍会导致集合与研究问题的相关性下降。
更多结构的模型
如果N个基因来自M个天然的分类，我们可以根据每种分类运用locfdr算法拟合，但是在小的分类中会引入评估问题。一种更好的做法是使用以下扩展模型：
$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$ 。它在保持了尾部特性同时，很好的兼容了不同均值和方差的类别。会在第10章讨论。
结合Fdr和fdr
其实没必要选择使用Fdr或fdr，它们可以合并使用。它们间是可以转换的。
贝叶斯的局限
经验贝叶斯推断的是 $\hat{fdr}(z_i)$ 即 $Pr\{ case\ i \ null | z_i\}$ ，不一定等同于
$Pr\{ case\ i \ null | z_1, z_2,...,z_N\}$ 特别是z值有相关性的情况下。会在第9章讨论。
假阳性和真阳性的期望
局部错误发现率控制下，对应的假阳性为“EFP”，对应的真阳性为“ETP”。
如果我们按个体来看，如果对 $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个的先验概率(可以取 $1/N$ )。
我们期望的是：通过调整阈值 $c=(c_1,c_2,...,c_N)$ ，在给定EFP前提下，最大化ETP。
由于 $\frac{\mathrm{d}EFP }{\mathrm{d} c_i} = \sum_{i=1}^Nw_i\pi_{i0}\frac{\mathrm{d}\alpha_i }{\mathrm{d} c_i} = -\sum_{i=1}^Nw_i\pi_{i0}f_{i0}(c_i)$
同样的 $\frac{\mathrm{d}ETP }{\mathrm{d} c_i} = -\sum_{i=1}^Nw_i\pi_{i1}f_{i1}(c_i)$ ，应用标准拉格朗日乘子法，对最佳的 $c$ ，存在常数 $\lambda$ 使得
$\pi_{i1}f_{i1}(c_i) = \lambda\pi_{i0}f_{i0}(c_i)$
由于知道先验知识时 $fdr_i(z_i) = \pi_0 f_{i0}(z_i) / f_i(z_i)$ ，则可推导出
$fdr_i(c_i) = 1 / (1 + \lambda)$
因此在给定EFP前提下最大化ETP：给定fdr下，z值等于阈值时。

5.4 power诊断

之前的讨论都主要专注于控制一类错误(正如fdr其名字)，本节主要讨论在局部错误发现率控制下的power诊断。
定义正确发现率：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个区间的基因，当然我们不能直接区分开null和non-null，但是可以进行估算
$y_{1k} = \hat{tdr}_k*y_k$
由于 $y_k$ 来自区间统计，会有较大波动(histogram noise)，一个更好的版本是结合之前的评估的概率密度函数：
$\hat{y}_{1k} = Nd \hat{tdr}_k * \hat{f}_k$
称为：smoothed non-null counts

上图是前面例子的对比。
此处有一个重要的区别，

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

中不再假设

f_0(z)

为

N(0,1)

，而是取empirical null

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

这会在第6章讨论。
在图中的105个smoothed non-null中，只有26.8个发生在

\hat{fdr}(z) \leq 0.2

的区域中，约占26%。也就是说这项研究的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}

图中还有一个模拟的高power示例（虚线）。

一个简单直接的关于power的统计量是 $\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大部分发生在fdr低的区域），反之为power低。

上表展示了一个模拟，可以发现增加z的个数并不会明显影响

\hat{Efdr_1}

，只是会让bias变小，真实

Efdr_1 = 0.329

，大部分bias会让

\hat{\pi}_0

评估偏大，从而降低

\hat{tdr}(z)

，这是因为采用了4.44中的

\hat{\pi}_0

估计法导致的。

在这种场景下，研究员经常会发现自己实验前认为相关的基因常常不会fdr拒绝域内，这可能是因为低power导致的。如前面所讲的，结合先验知识可能会有助改善此类问题。

最后值得注意的，所有本节的power诊断都是基于不需要先验知识的前提下，这是大规模研究的优点之一。

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

5 局部错误发现率

5.1 估计局部错误发现率

5.2 的泊松回归估计法

5.3 统计推断和局部错误发现率

5.4 power诊断

推荐阅读更多精彩内容

5.2 $f(z)$ 的泊松回归估计法