统计学
收集、整理、分析和展示数据的科学
后来回到有网络的世界，git push一下，就同步上网了，供各位客官公开下载（见于本站“作品”页面）。过了一阵子，领导的领导们看见了，觉得这样做有点不靠谱。再过了一阵子，又有好几个人建议我不要放出来公开下载。我呢，也绝非来自石器时代、对“世风日下人心不古”一无所知，佛曰“我不下地狱那哪个来下地狱呢？”佛敢下地狱，那是因为他是佛，是佛又咋地？佛有核心竞争力啊。现在这个忙碌年代，我们（别人）想尽办法折腾别人（我们），折腾的途径就是评价。网站不看内容要看备案，活人不看要看KPI，百姓生活不看要看GDP……而混学界的人就得靠著作来评价，于是使出浑身解数去发表，自己努力也好，坑蒙拐骗也好。我这样赤果果把作品放在网上供任何人下载，说出来真是一件很恐怖的事情，娘亲呐，多少豺狼虎豹盯着呢。
本来想继续扯到李一、马云等人，似乎有点扯远了，我是在想像马云这样伟大的宗教领袖式人物为何需要道教思想的支撑。俺们这个时代的人，真的越来越没底气了？打住打住。
最后我接受了腾飞的意见，把原本公开的书稿撤下来了。倒不是怕盗版抄袭，主要是一本不完整的书稿，对读者似乎不是一件好事，等正式完成之后，再集中发给大家提意见，这样效用会大一些；再想一想，出版之前公开下载也会对出版社带来一些麻烦，于是撤销下载。当然，目前我仍然内部开放书稿，我熟悉放心的人以及答应我在COS论坛帮别人回答图形方面100帖的人（鼓励为别人做贡献的人），皆可获得我的书稿。这本书的目录部分可以公开浏览：
各位客官若有任何建议，请随时与我联系，我会在完稿前尽量多采纳读者的意见。
—————————————
注1：平民老百姓。
这论文呢，自我评价可以给个80分。最大的问题在于没有花时间去整理文章的结构，所以构架上稍微有点散乱（俗称“意识流”）。内容上熟悉我的博客的客官一眼就能看出来，其实都是些博客文章的汇总，只不过用LaTeX让它们变得“人模狗样”一些而已，好在本小子平时也积攒了这些鸡零狗碎的东西，动过自己的脑子。我觉得群众大学的毕业论文，很多都是一个套路：经济/金融数据套一个神奇的模型，直到最后整个世界一片和谐，读者在最后一章都能隐约看到上帝老爷子在朝你挥手。其实也没啥，找工作不容易，地球人也都知道写论文就是忽悠——漫漫人生路上一道工序。
由于本小子是个小人（小小的活人），所以总关心小人关心的事情（俗称“人本主义”）。这论文嘛，窃以为也没什么上下高低之分，说出你怎么想的就可以了，而不要总说“他们”怎么想怎么做。一定要有数学上的创新？一定要有人家看不懂的公式才是好论文？一定要有综述？一定要有长长的参考文献列表才是好论文？一定要板起脸？不能写八卦？不准幽默？……嗨，作茧自缚。几年前看到一篇好文章，颇具恶搞性质，建议各位客官收藏：How to write Consistently Boring Scientific Literature。
言归正传：本文是厌倦八股文和数学理论的产物，从理论角度来说，几乎没什么价值，不过这篇文章是用Sweave写的，完全具有可重复性和100%透明度，对文中结果有怀疑的客官可以自行运行代码；其次，统计模拟和图形的声音在界内太微弱，大家都很忙，有人在忙着推公式，有人在忙着编数据，有人在忙着把公式用到不知道是不是编出来的数据上，本小子跟着瞎掺和了点别的东西，仅此而已。甭管有用没用，敬请拍砖。
－－－－－－－－－－外一篇：坛霸是怎样练成的－－－－－－－－－－
曾经有童鞋称呼在下为“坛霸”，这个……有时候确实有那么点意思，无图无真相（两个多月没怎么回帖了，一鼓作气）：
坛霸是怎样练成的
接下来我会陆续写第三届中国R语言会议、《现代统计图形》书稿和useR! 2010，若时间允许，我考虑一下电视剧《九阴真经》（93版）。
近日来收到邮件少了，但各个问题都不太好直接回答。比如这则关于McNemar检验的问题：McNemar检验可以有两种形式的统计量，一为(b – c)2/(b + c)，一为2b*log(2b/(b+c)) + 2c*log(2c/(b+c))，其中b和c是列联表非对角线上的频数。前者是McNemar检验本身的统计量，可以根据渐近正态分布得来（然后平方得到卡方），后者是似然比统计量（不带约束的似然除以带约束的，取对数，乘2）。McNemar检验看似复杂，实际上可以简化为检验b = c，或等价于检验一个n = b+c的二项分布中，是否p = 1/2（观察到X = b或c）。现在的问题是，这两种统计量有没有优劣之分？
作为一个懒得推公式的人，我向来喜欢用模拟回答问题，因为模拟的结果非常直截了当。我的考虑是，要看渐近统计量的优劣，那就看随着n增大，统计量和渐近分布有多接近好了。一个自然而然的想法当然是对若干统计量的观测值做分布检验了，比如KS检验。我们知道这两个统计量都是自由度为1的卡方分布，剩下的事情就是计算：
set.seed(123)
nmax = 1000
p = matrix(nrow = nmax, ncol = 2)
for (n in 2:nmax) {
# 生成服从二项分布的随机数，分别计算两种统计量并作KS检验、记录P值
b = rbinom(500, n, 0.5)
x1 = (b - (n - b))^2/n
x2 = 2 * b * log(2 * b/n) + 2 * (n - b) * log(2 * (n - b)/n)
p[n, 1] = ks.test(x1, "pchisq", df = 1)$p.value
p[n, 2] = ks.test(x2, "pchisq", df = 1)$p.value
}
# 调整一下数据格式，画图：随着n增大，P值如何变化？
library(ggplot2)
d = melt(p, varnames = c("n", "method"))
d$method = factor(d$method, labels = c("McNemar", "LRT"))
colnames(d)[3] = "p.value"
qplot(n, p.value, data = d, shape = method, geom = c("smooth", "point")) +
scale_shape_manual(values = c(2, 3))
McNemar检验统计量与卡方分布拟合的好坏
早上和其他两位评委Simon Urbanek以及Hadley Wickham进行了电话会议，我们将今年的Chambers奖授给Michael J. Kane和他的bigmemory包（剧透了剧透了）。通过看今年提交的参赛作品，我觉得拿下这个奖的困难并没有想象中那么大，国内的客官们努力努力，也是很有希望获奖的（比如我相信精于C++的颜大站长能独立写出bigmemory包的概率大于95%）。此前在COS论坛上呼吁大家踊跃参加，估计大家都觉得这是天方夜谭，明年我以95%的概率不会做评委了，不过这评奖过程给我几点感想可供后来人借鉴：
- 严格按照主办方的规则行事。主办方的评奖规则中怎么写，我们就对照这一条一条规则检查自己的作品是否都符合了要求。比如Chambers奖的规则描述是：
- 关于这些原则，如果参赛者能站在评委角度来考虑，肯定能为自己挣得不少分。也许有些作者软件写得很精妙，但缺少恰当的表现形式，所以就可能被埋没。记得有一位参赛者把自己三百多页的博士论文都发来了，满篇数学公式，想想如果自己是评委，看一篇博士论文的概率是多大呢？
- 包装很重要。由此又不得不说LaTeX与Word……呃，各位把本小子看作技术愤青的大人们，这真的不是技术问题，除非是Word高手，普通人用Word做出来的任何文档的排版质量跟LaTeX一比，评论只有两个字：垃圾。没得商量。LaTeX生产出来的论文，即使内容连垃圾都不如，其形式看起来也是正儿八经能唬住人的。拿着Word写的灰头土脸的PDF文档交上来，首先给人印象就是这童鞋以95%的概率不是高手，否则怎么连LaTeX都不会用呢。
The entries will be judged on a variety of dimensions, including the importance and relevance for statistical practice of the tasks performed by the software, ease of use, clarity of description, elegance and availability for use by the statistical community. Preference will be given to those entries that are grounded in software design rather than calculation.
最终评委的评分规则便根据三原则来：重要程度和与统计学的相关程度（多数作品都有很大的专业局限性，仅仅在自己的领域里针对某一特定模型写了软件包，不够通用，我也看不懂什么生物名词或天体物理名词）、创新和软件设计（想法是否足够新颖，没人做过当然最好，有人做过则要想想如何与众不同）、易用性和文档是否清楚（如果参赛者能多提供一些例子则会让评委更快了解你的软件，可以是录像、在线演示或动画、图形）。
再回头看本小子去年的申请，估计很大程度上得益于本小子的动画网站，以及每个动画函数下都有例子展示，评委不用动脑子去仔细研究函数的每一个参数怎么用，只需要端着咖啡看演示就可以了。另外，本小子处心积虑套用了John Chambers那句名言“To turn ideas into software, quickly and faithfully”（我把software换成了animations），这主要是为了体现软件包与统计学思想的联系，另一方面，一眼看去这和竞赛的主题切合得甚为紧密。
去年我在申请的时候还没接触到LyX，所以老老实实写LaTeX源代码然后老老实实编译，而且用的是和R News文章一样的字体（这也是评委之一Hadley常用的字体），呈上去给大佬们一看，嘿，眼熟，我看这小子和R有一定关系。
最后，这获奖者比其他选手还占一点优势，就是他在去年的JSM大会上做过Data Expo的poster，硕大的宣传板，我们三人都在那里看过，回头一想，脸熟啊。这一点呢，也和我去年类似，三位评委我见过两位，一位七分熟，一位三分熟（你煎牛排呢？），剩下一位未曾近距离接触，但我曾给他的一本书提过一处勘误。所以混圈子也是有用的。
这些不是告诉各位客官可以不劳而获或投机取巧，世上没那么多好事，而是用一个例子说明怎样小心地铺路，把自己能控制的因素都一步步做到最佳状态，剩下的事，或水到则渠成，或听天而由命……理想情况下，呔！手起刀落，砍他个人仰马翻。
又及：
只是老师布置的期末作业——选择合适的题目，收集相应的数据，建立统计模型，进行统计分析，直到你认为满意为止（咱们老师的原话）。是不是我的题目选的不好啊？再及：
如果直接访谈的话，我感觉暨很费时间，而且结果也不一定准确（其中也会受很多因素的影响），如果用回归的话，做问卷我可以改成匿名形式，只是变量的选取和设定变的有难度（对于我来说）原本这是一个如何保护受访者隐私的问题。这种问题有一些经典的解决办法，比如让受访者自己抛硬币，如果正面就回答“是否挂过科”，反面就回答“宿舍电话最后一位数字是否是奇数”，访员不干涉受访者填问卷的过程，最终我们也不知道受访者的硬币是正面还是反面（从而不知道他们究竟回答的是哪个问题），只知道他们回答了多少“是”和“否”。只要样本量充分大，我们就知道挂科的比例了。
至于变量间的相关性，这不是问题，因为世上完全独立的变量似乎还没生出来。即使再独立，不还存在所谓的“蝴蝶效应”么？退一步讲，回归也没有要求自变量相互独立，相关性太强也有解决办法。
本人怀着阴暗心理于今日公交车上翻IMS Bullutin的第39卷第1期，不幸看到第5页，又是一例统计结果不可重现的例子。一帮人，用可公开获得的数据获得了惊天的成功，到头来被人指责结果不可重复，而且不是一点半点不可重复，后人说的是“results are no better than chance”，嘿嘿，我心里冷笑着。你说说，他整一方法，号称威力无比，其实跟抛硬币得出来的结果没啥区别；这病人得没得癌症，抛个硬币决定吧。
难道这就是传说中的随机数发生器？又想起有些人用上百个变量做回归，这也是随机数发生器的一种，找出几个带显著性星号的系数不必欣喜若狂，要是找不出那才奇怪了呢。呜呼。
前几天看见这么一则报道，一直挂在我的浏览器中没有关掉：研究者称全国论文买卖去年销售额近10亿。初看这报道，心里弱弱地念了一句“骂了隔壁的”，你说说，这是谁在逼谁，这又是何苦要逼死这些“作者”们。难以理解。我觉得世上难以理解的事情只有两种，一种是纯粹的2，一种是精明之极。此处不展开。
之所以今天才写这事，主要是昨晚遇到了类似的事。有些老板要发论文，就逼学生分析数据，分析之前的结论都想好了，你就照着这个结论分析吧，还得人模狗样参考英文论文，论文三页纸，英文参考文献二三十篇。学生被逼急了只能造假，懂统计的可以高级造假（比如删掉几个数据使得检验显著），不懂统计的就低级造假（纯粹编假数）。老板可能也是被逼的，没论文没职称没钱没地位。经济方面的论文，编就编吧，反正大家都知道是假的，造个假数对大家都没影响；可这医学方面的论文，造数是不是不大好呢？如果论文跟治病救人没关系，那发论文就是堆垃圾了，何必要逼人发表；如果有关系，那这作者们良心何在？
回到我在统计之都新年构想中关于主站的目标一节：为什么期刊有存在的必要？为什么世上只有发表论文这一种指标来衡量人的工作和贡献？论文这个泥坑，学者有学者的痴狂，南郭先生有南郭先生的狡黠。跟买房一样，群体非理性，全然不顾是谁在背后蘸着口水数钱。
统计这玩意儿，一日不形成“reproducible”的规则，一日研究不成大器。
最后看个无关的短片，看什么叫“彪悍的人生不需要解释”：
对他这样的人，有没有必要用论文证明什么呢？
par(mar = c(4, 4, 0.5, 0.5), mfrow = c(1, 2))
set.seed(123)
plot(pval <- t(replicate(1000, {
x1 = rnorm(100, mean = 0, sd = runif(1, 0.5, 1))
x2 = rnorm(100, mean = 1, sd = runif(1, 2, 5))
c(t.test(x1, x2, var.equal = TRUE)$p.value, t.test(x1, x2,
var.equal = FALSE)$p.value)
})), xlab = "P-value: equal variance", ylab = "P-value: unequal variance",
pch = 20, asp = 1)
abline(0, 1)
plot(pval[, 1], pval[, 2] - pval[, 1], xlab = "P-value: equal variance",
ylab = "Diff of p-values (unequal var - equal var)", pch = 20)
过程是：从两个正态总体中生成样本，第一个总体均值为0，标准差随机取自U(0.5, 1)，第二个总体均值为1，标准差取自U(2, 5)，显然两个总体标准差不相等，那么在t检验时设定和不设定方差相等的选项对结果有多大影响？把两种情况的P值都画出来：左图是原始P值，可见基本在对角线上，说明大致相等，若眼神儿不好，可看右图，即P值的差异，可见方差不等时P值偏大（原因很简单，因为Welch校正的自由度小于等于不校正的自由度，样本量相等的时候统计量的分母即标准误一样，因此统计量完全一样，自由度越小，P值越大嘛），但大多少呢？其实也没大多少。
方差齐与不齐时t检验的结果对照
事情缘起于段炼同学9天前给我看的他的一篇博客：统计数字是不是拍脑袋出来的？87.53%。当时我在考试，没太仔细琢磨这件事情；现在邮件处理到了这一封，于是一层一层链接都打开来看，越看越摇头。这统计学在大家眼中敢情成了找借口的高级工具？抑或凡是有不正常的数字现象，都可以找到可能的“统计学”原因？这也太杯具了。
这个87.53%已经被证实只是个玩笑。在众多（只顾怀疑、相互抄袭、转载、或来路不明的）博客文章中，段炼的角度显然和所有人都不一样，他把所有的百分比数据的搜索频数都下载了下来，大家一看就知道，87.53这个数字本身并没有什么奇怪的，你去搜87.52或87.54都一样。众人纷纷解释这个0.53（100人中哪里来的0.53个人），不知道谁第一个提起了置信区间，总之我刚才看到的杯具有（考虑了一下，不是啥好事，就不给链接了）：
……在计算样本容量的时候要考虑一个置信区间的问题，也就是说调查了100个人，但是并不认为这100个人都是认真作答的，因此会在样本容量上再乘上一个置信度。
置信区间展现的是这个参数的真实值有一定概率落在测量结果的周围的程度。
第一种说法简直错了十万八千里，我闻所未闻，真是木有想到，置信度原来还有这种功效；第二种说法是对置信区间常见的误解；我正欲吐血时，竟然看见了维基百科的身影：置信区间。这下是真的杯具了，维基上赫然写着：
新华网关于学生冬季长跑的调查结果让人着实跌眼镜，一共调查了100人，报告中的结果都是xx.xx%形式的，例如“92.79%的学生认为强健了自己的身体”。这0.79个人是怎么来的？
咱们学统计的，应该对数字有一定的敏感性，比如当你看到小数位中含有667这样的数字（e.g. 0.291667）时就应该警觉：对方是否给出了样本量？如果没给的话，你就应该怀疑这个数字本来是0.29166666……如果你不知道这个比例是怎么来的，那么就拿一些整数去乘这个比例，看看哪个数字乘以这个比例能得到整数。最终你发现是24的倍数，样本量是7的倍数。然后你再想，7/24、14/48、28/96、……这一系列数字哪对更符合这个调查的背景。如：若你怀疑调查者很懒，那么不妨猜测他/她就调查了24个人。
以上只不过是低级的数字游戏，对统计来说根本没派上用场，现在很多人都琢磨着怎么建个模型整个P值去忽悠答辩委员会，而事实往往是，费尽千般心思，辛辛苦苦调查来的数据在建模之后根本没法用，要么系数是反的，要么不显著，或者有自相关，或有异方差，总之和初衷很不符，此时，离答辩往往只剩下几个星期，怎么办呢？只好眼睛一闭心一横，改数据吧！怎么改呢……【此处省略八千字】最后，王子和公主们过上了幸福生活。
我一般不相信经济学论文中的统计模型，原因之一就是数据。
近期评论