WinBUGS在统计分析中的应用Word格式.docx

1、下面我将以WinBUGS的统计应用为题，分几次来谈一谈WinBUGS这个软件。其中会涉及到空间数据的分析、GeoBUGS的使用、面向R及SPLUS的接口包R2WinBUGS的使用、GIS与统计分析等等衍生出的话题。如有问题，请大家留下评论，我会调整内容，择机给予回答。第一节 什么是WinBUGS?WinBUGS对于研究Bayesian统计分析的人来说，应该不会陌生。至少对于MCMC方法是不陌生的。WinBUGS （Bayesian inference Using Gibbs Sampling）就是一款通过MCMC方法来分析复杂统计模型的软件。其基本原理就是通过Gibbs sampling和Me

2、tropolis算法，从完全条件概率分布中抽样，从而生成马尔科夫链，通过迭代，最终估计出模型参数。引入Gibbs抽样与MCMC的好处是不言而喻的，就是想避免计算一个具有高维积分形式的完全联合后验概率公布，而代之以计算每个估计参数的单变量条件概率分布。具体的算法思想，在讲到具体问题的时候再加以叙述，在此不过多论述。就不拿公式出来吓人了（毕竟打公式也挺费劲啊）。第二节 为什么要用WinBUGS?第一、因为同类分析软件中它做得最好。同类的软件：OpenBUGS、JAGS等在成熟度、灵活性以及兼容性方面和它相比还有一定距离。在处理spatial data的方面，它采用了Gibbs抽样和MCMC的方法，

3、在模型支持方面又具有极大的灵活性，较之名声大噪的GeoR包，虽然也实现了bayesian的手法，但是灵活性还是不及WinBUGS。第二、因为它免费。免费的东西总有吸引人之处。第三、有各色的R包为WinBUGS实现了针对R的、SPLUS的、Matlab的软件接口。只要你喜欢，就直接在R下调用WinBUGS吧，无非是装个R2WinBUGS包，简简单单。第四、详细的文档、帮助、指导、范例。当然没有中文版的，小小一点遗憾。第三节 如何得到WinBUGS?WinBUGS目前是一款免费的软件，去http:/www.mrc-bsu.cam.ac.uk/bugs/下载就好了。不过要用高级功能（如GeoBUGS

4、）的话，还是去http:/www.mrc-bsu.cam.ac.uk/bugs/winbugs/contents.shtml注册一下好了，挺方便的。系统会立即把注册码发到你邮箱（真是好人啊）。不过只可以用一个月。这倒无妨，到时在注册一下就好了。第四节 初试WinBUGSWinBUGS-GUI我们先找一个例子来实际地运行一下WinBUGS，感受一下基本的操作流程，然后再考虑高级的操作。第一步，打开WinBUGS，通过菜单File - New新建一个空白的窗口第二步，在第一步中新建的空白窗口中输入三部分内容：模型定义、数据定义、初始值定义（代码见附录）第三步，点击菜单Model - Specifi

5、cation，弹出一个Specification Tool面板。第四步，在第二步中的提到的那个窗口中，将model这个关键字高亮起来，点击check model。你会看到WinBUGS的左下角状态栏上显示”model is syntactically correct.”第五步，把定义的data前的关键字list也高亮起来，点Specification Tool面板上的load data第六步，改Specification Tool面板上的马尔科夫链的数目，默认为1就好了第七步，点击Specification Tool面板上的compile第八步，把定义的初始值中的list关键字也高亮起来，再点

6、击Specification Tool面板上的load inits第九步，关了Specification Tool面板第十步，点击菜单Inference - Samples，弹出一个Sample Monitor Tool面板。第十一步，在Sample Monitor Tool面板的node中填要估计的参数名，这里可以是tau, alpha0, alpha1, b, 把它们一个一个填在node中，逐一点set。第十二步，关了Sample Monitor Tool面板第十三步，点击菜单Model - Update，弹出一个Update Tool面板。第十四步，将Update Tool面板中的upd

7、ates改大点，比如50000，点update按钮。第十五步，运行完后，关了Update Tool面板第十六步，点击菜单Inference - Samples第十七步，在弹出的Sample Monitor Tool面板上选一个node第十八步，点history看所有迭代的时间序列图，点trace看最后一次迭代的时间序列图，点auto cor看correlogram时间序列图，点stat看参数估计的结果Estimation results by WinBUGS 1.4附第二步中的代码如下：#MODELmodel for （i in 1:N） Oi dpois（mui） log（mui） - lo

8、g（Ei） + alpha0 + alpha1 * Xi/10 + bi # Area-specific relative risk （for maps） RRi - exp（alpha0 + alpha1 * Xi/10 + bi） # CAR prior distribution for random effects: b1:N car.normal（adj, weights, num, tau） for （k in 1:sumNumNeigh） weightsk - 1 # Other priors: alpha0 dflat（） alpha1 dnorm（0, 1e-05） tau d

9、gamma（0.5, 5e-04） # prior on precision sigma - sqrt（1/tau） # standard deviation#DATAlist（N = 56, O = c（9, 39, 11, 9, 15, 8, 26, 7, 6, 20, 13, 5, 3, 8, 17, 9, 2, 7, 9, 7, 16, 31, 11, 7, 19, 15, 7, 10, 16, 11, 5, 3, 7, 8, 11, 9, 11, 8, 6, 4, 10, 8, 2, 6, 19, 3, 2, 3, 28, 6, 1, 1, 1, 1, 0, 0）, E = c（1.

10、4, 8.7, 3, 2.5, 4.3, 2.4, 8.1, 2.3, 2, 6.6, 4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4, 10.5, 22.7, 8.8, 5.6, 15.5, 12.5, 6, 9, 14.4, 10.2, 4.8, 2.9, 7, 8.5, 12.3, 10.1, 12.7, 9.4, 7.2, 5.3, 18.8, 15.8, 4.3, 14.6, 50.7, 8.2, 5.6, 9.3, 88.7, 19.6, 3.4, 3.6, 5.7, 7, 4.2, 1.8）, X = c（16, 16, 10,

11、24, 10, 24, 10, 7, 7, 16, 7, 16, 10, 24, 7, 16, 10, 7, 7, 10, 7, 16, 10, 7, 1, 1, 7, 7, 10, 10, 7, 24, 10, 7, 7, 0, 10, 1, 16, 0, 1, 16, 16, 0, 1, 7, 1, 1, 0, 1, 1, 0, 1, 1, 16, 10）, num = c（3, 2, 1, 3, 3, 0, 5, 0, 5, 4, 0, 2, 3, 3, 2, 6, 6, 6, 5, 3, 3, 2, 4, 8, 3, 3, 4, 4, 11, 6, 7, 3, 4, 9, 4, 2,

12、4, 6, 3, 4, 5, 5, 4, 5, 4, 6, 6, 4, 9, 2, 4, 4, 4, 5, 6, 5）, adj = c（19, 9, 5, 10, 7, 12, 28, 20, 18, 19, 12, 1, 17, 16, 13, 10, 2, 29, 23, 19, 17, 1, 22, 16, 7, 2, 5, 3, 19, 17, 7, 35, 32, 31, 29, 25, 29, 22, 21, 17, 10, 7, 29, 19, 16, 13, 9, 7, 56, 55, 33, 28, 20, 4, 17, 13, 9, 5, 1, 56, 18, 4, 50

13、, 29, 16, 16, 10, 39, 34, 29, 9, 56, 55, 48, 47, 44, 31, 30, 27, 29, 26, 15, 43, 29, 25, 56, 32, 31, 24, 45, 33, 18, 4, 50, 43, 34, 26, 25, 23, 21, 17, 16, 15, 9, 55, 45, 44, 42, 38, 24, 47, 46, 35, 32, 27, 24, 14, 31, 27, 14, 55, 45, 28, 18, 54, 52, 51, 43, 42, 40, 39, 29, 23, 46, 37, 31, 14, 41, 3

14、7, 46, 41, 36, 35, 54, 51, 49, 44, 42, 30, 40, 34, 23, 52, 49, 39, 34, 53, 49, 46, 37, 36, 51, 43, 38, 34, 30, 42, 34, 29, 26, 49, 48, 38, 30, 24, 55, 33, 30, 28, 53, 47, 41, 37, 35, 31, 53, 49, 48, 46, 31, 24, 49, 47, 44, 24, 54, 53, 52, 48, 47, 44, 41, 40, 38, 29, 21, 54, 42, 38, 34, 54, 49, 40, 3

15、4, 49, 47, 46, 41, 52, 51, 49, 38, 34, 56, 45, 33, 30, 24, 18, 55, 27, 24, 20, 18）, sumNumNeigh = 234）#INITIAL VALUESlist（tau = 1, alpha0 = 0, alpha1 = 0, b = c（0, 0, 0, 0, 0, NA, 0, NA, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0

16、, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0）WinBUGS在统计分析中的应用 第一部分完1. 胡江堂 on 2008/12/09 at 18:03 WinBugs，头一回听说啊。老齐在annpro，有空分享一下工作心得吧。回复 2. 刘思喆 on 2008/12/10 at 11:13 土了，还以为是 windows 下的 debug 工具。原来是这 ！3. 谢益辉 on 2008/12/15 at 20:14 呵呵，能写点公式还是写点吧，看代码不知道做的问题是什么4. 谢益辉: 统计之都本周导读第三辑 | 统计之都 on 2008/12/15 at 20:23 . WinB

17、UGS是贝叶斯统计的有力工具（而不是Windows下面的Debug工具），齐韬加入COS主站之后发表了第一篇文章WinBUGS在统计分析中的应用（第一部分），为我们讲述了WinBUGS的一些基本操作。 .5. 齐韬 on 2008/12/16 at 11:45 我会在第二部分中将相关的理论部分加上，看来得要学习一下怎么嵌入LaTex了.6. 王化儒 on 2008/12/16 at 16:57 嗯，楼主辛苦了！跟着一点一点学，期待着空间数据分析。7. 胡江堂 on 2008/12/18 at 11:09 贴一首诗吧，刚BUGS team发过来的，非常有意思：Each year you wait

18、 with bated breath,The old WinBUGS key nearing death,And will the brand new key appearIn time to join the festive cheer?The waitings over raise your glass,And drink to rituals that pass.Relax, sit back and have a chortle;This time your WinBUGS keys immortal.8. xjx on 2008/12/27 at 18:18 不知道在linux下能不

19、能用?9. hong on 2009/04/25 at 21:55 大家好麻烦问一下这个迭代次数如何选择.谢谢回答10. DJ on 2009/05/23 at 17:50 很好的教程，谢谢了。MCMC万岁！11. 左伊秩訾 on 2009/05/26 at 22:35 前几天上课时听说了这个软件，真是及时雨！谢谢了！12. icwei on 2009/06/01 at 21:小弟现正在剑桥mrc-bsu做postdoc，具体的项目就是BUGS的开发以及在生物及医学方面应用。WinBUGS这个软件是我两个老板David Spiegehalter,Dave Lunn 和其它一些牛人共同开发的。我

20、们现在正在从WinBUGS 转向openBUGS,目的是将它做成open source的软件以应用在更广的领域。 我现在正在开发BUGS中的WBDiff部分并将它应用在二型糖尿病的动态系统的数据分析中。有兴趣或有问题的同学可以和我联系：chen.weimrc-bsu.cam.ac.uk还有我们这里每年会举办3-4次BUGS的培训，2天的课程，在英国或能到英国出差的同学有兴趣的话可以参加，主讲人是David speigehalter 和 Dave Lunn。o DJ on 2009/06/07 at 22:08 楼上的大牛人啊！o maple on 2009/11/09 at 17:这位仁兄真牛

21、o 海涛 on 2009/11/14 at 10:53 请问在winbugs中我编译FRAILTY模型时提示educational version cannot do this model 难道这个软件还有别的版本？我用的版本是1.4我编译的模型是帮助文件Examples Volume 1中的最后一个文件Cox regression with random effects13. zyn on 2009/06/16 at 15:42 希望能和您联系,我这边看了gibbs sampling有不少问题，不知道您可否提供帮助？14. phlissia on 2009/09/18 at 11:10 弱弱

22、的问一句，执行到第四步check model的时候，winbugs坐下角显示的不是“model is syntactically correct”，而是在alpha1 dnorm（0, 1e-05）一句1e处显示e处应为”expected right parenthesis，在代码最后0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0）指示”invalid or unexpected token scanned”，这要怎么弄呢？ 第一次接触这个软件，还不会使，请高手指点。WinBUGS在统计分析中的应用（第二部分）By 齐韬 2008/12/18

23、 Bayesian Analysis, SAS, WinBUGS 分类：数据分析, 统计软件, 贝叶斯统计齐韬. WinBUGS在统计分析中的应用（第二部分）. 统计之都, 2008.12. URL:/cos.name/2008/12/statistical-analysis-and-winbugs-part-2/. TITLE = WinBUGS在统计分析中的应用（第二部分）,/cos.name/2008/12/statistical-analysis-and-winbugs-part-2/,第一节 WinBUGS数据分析案例在这一节中，我将拿一个经典的研究数据，利用WinBUGS给出简单的

24、分析。首先介绍一下这个数据：Seedsseed O. aegyptiaco 75 seed O. aegyptiaco 73 Bean Cucumber Bean Cucumberr n r/n r n r/n r n r/n r n r/n10 39 0.26 5 6 0.83 8 16 0.5 3 12 0.2523 62 0.37 53 74 0.72 10 30 0.33 22 41 0.5423 81 0.28 55 72 0.76 8 28 0.29 15 30 0.526 51 0.51 32 51 0.63 23 45 0.51 32 51 0.6317 39 0.44 46

25、79 0.58 0 4 0 3 7 0.4310 13 0.77这个数据来自Crowder （1978）。之后Breslow and Clayton （1993） 作为例子，也分析过这个数据。数据反映的是某一品种的豆类种子和某一品种的黄瓜种子，分别放在21个培养皿（plates）中分别培养，在根提取物aegyptiaco 75和aegyptiaco 73的作用下的出芽率的差异。其中r是出芽的个数，n是种子的个数，而r/n是出芽率。我们用random effect logistic regression模型来进行分析（注意，在Bayesian分析中，通常是将covariates看做是服从某一个分

26、布的随机变量，这和一般意义上的GLM, GLME, LME中对于covariates解释是不同的）：其中是种子的类型，是根提取物的类型，是交互项, 是给定的独立的 “noninformative” 先验参数。在Bayesian分析中，通常我们会定义一个DAG图（即Directed Acyclic Graph有向无圈图） 。我们可以在WinBUGS中通过设计DAG图来定义模型。不过这一节中我们还是用WinBUGS中的BUGS语言来定义模型，如何在WinBUGS中通过设计DAG图来定义模型我将在下一节中详细介绍，但是必须要说明的是BUGS语言比DAG图灵活，不过直观性不如后者。模型 for（ i

27、in 1 : N ） ri dbin（pi,ni） bi dnorm（0.0,tau） logit（pi） - alpha0 + alpha1 * x1i + alpha2 * x2i + alpha12 * x1i * x2i + bi alpha0 dnorm（0.0,1.0E-6） alpha1 dnorm（0.0,1.0E-6） alpha2 dnorm（0.0,1.0E-6） alpha12 dnorm（0.0,1.0E-6） tau dgamma（0.001,0.001）- 1 / sqrt（tau）WinBUGS doodle模型数据list（r = c（10, 23, 23, 26, 17, 5, 53, 55, 32, 46, 10, 8, 10, 8, 23, 0, 3, 22, 15, 32, 3）, n = c（39, 62, 81, 51, 39, 6, 74, 72, 51, 79, 13, 16, 30, 28, 45, 4, 12, 41, 30, 51, 7）, x1 = c（0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1）, x2 = c（0, 0, 0, 0, 0