数据分析
数据分析示例
近日来收到邮件少了，但各个问题都不太好直接回答。比如这则关于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检验统计量与卡方分布拟合的好坏
本人怀着阴暗心理于今日公交车上翻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个人都是认真作答的，因此会在样本容量上再乘上一个置信度。
置信区间展现的是这个参数的真实值有一定概率落在测量结果的周围的程度。
第一种说法简直错了十万八千里，我闻所未闻，真是木有想到，置信度原来还有这种功效；第二种说法是对置信区间常见的误解；我正欲吐血时，竟然看见了维基百科的身影：置信区间。这下是真的杯具了，维基上赫然写着：
如果你想掩盖数据，那么就把它们离散化吧！
不知道为什么这么多人钟爱于将连续数据离散化，例如明明有年龄数据，在分析的时候非要分成老幼青壮这样的分类变量；明明有原始的计数数据，非要搞成“0-5、6-10、……”这样的频数表。大概是数据得来不花钱吧，这样毁灭信息一点都不心疼。
某年我在某医学统计会议上专门强调了这个愚蠢的问题，结果后面还有某小师妹没理解我的意思，把我批驳了一番，依然支持离散化，我无语，只能摇摇头叹口气。去年useR! 2008会议上，Frank Harrell也提到了这个问题，他也想不通，为什么人们喜欢离散化。
如果你问一位lady：请问姑娘芳龄多少哇？姑娘回答：臣妾属于0~100岁这一组的。我想，此时这些人该能理解离散化的毛病所在了吧。
刚看到这么一则报导：Recipe for Disaster: The Formula That Killed Wall Street
话说David X. Li的高斯Copula模型毁了华尔街，进而导致全面金融危机，真是好笑，老夫甚为幸灾乐祸。看看Li的维基页面，其中有这么一句旁白：
People got very excited about the Gaussian copula because of its mathematical elegance, but the thing never worked. Co-association between securities is not measurable using correlation.
窃以为数学的追求是多种角度的美感，包括理性、简洁、严谨等，上面第一句话的前后两半句的鲜明对比让我狂笑不已。昨晚我刚刚提到John W. Tukey在四十多年前的文章“The Future of Data Analysis”，就说了他关于“数学不是科学”的观点，原因是数学不需要实践检验。
我一直对纯数学理论的研究不感兴趣，而很多人则认为把这些东西装在脑子里才像“贵族”一些。传说古希腊某哲人给学生上课时，一位学生提出异议问这些理论有什么用，该哲人叫人去拿两个铜板给他，说你这么看重实用，那就给你两个铜板好了。
这件事情Li本人是无辜的，他知道数学模型的局限性，只不过“一小撮”人过于崇拜（他们看不懂的）模型了。
想当年，老夫对Copula也是非常感兴趣，也被它的简洁所吸引。各位客官若有兴趣，可以任意用百度或者Google搜索一下“Copula+论文”，第一页就可以看见我在COS论坛发的帖子。
好像我越是打击结构方程模型，找我咨询的人越多，不知怎么搞的（注1）。莫非搜索引擎不知道怎么分析我的语义，而仅仅只是看我的网站中有没有“结构方程模型”这六个字？
话说前两天发生的一件事情：某高校一位教授找到我，希望我帮忙算一算SEM，我本不太情愿干这种和SEM有关的事，有碍于是导师大人介绍的，也就答应下来，哪知前天突然得知，他们圈内有规定，SEM必须用LISREL求解，不能用AMOS，论文必须附注LISREL代码。
娘希匹，迷信SEM建模本身的方法论就忍了，居然对软件都有宗教般的信仰，你们说说，这帮人是不是脑子让猪栏门给夹了？（注2）
注1：本忘了写这件事，刚才又有人发邮件问，“您在验证模型出现过五个拟合指数只出现其中的两三个不全的情况吗？”我很想知道，这“五个拟合指数”是什么意思，难道又是某大人规定了SEM的五大金律？也许是我抵触情绪太大，想歪了吧。
注2：“脑子让门夹了”是很久以前网络上流行的一句话，我前年寒假回家，在家听到一个农村版本，就是“脑子让猪栏门夹了”，当即忍俊不禁。
根据数据检验总体的分布在我看来几乎没有什么用处，不过历史上已经出现了无数种关于分布的检验，例如Chi-square检验、KS检验、Shapiro正态性检验等等。我觉得检验没有实际用处的原因有二：
一、若拒绝零假设，即数据不服从某种分布，那么往往会使得下面要做的工作的前提假设不成立——这显然会很惨；
二、若不拒绝零假设——这几乎是无用的结论，因为不拒绝这个零假设，不代表能拒绝其它零假设，因此你仍然不知道数据是什么分布——这显然更惨；
所以我们要把自己的眼睛捂上，假装看不见，像数理统计学家那样，我们假定X服从帕累托分布，然后咋地咋地。
附1：本文是回答Renxiang Yan同志的邮件，因为我写了好半天，然后点“Reply”，发现从163返回了错误信息，说收件人拒绝收信，真是气人。我也不知道是Yan同学自己邮箱设置有问题，还是163有时候会拒绝Gmail的邮件。估计后者可能性大一些。
附2：还要补充说明一点，关于分布的假设检验中，零假设往往是确定的分布，而不是带有未知参数的笼统的分布，即分布的参数都是确定的值。只有少数几个关于正态分布的检验除外，因为它们有渐近性质。因此，提问时最好不要抽象地问怎么检验样本是否是广义极值分布。
从论坛上最近的几个帖子来看，貌似相当一部分人的脑子已经被正态分布严重“毒化”了。例如，有人问，若总体不是正态分布能否求均值和标准差（作为描述统计量）？还有人问，新闻标题的字符长度服从什么分布（有意用正态分布），下面回复中也几乎都是正态分布的天下。为什么我们如此教条地对待正态分布？
我猜想，一个致命的原因就是数理统计的理论发展有相当大一部分都是建立在正态分布的根基上，这与实际应用中的统计需求存在明显的矛盾。统计分析并不同于找对象（在心里想好了要相亲的对象的性格、外貌、喜不喜欢猫、是否爱吃辣等“理论框架”之后再赴约），而是带有探索的意味。以我愚见，这可能是John Tukey在”The Future of Data Analysis”一文中强调的重点之一，不幸的是，我们至今仍然把数理统计高高供奉在统计学的神坛上，甘愿成为“正态教”的信徒。
作为实例，下面是数理统计版中100个帖子标题的字符长度，感兴趣的客官不妨琢磨琢磨它的分布：
20 20 13 15 2 11 31 10 12 20 13 56 7 13 19 46 16 19 14 9 20 10 22 13 2 43 11 15 20 14 26 10 19 33 15 15 65 7 16 18 10 32 14 17 14 24 19 60 13 17 27 7 12 7 11 70 50 8 13 8 15 2 20 27 39 7 7 26 21 19 22 8 26 42 8 17 37 17 5 14 21 8 28 18 69 12 23 12 17 14 17 8 20 31 36 25 20 6 6 11
近期评论