6 Theoretical, Permutation, and

2022-02-16 本文已影响0人老姚记事本

(不知道怎么翻译标题好了……)
在经典的假设检验中，零假设的分布是魔鬼的代言人：观察值必须超过的标准，以使科学界相信发生了一些有趣的事情（比如在零假设分布中超过1.96的中心距）。

第3章可以看出，学术界之前已经付出了很多努力使经典模型可适用于大规模推断场景。但当N特别大时，一些不同点让零假设分布的角色发生了变化：

例如当N = 10000时…显然基于经典零假设分布的模型不适合当前场景；
单假设检验场景下，常常希望拒绝零假设的power是80%。大规模推断场景下，研究者希望绝大多数假设接受零假设，仅仅保留少数有意思的点；
零假设分布的形状（比如 $\mu\sim N(0, \sigma^2)$ ），在大规模推断中相对没那么重要。N中大部分case会有一个很小的非零 $\mu$ 。这扩宽了经典零假设；
大规模的研究中允许经验贝叶斯分析，这会让零假设和备择假设混在一起；
大规模的研究中estimation和testing间的线变得模糊。

在之前的例子中，理论上的零假设分布表现的不错，这并不是常见的。例如在下面四个案例中就有严重问题，后续会基于这四个案例讨论。

6.1 四个案例

下面的图示中展示了以下信息：

N个case的直方图；
$\pi_0$ 的估计值 $\hat{\pi}_{00}$ （按第四章的方法，取 $\alpha_0$ 为0.5），理论零分布取 $f_0(z) \sim N(0,1)$ ；
按后面6.3方法估计的 $(\hat{\delta}_0, \hat{\sigma}_0, \hat{\pi}_0)$ ，零分布模型是 $\hat{f}_0(z) \sim N(\hat{\delta}_0, \hat{\pi}_0)$ ；
粗实线代表经验分布；
一条轻的点曲线成比例于 $\hat{\pi}_{00} * \varphi(z)$ ，其中 $\varphi(z) = exp\{ -\frac{1}{2}*z^2 \} / \sqrt{2 * \pi}$ ；
x轴上的小三角标注出按经验分布评估的fdr小于等于0.2的阈值，混合分布评估采用的是前一章的方法；
x轴下的小红线标出通过fdr过滤的z值；
通过过滤的z值数量。

A. 白血病的研究

高密度寡核苷酸微阵列：N=7128对基因，供72个病人参与研究，其中 $n_1 = 45$ 个ALL（急性淋巴细胞白血病），其中 $n_2 = 27$ 个AML（急性髓性白血病），后者更加严重。
原始阵列已经转换为了一个normal score：
$x_{ij} = \Phi^{-1}(\frac{rank(X_{ij})-0.5}{N})$
其中 $X_{ij}$ 代表j个病人的i基因， $rank(X_{ij})$ 代表 $X_{ij}$ 在N中的排名。Z值来自ALL和AML的双样本t检验。
图中经验零分布为 $N(0.09, 1.68*2)$ 且 $\hat{\pi}_0 = 0.937$ ，其中173个基因的 $\hat{fdr}(z) \leq 20%$ 。
理论零分布为 $N(0,1)$ ，对应的 $\pi_{00}$ 是0.654，有1548的 $\hat{fdr}(z) \leq 20%$ 。也许 $(1 - \pi_{00}) * N = 2464$ 个基因会不同，但更大可能是理论上的零分布不合时宜。

B. 卡方数据

本实验研究了N = 16 882个基因中某些化学标签在位点的结合。每个基因的 K 位点数从三个到几百个不等，中位数K=12。在每个基因内的每个位点，对结合标签的数量进行计数。计数是在两种不同的实验条件下进行的，研究的目的是识别两种条件下标签比例不同的基因。
下表中统计了两组方案中第一个K点数量的分布。

i基因对应的

z_i

从

table_i

中算出：
(i) 按将表里每个cell进行计算；
(ii) 对表计算

S_i

——一个独立情况下的卡方检验统计值；
(iii) 通过卡方统计值计算p值；
(iv) 转化为

z_i

：

z_i = \Phi^{-1}(1- p_i)

方法中不需要标准卡方定义等检验统计量的经典形式，但它们确实依赖于能够用正态曲线逼近 z 值直方图的中心。这导致了同时推理中的可比性和相关性问题，会在第10章讨论。

其经验零分布为 $N(0.32, 1.25^2)$ ，只有10个基因的预估fdr小于0.2。

C. 警方数据

2006年进行了一项关于纽约警察要求行人停止是否有种族歧视的研究。计算了N = 2749个警察的代表偏见程度的数值 $z_i$ 。
定义 $x_{ij}$ 是警官i在对j次停止的协变量。一个简化的逻辑回归模型是
$logit\{ Pr(y_{ij} = 1)\} = \beta_i + \gamma'x_{ij}$
其中 $y_{ij}$ 代表被停止的人是不是少数族裔， $\beta_i$ 代表“警官效果”， $\gamma$ 是协变量的回归系数向量。警察i的z值是
$z_i = \hat{\beta}_i / se(\hat{\beta}_i)$
其中 $\hat{\beta}_i$ 是评估值， $se(\hat{\beta}_i)$ 是 $\beta_i$ 的标准误。

这个例子中理论零假设与经验零假设有的结果有巨大差异。

D. HIV数据

对 $n_1=4$ 个健康人与 $n_2=4$ 个HIV阳性者的 $N=7680$ 个基因进行研究。通过双样本T检验计算p值，在转换为正态情况下对应的z值。

这个例子的经验零假设比较接近理论零假设，

\hat{f}_0(z) = N(0.12, 0.77^2)

，

\hat{\pi}_0=0.949

。
下面展示了一个人造的例子。
基于贝叶斯层次模型

\mu \sim g(.)

且

z|\mu \sim N(\mu, 1)

，此时对

g

选择为

g(\mu) = 0.9*\varphi_{0,0.05}(\mu) + 0.1 * \varphi_{2.5,0.5}(\mu)

公式中

\varphi_{a,b}

代表

N(a, b^2)

分布的密度函数。
然后此时混合后z值得密度函数不是单峰的，如下图

此时真正感兴趣的部分应该是

\mu>1.5

的区间。
假设不知道它们的真实先验，通过模拟数据估计会得到

\hat{\pi}_0 = 0.93

，

\hat{f}_0(z) = N(0.02, 1.14^2)

。可以即使不知道确切先验，通过观测值估计，仍然能很好的找出感兴趣的区域。

6.2 评估经验零假设

上面四个例子展示了理论零假设不太合理。经验零假设通过数据评估一个合适的零分布。零假设占比一般很高，设 $\pi_0 \leq 0.9$ ，给了我们零分布的可能。
定义
$f_{\pi_0}(z) = \pi_0f_0(z)$
则
$fdr(z) = f_{\pi_0}(z) / f(z)$
如果设 $f_0(z)$ 是正态分布，但不是标准正态分布：
$f_0(z) \sim N(\delta_0, \sigma_0^2)$
这会得到关于z值的二次函数：
$log(f_{\pi_0}(z)) = [log(\pi_0) - \frac{1}{2}\{ \frac{\delta_0^2}{\sigma^2_0} + log(2\pi\sigma^2_0)\}] + \frac{\delta_0}{\sigma_0^2} z- \frac{1}{2\sigma_0^2}z^2$
通过Central matching法假设 $log(f(z))$ 在 $z=0$ 附近是一个二次函数来评估 $f_0(z)$ 和 $\pi_0$ ：
$log(f(z)) \doteq \beta_0 + \beta_1z + \beta_2 z^2$
通过在 $z=0$ 另附近的数量 $y_k$ 来评估 $(\beta_0,\beta_1,\beta_2)$ 并与原公式进行匹配：比如 $\sigma^2 = -1 / (2\beta_2)$

上图展示了HIV数据的计算。用前面5.2中的方法，通过中心附近的z值拟合（由于这区间内

\pi_1

非常小，可以近似）。
显然这个评估是有偏的，但是以下模型时（第二章提到的模型），效果近似于无偏：

\mu \sim g(.) \ and \ z|\mu \sim N(\mu, 1)

g(\mu) = \pi_0I_0(\mu) + \pi_1g_1(\mu)

通过以下模拟可以证明在

\pi_0

比较大时，通过上述方法可以得到很好的估计效果：
取

f_0(z) = \varphi(z)

为标准正态分布，

f_1(z) = \int_{-\infty }^{\infty}g(\mu)\varphi(z-\mu)d\mu

和固定的

\pi_0

来模拟，根据观测值评估的

(\delta_0, \sigma_0)

，

(\delta_g,\sigma_g)

是central matching评估结果：

\delta_g = argmax\{ f(z) \}\ and \ \sigma_g = [-\frac{d^2}{dz^2}logf(z)]_{\delta_g} ^ {-\frac{1}{2}}

可以对比

(\delta_g,\sigma_g)

与真实值

(0,1)

差多少。对指定

\pi_0

，定义评估最差的情况：

\delta_{max} = max\{ |\delta_g| \}\ and\ \sigma_{max} = max\{ \sigma_g\}

根据下表结果可知在

\pi_0>=0.9

的情况下，central matching评估的偏差不严重。

上图以

\delta_{max}

和

\sigma_{max}

作为

\pi_1 = 1 - \pi_0

的函数进行了展示。图中还画出了限制了

g_1

在0处对称、对称且正态时的情况。
locfdr包中默认使用的是MLE方法，而不是central matching。因为中心直方图中的轻微不规则性，可能会破坏中心匹配。MLE更稳定，但是可能增大bias。

6.3 MLE经验零分布

MLE是一种更直接的方式。基于认为落在中心几乎全是零假设的z值集合评估 $(\hat{\delta}_0,\hat{\sigma}_0,\hat{\pi}_0 )$ 。相比上一节的方法，波动性更小但更容易偏差。
全集为 $z = (z_1,z_2,...,z_N)$ ， $N_0$ 是选中的集合， $I_0$ 是他们的索引：
$I_0 = \{ i:z_i \in \mathcal A_0\}\ \ and \ N_0 = \#I_0\ and \ z_0 = \{ z_i, i \in I_0\}$
并且 $\varphi_{\delta_0, \sigma_0}(z)$ 是 $N(\delta_0, \sigma_0)$ 的密度函数，则落入区域的概率为
$H_0(\delta_0, \sigma_0) \equiv \int_{ \mathcal A_0}\varphi_{\delta_0, \sigma_0}(z)dz$
假设z值独立且来自wo-groups模型： $f_0 \sim N(\delta_0, \sigma_0)$ ， $f_1(z) = 0\ for\ z \in \mathcal A_0$ 。
则 $z_0$ 的似然函数为：
$f_{\delta_0, \sigma_0,\pi_0}(z_0) = [\binom{N}{N_0}\theta^{N_0}(1 - \theta)^{N-N0}][\prod _{I_0}\frac{\varphi_{\delta_0, \sigma_0}(z)}{H_0(\delta_0, \sigma_0)}]$
其中 $\theta = \pi_0H_0(\delta_0, \sigma_0) = Pr\{ z_i \in \mathcal A_0\}$

下表是一个蒙特卡洛模拟结果

通过上表可知，MLE方法相比CM方法有更小的标准差，但是偏差更大（模拟的

\delta_0=-0.125

），特别是

\pi_0

。
其中

\theta

的估计值是

N_0/N

，

\delta

和

\sigma

的估计值通过MLE得到，因此可计算得到

\hat{\pi}_0 = \hat{\theta} / H_0(\hat{\delta}_0, \hat{\sigma}_0)

6.4 为何理论零假设失效

控制Fdr的关注点，是找到可以控制的中心距，而不是 $N(0,1)$ 相对的距离。以下是无法使用 $N(0,1)$ 的常见原因：

(I) 违背了数学假设

比如双样本t检验，常常假设样本来自独立同分布的正态分布；

(II) 随机单元之间的联系

不像双样本随机实验可以保证随机采样，很多情况下是自然实验；

(III) 检验结果之间的相关性

即使每个z值服从标准正态分布，但是z值之间的相关性导致理想零分布无法控制错误发现率。
下面是一个模拟的例子，零分布的z值具有相关性，导致通过理论零假设控制Fdr时的效果较差。

(IV) 未观测到的协变量

比如白血病的研究并不是一个随机实验， AML/ALL是通过观测区分的，还有其它未观测的协变量比如年龄、性别、健康程度等等，它们也会影响结果。

6.5 置换零分布

置换技术介于理论零分布与经验零分布之间，但更偏向于前者。
将原来的两组打散，再随机分组，并产生B组结果，得到一个N*B的矩阵：
$Z^* = (z^{*1}, z^{*2},...,z^{*B})$
此时一般的置换零分布为：
$\hat{f}^{perm}_0 = N · B个z_i^{*b}值的经验分布$

置换的零分布也会出现失效的情况。考虑以下几点：

对于上一节的4种理论零假设失效的场景，置换零分布可以很好的适用于(I)，因为它是基于随机排列模拟的
无法解决(II)，重排列基于假设样本间独立
无法解决(IV)；
置换方法的一个优点是它们保留了案例之间的相关性，然而无助于场景(III)；
事实上置换零分布会分非常接近 $N(0,1)$ ；
置换方法和经验方法可以结合；
置换方法不仅局限于两个分组的场景。