标签归档:P值

美国统计协会开始正式吐槽(错用)P值啦

(图片来源:https://xkcd.com/1478,一幅讽刺滥用P值的漫画)

今天美国统计协会(ASA)正式发布了一条关于P值的声(吐)明(槽),算起来可以说是近期统计学界的一件大事了。为什么这么说呢?首先,P值的应用太广,所以对P值进行一些解释和声明非常有必要。其次,对P值的吐槽历来有之,但今天是第一次被一个大型的专业协会以非常正式的形式进行澄清,多少带有一些官方的意思。声明的全文可以在这个页面中下载。

那么这则声明里面都说了什么呢?小编整体读了一遍,把我认为重要的信息概括在这篇文章之中。

首先,ASA介绍了一下这则声明诞生的背景。2014年,ASA论坛上出现了一段如下的讨论:

问:为什么那么多学校都在教 p = 0.05?

答:因为那是科学团体和期刊编辑仍然在用的标准。

问:为什么那么多人还在用 p = 0.05?

答:因为学校里还在这么教。

看上去多少有点讽刺的味道,但事实却也摆在眼前。从舆论上看,许许多多的文章都在讨论P值的弊端,小编摘录了几条言辞比较激烈的:

这是科学中最肮脏的秘密:使用统计假设检验的“科学方法”建立在一个脆弱的基础之上。——ScienceNews(Siegfried, 2010)

假设检验中用到的统计方法……比Facebook隐私条款的缺陷还多。——ScienceNews(Siegfried, 2014)

针对这些对P值的批评,ASA于是决定起草一份声明,一方面是对这些批评和讨论作一个回应,另一方面是唤起大家对科学结论可重复性问题的重视,力图改变长久以来一些已经过时的关于统计推断的科学实践。经过长时间众多统计学家的研讨和整理,这篇声明今天终于出现在了我们面前。

P值是什么

这份声明首先给出了P值一般的解释:P值指的是在一个特定的统计模型下,数据的某个汇总指标(例如两样本的均值之差)等于观测值或比观测值更为极端的概率。

这段描述是我们通常能从教科书中找到的P值定义,但在实际问题中,它却经常要么被神话,要么被妖魔化。鉴于此,声明中提出了六条关于P值的准则,作为ASA对P值的“官方”态度。这六条准则算是这条声明中最重要的部分了。

继续阅读美国统计协会开始正式吐槽(错用)P值啦

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

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

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

假设检验初步

准备尝试一下,用大白话叙述一遍统计推断中最基础的东西(假设检验、P值、……),算是把这段时间的阅读和思考做个梳理(东西不难,思考侧重在如何表述和展示)。这次打算用一种“迂回的”表达方式,比如,本文从我们的日常逻辑推理开始说起。

0.普通逻辑

复习一下普通逻辑的基本思路。假设以下陈述为真:

你打了某种疫苗P,就不会得某种流行病Q。

我们把这个先决条件表述如下:

如果P 则非Q

其中,

P表示打了疫苗P,

Q表示得流行病Q

或者,更形式化一点:

if P then NOT Q

然后,如果观察到你得了流行病Q,那么就可以推出你没有打疫苗P——这个推断只不过是上述前提条件的逆反命题而已。我们把以上推理过程表述如下:

if P then NOT Q (先决条件)

Q (前提)

———————–

then NOT P (结论)

还有,如果你没有得流行病Q,就能推断出你打了疫苗P吗?显然不能。打疫苗P是不得流行病Q的充分条件,但非必要条件:你没有得流行病Q,可能是因为打了疫苗P,也可以是因为其他任何原因。即,if P then NOT Q,不能够推出if NOT Q then P。

到此为止没有任何令人惊奇的地方。下面将表明,假设检验背后的统计推断规则也只不过是我们以上日常逻辑推理的一个衍生而已。这只需要思维的一次小小的“跳跃”。

1.假设检验

在统计推断中,我们不说“你打了疫苗P,就不会得流行病Q”,而是说,比如,“你打了疫苗P,就有95%的把握不会得流行病Q”,即if P then probably NOT Q。把上面的逻辑推理规则改写成统计推断规则:

if P then probably NOT Q    (先决条件)

Q                                                     (前提)

———————–

then probably NOT P         (结论)

回到以前“万能”的硬币实验,我们做实验来考察一枚硬币是不是均匀的。改写成现在我们熟悉的形式:

P:硬币是均匀的。

Q:在100次投掷中,得到90次正面,10次反面。

我们说,如果是一个均匀的硬币,就不太可能发生这样的情形:投100次,出现90次正面,10次反面(if P then probably NOT Q)。现在如果在100次投掷实验中,观察到出现90次正面,10次反面(Q),那就可以有把握地说,这个硬币不是均匀的(NOT P)。这个推理可以写成与上面一致的统计推断的形式,其中,P是原假设H0,NOT P是备择假设Ha:

H0:硬币是均匀的  (P

Ha:硬币是有偏的 (NOT P

如果原假设为真,即硬币是均匀的,就不太可能发生这样极端的事情,比如:在100次投掷实验中,观察到出现90次正面,10次反面(Q)。如果真的观察到这样极端的事情,你就有把握认为硬币不是均匀的,即拒绝原假设(P),接受备择假设(NOT P)。

另外,如果在100次投掷实验中,观察到60个正面,40个反面(NOT Q)。这时你就不好下结论了,因为一个均匀的硬币可能投出这样的结果,一个有偏的硬币也可能投出这样的结果。最后,你只能说,如果实验结果是这样的,那就没有把握拒绝原假设。这枚硬币是否有偏,需要更多的证据来证明(这通常意味着更多的实验,比如,再投1000次)。

总结一下。在搜集数据之前,我们把想证明的结论写成备择假设,把想拒绝的结论写成原假设。之所以写成这个形式,因为从上面不厌其烦的讨论中得知,这是方便逻辑/统计推断的形式:当我们难以拒绝原假设时,只能得到结论,原假设也许是真的,现在还不能拒绝它;而当我们能够拒绝原假设时,结论是:它就很有把握是不真的。注意,在看到数据之前,我们不知道自己想证明的结论是否能够被证据所支持。

在确定假设检验的形式的同时,我们对之前一直随意说的“把握”、“可能”也做一个限定,即指定一个显著性水平α(significance level),也叫犯第一类错误的概率(type I error,在上面的硬币实验中,就是否定一个均匀硬币的错误,也叫“弃真”错误)。

根据某些保守或稳健的原则(比如,我们认为,把一个无辜的人判决为有罪,比放掉一个有罪的人,后果更为严重),我们要尽量把犯“弃真”错误的概率控制在一个很小的水平里。通常α=0.05,这时候就是说,如果拒绝了原假设,你就有95%的把握说原假设是不真的。这里,95%(=1-α)就是置信水平(confidence level)。

又,放掉一个有罪的人,即把一个有罪的人判为无罪,这犯的是第二类错误β(type II error,在硬币实验中,就是把一个有偏的硬币当成均匀硬币的错误,也叫“取伪”错误)。关于第一类和第二类错误之间的权衡取舍(trade off),详见《决策与风险》。在我们的假设检验里,我们认为犯一类错误的后果比犯第二类错误的后果更为严重。

需要注意的是,在这里,我强调的是先提出需要检验的假设,然后再搜集收据。这是统计推断的原则之一。如果看到了数据之后再提出假设,你几乎可以得到所有你想要的结果,这是不好的机会主义的倾向。强调这些,是因为在学校里,我们大多是看了别人搜集好的数据之后再做统计练习。

事先确定好你想拒绝/证明的假设,在看到数据之前,你不知道结果如何。

2.P值(P Value)

上面提到“极端”事件,比如,在100次硬币投掷实验中,观察到出现90次正面,10次反面(Q)。怎么样的事件才是“极端的”?简单地说,一个事件很极端,那么少比它本身“更极端”的事件就非常少(比如,只有“91次正面,9次反面”、“91次反面,9次正面”等情况才比它更极端)。

但这个Q只是从一次实验中得出的。我们可以重复做这个实验,比如100次,每次都投掷100次,记录下的正面数X,它构成一个二项分布,X~B(n,p),其中,n=100,p=0.5。根据某个中心极限定理,正态分布是二项分布的极限分布,上面的二项分布可以由均值为np=50,方差为np(1-p)=25的正态分布来近似。我们在这个近似的正态分布的两端来考察所谓“更极端”的事件,那就是正面数大于90或者小于10。

重复一遍,“P值就是当原假设为真时,所得到的样本观察结果更极端的结果出现的概率”。如果P值很小,就表明,在原假设为真的情况下出现的那个分布里面,只有很小的部分,比出现的这个事件(比如,Q)更为极端。没多少事件比Q更极端,那就很有把握说原假设不对了。

在上述近似的正态分布中,P值就等于X>90 或 X<10的概率值(记做,P{X>90 or X<10})。根据对称性,这个概率值等于2*P{X<10}=1.2442E-15。

上面我们的确求出了一个非常小的P值,但如何不含糊地确定它就是很“极端”呢? 事先确定的显著性水平α,本身就是一个判定法则。只要P值小于显著性水平α,我们就认为,在认为原假设为真的情况下出现的事件Q,是如此地极端,以至于我们不再相信原假设本身。一句话,我们的判定法则是:

P值小于显著性水平α,拒绝原假设。

3.一个手算示例

用一个双侧的单样本T检验做例子。假设我们想知道,螃蟹的平均温度,跟空气的温度(24.3)有没有统计差别(α=0.05)。事先确定的假设检验的形式表达如下:

零假设(H0):   μ=24.3°C

备择假设(Ha):  μ≠24.3°C

以下是25只螃蟹在温度为24.3°C下的体温(单位:°C):

25.8    24.6    26.1    22.9    25.1
27.3    24        24.5    23.9    26.2
24.3    24.6    23.3    25.5    28.1
24.8    23.5    26.3    25.4    25.5
23.9    27        24.8    22.9    25.4

一些基本的算术结果:

样本均值:$\bar{X}=25.3$

样本量:n=25

样本方差:$s^2$=1.8

样本均值的标准误差:$s(\bar{X})=\sqrt{s^2/n}=0.27$

这里T检验的思路如下:

  1. 我们先假设H0为真,即认为螃蟹的平均温度跟空气温度没有差异(P),  μ=24.3°C。有一个极端事件Q,如果原假设H0成立,Q就不成立(if H0 then probably NOT Q);但如果在原假设为真的情况下,出现了这么一个Q,那我们就有把握拒绝原假设。
  2. 样本均值:$\bar{X}$是总体均值μ的最好的估计,在本例中,$\bar{X}=25.03$。
  3. 这个样本均值只是一个估计值。它只是从总体的一个随机样本中得到的(样本是上述25只螃蟹)。我们不知道这次实验结果是不是“极端”事件。而判断一个事件是不是极端事件,根据第二节的讨论,我们可以重复做上述实验,比如100次,每次都抓25只螃蟹,都在空气温度24.3的状态下测量其体温,然后也各自求出一个样本均值来。
  4. 容易得出,这种实验出来样本均值,辅以适当的数学形式,就服从一个自由度为24(=25-1)的t分布,即$(\bar{X}-\mu)/s(\bar{X})\sim t(24)$。
  5. 样本均值$\bar{X}=25.03$,在这个自由度为24的t分布下,有一个对应的t值,t=25.03-24.3/0.27=2.704。现在我们可以在整个分布里考察这个t值。在这个自由度为24的t分布里,我们看 t=2.704是不是一个“极端”事件Q。根据对称性,比Q更极端的是那些大于2.704或者小于-2.704的点。

t

从上图可以看到,在这个t分布里,比t=2.704更“极端”的点占整个分布的0.0124。这个0.0124就是我们要求的P值。这个P值小于我们事先选定的显著性水平α=0.05,因此我们可以拒绝原假设,认为这批螃蟹的平均体温不等于空气温度。

这个双侧P值可以手算如下:

在SAS里,P=2*(1-probt(t,df))=2*(1-probt(2.704,24))=0.012392

在R里,     P=2*(1-pt(t,df))=2*(1-pt(2.704,24))=0.012392

———-

以上是用P值作为判定条件。一个等价的做法是用临界值来判断。我们事先给定的显著性水平α=0.05,在这个自由度为24的t分布里,就对应着一个临界t值2.064。下图的阴影部分,也称作拒绝区域。上面求出的跟样本均值$\bar{X}=25.03$对应的t值=2.704,处在这个拒绝区域内(2.704>2.064),于是我们一样拒绝原假设。

t2

又,上述临界值可以手算(或查表)如下:

在SAS里,tCritic=tinv(1-alpha/tail,df)=2.06390

其中,alpha=0.05,tail=2表示双侧检验,df=24.

在R里,tCritic=qt(1-alpha/tail,df)=2.063899

4.注

本文是对近期阅读做的一个笔记。作为一个非统计科班出身的程序员,我一直在思考,如何来理解统计概念,以及如何把自己的理解向同行传达。关于用日常逻辑推理来理解假设检验的思路,来自

Common Statistical Methods for Clinical Research with SAS Examples(2nd edition, SAS Inc., 2002, by Glenn A. Walker)

关于决策与风险的讨论,参考了

维恩堡《数理统计初级教程》(常学将等译,太原:山西人民出版社,1986,Statistics: An Intuitive Approach By George H. Weinberg and John Abraham Schumaker)

第三节示例的数据,来自

Biostatistical Analysis (5th Edition) by Jerrold H. Zar, Prentice Hall, 2009

第三节的t分布图,来自一个在线的t分布生成器(很好用):

http://onlinestatbook.com/analysis_lab/t_dist.html

附录: 用SAS来计算

上面的文字尽量做到“平台无关”。这里附出SAS例子,是想把以上的手算结果跟机器结果做个对照,让读者更有信心一些。 欢迎读者贴出自己趁手的工具得出的结果。

/*data*/
data body;
input temp @@;
h0=24.3;
diff=temp-h0;
datalines;
25.8    24.6    26.1    22.9    25.1
27.3    24      24.5    23.9    26.2
24.3    24.6    23.3    25.5    28.1
24.8    23.5    26.3    25.4    25.5
23.9    27      24.8    22.9    25.4
;

/*method 1: use proc means*/
proc means data=body T PRT;
var diff ;
run;

结果是:

t Value    Pr > |t|
——————-
2.71      0.0121
——————-

上面的t Value 就是计算出来的t值,Pr > |t| 就是P值(这里的|t|就是上面计算出来的t值2.704,Pr > |t|求的是比t值更极端的概率,即P值)。proc means没有提供临界t值(即通常说的查表得出的t值),下同。

/*method 2 (prefered): use proc ttest*/
proc ttest data=body h0=24.3 alpha=0.05;
var temp;
run;

proc ttest的结果更为丰富:

N      Mean     Std Dev  Std Err    Minimum   Maximum

25     25.0280      1.3418      0.2684 22.9000        28.1000

Mean     95% CL     Mean       Std Dev     95% CL   Std Dev

25.0280 24.4741  25.5819            1.3418             1.0477   1.8667

DF    t Value    Pr > |t|

24       2.71           0.0121

也谈提高R语言的运算效率

用过底层语言做计算的人转入R语言的时候一般都会觉得R语言的运算太慢,这是一个常见的对R的误解或者对R的设计的不理解。在二三十年前Chambers等一群人在贝尔实验室设计S语言之前,统计计算主要是通过那些底层语言实现的,典型的如Fortran。当时有一个基于Fortran的统计分析库,用它的麻烦就在于无论做什么样的数据分析,都涉及到一大摞繁琐的底层代码,这让数据分析变得很没劲,统计学家有统计学家的事情,天天陷在计算机程序代码中也不是个办法,要摆脱底层语言,那就只能设计高层语言了。有所得必有所失,众所周知,高层语言一般来说比底层语言低效,但对用户来说更友好。举个简单的例子,用C语言计算均值时,我们得对一个向量(数组)从头循环到尾把每个值累加起来,得到总和,然后除以向量的长度,而均值在统计当中简直是再家常便饭不过了,要是每次都要来这么个循环,大家也都甭干活儿了,天天写循环好了。

前两天COS论坛上有个帖子提到“R语言的执行效率问题”,大意如下:

有120000行数据,用perl写成12万行R命令做t.test,然后执行,大概几分钟就算完了。如果只用R语言,把所有数据先读入,然后用循环,每一行做t.test,花了几个小时,到现在还没有算完。这说明一个问题,在R里执行单行命令要比用循环快,涉及到循环的问题,最好写成单行命令执行。

我不清楚作者是如何写这“12万行”R命令的,本文假设是把t.test(dat[i, ]), i = 1, 2, ..., 120000写入一个文件,然后用source()执行之。面对这样一个问题,我们有若干种改进计算的可能性。首先我们看“硬”写入程序代码的方法:

## 为了使计算可重复,设定随机数种子
set.seed(123)
## 生成数据,随机数,120000行 x 100列矩阵
dat = matrix(rnorm(120000 * 100), 120000)
nr = nrow(dat)
nc = ncol(dat)
## 六种方法的P值向量
p1 = p2 = p3 = p4 = p5 = p6 = numeric(nr)

## via source()
f = file("test.t")
writeLines(sprintf("p1[%d] = t.test(dat[%d, ])$p.value",
    seq(nr), seq(nr)), f)
system.time({
    source(f)
})
#   user  system elapsed
#  95.36    0.19   95.86
close(f)
unlink("test.t")

1、向量式计算:apply()以及*apply()

当我们需要对矩阵的行或者列逐一计算时,apply()系列函数可能会提高效率。本例是对矩阵的行做t检验,那么可以这样计算:

## via apply()
system.time({
    p2 = apply(dat, 1, function(x) {
        t.test(x)$p.value
    })
})
#   user  system elapsed
#  63.12    0.06   63.50
identical(p1, p2)
# [1] TRUE

结果比第一种方法快了大约半分钟,而且计算结果完全一样。apply()本质上仍然是循环,但它在某些情况下比直接用显式循环要快:

## via for-loop
system.time({
    for (i in seq(nr)) p3[i] = t.test(dat[i, ])$p.value
})
#   user  system elapsed
#  69.88    0.03   70.05
identical(p2, p3)
# [1] TRUE

不过apply()系列函数在提高运算速度上优势并不会太明显,提倡用它的原因是它和统计中的矩阵运算相似,可以简化代码,相比起$\sum_{i=1}^n x_i/n$,我们可能更愿意看$\bar{x}$这样的表达式。很多R内置函数都是用底层语言写的,我们需要做的就是把一个对象作为整体传给函数去做计算,而不要自行把对象分解为一个个小部分计算,这个例子可能更能体现向量式计算的思想:

system.time(apply(dat, 1, mean))
#   user  system elapsed
#   5.28    0.04    5.25
system.time({
    x = numeric(nr)
    for (i in 1:nr) x[i] = mean(dat[i, ])
})
#   user  system elapsed
#   4.44    0.02    4.42
system.time(rowMeans(dat))
#   user  system elapsed
#   0.11    0.00    0.13

2、明确计算的目的

很多情况下,R函数返回的不仅仅是一个数字作为结果,而是会得到一系列诸如统计量、P值、各种系数等对象,在我们调用R函数之前如果能想清楚我们究竟需要什么,也许对计算的速度提升有帮助。比如本例中,也许我们仅需要12万个双边P值,其它数字我们都用不着,那么可以考虑“手工”计算P值:

## "hand" calculation in R
system.time({
    p4 = 2 * pt(apply(dat, 1, function(x, mu = 0) -abs((mean(x) -
        mu)/sqrt(var(x)/nc))), nc - 1)
})
#   user  system elapsed
#  15.97    0.07   16.08
identical(p3, p4)
# [1] TRUE

上面的计算更“纯净”,所以计算速度有了本质的提升,而且计算的结果和前面毫无差异。在做计算之前,人的脑子多思考一分钟,也许计算机的“脑子”会少转一个小时。

3、把四则运算交给底层语言

R是高层语言,把它拿来做简单的四则运算是很不划算的,而且容易导致程序低效。加加减减的活儿是C和Fortran等底层语言的强项,所以可以交给它们去做。以下我们用一段C程序来计算t统计量,然后用R CMD SHLIB将它编译为dll(Windows)或so(Linux)文件,并加载到R中,用.C()调用,最终用R函数pt()计算P值:

## "hand" calculation in C for t-statistic
writeLines("#include <math.h>

void calc_tstat(double *x, int *nr, int *nc, double *mu, double *tstat)
{
    int i, j;
    double sum = 0.0, sum2 = 0.0, mean, var;
    for (i = 0; i < *nr; i++) {
        for (j = 0; j < *nc; j++) {
            sum += x[i + j * *nr];
        }
        mean = sum / (double) *nc;
        sum = 0.0;
        for (j = 0; j < *nc; j++) {
            sum2 += (x[i + j * *nr] - mean) * (x[i + j * *nr] - mean);
        }
        var = sum2 / (double) (*nc - 1);
        sum2 = 0.0;
        tstat[i] = (mean - *mu) / sqrt(var / (*nc - 1));
    }
}", "calc_tstat.c")
system("R CMD SHLIB calc_tstat.c")
dyn.load(sprintf("calc_tstat%s", .Platform$dynlib.ext))
system.time({
    p5 = 2 * pt(-abs(.C("calc_tstat", as.double(dat), nr, nc,
        0, double(nrow(dat)))[[5]]), nc - 1)
})
#   user  system elapsed
#   0.86    0.06    0.92
dyn.unload(sprintf("calc_tstat%s", .Platform$dynlib.ext))

因为R可以用system()调用系统命令,所以整个过程全都可以用R完成,Windows用户需要安装Rtools并设置系统环境变量PATH才能使用R CMD SHLIB

更进一步,因为R自身的一些C程序也是可供用户的C程序调用的,所以我们可以把整个P值的计算过程全都扔进C代码中,一步完成:

## "hand" calculation in C calling Rmath.h
writeLines("#include <Rmath.h>
void calc_pvalue(double *x, int *nr, int *nc, double *mu, double *pval)
{
    int i, j;
    double sum = 0.0, sum2 = 0.0, mean, var;
    for (i = 0; i < *nr; i++) {
        for (j = 0; j < *nc; j++) {
            sum += x[i + j * *nr];
        }
        mean = sum / (double) *nc;
        sum = 0.0;
        for (j = 0; j < *nc; j++) {
            sum2 += (x[i + j * *nr] - mean) * (x[i + j * *nr] - mean);
        }
        var = sum2 / (double) (*nc - 1);
        sum2 = 0.0;
        pval[i] = 2 * pt(-fabs((mean - *mu) / sqrt(var / (*nc - 1))),
                      (double) (*nc - 1), 1, 0);
    }
}", "calc_pvalue.c")
system("R CMD SHLIB calc_pvalue.c")
dyn.load(sprintf("calc_pvalue%s", .Platform$dynlib.ext))
system.time({
    p6 = .C("calc_pvalue", as.double(dat), nr, nc, as.double(0),
        double(nrow(dat)))[[5]]
})
#   user  system elapsed
#   0.83    0.07    0.91
dyn.unload(sprintf("calc_pvalue%s", .Platform$dynlib.ext))

头文件Rmath.h的引入使得我们可以调用很多基于C程序的R函数,详情参考手册Writing R Extensions。通过C计算出来的P值和前面用R算的略有差异,下面画出p6 - p1 vs p1以及p6 - p5 vs p5的图:

P值的差异
P值的差异

导致差异的原因此处不细究,感兴趣的读者可以帮忙检查一下。

小结

  1. 若你熟悉底层语言,计算又不太复杂,那么可用底层语言写,然后用R调之;
  2. 否则把R对象当做整体去计算,能做x + 1就不要做for (i in length(x)) x[i] + 1
  3. 不要低估R core们的编程水平,他们已经做了很多工作让用户脱离底层编程

注:本文中的运算时间可能不可重复,这与计算机所处的状态有关,但大体来说,运算速度的快慢是可以肯定的。本文仅仅是关于统计计算的一个微小的例子,以后若有更好的例子,可能会更新本文;也欢迎各位提供更多示例。

不拒绝零假设意味着什么

由一道试题引发的一点思考

2008年统计学考研真题第四题“食品厂家说:净含量是每袋不低于250g。但有消费者向消协反映不是250g,消协据此要求厂家自检,同时消协也从中随机抽取20袋检验”

(1)如果厂家自己检验,你认为提出什么样的原假设和备则假设?并说明理由。
(2)如果从消费者利益出发,你认为应该提出什么样的原假设和备则假设?并说明理由。

……

作为统计专业的学生来说,熟悉得不能再熟悉了。但是,通过做上面的题目,我发现自己在理解假设检验的问题上犯了一个十分严重的错误。这个问题主要是由于我们学的教材上面写着:“假设检验要么P-value小于a拒绝原假设,P-value大于a接受原假设……”。后来再看看其他教材,发现绝大多数都是这样写的。其实“P-value大于a接受原假设”这种说法是错误的。

P-value大于a的时候,结论到底是什么呢?最早提出这个问题的是:

E·皮尔逊问耶日·奈曼,在检验一组数据是否为正态分布时,如果没能得到一个显著性的 P值,那么怎样才能看这组数据是正态分布的呢?

费歇尔其实已经间接地回答了这个问题。费歇尔把比较大的 P 值(代表没有找到显著性证据)解释为:根据该组数据不能做出充分的判断。依据费歇尔的解释,我们绝对不会得出这样的推理,即没有找到显著性的证据,就意味着待检验的假设为真。这里
引用费歇尔的原话:

“相信一个假设已经被证明是真的,仅仅是由于该假设与已知的事实没有发生相互矛盾,这种逻辑上的误解,在统计推断上是缺乏坚实根基的,在其它类型的科学推理中也是如此。当显著性检验被准确使用时,只要显著性检验与数据相矛盾,这个显著性检验就能够拒绝或否定这些假设,但该显著性检验永远不能确认这些假设一定是真的,……”

所以假设检验的目的在于试图找到证据拒绝原假设,而不在于证明什么是正确的。当没有足够证据拒绝原假设时,不采用“接受原假设”的表述,而采用“不拒绝原假设”的表述。“不拒绝”的表述实际上意味着并未给出明确的结论,我们没有说原假设正确,也没有说它不正确。

举个例子来说:比如原假设为H0: m =10,从该总体中抽出一个随机样本,得到`x=9.8,在a=0.05的水平上,样本提供的证据没有推翻这一假设,我们说“接受”原假设,这意味着样本提供的证据已经证明m=10是正确的。如果我们将原假设改为H0: m =10.5,同样,在a=0.05的水平上,样本提供的证据也没有推翻这一假设,我们又说“接受”原假设。但这两个原假设究竟哪一个是“真实的”呢?我们不知道。

总之,假设检验的主要目的是为了拒绝而不是接受。

由一道试题引发的另一道试题

这让我想到2007年统计学考研真题第一题问:“正态分布的假定能不能用数据证明?”

其实也是关于假设检验的问题,具体内容可以参看: