WinBUGS在统计分析中的应用.docx
《WinBUGS在统计分析中的应用.docx》由会员分享，可在线阅读，更多相关《WinBUGS在统计分析中的应用.docx（34页珍藏版）》请在冰豆网上搜索。
WinBUGS在统计分析中的应用
WinBUGS在统计分析中的应用（第一部分）
By齐韬@2008/12/08
关键词：
MCMC,R,WinBUGS,空间统计分类：
统计软件
作者信息：
ComputationalMathematicianinAnnproAnalyticTechnologies,Inc.
版权声明：
本文版权归原作者所有，未经许可不得转载。
原文可能随时需要修改纰漏，全文复制转载会带来不必要的误导，若您想推荐给朋友阅读，敬请以负责的态度提供原文链接；点此查看如何在学术刊物中引用本文
常规引用方式
齐韬.WinBUGS在统计分析中的应用（第一部分）.统计之都,2008.12.URL:
http:
//cos.name/2008/12/statistical-analysis-and-winbugs-part-1/.
BibTeX引用
@ARTICLE{,
AUTHOR={齐韬},
TITLE={WinBUGS在统计分析中的应用（第一部分）},
JOURNAL={统计之都},
YEAR={2008},
month={12},
URL={http:
//cos.name/2008/12/statistical-analysis-and-winbugs-part-1/},
}
开篇词
首先非常感谢COS论坛提供了这样一个良好的平台，敝人心存感激之余，也打算把一些学习心得拿出来供大家分享，文中纰漏之处还请各位老师指正。
下面我将以WinBUGS的统计应用为题，分几次来谈一谈WinBUGS这个软件。
其中会涉及到空间数据的分析、GeoBUGS的使用、面向R及SPLUS的接口包R2WinBUGS的使用、GIS与统计分析等等衍生出的话题。
如有问题，请大家留下评论，我会调整内容，择机给予回答。
第一节什么是WinBUGS?
WinBUGS对于研究Bayesian统计分析的人来说，应该不会陌生。
至少对于MCMC方法是不陌生的。
WinBUGS（BayesianinferenceUsingGibbsSampling）就是一款通过MCMC方法来分析复杂统计模型的软件。
其基本原理就是通过Gibbssampling和Metropolis算法，从完全条件概率分布中抽样，从而生成马尔科夫链，通过迭代，最终估计出模型参数。
引入Gibbs抽样与MCMC的好处是不言而喻的，就是想避免计算一个具有高维积分形式的完全联合后验概率公布，而代之以计算每个估计参数的单变量条件概率分布。
具体的算法思想，在讲到具体问题的时候再加以叙述，在此不过多论述。
就不拿公式出来吓人了（毕竟打公式也挺费劲啊）。
第二节为什么要用WinBUGS?
第一、因为同类分析软件中它做得最好。
同类的软件：
OpenBUGS、JAGS等在成熟度、灵活性以及兼容性方面和它相比还有一定距离。
在处理spatialdata的方面，它采用了Gibbs抽样和MCMC的方法，在模型支持方面又具有极大的灵活性，较之名声大噪的GeoR包，虽然也实现了bayesian的手法，但是灵活性还是不及WinBUGS。
第二、因为它免费。
免费的东西总有吸引人之处。
第三、有各色的R包为WinBUGS实现了针对R的、SPLUS的、Matlab的软件接口。
只要你喜欢，就直接在R下调用WinBUGS吧，无非是装个R2WinBUGS包，简简单单。
第四、详细的文档、帮助、指导、范例。
当然没有中文版的，小小一点遗憾。
第三节如何得到WinBUGS?
WinBUGS目前是一款免费的软件，去http:
//www.mrc-bsu.cam.ac.uk/bugs/下载就好了。
不过要用高级功能（如GeoBUGS）的话，还是去http:
//www.mrc-bsu.cam.ac.uk/bugs/winbugs/contents.shtml注册一下好了，挺方便的。
系统会立即把注册码发到你邮箱（真是好人啊）。
不过只可以用一个月。
这倒无妨，到时在注册一下就好了。
第四节初试WinBUGS
WinBUGS-GUI
我们先找一个例子来实际地运行一下WinBUGS，感受一下基本的操作流程，然后再考虑高级的操作。
第一步，打开WinBUGS，通过菜单File->New新建一个空白的窗口
第二步，在第一步中新建的空白窗口中输入三部分内容：
模型定义、数据定义、初始值定义（代码见附录）
第三步，点击菜单Model->Specification，弹出一个SpecificationTool面板。
第四步，在第二步中的提到的那个窗口中，将model这个关键字高亮起来，点击checkmodel。
你会看到WinBUGS的左下角状态栏上显示”modelissyntacticallycorrect.”
第五步，把定义的data前的关键字list也高亮起来，点SpecificationTool面板上的loaddata
第六步，改SpecificationTool面板上的马尔科夫链的数目，默认为1就好了
第七步，点击SpecificationTool面板上的compile
第八步，把定义的初始值中的list关键字也高亮起来，再点击SpecificationTool面板上的loadinits
第九步，关了SpecificationTool面板
第十步，点击菜单Inference->Samples，弹出一个SampleMonitorTool面板。
第十一步，在SampleMonitorTool面板的node中填要估计的参数名，这里可以是tau,alpha0,alpha1,b,把它们一个一个填在node中，逐一点set。
第十二步，关了SampleMonitorTool面板
第十三步，点击菜单Model->Update，弹出一个UpdateTool面板。
第十四步，将UpdateTool面板中的updates改大点，比如50000，点update按钮。
第十五步，运行完后，关了UpdateTool面板
第十六步，点击菜单Inference->Samples
第十七步，在弹出的SampleMonitorTool面板上选一个node
第十八步，点history看所有迭代的时间序列图，点trace看最后一次迭代的时间序列图，点autocor看correlogram时间序列图，点stat看参数估计的结果
EstimationresultsbyWinBUGS1.4
附第二步中的代码如下：
#MODEL
model
{
for（iin1:
N）{
O[i]~dpois（mu[i]）
log（mu[i]）<-log（E[i]）+alpha0+alpha1*X[i]/10+
b[i]
#Area-specificrelativerisk（formaps）
RR[i]<-exp（alpha0+alpha1*X[i]/10+b[i]）
}
#CARpriordistributionforrandomeffects:
b[1:
N]~car.normal（adj[],weights[],num[],tau）
for（kin1:
sumNumNeigh）{
weights[k]<-1
}
#Otherpriors:
alpha0~dflat（）
alpha1~dnorm（0,1e-05）
tau~dgamma（0.5,5e-04）
#prioronprecision
sigma<-sqrt（1/tau）
#standarddeviation
}
#DATA
list（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.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,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,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,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,37,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,34,49,47,46,41,52,51,49,38,34,
56,45,33,30,24,18,55,27,24,20,18）,sumNumNeigh=234）
#INITIALVALUES
list（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,0,0,0,0,0,0,0,0,0,0））
WinBUGS在统计分析中的应用第一部分完
1.胡江堂on2008/12/09at18:
03
WinBugs，头一回听说啊。
老齐在annpro，有空分享一下工作心得吧。
回复
2.刘思喆on2008/12/10at11:
13
土了，还以为是windows下的debug工具。
原来是这！
！
回复
3.谢益辉on2008/12/15at20:
14
呵呵，能写点公式还是写点吧，看代码不知道做的问题是什么……
回复
4.谢益辉:
统计之都《本周导读》第三辑|统计之都on2008/12/15at20:
23
[...]WinBUGS是贝叶斯统计的有力工具（而不是Windows下面的Debug工具），齐韬加入COS主站之后发表了第一篇文章WinBUGS在统计分析中的应用（第一部分），为我们讲述了WinBUGS的一些基本操作。
[...]
回复
5.齐韬on2008/12/16at11:
45
我会在第二部分中将相关的理论部分加上，看来得要学习一下怎么嵌入LaTex了.
回复
6.王化儒on2008/12/16at16:
57
嗯，楼主辛苦了！
跟着一点一点学，期待着空间数据分析。
。
。
回复
7.胡江堂on2008/12/18at11:
09
贴一首诗吧，刚BUGSteam发过来的，非常有意思：
Eachyearyouwaitwithbatedbreath,
TheoldWinBUGSkeynearingdeath,
Andwillthebrandnewkeyappear
Intimetojointhefestivecheer?
Thewaiting’sover–raiseyourglass,
Anddrinktoritualsthatpass.
Relax,sitbackandhaveachortle;
ThistimeyourWinBUGSkey’simmortal.
回复
8.xjxon2008/12/27at18:
18
不知道在linux下能不能用?
回复
9.hongon2009/04/25at21:
55
大家好
麻烦问一下这个迭代次数如何选择.谢谢回答
回复
10.DJon2009/05/23at17:
50
很好的教程，谢谢了。
MCMC万岁！
回复
11.左伊秩訾on2009/05/26at22:
35
前几天上课时听说了这个软件，真是及时雨！
谢谢了！
回复
12.icweion2009/06/01at21:
35
小弟现正在剑桥mrc-bsu做postdoc，具体的项目就是BUGS的开发以及在生物及医学方面应用。
WinBUGS这个软件是我两个老板DavidSpiegehalter,DaveLunn和其它一些牛人共同开发的。
我们现在正在从WinBUGS转向openBUGS,目的是将它做成opensource的软件以应用在更广的领域。
我现在正在开发BUGS中的WBDiff部分并将它应用在二型糖尿病的动态系统的数据分析中。
有兴趣或有问题的同学可以和我联系：
chen.wei@mrc-bsu.cam.ac.uk
还有我们这里每年会举办3-4次BUGS的培训，2天的课程，在英国或能到英国出差的同学有兴趣的话可以参加，主讲人是Davidspeigehalter和DaveLunn。
回复
oDJon2009/06/07at22:
08
楼上的大牛人啊！
回复
omapleon2009/11/09at17:
14
这位仁兄真牛
回复
o海涛on2009/11/14at10:
53
请问在
winbugs中我编译FRAILTY模型时提示
educationalversioncannotdothismodel难道这个软件还有别的版本？
？
？
我用的版本是1.4
我编译的模型是帮助文件ExamplesVolume1中的最后一个文件Coxregressionwithrandomeffects
回复
13.zynon2009/06/16at15:
42
希望能和您联系,我这边看了gibbssampling有不少问题，不知道您可否提供帮助？
回复
14.phlissiaon2009/09/18at11:
10
弱弱的问一句，执行到第四步checkmodel的时候，winbugs坐下角显示的不是“modelissyntacticallycorrect”，而是在alpha1~dnorm（0,1e-05）一句1e处显示e处应为”expectedrightparenthesis，在代码最后0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0））指示”invalidorunexpectedtokenscanned”，
这要怎么弄呢？
第一次接触这个软件，还不会使，请高手指点。
WinBUGS在统计分析中的应用（第二部分）
By齐韬@2008/12/18
关键词：
BayesianAnalysis,SAS,WinBUGS分类：
数据分析,统计软件,贝叶斯统计
作者信息：
ComputationalMathematicianinAnnproAnalyticTechnologies,Inc.
版权声明：
本文版权归原作者所有，未经许可不得转载。
原文可能随时需要修改纰漏，全文复制转载会带来不必要的误导，若您想推荐给朋友阅读，敬请以负责的态度提供原文链接；点此查看如何在学术刊物中引用本文
常规引用方式
齐韬.WinBUGS在统计分析中的应用（第二部分）.统计之都,2008.12.URL:
http:
//cos.name/2008/12/statistical-analysis-and-winbugs-part-2/.
BibTeX引用
@ARTICLE{,
AUTHOR={齐韬},
TITLE={WinBUGS在统计分析中的应用（第二部分）},
JOURNAL={统计之都},
YEAR={2008},
month={12},
URL={http:
//cos.name/2008/12/statistical-analysis-and-winbugs-part-2/},
}
第一节WinBUGS数据分析案例
在这一节中，我将拿一个经典的研究数据，利用WinBUGS给出简单的分析。
首先介绍一下这个数据：
Seeds
seedO.aegyptiaco75seedO.aegyptiaco73BeanCucumberBeanCucumber
rnr/nrnr/nrnr/nrnr/n
10390.26560.838160.53120.25
23620.3753740.7210300.3322410.54
23810.2855720.768280.2915300.5
26510.5132510.6323450.5132510.63
17390.4446790.58040370.43
10130.77
这个数据来自Crowder（1978）。
之后BreslowandClayton（1993）作为例子，也分析过这个数据。
数据反映的是某一品种的豆类种子和某一品种的黄瓜种子，分别放在21个培养皿（plates）中分别培养，在根提取物aegyptiaco75和aegyptiaco73的作用下的出芽率的差异。
其中r是出芽的个数，n是种子的个数，而r/n是出芽率。
我们用randomeffectlogisticregression模型来进行分析（注意，在Bayesian分析中，通常是将covariates看做是服从某一个分布的随机变量，这和一般意义上的GLM,GLME,LME中对于covariates解释是不同的）：
其中
是种子的类型，
是根提取物的类型，
是交互项,
是给定的独立的“noninformative”先验参数。
在Bayesian分析中，通常我们会定义一个DAG图（即DirectedAcyclicGraph有向无圈图）。
我们可以在WinBUGS中通过设计DAG图来定义模型。
不过这一节中我们还是用WinBUGS中的BUGS语言来定义模型，如何在WinBUGS中通过设计DAG图来定义模型我将在下一节中详细介绍，但是必须要说明的是BUGS语言比DAG图灵活，不过直观性不如后者。
模型
model
{
for（iin1:
N）{
r[i]~dbin（p[i],n[i]）
b[i]~dnorm（0.0,tau）
logit（p[i]）<-alpha0+alpha1*x1[i]+alpha2*x2[i]+
alpha12*x1[i]*x2[i]+b[i]
}
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）
sigma<-1/sqrt（tau）
}
WinBUGSdoodle模型
数据
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
升级会员 