分类目录归档:统计推断

数理统计(传统参数估计、假设检验等理论)

因果推断简介之四:观察性研究,可忽略性和倾向得分

这节采用和前面相同的记号。$Z$ 表示处理变量($1$ 是处理,$0$ 是对照),$Y$ 表示结果,$X$ 表示处理前的协变量。在完全随机化试验中,可忽略性 $Z \bot \{Y(1), Y(0)\} $ 成立,这保证了平均因果作用 $ACE(Z\rightarrow Y) = E\{Y(1) – Y(0)\} = E\{Y\mid Z=1\} – E\{Y\mid Z=0\}$ 可以表示成观测数据的函数,因此可以识别。在某些试验中,我们“先验的”知道某些变量与结果强相关,因此要在试验中控制他们,以减少试验的方差。在一般的有区组(blocking)的随机化试验中,更一般的可忽略性 $Z \bot \{Y(1), Y(0)\} | X$ 成立,因为只有在给定协变量 $ X $ 后,处理的分配机制才是完全随机化的。比如,男性和女性中,接受处理的比例不同,但是这个比例是事先给定的。

在传统的农业和工业试验中,由于随机化,可忽略性一般是能够得到保证的;因此在这些领域谈论因果推断是没有太大问题的。Jerzy Neyman 最早的博士论文,就研究的是农业试验。但是,这篇写于 1923 年的重要统计学文章,迟迟没有得到统计学界的重视,也没有人将相关方法用到社会科学的研究中。1970 年代,Donald Rubin 访问 UC Berkeley 统计系,已退休的 Jerzy Neyman 曾问起:为什么没有人将潜在结果的记号用到试验设计之外?正如 Jerzy Neyman 本人所说 “without randomization an experiment has little value irrespective of the subsequent treatment(没有随机化的试验价值很小)”,人们对于观察性研究中的因果推断总是抱着强烈的怀疑态度。我们经常听到这样的声音:统计就不是用来研究因果关系的!

在第一讲 Yule-Simpson 悖论的评论中,有人提到了哲学(史)上的休谟问题(我的转述):人类是否能从有限的经验中得到因果律?这的确是一个问题,这个问题最后促使德国哲学家康德为调和英国经验派(休谟)和大陆理性派(莱布尼兹-沃尔夫)而写了巨著《纯粹理性批判》。其实,如果一个人是绝对的怀疑论者(如休谟),他可能怀疑一切,甚至包括因果律,所以,康德的理论也不能完全“解决”休谟问题。怀疑论者是无法反驳的,他们的问题也是无法回答的。他们存在的价值是为现行一切理论起到警示作用。一般来说,统计学家不会从过度哲学的角度谈论问题。从前面的说明中可以看出,统计中所谓的“因果”是“某种”意义的“因果”,即统计学只讨论“原因的结果”,而不讨论“结果的原因”。前者是可以用数据证明或者证伪的;后者是属于科学研究所探索的。用科学哲学家卡尔·波普的话来说,科学知识的积累是“猜想与反驳”的过程:“猜想”结果的原因,再“证伪”原因的结果;如此循环即科学。

继续阅读因果推断简介之四:观察性研究,可忽略性和倾向得分

因果推断简介之三:R. A. Fisher 和 J. Neyman 的分歧

R.A. Fisher
这部分谈到的问题非常微妙:完全随机化试验下的 Fisher randomization test 和 Neyman repeated sampling procedure。简单地说,前者是随机化检验,或者如很多教科书讲的Fisher 精确检验 (Fisher exact test);后者是 Neyman 提出的置信区间 (confidence interval)理论。

我初学因果推断的时候,并没有细致的追求这些微妙的区别,觉得了解到简介之二的层次就够了。不过在 Guido Imbens 和 Donald Rubin 所写的因果推断教科书(还未出版)中,这两点内容放在了全书的开端,作为因果推断的引子。在其他的教科书中,是看不到这样的讲法的。平日里常常听到 Donald Rubin 老爷子对 Fisher randomization test 的推崇,我渐渐地也被他洗脑了。

Fisher 的随机化检验,针对的是如下的零假设,又被称为 sharp null:  $$H_0 : Y_i(1) = Y_i(0), \forall i = 1,\cdots,n.$$ 坦白地说,这个零假设是我见过的最奇怪的零假设,没有之一。现行的统计教科书中,讲到假设检验,零假设都是针对某些参数的,而 Fisher 的 sharp null 看起来却像是针对随机变量的。这里需要讲明白的是,当我们关心有限样本 (finite sample)的因果作用时,每个个体的潜在结果 $ \{Y_i(1), Y_i(0)\} $ 都是固定的,观测变量 $Y_i = Z_i Y_i(1) + (1 – Z_i)Y_i(0) $ 的随机性仅仅由于“随机化” $Z_i $ 本身导致的。理解清楚这点,才能理解 Fisher randomization test 和后面的 Neyman repeated sampling procedure。如果读者对于这种有限样本的思考方式不习惯,可以先阅读一下经典的抽样调查教科书,那里几乎全是有限样本的理论,所有的随机性都来自于随机采样的过程。

如果认为潜在结果是固定的数,那么 Fisher sharp null 就和现行的假设检验理论不相悖。这个 null 之所以“sharp”的原因是,在这个零假设下,所有个体的潜在结果都固定了,个体的因果作用为零,唯一的随机性来自于随机化的“物理”特性。定义处理分配机制的向量为 $$ \overrightarrow{Z} = (Z_1, \cdots, Z_n).$$  结果向量为 $$\overrightarrow{Y} = (Y_1, \cdots, Y_n).$$

此时有限样本下的随机化分配机制如下定义:

$$P(  \overrightarrow{Z} |   \overrightarrow{Y}  ) = \binom{n}{m}^{-1}, \forall  \overrightarrow{Y} ,$$

其中, $m = \sum\limits_{i=1}^n Z_i $ 为处理组中的总数。这里的“条件期望”并不是说  $ \overrightarrow{Y}  $ 是随机变量,而是强调处理的分配机制不依赖于潜在结果。比如,我们选择统计量 $$T = T(\overrightarrow{Z}, \overrightarrow{Y}) = \frac{1}{m} \sum\limits_{i=1}^n Z_i Y_i –  \frac{1}{n-m }\sum\limits_{i=1}^n (1 – Z_i) Y_i $$
J. Neyman来检验零假设,问题在于这个统计量的分布不易求出。但是,我们又知道,这个统计量的分布完全来自随机化。因此,我们可以用如下的“随机化”方法 (Monte Carlo 方法模拟统计量的分布):将处理分配机制的向量 $ \overrightarrow{Z}  $ 进行随机置换得到$ \overrightarrow{Z}^1 = (Z_1^1, \cdots, Z_n^1) $,计算此时的检验统计量 $ T^1 = T(\overrightarrow{Z}^1, \overrightarrow{Y}) $;如此重复多次($n$ 不大时,可以穷尽所有的置换),便可以模拟出统计量在零假设下的分布,计算出 p 值。

有人说,Fisher randomization test 已经蕴含了 bootstrap 的思想,似乎也有一定的道理。不过,这里随机化的方法是针对一个特例提出来的。

下面要介绍的 Neyman 的方法,其实早于 Fisher 的方法。这种方法在 Neyman 1923 年的博士论文中,正式提出了。这种方法假定 $n $ 个个体中有 $m $个随机的接受处理,目的是估计(有限)总体的平均因果作用:$$ \tau = \frac{1}{n} \sum\limits_{i=1}^n \{ Y_i(1) – Y_i(0) \} .$$ 一个显然的无偏估计量是  $$\hat{\tau}  = \bar{y}_1 – \bar{y}_0 =  \frac{1}{m} \sum\limits_{i=1}^n Z_i Y_i –  \frac{1}{n-m} \sum\limits_{i=1}^n (1 – Z_i) Y_i .$$ 但是,通常的方差估计量,

$\hat{\text{Var}}(\hat{\tau})  =  \sum\limits_{Z_i=1} (Y_i – \bar{y}_1)^2 /(m-1)m+  \sum\limits_{Z_i=0} (Y_i – \bar{y}_0)^2/(n-m-1)(n-m) $

高估了方差,构造出来的置信区间在 Neyman – Pearson 意义下太“保守”。可以证明,在个体处理作用是常数的假定下,上面的方差估计是无偏的。

通常的教科书讲假设检验,都是从正态均值的检验开始。Neyman 的方法给出了 $ \tau $ 的点估计和区间估计,也可以用来检验如下的零假设:$$H_0:  \tau = 0.$$

实际中,到底是 Fisher 和零假设合理还是 Neyman 的零假设合理,取决于具体的问题。比如,我们想研究某项政策对于中国三十多个省的影响,这是一个有限样本的问题,因为我们很难想象中国的省是来自某个“超总体”。但是社会科学中的很多问题,我们不光需要回答处理或者政策对于观测到的有限样本的作用,我们更关心这种处理或者政策对于一个更大总体的影响。前者,Fisher 的零假设更合适,后者 Neyman 的零假设更合适。

关于这两种角度的争论,可以上述到 Fisher 和 Neyman 两人。1935 年,Neyman 向英国皇家统计学会提交了一篇论文“Statistical problems in agricultural experimentation”,Fisher 和 Neyman 在讨论文章时发生了激烈的争执。不过,从今天的统计教育来看,Neyman 似乎占了上风。

用下面的问题结束:

  1. 在 sharp null下,Neyman 方法下构造的 T 统计量,是否和 Fisher randomization test 构造的统计量相同?分布是否相同?
  2. Fisher randomization test 中的统计量可以有其他选择,比如 Wilcoxon 秩和统计量等,推断的方法类似。
  3. 当 $ Y $ 是二值变量时,上面 Fisher 的方法就是教科书中的 Fisher exact test。在没有学习 potential outcome 这套语言之前,理解 Fisher exact test 是有些困难的。
  4. 证明 $ E \{  \hat{\text{Var}}(\hat{\tau})  \}  \geq \text{Var}(\hat{\tau})  $。
  5. 假定 $n $ 个个体是一个超总体(super-population)的随机样本,超总体的平均因果作用定义为 $$ \tau_{SP} = E\{ Y(1) – Y(0) \}.$$ 那么 Neyman 的方法得到估计量是超总体平均因果作用的无偏估计,且方差的表达式是精确的;而 sharp null 在超总体的情形下不太适合。
原始的参考文献是:
  • Neyman, J. (1923) On the application of probability theory to agricultural experiments. Essay on principles. Section 9. reprint in Statistical Science. 5, 465-472. with discussion by Donald Rubin.

最近的理论讨论是:

  • Ding, P. (2014). A Paradox from Randomization-Based Causal Inference. On Arxiv.

 

因果推断简介之二:Rubin Causal Model (RCM)和随机化试验

因果推断用的最多的模型是 Rubin Causal Model (RCM; Rubin 1978) 和 Causal Diagram (Pearl 1995)。Pearl (2000) 中介绍了这两个模型的等价性,但是就应用来看,RCM 更加精确,而 Causal Diagram 更加直观,后者深受计算机专家们的推崇。这部分主要讲 RCM。

设  $Z_i$ 表示个体 $i$ 接受处理与否,处理取 $1$,对照取$0$ (这部分的处理变量都讨论二值的,多值的可以做相应的推广);$Y_i$ 表示个体 $i$ 的结果变量。另外记 $ \{  Y_i(1),Y_i(0)\}  $ 表示个体 $i$ 接受处理或者对照的潜在结果 (potential outcome),那么 $Y_i(1) -Y_i(0)$  表示个体 $i$  接受治疗的个体因果作用。不幸的是,每个个体要么接受处理,要么接受对照 $\{Y_i(1),Y_i(0)\} $ 中必然缺失一半,个体的因果作用是不可识别的。观测的结果是 $ Y_i = Z_i Y_i(1) + (1 – Z_i) Y_i(0)  $。 但是,在 $Z$ 做随机化的前提下,我们可以识别总体的平均因果作用 (Average Causal Effect; ACE):

$$
ACE(Z\rightarrow Y) = E\{Y_i(1) – Y_i(0)\}
$$

这是因为
$$
\begin{eqnarray*}
ACE(Z \rightarrow Y) & = & E\{Y_i(1)\} -E\{Y_i(0)\} \\
& = & E\{Y_i(1) \mid Z_i =1\}  -E\{Y_i(0)\mid Z_i=0\} \\
& = & E\{Y_i \mid Z_i =1\} – E\{Y_i \mid Z_i=0\}
\end{eqnarray*}
$$
最后一个等式表明 $ACE$ 可以由观测的数据估计出来。其中第一个等式用到了期望算子的线性性(非线性的算子导出的因果度量很难被识别!);第二个式子用到了随机化,即 $$Z\bot\{Y(1),Y(0)\} $$ 其中, $\bot$ 表示独立性。由此可见,随机化试验对于平均因果作用的识别起着至关重要的作用。

当 $ Y $ 是二值的时候,平均因果作用是流行病学中常用的“风险差”(risk difference; RD):
$$
\begin{eqnarray*}
CRD(Z\rightarrow Y)  & = & P(Y(1) = 1)  –  P(Y(0)=1) \\
& = & P(Y=1\mid Z=1) – P(Y=1\mid Z=0)
\end{eqnarray*}
$$

当然,流行病学还常用“风险比”(risk ratio; RR):
$$
\begin{eqnarray*}
CRR(Z \rightarrow Y) & = & \frac{P(Y(1) = 1)}{P(Y(0)=1)}\\
& = & \frac{P(Y=1\mid Z=1)}{P(Y=1\mid Z=0)}
\end{eqnarray*}
$$

和“优势比”(odds ratio; OR):
$$
\begin{eqnarray*}
COR(Z \rightarrow Y) & = & \frac{P(Y(1) = 1)P(Y(0)=0) }{P(Y(0)=1)P(Y(1)=0) } \\
& = & \frac{P(Y=1\mid Z=1)P(Y=0\mid Z=0)}{P(Y=1\mid Z=0) P(Y=0\mid Z=1)}
\end{eqnarray*}
$$

上面的记号都带着“C”,是为了强调“causal”。细心的读者会发现,定义 CRR 和 COR 的出发点和 ACE 不太一样。ACE 是通过对个体因果作用求期望得到的,但是 CRR 和 COR 是直接在总体上定义的。这点微妙的区别还引起了不少人的研究兴趣。比如,经济学中的某些问题,受到经济理论的启示,处理的作用可能是非常数的,仅仅研究平均因果作用不能满足实际问题的需要。这时候,计量经济学家提出了“分位数处理作用”(quantile treatment effect: QTE):
$$
QTE(\tau) = F^{-1}_{Y(1)}(\tau) – F^{-1}_{Y(0)}(\tau)
$$

在随机化下,这个量也是可以识别的。但是,其实这个量并不能回答处理作用异质性(heterogenous treatment effects)的问题,因为处理作用非常数,最好用如下的量刻画:
$$
\Delta(\delta) = P(Y(1) – Y(0) \leq \delta )
$$

这个量刻画的是处理作用的分布。不幸的是,估计 $\Delta(\delta) $ 需要非常强的假定,通常不具有可行性。

作为结束,留下如下的问题:

  1. “可识别性”(identifiability)在统计中是怎么定义的?
  2. 医学研究者通常认为,随机对照试验(randomized controlled experiment)是研究处理有效性的黄金标准,原因是什么呢?随机化试验为什么能够消除 Yule-Simpson 悖论?
  3. $QTE(\tau) $ 在随机化下是可识别的。另外一个和它“对偶”的量是 Ju and Geng (2010) 提出的分布因果作用(distributional causal effect: DCE):$  DCE(y) = P(Y(1) \geq y) – P(Y(0) \geq y) $ ,在随机化下也可以识别。
  4. 即使完全随机化,$\Delta(\delta)$ 也不可识别。也就是说,经济学家提出的具有“经济学意义”的量,很难用观测数据来估计。这种现象在实际中常常发生:关心实际问题的人向统计学家索取的太多,而他们提供的数据又很有限。

*关于 RCM 的版权,需要做一些说明。目前可以看到的文献,最早的是 Jerzy Neyman 于 1923 年用波兰语写的博士论文,第一个在试验设计中提出了“潜在结果”(potential outcome)的概念。后来 Donald Rubin 在观察性研究中重新(独立地)提出了这个概念,并进行了广泛的研究。Donald Rubin 早期的文章并没有引用 Jerzy Neyman 的文章,Jerzy Neyman 的文章也不为人所知。一直到 1990 年,D. M. Dabrowska 和 T. P. Speed 将 Jerzy Neyman 的文章翻译成英文发表在 Statistical Science 上,大家才知道 Jerzy Neyman 早期的重要贡献。今天的文献中,有人称 Neyman-Rubin Model,其实就是潜在结果模型。计量经济学家,如 James Heckman 称,经济学中的 Roy Model 是潜在结果模型的更早提出者。在 Donald Rubin 2004 年的 Fisher Lecture 中,他非常不满地批评计量经济学家,因为 Roy 最早的论文中,全文没有一个数学符号,确实没有明确的提出这个模型。详情请见,Donald Rubin 的 Fisher Lecture,发表在 2005 年的 Journal of the American Statistical Association 上。研究 Causal Diagram 的学者,大多比较认可 Donald Rubin 的贡献。但是 Donald Rubin 却是 Causal Diagram 的坚定反对者,他认为 Causal Diagram 具有误导性,且没有他的模型清楚。他与James Heckman (诺贝尔经济学奖), Judea Pearl (图灵奖) 和 James Robins 之间的激烈争论,成为了广为流传的趣闻。

参考文献:

  1. Neyman, J. (1923) On the application of probability theory to agricultural experiments. Essay on principles. Section 9. reprint in Statistical Science. 5, 465-472.
  2. Pearl, J. (1995) Causal diagrams for empirical research. Biometrika, 82, 669-688.
  3. Pearl, J. (2000) Causality: models, reasoning, and inference. Cambridge University Press。
  4. Rubin, D.B. (1978) Bayesian inference for causal effects: The role of randomization. The Annals of Statistics, 6, 34-58.

因果推断简介之一:从 Yule-Simpson’s Paradox 讲起

在国内的时候,向别人介绍自己是研究因果推断(causal inference)的,多半的反应是:什么?统计还能研究因果?这确实是一个问题:统计研究因果,能,还是不能?直接给出回答,比较冒险;如果有可能,我需要花一些篇幅来阐述这个问题。

目前市面上能够买到的相关教科书仅有 2011 年图灵奖得主 Judea Pearl 的 Causality: Models, Reasoning, and Inference。Harvard 的统计学家 Donald Rubin 和 计量经济学家 Guido Imbens 合著的教科书历时多年仍尚未完成;Harvard 的流行病学家 James Robins 和他的同事也在写一本因果推断的教科书,本书目前只完成了第一部分,还未出版(见此处)。我本人学习因果推断是从 Judea Pearl 的教科书入手的,不过这本书晦涩难懂,实在不适合作为入门的教科书。Donald Rubin 对 Judea Pearl 提出的因果图模型(causal diagram)非常反对,他的教科书中杜绝使用因果图模型。我本人虽然脑中习惯用图模型进行思考,但是还是更偏好 Donald Rubin 的风格,因为这对于入门者,可能更容易。不过这一节,先从一个例子出发,不引进新的统计符号和概念。

天才的高斯在研究天文学时,首次引进了最大似然和最小二乘的思想,并且导出了正态分布(或称高斯分布)。其中最大似然有些争议,比如 Arthur Dempster 教授说,其实高斯那里的似然,有贝叶斯或者信仰推断(fiducial inference)的成分。高斯那里的“统计”是关于“误差”的理论,因为他研究的对象是“物理模型”加“随机误差”。大约在 100 多年前,Francis Galton 研究了父母身高和子女身高的“关系”,提出了“(向均值)回归”的概念。众所周知,他用的是线性回归模型。此时的模型不再是严格意义的“物理模型”,而是“统计模型” — 用于刻画变量之间的关系,而不一定是物理机制。之后,Karl Pearson 提出了“相关系数”(correlation coefficient)。后世研究的统计,大多是关于“相关关系”的理论。但是关于“因果关系”的统计理论,非常稀少。据 Judea Pearl 说,Karl Pearson 明确的反对用统计研究因果关系;有意思的是,后来因果推断为数不多的重要文章(如 Rosenbaum and Rubin 1983; Pearl 1995)都发表在由 Karl Pearson 创刊的 Biometrika 上。下面讲到的悖论,可以说是困扰统计的根本问题,我学习因果推断便是由此入门的。

在高维列联表分析中, 有一个很有名的例子,叫做 Yule-Simpson’s Paradox。有文献称,Karl Pearson 很早就发现了这个悖论 — 也许这正是他反对统计因果推断的原因。此悖论表明,存在如下的可能性:$X$ 和 $Y$ 在边缘上正相关;但是给定另外一个变量 $Z$ 后,在 $Z$ 的每一个水平上,$X$ 和 $Y$ 都负相关。Table 1 是一个数值的例子,取自 Pearl (2000)。

simpson

Table 1 中,第一个表是整个人群的数据:接受处理和对照的人都是 40 人,处理有较高的存活率,因此处理对整个人群有“正作用”。第二个表和第三个表是将整个人群用性别分层得到的,因为第一个表的四个格子数,分别是下面两个表对应格子数的和: $$20 = 18+2, 20 = 12+8, 16 = 7+9, 24 = 3+21.$$ 奇怪的是,处理对男性有“负作用”,对女性也有“负作用”。一个处理对男性和女性都有“负作用”,但是他对整个人群却有“正作用”:悖论产生了!

有人可能会认为这种现象是由于随机性或者小样本的误差导致的。但是这个现象与样本量无关,与统计的误差也无关。比如,将上面的每个格子数乘以一个巨大的正数,上面的悖论依然存在。

纯数学的角度,上面的悖论可以写成初等数学 $ \frac{a}{b} < \frac{c}{d}, \frac{a’}{b’} < \frac{c’}{d’} , \frac{a  + a’}{b + b’} > \frac{c + c’}{d + d’} $;这并无新奇之处。但是在统计上,这具有重要的意义 — 变量之间的相关关系可以完全的被第三个变量“扭曲”。更严重的问题是,我们的收集的数据可能存在局限性,忽略潜在的“第三个变量”可能改变已有的结论,而我们常常却一无所知。鉴于 Yule-Simpson 悖论的潜在可能,不少人认为,统计不可能用来研究因果关系。

上面的例子是人工构造的,在现实中,也存在不少的实例正是 Yule-Simpson’s Paradox。比如,UC Berkeley 的著名统计学家 Peter Bickel 教授 1975 年在 Science 上发表文章,报告了 Berkeley 研究生院男女录取率的差异。他发现,总体上,男性的录取率高于女性,然而按照专业分层后,女性的录取率却高于男性 (Bickel 等 1975)。在流行病学的教科书 (如 Rothman 等 2008) 中,都会讲到“混杂偏倚”(confounding bias),其实就是 Yule-Simpson’s Paradox,书中列举了很多流行病学的实际例子。

由于有 Yule-Simpson’s Paradox 的存在,观察性研究中很难得到有关因果的结论,除非加上很强的假定,这在后面会谈到。比如,一个很经典的问题:吸烟是否导致肺癌?由于我们不可能对人群是否吸烟做随机化试验,我们得到的数据都是观察性的数据:即吸烟和肺癌之间的相关性 (正如 Table 1 的合并表)。此时,即使我们得到了吸烟与肺癌正相关,也不能断言“吸烟导致肺癌”。这是因为可能存在一些未观测的因素,他既影响个体是否吸烟,同时影响个体是否得癌症。比如,某些基因可能使得人更容易吸烟,同时容易得肺癌;存在这样基因的人不吸烟,也同样得肺癌。此时,吸烟和肺癌之间相关,却没有因果作用。

相反的,我们知道放射性物质对人体的健康有很大的伤害,但是铀矿的工人平均寿命却不比常人短;这是流行病学中有名的“健康工人效应”(healthy worker effect)。这样一来,似乎是说铀矿工作对健康没有影响。但是,事实上,铀矿的工人通常都是身强力壮的人,不在铀矿工作寿命会更长。此时,在铀矿工作与否与寿命不相关,但是放射性物质对人的健康是有因果作用的。

这里举了一个悖论,但没有深入的阐释原因。阐释清楚这个问题的根本原因,其实就讲清楚了什么是因果推断。这在后面会讲到。作为结束,留下如下思考的问题:

  1. Table 1 中,处理组和对照组中,男性的比例分别为多少?这对悖论的产生有什么样的影响?反过来考虑处理的“分配机制”(assignment mechanism),计算 $ P(\text{Treatment} \mid \text{Male})$ 和 $P(\text{Treatment} \mid \text{Female})  $。
  2. 假如 $ (X, Y, Z) $ 服从三元正态分布,$ X $ 和 $ Y$  正相关,$ Y$ 和 $ Z$ 正相关,那么 $ X $ 和 $ Z$ 是否正相关?(北京大学概率统计系 09 年《应用多元统计分析》期末第一题)
  3. 流行病学的教科书常常会讲各种悖论,比如混杂偏倚 (confounding bias)和入院率偏倚(Berkson’s bias)等,本质上是否与因果推断有关?
  4. 计量经济学中的“内生性”(endogeneity)怎么定义?它和 Yule-Simpson 悖论有什么联系?

参考文献:

  1. Bickel, P. J. and Hammel, E. A. and O’Connell, J. W. (1975) Sex bias in graduate admissions: Data from Berkeley. Science, 187, 398-404.
  2. Pearl, J. (2000) Causality: models, reasoning, and inference. Cambridge University Press。
  3. Rosenbaum, P.R. and Rubin, D.B. (1983) The central role of the propensity score in observational studies for causal effects. Biometrika, 70, 41-55.
  4. Rothman, K., Greenland, S. and Lash, T. L.  (2008) Modern Epidemiology. Lippincott Williams & Wilkins.

Diovan

真理在缩水,还是上帝在掷骰子?

最近在Google Reader中看见科学松鼠会有两篇文章被频繁分享,名为《真理在缩水——现代科学研究方法并不尽善尽美?》()与(),下文简称《缩水》。文章很有意思,而实际上说的是我们的老本行——统计学,因此我在这里也发表一些我的想法和理解,包括这两年我在美帝学习的一些思考,部分内容受益于两位老师Kaiser和Nettleton教授,先向他们致谢(尽管他们永远都不会看到这篇文章)。同时我也要先说明一下,读这篇文章可能会很花时间(至少我花了大约二十小时写这篇文章),即使我的观点没有价值,我相信里面的引用文献是有价值的。

初读文章,我脑子里冒出的一句话是“上帝在跟我们掷骰子”,文中给出了大量的不可重复的试验,仿佛就像那些号称“具有统计学意义”(下文我再说这个所谓的“意义”)的试验结果只是若干次骰子中的一次运气好的结果而已。读完文章,我们可能不禁要问,到底是真理在缩水,还是它根本就不曾存在?下面我从四个方面来展开,分别说明人对随机性的认识、统计推断的基石、让无数英雄折腰的P值、以及可重复的统计研究。

一、感知随机

随机变量在统计分析中占据中心地位,数学上关于随机变量的定义只是一个“干巴巴的函数”,从样本空间映射到实数集,保证从实数集上的Borel域逆回去的集合仍然在原来的sigma域中即可。随机变量的性质由其分布函数刻画。写这段话的目的不是为了吓唬你,也不是为了作八股文,而是来说明我为什么不喜欢数学的理由,对我而言,我觉得有些数学工具只是为了让自己不要太心虚,相信某时某刻某个角落有个理论在支撑你,但后果是弱化了人的感知,当然,也有很多数学工具有很强的直觉性(如果可能,我想在未来下一篇文章里面总结这些问题)。我一直认为很多人对随机性的感知是有偏差的,对概率的解释也容易掉进陷阱(参见Casella & Berger的Statistical Inference第一章,例如条件概率的三囚徒问题)。

《缩水》一文发表了很多不可重复的试验案例,我们应该吃惊吗?我的回答是,未必。举两个简单的例子:

第一个例子:很多数据分析人员都很在意所谓的“离群点”,论坛上也隔三差五有人问到这样的问题(如何判断和处理离群点),而且也有很多人的做法就是粗暴地删掉,我从来都反对这种做法。除了基于“数据是宝贵的”这样简单的想法之外,另一点原因是,离群点也许并非“异类”。离群点是否真的不容易出现?请打开R或其它统计软件,生成30个标准正态分布N(0, 1)随机数看看结果,比如R中输入rnorm(30),这是我运行一次的结果:

> rnorm(30)
 [1]  1.19062761 -0.85917341  2.90110515  0.59532402 -0.05081508 -0.06814796
 [7]  2.08899701  0.76423007  0.92587075 -1.16232929 -0.68074378 -1.40437532
[13] -0.17932604 -0.72980545 -0.53850923  0.21685537 -0.35650714 -1.32591808
[19] -0.88071526 -1.25832441  0.24001498 -0.41682799 -0.09576492 -0.17059052
[25] -0.99947485  0.25108253 -0.47566842 -0.28028786  0.79856649 -0.13250974

30在现实中是一个比较小的样本量,你看到了什么?2.901?它接近3倍标准差的位置了。还有2.089?……如果你不知道这批数据真的是从标准正态分布中生成出来的,现在你会有什么反应?把2.9删掉?标准正态分布是一个在我们眼中很“正常”的分布,而一个不太大的样本量一次试验足以生成几个“离群点”,那么要是成千上万的试验中没能产生几项震惊世界的结果,你会怎样想?(上帝的骰子坏掉了)

另一个例子和统计学结合紧密一点,我们谈谈残差的QQ图。QQ图是用来检查数据正态性的一种统计图形,与腾讯无关,细节此处略去,大意是图中的点若呈直线状(大致分布在对角线上),那么可以说明数据的正态性比较好,因此QQ图经常被用在对回归模型残差的正态性诊断上。我的问题是,即使数据真的是正态分布,你是否真的会看见一些分布在直线上的点?若答案是否定的,那么我们就得重新审视我们对分布和随机的认识了。下图是一幅教科书式的QQ图(仍然基于30个正态分布随机数):

“正常的”QQ图
“正常的”QQ图(来自R语言qqnorm(rnorm(30)))

随机性并没有这么美好,即使数据真的来自正态分布,你也有可能很容易观察到歪歪扭扭的QQ图,尤其是小样本的情况下。比如下图是50次重复抽样的正态数据QQ图,它和你想象的QQ图本来的样子差多远?

library(animation)
set.seed(710)
ani.options(interval = 0.1, nmax = 50)
par(mar = c(3, 3, 2, 0.5), mgp = c(1.5, 0.5, 0), tcl = -0.3)
sim.qqnorm(n = 30, pch = 19, col = "red", last.plot = expression(abline(0, 1)))
真实的正态分布QQ图
真实的正态分布QQ图(图中直线为y = x)

正态分布是统计学中比较“正常”的一类分布(台湾学者甚至译为“常态分布”),我们尚不能很好感知它的随机性,就不必说其它“怪异”的分布了。

这是我想说的第一点,作为人类,我们对上帝的骰子至少在感知上就可能有问题(别误会,我不信教),接下来从理性角度分析一下。

二、统计推断

《缩水》一文中提到的基本上都是统计推断方法带来的结果,为了理解这些结果,我们必须三思统计推断本身的过程。一般说来,统计推断有两种途径:随机试验和(概率)统计模型,或者说基于试验的推断和基于模型的推断。前者的代表性方法为置换检验(Permutation test),不过它似乎被大多数人遗忘了,更多的人拿到数据首先想的就是用哪个统计模型和分布;实际上,置换检验是一种极具代表性的统计推理方法,可以用典型的“三段论”来说明它(参见去年江堂的文章):

  1. 要么A,要么B
  2. 若有A,则有C
  3. 若非C,则非A,于是B

置换检验的场景是,一次试验中,我们为试验单元随机分配不同的处理(treatment),为了简单起见,假设这里只有两种处理水平A和B,我们想知道两种处理在试验单元的某项指标上是否有显著差异。逻辑是这样:假设处理毫无效果,那么某一个试验对象的指标将不受处理影响,不管我们给老鼠嗑的是A药还是B药,结果都一样,那么我们可以把处理的标签随机打乱(某些A、B随机互换),打乱后A组和B组的差异不应该会和原来的差异很不一样(因为药不管用嘛),否则,我们恐怕得说,药还是管用的。就这样,我们随机打乱标签很多次,形成一个“人工生成”的AB差异分布(因为我们生成了很多个差异值),看原来的差异在这个分布的什么位置上,如果在很靠近尾巴的位置上,那么就认为P值很小。当了个当,当了个当,P值出场了。对置换检验熟悉的人是否有想过,好像我们一直没谈到什么分布的假设,那这个概率值(P值)是从哪里生出来的?答案是最初的“随机分配处理到试验单元上”。这就涉及到试验设计的一大原则:随机化。为什么要随机化?因为试验单元之间可能本来就有差异,换句话说,如果你不随机化,那么你观察到的差异有可能来自试验单元本身。比如,你从笼子里抓前10只老鼠给嗑A药,后10只老鼠给B药,这就没有随机化,前10只老鼠可能很笨或是老弱病残,容易被你抓住,而后10只老鼠跑得贼快。随机化将这些个体差异转变为了随机的误差,例如两只老鼠之间的确本身有差异,但我要是把它们随机分配给处理,那么这种个体差异就会随机进入处理组,保证统计推断有根基。如果这一点没有理解透彻,试验人员很容易在数据收集阶段就已经收集了错误的数据。《缩水》一文中的试验都是怎么做的,我没空去细究。

基于模型的推断的一大特征就是开始对随机变量做一些分布上的假设,例如两样本t检验,假设样本来自独立同方差的正态分布。仔细想想这里面的问题,对建模和理解模型结果有至关重要的作用。一个最直接的问题就是,假设条件是否可靠?George Box大人很久很久以前就说了一句被引用了无数遍的话:所有的模型都是错的,但有些是有用的。统计学方法很“滑”(用麦兜的话说),它的科学与艺术双重身份,总让人感觉拿捏不准。学数学的人可能觉得艺术太不靠谱,其它外行人士可能觉得科学太神秘。这个问题我不想作答,也无法作答,搁在一边,先说一些我曾经考虑过的问题,它们与《缩水》一文可能并没有直接联系,但应该能或多或少启发我们从源头考虑统计模型,直接上手统计模型的隐患就在于你认为一切都理所当然,随着时间推移,假设渐渐变成了“公理”和“常识”,我觉得这是很可怕的。

第一个问题是似然函数(likelihood function),它是频率学派的命脉,而我们大多数人仍然都是频率学派的“教徒”。对于离散变量来说,基于似然函数的方法如极大似然估计很容易理解:我们要找一个参数值,使得观察到的数据发生的概率最大。这里的“概率”二字应该被重重划上记号!接下来我要提一个让我曾经觉得后背发凉的问题:

为什么对连续变量来说,似然函数是密度函数的乘积?

你是否想过这个问题?我们知道连续变量取任何一个值的概率都是0,也就是说,如果按照概率的方式解释似然函数,那么连续变量的似然函数应该是0才对,换句话说,你观察到的数据发生的概率是0。现在你觉得似然函数还是一个理所当然的统计工具吗?

有一位统计学家叫J. K. Lindsey,1998年在(英国)皇家统计学会宣读了一篇论文,叫Some statistical heresies(一些统计异端邪说),如果你没见过统计学家打仗,那么这篇论文将会让你看到一场超大规模的战争,后面的讨论者中包括Murray Aitkin、D. R. Cox和J. A. Nelder等老江湖。Lindsey的文章几乎是被大家围攻了(期待我这篇文章也能被围攻),不过他对于似然函数的解释倒是让我有点茅塞顿开。细节我不展开,大意是,似然函数也是一种近似(用积分中值定理去想)。

第二个问题是渐近统计(asymptotic statistics),同样,这也是统计学家的日常工具之一,因为太常见,所以也有一种理所当然的味道。比如我们看到列联表就想到卡方检验(检验行列变量的独立性),殊不知卡方检验只是一种大样本下的近似方法。渐近统计的基石是什么?如果你要挖,我想挖到头你一定会挖到泰勒展开。至少目前我认为渐近统计“基本上就是泰勒展开的第二项(或少数情况下第三项)”。理论上泰勒展开有无穷多项,我们往往根据一些假设条件,把后面的高阶项都消灭掉,剩下一次项或二次项。比如你展开似然函数,就得到了似然比检验的卡方分布;你展开极大似然估计的一个连续函数,你就有了Delta方法(当然,需要依分布收敛的前提);就这样左展右展,展出了中心极限定理(对特征函数或母函数展开),展出了拉普拉斯近似(对对数密度函数展开)。如果你能看到这一点,就不会奇怪为什么正态分布在统计推断中有如此中心地位(谁叫正态分布的对数密度函数是个二次函数呢)。

第三个问题是,贝叶斯方法又如何?既然频率学派中几乎处处是近似,贝叶斯学派是否会好一些?我的回答是好不到哪儿去。贝叶斯的逻辑很简单,但由于灵活性太强,应用时非常摧残人的脑力,导致争议不断(先验分布怎么取、MCMC是否收敛等)。在《缩水》一文中,恐怕是没有基于贝叶斯方法的试验,因为在所谓的科学试验中,人们往往排斥“先验”这种想法,就好像先验一定代表主观、而客观一定正确一样。逻辑上,这是荒谬的。关于这些问题的重磅讨论,可参考Efron去年发表的The Future of Indirect Evidence以及文章后面Gelman他们三个人的讨论。我总感觉这不是我这个年龄应该看的东西,太哲学了。

我提到这些问题,本意是给统计学家看的,如果你是一个合格的统计学家,你自己首先应该意识到统计学的局限性,而不是拿到数据就分析。

三、万能的P值?

早些年当我还是个无知轻狂小子的时候,我曾戏称“统计软件就是为了放个P”,这里的P指的是P值,不是粗话。这话好像也不全然轻狂无知。使用统计方法的人,难道不是在追逐一个小于0.05的P值吗?如果你的结果不显著,那么肯定发表不了。换句话说,发表的结果很有可能是我们在自欺欺人。下面的漫画生动刻画了人们寻找P值的过程(来自xkcd):

Significant
Significant

重大科学发现!吃绿色的软糖会让你长痘痘!95%置信度!

当你看到95%的时候,你看不到红色的、灰色的、褐色的、橙色的、黄色的软糖……这便是《缩水》一文中说的“发表偏见”(publication bias,“偏见”翻译似乎不妥),即发表的结果是经过人工选择的,你已经不能用正常的概率意义去解读它,或者说,概率已经变了样。

插一句“统计学意义”:这个概念本来的英文是statistical significance,但是被很多专业的人翻译为“统计学意义”,我一直认为这很不妥,给人一种错觉,仿佛统计学保证了什么东西有意义一样,我提倡的译法是“统计显著性”。尤其是“由于P值小于0.05,所以具有统计学意义”这种话,我觉得应该见一句删一句。

上面的软糖问题涉及到传统的多重比较(multiple comparison)与P值调整,这里“传统”指的是要控制族错误率(Familywise Error Rate,下文简称FWER)。所谓控制FWER,也就是要使得一族(多个)统计检验中,一个第一类错误都不犯的概率控制在一定水平之下,比如0.05。让多个检验全都不犯错和单个检验不犯错(指第一类错误)显然是有区别的,比如假设所有的检验都是独立的,一个检验不犯错的概率是95%,两个都不犯错的概率就变成了95% * 95% = 90.25%,检验越多,不犯错的概率就越小。把整体的第一类错误率控制在某个alpha值之下,就意味着单个检验必须更“严格”,因此我们不能再以0.05去衡量每个检验是否显著,而要以更小的值去衡量,这就是P值的调整,老办法有Bonferroni等方法。

控制FWER显得有些苛刻,比如有10000个基因都需要检验在不同处理下的表达差异,那么要是按照传统的P值调整方法,恐怕是很难得出某个基因显著的结论(因为alpha值会被调得极其小)。FWER的目标是一个错误都不能犯,但实际上我们也许可以容忍在那些我们宣称显著的结果中,有少数其实是犯错的,就看你是不是“宁愿错杀三千,也不放过一个”了。

Efron在前面我提到的文章中把Benjamini和Hochberg在1995年的论文称为“二战以来统计界第二伟大的成果”(他把第一名给了James & Stein的有偏估计),那么B&H做了什么?答案就是虚假发现率(False Discovery Rate,下文简称FDR)。FDR要控制的是在宣称显著的结论中错误的结论的比例,比如10000个基因中有100个基因显著,但实际上有5个是虚假的发现(即本来这5个基因分别在两种处理下的表达并没有差异)。尽管有错误存在,但我们认为可以承受。

初学者到这里应该可以意识到了,通过FDR选出来的结果在理论上就已经存在错误了,当然这只是小问题,更大的问题在于,FDR的定义实际上是一个期望的形式:真实的零假设个数除以被拒绝的零假设个数的期望(零假设是没有差异)。凡是涉及到期望的东西,我们都应该冷静三秒想一下,期望意味着什么?

假设有一个游戏,你获胜的概率是70%,要是赢了,你得到一百万,要是输了,你付出一百万;获利的期望是40万。现在我请你去玩,一把定输赢,你玩不玩?我相信大多数人不会玩(除非你实在太有钱了),为什么期望是赚40万你也不玩?因为期望往往是“样本量无穷大”你才能隐约看到的东西,玩一次游戏,输掉一百万的概率是30%,恐怕你不敢。

FDR是个期望,也就是你做很多次试验,平均来看,FDR在某个数值附近。一次试验中它究竟在哪儿,谁都不知道。就像(频率学派的)置信区间一样,我给你一个区间,其实你并不知道参数在不在区间里,你也无法用概率的方式去描述单个区间,比如你不能说“参数落在这个区间内的概率是95%”(只有无数次抽样重新计算区间,这无数个区间覆盖真实参数的概率才是95%)。

所以,某种意义上,概率论源于赌博,而统计学在骨子里一直都是赌博。旁观者看赌徒,总觉得他在赚钱。当然,统计学家是“高级赌徒”,他们不是随机乱赌。

四、可重复的统计研究

看到这里,你大概也脑子有点抽筋了(如果你把我提到的Lindsey和Efron的文章都看过了的话),我们换个轻松点的话题:可重复的统计研究。这不是我第一次提它,我们一直都在号召可重复的统计研究(如《Sweave:打造一个可重复的统计研究流程》)。还是老话一句,不谈道德问题,因为这是不可控的因素,我们可控的只有制度和工具。源代码不会撒谎。

我们期待学术研究过程的透明化,至少统计之都在努力。