N岸N地雾霾之争已经持续了一个星期。相关的博文我写了好几篇发在珍珠湾与科学网,按理说应该是非常透彻了。。。。所以我来做一个简单的回顾。。。
在这篇博文中,我写到【那张图显示的是每年死亡率增加的百分比 (Percent INCREASE),而不是死亡率本身,图中数据只要是大于零,就说明死亡率在增加。因此,原图表明,死亡率确实是在逐年递增。】,【PM2.5 吸入剂量 Z 随时间线性增长,而死亡率在2006-2009期间也几乎是线性增长。】
实际死亡的数据证实了2006-2009年心血管疾病死亡率几乎线性增长的判断。以北京城区为例,2008的循环系统病死亡率比2007年增长了 9%。
参见下图(北京数据与全国数据)
2. 对 《方舟子妄批柴静捏造数据的错误》的补充
以上结论是在我没有看到论文,而只是看到上图6的情况下做出的。实际上论文图6的文字描述有误(或者疏漏),在科学网的讨论中,网友熊兵提供了论文的摘要,根据这个摘要,我进行了 PS2 补充与纠正:【原论文图的正确描述应该是【The inter-annual variability of the estimated percent increases in daily mortality ASSOCIATED WITH 10 ug/m^3 INCREASE of PM2.5 in recent years. 】
原论文图下的描述少了大写的部分。也就是说,原图是 PM2.5 值每增加10,死亡率增加的百分比。
3. 方舟子的回应
方舟子在我的博文之后,改变了原来认为曲线是描述死亡率的说法。他改口说:”与基数相比,2007年北京循环系统疾病每日死亡率增加2%,而 2008年增加1.5%,那么2008年与2007年相比,死亡率实际上是下降的,绝非像岳东晓理解的那样每年死亡率都在直线上升,没过几年北京人都死绝 了。“
首先,方舟子的理解与实际数据完全不符合。上面的实际数据显示,2008年的心血管死亡率较2007年增长9 %,而不是像方舟子想象的那样下降。而中科院这篇论文是根据实际死亡数据做的分析。
方舟子没有说他理解的基数是什么,但根据他误以为2008年死亡率下降的情况,应该是以为“基数”是某个固定的数字。这显然是错误的(参见我写的《详解方舟子没看懂的“柴静曲线"》。
至于方舟子说【每年死亡率都在直线上升,没过几年北京人都死绝 了】,缺乏概念。直线上升不等于垂直上升。直线上升的快慢取决于曲线的坡度。
目前CM死亡率为 1/500,如果按每年10%的增长率直线(或者指数)增长,几十年后才能达到 1 %, 而中国出生率 是 1.2%。即使CM疾病死亡率持续保持增长,几百年内是死不光的。具体计算留给大家作为练习了。
另外,我们说的是2006-2009年死亡率直线增长(这也是实际数据),显然这个直线增长不可能无限的保持。
再者,我补充一句,虽然我在前文中论证了方舟子没脑子,但是方舟子好歹是福建云霄的高材生,考入中科大的,也是正规的 Ph.D., 虽然学生物不要多少数学,但这点概念应该有,更可能是没过脑子。
中科院
论文使用的模型公式是
其中 E(Y )代表死亡人数期待值,X代表PM指数 (如PM 2.5、PM10),PM指数X 前β是需要通过分析确定的系数。
通过简单的运算,我们看出 β=1Y∂Y∂X。
具体的方法称为 Generalized Additive Model (GAM),这是1990年代之后才开始广泛使用的统计分析方法。我们大多数人应该没有学过。这次为了完全理解这篇论文,我阅读了一本这方面的专著、以及相关的编程,算是掌握了一些 working knowledge。这也是进行这种讨论的好处。可以学到新东西。网上那些看了本人的博文之后,只能学舌的,显然不知道这些深度知识,也就只能装了。
5. 中科院论文图6 的通俗讲解
β=1Y∂Y∂X的确定不是以什么基数算百分比。论文依据的数据仅仅是每天的死亡人数与PM指数,根本不可能知道哪些人是因为雾霾而死,哪些人不是因为雾霾而死。
而且
这个死亡人数也是每天都在变化的,而不是固定值。 具体怎么计算可以说相当复杂,有兴趣的应该学习
Generalized Additive Model 的知识。这里就不再详述了。
如果用通俗的语言来讲,这张曲线表达的是什么呢? 科学网我的博客下一位读者提出毒力的概念,很生动。
用通俗的话说,论文图6的纵轴是雾霾的毒性大小。毒性大小是什么意思? 我们说氰化钾与砒霜毒性大,是说同样重量的氰化钾能够毒死的人比砒霜毒死的人多。这里,雾霾的毒性是指单位 PM2.5 能造成的杀伤力(这里是死亡率的增率)。这张图的横轴是时间。
因此,这张图描述的是雾霾毒性随时间的变化。从图中的趋势可以看出,除了2008年雾霾毒性下降外,雾霾的毒性总的来说是在增长的。这也是该论文的一个结论。
至于为什么雾霾毒性还在增强,这当然不是这个统计分析能够解答的。
我在之前的博文中论述了雾霾的积累效应。有时间,我把相关的数学模型写一下。
6. 柴静(或者其团队)对图6的理解正确
论文结论是:【2005-2009的时序分析表明PM2.5值每增加10,非事故死亡率增加 0.65%, 呼吸系统疾病死亡率增加0.63%, 心血管病死亡率增加 1.38%】
柴静在视频中说道:【当pm2.5值升高的时候,人群的死亡率是随之上升的。】这完全符合上述论文的结论 --- 该论文都算了,PM值升高10,死亡率增加多少。
柴静在另一处说: 【这是中科院做的一个测算,他们有一个简单的结论,pm2.5上升到200微克每立方米的时候,我们各种呼吸系统疾病,和心血管系统的疾病死亡率会 增加14%到26%】
这应该是根据这篇论文图6的计算结果,1.38% * 190/10 = 26.2%, 0.63* 190 /10 = 12 %。与柴静说的数字接近---误差可能是近似。这可以说明柴静(或者其团队)不但用了图6,还用图6 的数据进行了假设性死亡率计算。
到此,大家应该明白,原图中最上面的黑色柱子只是记载PM值的参考信息,完全可以去掉,对于【当pm2.5值升高,人群的死亡率是随之上升的】的结论没有任何影响。这些年中,平均PM2.5指数在 65到84之间,相差30%。柴静陈述的观点中没有任何一点依赖于这个PM2.5指数的变化,也没有任何一点加入这个变化就不成立。有人认为这些PM2.5平均值是【模式计算疾病死亡率的初始值】因而很重要,完全是一窍不通。
整个计算是根据北京每天的PM2.5值,整个计算是直接输入全部数据以及拟合函数(包括统计分布类型-- 此处为 Poisson),调用 mgcv中的gam 函数。如果做过这样的练习就知道,根本没有要先算出年平均 PM2.5 。
It can be seen that the PM2.5 levels generally declines in recent years in Beijing, whereas the relative risk shows different change trends. Overall, the level of the estimated percentage increase assumes an escalating tendency
during the study period... it bounced off the bottom and started a new uptrend after 2008. 】 我拨冗翻译如下:【从图中(黑柱)看出,近年来北京的 PM2.5 指数一般呈下降趋势,但是雾霾的相对危险度的趋势却相反。总的来说,在这段时间单位PM2.5造成的死亡率增加呈加速上升趋势。。。这个危险度从2008年的低谷反弹从那以后开始一个新的上升趋势】。
由上可见,黑色柱子是为了说明这一点,虽然PM2.5下降了,但是雾霾的毒性却加强了。除非柴静要详细讲解雾霾毒性随时间加强这一点,这个黑色的柱子是多余的,只会让人把本来就已经相当费解(方舟子都完全看不懂)的图搞得更加难懂。
后来我发现看到柴静视频上实际上列出了这条曲线的系数,其实这也是一条三次曲线,但三次项系数很小。柴静的视频还给出了 R^2 值。这个的R^2值明显偏小。有统计学知识的也许可以在这个问题上进行动机判断。】
R^2 = 1 - RSS/TSS
柴静这条曲线的 R^2 只有 0.19 ,这意味这条曲线并未与四个数据点得到很好的拟合(而另外两条曲线的R^2 都接近1)。这一点也是显然的,三次曲线有四个独立参数,应该被四个数据点完全确定,而红色曲线没有通过四个点。
R^2 值小固然是一个疑问,但是由于拟合并不影响结论,我个人认为很难做什么判断。】
曲线是显示雾霾毒性随时间的变化,柴静的演讲中并没有涉及这个毒性变化的问题,参见上面的 6、7。
CONCLUSION
上面我就几点进行的简单陈述,足以证明科学来不得虚假。