所有由谢益辉发布的文章

关于谢益辉

RStudio码了个工,Iowa State University统计系博了个士。统计之都网站创办者;研究兴趣为统计图形及数据可视化,对统计模型方法的发展感兴趣但不喜欢纯粹抽象的数学理论,以直观、实用为学习标准;偏好以R语言为工具;Email:[email protected];个人主页:http://yihui.name

用交互式图形探索一个五百年前的脑洞

按惯例先跑几段火车,赶时间的请直接从下面油画处开读。我很少看电影,欠的稿子都写不完还看毛线电影,不过前段时间《大鱼海棠》的精美海报画面还是吸引了我的注意力(又是从涛妹的票圈看到的),深为赞叹现在国内的动画制作技术。然而过了几天,好像评论的风向就变了。可惜了情怀这个词,现在也成了为人不齿的陈词滥调了:情怀,情你个锤子的怀,你才情怀,你全家都情怀。遥想当年,萌主(周扬)在明德楼地下咖啡厅的小房间里给我们展示 R/ECharts/Shiny 的时候,第一次提到情怀一词,小板凳上的我们都感受到了内心的一团火。“厉害啊!”萌主洋洋自得。

据说《大鱼海棠》可惜在用了辣么精良的画面,却愣是没讲好一个故事(重申一遍:我没看,只是据说);相比之下,人家徐克老爷子二十年前用简陋的技术却做出动画片《小倩》,同样是用中国传统故事素材,但比《大鱼海棠》不知道高到哪里去了。

“要来找我哦,我家就在北村,门口有棵好大好大的桃花树,记得一定要来找我哦!”

那我们就来谈谈讲故事的事。在统计之都十周年感言中我曾经提到过“精致的脑洞”,今天我就给大家解说一个五百年前的脑洞,这个脑洞是我压箱底的货,一般人我不告诉他,我第一次讲它是在我的博士论文答辩会上,后来便很少提它了。讲这个洞有两个目的,一是谈谈我对讲故事本身的一些想法(讲故事本不是我擅长的,但这个洞很适合讲故事),二是演示一下交互式图形的基本概念。 继续阅读用交互式图形探索一个五百年前的脑洞

统计之都十周年感言

呐,统计之都已经创建十周年。作为所谓的创始人,自然也是时候卷起袖子跟大家一起干一大碗鸡汤,毕竟十年这个时间长度听起来好像还蛮厉害的。不巧最近这些天挺忙,加上我其实并不太喜欢专门写文字给乌泱乌泱的客官们看,年龄越大,就越不愿意去安利别人。这篇十周年感言如何写,每天晚上苦苦思考三分钟之后就睡得特别香。想全面概括这十年的发展是不可能的,只能用我最擅长的意识流方式想到哪儿说到哪儿了,这篇文章基本上是纯个人视角,无意借机强行输出价值观。10th

想当初创建统计之都那会儿,朕是天不怕地不怕,心高气傲,不懂就放狗搜了若无其事地回来装懂,在满论坛的点赞声中深藏功与名。有时候在论坛回帖回到半夜一两点,就是那种“扶朕起来朕还能回”的感觉。时间长了,偶尔会有不知情的热心网友称谢教授,嗨别介,我就一逗比本科生而已,然而心里自然窃喜不已,三条腿的蛤蟆好找,二十二岁的教授不常见啊。现在想想人不轻狂枉少年,虚荣心也不是什么坏事,我的 R 语言技能,大致就是在三天两头帮人看函数文档中学出来的。话说那时候 R 帮助文档还是 CHM 格式的,现在的娃估计都不知道咩是 CHM 了。这话题切换到 R 也忒快了。 继续阅读统计之都十周年感言

数据科学家的崛起

美国2012总统大选是奥巴马的胜利,但实际上也是统计学家的胜利。奥巴马当选之夜,我看见推特上有一条消息被疯狂转载:

NATE SILVER ELECTED 44TH PRESIDENT OF UNITED STATES

当然这是一句玩笑话,但Nate Silver是谁?他号称“竞选预测之神谕”:2008年的总统大选他预测对了最终结果,而且美国50州的投票结果他预测对了49个;今年的大选他又预测对了,并且是50州全对。Silver是一名统计学家,毕业于芝加哥大学,随后在毕马威会计师事务所“度过了令自己后悔的四年时间”(不喜欢那里的工作),后来转向预测棒球选手的成绩,再后来转向政治方面的数据分析和预测。总统大选的预测是一件噪声很大的工作,各家有各家的预测和分析,各种突发事件可能会导致某位候选人的支持短期内大幅变动。Silver的工作就像机器学习中的“集成学习”(他自己的描述是“贝叶斯统计”,用自己的先验信息和数据得到后验),集合众多民意调查结果,根据自己的经验判断去平均它们(具体过程我不清楚)。 继续阅读数据科学家的崛起

knitr与可重复的统计研究(花絮篇)

2010年年底我写了章,关于Sweave/LyX/pgfSweave,顺便引出可重复研究(Reproducible Research)的概念。一年过后,我逐渐意识到这一系列基于Sweave的工具都有致命的设计缺陷,束缚感越来越强,屡屡冒出要重复造轮子的想法。于是就在“造乎?不造乎?”的犹豫中最终痛下决心全盘重造,knitr包就诞生了。在第五届中国R语言会议上魏太云已经对它作了初步介绍,我会在统计之都以系列文章全面介绍它,本篇先以各种花絮开头。过去几天里我和RStudio的作者先后在我们Ames村办大学、明尼苏达R用户组和纽约R用户组分别做了knitr与RStudio的报告,下周R官方会议useR! 2012在田纳西州举办,我们也有幸得到了在会上做邀请报告的机会。在这个报告里,我要谈的就是一些开发中的思考,本文先给出这些思考的一个预览。如果你之前不熟悉Sweave,下面的内容可能不太容易理解,但没关系,一来很多东西你已经没有理解的必要了(旧世界的糟粕),二来今后我还会详细介绍knitr的功能。

我自从09年来美帝开始,所有的作业和报告都是用Sweave写的(纯数学的除外),因此Sweave里面的边边角角我都比较熟悉,源代码也是看了一遍又一遍,包括后来基于Sweave扩展的pgfSweave包,我也是翻了很多遍源代码。最终结论是,Sweave继承了一个伟大的想法,但在具体实现上走入了一个死角,默认功能不强,扩展性又太差。随后在我给一门R课程做助教的时候,每次看学生用Word文档交来的作业都觉得丑陋不堪(少数人会精心调整排版,但你懂的),要重跑他们的代码实在太麻烦了。没有可重复的作业,何谈可重复的科学研究?在knitr的各种反馈中,我看到一条推特消息最令我欣慰:

学knitr吧!

他描述的是一个普遍事实:大多人都还在复制粘贴时代。然而,表面上看起来最直接的办法往往深藏隐患。复制粘贴不仅麻烦,而且将结果置于难以重复的境地。要是别人想重复你的分析,你得详细交待每一个操作步骤。万一一个步骤出错,可能会导致后面的全都错掉,并且修改起来也麻烦。代码能很好避免这些问题,一处代码改动,可以让后续结果全都自动更新。

为了让多数人走上正确的道路,我们只有一个选择,那就是:让正确的路比错误的路更容易走。如果你做不到比复制粘贴更快更简单,那么任何说教都是无效的。基于这个想法,我列举knitr的九条设计原则如下。

1、默认美观

软件默认设定非常重要,它决定了用户的第一印象。knitr默认代码高亮(无论什么输出格式)以及代码重整理,这都是为了增强代码和结果的可读性,面对一堆毫无生气的代码,谁都觉得累。为了设计默认的高亮主题,我专门请教了我们颜林林大站长和李龑大设计师;如果对默认主题不满意,knitr自带上百个高亮颜色主题,很方便切换。代码重整理的意思是,无论你的源代码多乱,我都给你自动重新整理整齐,熟悉我的工作的人可能能猜出来这是formatR包的功能。当初我向Sweave作者进谏重整理代码的功能被谢绝了,后来pgfSweave作者采纳了我的建议,现在这功能回到我自己的包中了。

knitr代码高亮

2、自然输出

就像德鲁克说(管理方面)好的企业看起来平淡无奇一样,好的软件也不应该有太多“惊喜”。Sweave有很多让用户感到意外的特征,比如基于grid的图形(如lattice和ggplot2)必须要print()才能被画出来,一个代码段中最多只能产生一幅图,要让输出中有图形,必须专门设置选项,等等。knitr秉承的设计理念是,同样的代码粘贴在R中看到的结果全部都会在knitr的默认输出中看到,有图出图,有表出表,不需要设置任何选项,一切自然而然。让用户必须记忆选项的软件不是好软件。

3、以分析为中心

在Sweave旧社会我们经常看到诸如cat('\\includegraphics{}')之类的代码,这样的代码往往是设计缺陷的症状,因为设计中缺乏某些功能,导致用户必须在R的层面上去弥补那些缺陷,这样数据分析代码和那些暗黑代码就混在了一起,数据分析者一会儿考虑统计方面的东西,一会儿考虑LaTeX方面的问题,精力难免分散。knitr去掉了所有需要用黑客方式去解决的问题,比如过去每幅图形的输出宽度设置很麻烦,Sweave引进了一项非常暗黑的LaTeX技巧,叫\setkeys{Gin},如果你不知道这个东西,建议你永远不要知道它。knitr解放了图形大小设置的问题,你可以对每一幅图形设置输出宽度(out.width选项)。

4、可重用的输出

这个想法很简单,就是让那些提示符>和续行符+有多远滚多远。我们常常看到这样的输出:

> if (TRUE) {
+ 1}
[1] 1

提示符对我来说毫无意义,它唯一的作用就是糟蹋源代码,让我没办法复制粘贴代码去运行。knitr默认输出是这样的:

if (TRUE) {
  1
}
## [1] 1

这是我这两年做助教恨得咬牙切齿的问题之一,现在我终于可以把提示符去掉了,输出也被注释掉,不影响复制代码运行。

5、功能模块化

道理说起来谁都懂,可到了现实世界中,无数码农仍然是一个五百行的函数打天下,各种功能拉不开扯不散,混在一起,难维护、难测试、难扩展。knitr的设计主线也就三部分:文档解析器、代码运行器、输出生成器。一个文档拿来,先抽代码出来,运行它,再根据运行结果写入输出文档。knitr在生成器上解放了生产力,引进了输出钩子函数的概念,让用户可以自定义结果输出方式,比如1 + 1在R里面会打印一个字符串[1] 2,利用knitr的钩子函数,你可以决定如何装裱这个字符串,可以是特殊的LaTeX环境:

\begin{mySource}
[1] 2
\end{mySource}

也可以是HTML代码:

<div class="mySource">
[1] 2
</div>

又或者是Markdown:

```
[1] 2
```

总之,你愿意怎么安排就怎么安排,knitr把运行过的代码和结果都给你。

6、好的功能照单全收

过去大家对扩展Sweave做了各种尝试,如pgfSweave、cacheSweave和weaver等包。你仔细看看这些包就会觉得无奈,每个包都先把Sweave那上千行源代码先复制一遍,再在局部进行一些修改,以实现增加新功能的目的。随着R自身的更新,这些被复制的源代码逐渐也落后于R,于是包的维护渐渐就成了问题,我基本上亲眼目睹了pgfSweave的兴衰过程。knitr收录了大多数跟Sweave有关的包的功能,这些功能基本上都以更简单的代码重写了,并且不需要复制八百行代码。其中我个人比较喜欢的是tikz图形、缓存和动画功能。

7、照顾初学者

每当我说LaTeX可能是壁垒时,总有人怀疑我(会R的人怎么会不会LaTeX)。knitr自今年初出道一来,让我感觉推广阻力最大的人群是org-mode的人。Emacs是万能的,嗯。JSS上今年出的org-babel论文四个月下载九千次,我关于knitr的一篇日志四个星期浏览九千次。最可怕的开发者就是认为用户应该懂这懂那,最好是通读自己的源代码。有时候这种高期望是对的,比如统计学,你要是不懂统计方法最好不要乱用函数,但有时候用户即使无知也无害,比如怎么把Markdown转化为HTML,这种事情他知道与不知道又有什么关系呢?如果点一下按钮就能生成结果,那么让用户点就是了,不必非得了解背后是怎么回事。

为了让初学者尽快入门,我最初在LyX 2.0.3中加入了knitr模块支持,让一键生成PDF变得可能,但LyX背后仍然是LaTeX,所以我需要一个不是非用LaTeX不可的编辑器支持。大约两个月前,RStudio的开发者联系到我,我们首先对LaTeX文档添加了knitr的支持,后来在我的建议下,又陆续添加了HTML和Markdown的支持。最近各种R Markdown的应用风生水起,与RStudio的支持密不可分。我选择Markdown作为给初学者入门的媒介,原因就是它超级简单,你可以在五分钟之内基本学会它的用法,若再多花点时间,完全有可能学完它的用法,注意是“学完”。这世上能被学完的语言不多,因为大多数语言都想让自己功能多,而Markdown是为了让功能少。

8、开放源代码需要开放

knitr是一个R包,当然也是开放源代码的,但对“开源”二字来说,存在一个“到底有多开放”的问题。有些开源产品有很好的API设计(如Wordpress),但有些则未必。knitr里除了核心的运行代码部分,其它几乎处处开放,举一个小例子:尽管knitr基于R,但它不一定非得运行R代码,如果你乐意,你可以嵌入Python或AWK或其它语言代码,这体现在engine参数上。

9、文学化编程也是编程

文学化编程(Literate Programming)是整个设计的核心思想,但过去的模式局限在“代码+文档”的简单模型上,knitr使得一份文档变得可编程。为了说明这个可编程的特性,举一个钩子函数例子(伪代码):

```{r tweet-hook, cache=FALSE, include=FALSE}
knit_hooks$set(tweet = function(before, options, envir) {
  library(twitteR)
  # Authentication with OAuth here, then
  if (!before) {
    msg = paste('I have finished the chunk',
                options$label, ', my Lord!')
    tweet(msg)
  }
})
# enable the chunk hook
opts_chunk$set(tweet = TRUE)
```

所谓钩子函数就是挂在代码段选项上的函数,当选项不为空(NULL)的时候,这个函数就会被执行。上面的tweet钩子的大意就是用twitteR包发推特消息,每当一个代码段运行完之后,就把该代码段的标签写入一个消息,然后发推特,这样随着整个文档被编译,推特上就会逐渐显示编译进度。钩子函数让一份文档超越了仅仅运行代码段的功能,你还可以用它执行一些附加任务。顺便再说一则花絮,6月12日第8届国际R语言会议上有一位讲师最近在准备培训材料,突然冒出一个想法,试探性问我有没有可能明年的R会议做出来,结果是不用等一年,用钩子函数5分钟就够实现了。这样的事情在Sweave的世界里几乎不可能完成。

其实关于knitr这个包我早已经写完一份中文介绍,感兴趣的可以先下手了。大多数文档仍然处于英文状态,但除非你是高级忍者,否则所有的英文文档只需要看选项文档基本就够了。

最后向大家介绍两个应用的例子:

  1. 云端的报告生成器:你什么都不需要安装,只需要一个浏览器,就够你生成报告了,后台基于OpenCPU(一个年轻但相当猛的REST架构云端R);
  2. RPubs.com:这又是一个基于knitr的云端服务,但需要你在本地RStudio中事先生成报告,再上传过去,相信在不久的将来,我们的作业和报告会变得漂亮,彻底告别那恶心的Word文档;

还有其它诸如HTML5幻灯片的例子在此就先不介绍了。如果你要学习knitr,建议从RStudio和Markdown起步(示例)。到目前为止从knitr的反馈来看,大家对Markdown都比较感兴趣,它可能的确迎合了初学者的需要:简单、可用。2012都来了,抓紧学点儿基本网页知识,相信不久的将来(如果还有将来的话)你一定会意识到它无穷的回报。

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

最近在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:打造一个可重复的统计研究流程》)。还是老话一句,不谈道德问题,因为这是不可控的因素,我们可控的只有制度和工具。源代码不会撒谎。

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