GeoBUGS疾病制图法在条件自回归模型中的应用*

兰州大学公共卫生学院流行病与卫生统计学研究所(730000)

高文龙# 张继巍# 拉扎提·木拉提 李学朝 秦天燕 李娟生△

教育部人文社科项目(项目号：15XJC910001)；中央高校基本科研业务专项资金(项目号：LZUjbky-2016-025)

#共同第一作者

△通信作者：李娟生，E-mail：lijsh@lzu.edu.cn

GeoBUGS疾病制图法在条件自回归模型中的应用*

兰州大学公共卫生学院流行病与卫生统计学研究所(730000)

高文龙#张继巍#拉扎提·木拉提 李学朝 秦天燕 李娟生△

贝叶斯统计起源于英国学者贝叶斯在1763年的一篇题为“机遇理论中一个问题的解”的论文,他提出了著名的贝叶斯公式[1]。贝叶斯统计方法与经典统计方法最根本的区别在于不仅利用总体信息和样本信息进行统计推断,而且充分利用了参数的先验信息,它将每一个不确定的参数都看成一个随机变量,通过给予先验分布,结合马尔科夫链蒙特卡洛(markov chain monte carlo,MCMC)法进行Gibbs抽样,得出参数的后验分布,因此可以提高统计推断的效果。其广泛应用于经济、金融、医学、生物统计、自然科学和社会科学等各个领域[2-4]。随着OpenBUGS软件[5]的成功开发,对GeoBUGS模块的功能和界面做了相应的调整和优化,与其在WinBUGS中相比,增加了新的贝叶斯地图模板和应用案例,也为ArcView格式文件的导入提供了接口,促进了其在疾病空间模型的构造和空间地图的绘制方面的发展[6-7]。

疾病制图是空间流行病学研究的主要任务,其目的在于将疾病危险的空间变异可视化在地图上,确定病例聚集地点和空间分布轮廓,揭示疾病空间分布联系模式,为进一步疾病病因和危险因素的研究提供线索[8]。而贝叶斯统计分析软件OpenBUGS中的模块GeoBUGS就是用于空间数据的分析和空间地图的绘制。但目前国内有关GeoBUGS在空间地图中的应用报道的比较少,并且没有一个详细的介绍,使广大读者不能真正去掌握和应用。本文主要通过实例详细介绍如何实现GeoBUGS在疾病制图中的应用,希望能够对广大读者起到抛砖引玉的作用。

GeoBUGS软件介绍

1.GeoBUGS的功能:GeoBUGS是贝叶斯统计分析软件OpenBUGS的一个附加模块[9],专门用来分析空间数据并生成空间地图。GeoBUGS通过Map下拉菜单为我们提供了一个窗口界面,我们可以通过鼠标指令完成疾病地图的绘制。

2.GeoBUGS包含的现成地图有:Belgium、Elevation、Foreset、France、GB_Counties、GreeceNomoi、grid、gridepimap、gridsplus、gridarcinfo、HalfRongelap、Huddersfield_750m_grid、LHA、Munich、Rongelap、Sardinia、Scotland、WestYorkshire。

3.基本操作:GeoBUGS基本操作包括Mapping Tool、Adjacency Tool、Import ArcInfo、Import Epimap、Import Splus和Export Splus命令。

(1)Mapping Tool:用于生成空间地图。

(2)Adjacency Tool:用于生成空间邻接矩阵。

(3)Import ArcInfo:用于导入用户在ArcInfo软件中自定义的多边形地图。

(4)Import Epimap:用于导入用户在Epimap软件中自定义的多边形地图。

(5)Import Splus:用于导入用户在Splus中自定义的多边形地图。

(6)Export Splus:用于将GeoBUGS生成的空间地图导出为Splus格式。

4.GeoBUGS操作步骤:

(1) 生成GeoBUGS格式的地图:目前分析用的地图普遍采用.shp格式,而GeoBUGS只能识别自带的地图,对于ArcInfo、EpiMap和Splus的地图格式,需要进行适当的转化才能识别。因此我们首先需生成一个GeoBUGS格式的目标地图,可以通过ArcInfo或R软件来实现。

(2)生成空间邻接矩阵:在OpenBUGS窗口界面,点击“Map->Adjacency Tool”启动 Adjacency Tool对话框,先在map标签中选择自己将要分析的地图,点击adj map按钮生成此地图,此时adj matrix 按钮和 show region按钮被激活,然后点击adj matrix 按钮生成邻接矩阵(注意:GeoBUGS在计算邻接矩阵时,存在0.1米的误差[9];生成的矩阵在下一步的数据加载中将会用到);通过show region 按钮和标签中的数字凸显地图中指定的区域。

(3)进行Gibbs抽样:具体过程参见OpenBUGS用户指南[10]。

(4)生成疾病地图:打开Map Tool对话框,如图1所示:

①在Map标签的下拉菜单中选择我们想要绘制的地图;

②在Variable标签的空白框中输入希望在地图中呈现的模型参数;

③根据Variable中变量的类型,在Quantity下拉菜单中选择相应的值。如果变量是data(例如SMR、期望值E或者协变量),则选择Value;如果变量是stochastic quantity(如相对危险度),我们可以选择Quanntity菜单中的任意一个;

a)当为变量设置了summary monitor时,我们只能选择mean(summary);

图1 Map Tool对话框

b)当为变量设置了samples monitor,我们可以选择Quantity中其他类型;(当选择Percentile时,quantile标签会被激活,输入适当的百分位数,就会绘制出变量后验分位数地图;当选择prob greater或者prob less时,threshold标签就会被激活,输入相应的值,将会绘制变量值大于等于或者小于等于指定值的后验概率地图)

④设置分割点:在GeoBUGS中有两种地图分割点,分别为绝对值分割点(abs value)和百分位数分割点(percentile),对于绝对值分割点,GeoBUGS选择了一组基于变量绝对值的默认间隔来绘图(这些间隔一般为等距间隔);对于百分位数分割点,GeoBUGS选择变量先验分布的第10、第50和第90百分位数去绘制地图。GeoBUGS也允许用户自定义切割点的值。

⑤设置颜色:我们可以通过Palette下拉菜单来编辑地图的颜色,也可以单击地图上注释的小方块来改变颜色。

⑥我们可以通过 set cuts按钮去更新目前选择的地图,也可以通过plot按钮去生成一个新的地图。

⑦地图的输出和保存:通过GeoBUGS工具生成的疾病地图,通过“File<-Save As”将其保存为OpenBUGS能直接调用的.odc格式,以便我们重新编辑地图的切割点和地图颜色;另外,我们可以通过同时按住“Ctrl”和“Space”键选中地图,将其复制、粘贴到Word或PowerPoint中。

GeoBUGS应用实例

1.资料来源:资料来源于国家卫计委网站《2009年中国卫生统计年鉴》中女性乳腺癌的患病数据[11],利用乳腺癌的患病率(pi)和实查人数(ni)计算实际发病数(Yi=pi*ni),根据年龄分布计算其期望发病数Ei=∑(Nij*Pij),其中Nij和Pij分别是各年龄组的总人数和乳腺癌发病率;.shp格式的中国地图来源于国家地理信息系统。目的是基于贝叶斯统计方法实现GeoBUGS在疾病制图中的应用。

2.研究方法与结果:首先利用R软件将.shp格式的地图转换成GeoBUGS格式;然后依据贝叶斯统计推断的基本原理,结合OpenBUGS软件,进行MCMC模拟,构建贝叶斯条件自回归模型(conditional autoregressive,CAR)如下:Log(mu[i])=Log(E[i])+alpha0+b[i]

式中alpha0反映的是各个区域间患病的相对基准风险;b[i]反应的是与地域相关的潜在的患病风险因子。上述贝叶斯CAR模型的OpenBUGS代码如下:

model { for (i in 1:N)

{O[i] ～ dpois(mu[i]) #观察病例数服从泊松分布

log(mu[i]) <- log(E[i]) + alpha0 + b[i] #条件自回归模型

RR[i] <- exp(alpha0 + b[i])} #地图中第i个区域的相对危险度

b[1:N] ～ car.normal(adj[],weights[],num[],tau) #设定b[i]通过car.normal先验分布来描述

for(k in 1:sumNumNeigh) {weights[k] <- 1}

alpha0 ～ dnorm(0,0.0001) #设定截距的无信息先验服从正态分布

tau ～ dgamma(0.5,0.0005) #设定参数精度的无信息先验服从伽马分布

sigma <- sqrt(1 / tau) #通过方差求解标准差

b.mean <- sum(b[])} #通过sum求解随机效应b[i]的均数

其中:adj[]表示每个区域邻接区域的编号;weights[]表示各个区域间的权重因子;num[]表示每个区域相邻区域的个数;tau表示条件自回归模型先验参数的精度;sumNumNeigh表示每个区域相邻区域个数的合计。

本例采用两条链,抛去前1000次迭代,以提高迭代的稳定性和模型的收敛性。待模型收敛后得到参数后验分布的均数、标准差、中位数等信息如图2和图3。

最后利用GeoBUGS绘制的疾病地图如图4:

从图中我们能直观的看出,2008年女性乳腺癌发病RR≥6的有吉林、江苏、浙江和贵州4省;4≤RR<5的只有广西省;3≤RR<4的有山东、山西、湖北、青海、四川、云南和安徽7个省;而1≤RR<2的省(市)、自治区最多,总计11个。除此之外,我们能从图中清楚的看出2008年全国各省(市)、自治区女性乳腺癌的空间分布轮廓和空间聚集现象。对那些发病相对危险度高的地区我们应该加强检查力度,投入更多的卫生资源,做到早发现、早诊断、早治疗,降低疾病的发病风险,减少疾病的发生。因此,依据地图提供的信息,不仅有利于我们合理优化配置卫生资源,也有利于我们发现疾病的聚集倾向和分布轮廓,为进一步探索疾病病因和危险因素提供线索。

图2 Gibbs抽样结果

图3 模型迭代结果

图4 GeoBUGS输出结果

除此之外,GeoBUGS也可以绘制出疾病期望发病数(Ei)、实际发病数(Oi)、标准化死亡比(SMR)和空间模型协变量等参数的地图。当然,GeoBUGS不仅适用于CAR模型,也实现了在Multivariate CAR模型、Gaussian kriging模型(Diggle等人)、Poisson-gamma convolution 模型(Best等人)和Shared comp-onent模型(Knorr-Held等人)等空间贝叶斯模型中的应

用[9-10]。因此,我们可以根据自己的需要,绘制不同的疾病地图。

讨 论

本文着重介绍了贝叶斯空间分析工具GeoBUGS的功能及具体操作步骤,并通过实例分析,实现了GeoBUGS在疾病制图中的应用,使读者对该工具有一个初步的认识,并能掌握其具体操作步骤和基本要求。相比于其他疾病制图软件(ArcGIS[12]、EpiMap[13]、ArcInfo),GeoBUGS不仅容易获得,更重要的是其基于贝叶斯统计方法充分利用了先验信息[1],避免了频率统计学方法中要求样本之间相互独立和无法利用先验信息的缺点,使绘制的疾病地图更加精确。

当然,GeoBUGS工具也有自己的限制性缺点,因为其只能识别自带的地图格式,因而对很多我们需要分析的地图(中国各省市地图),必须通过其他软件(ArcInfo、R、EpiMap)转化成GeoBUGS格式才能使用。但随着贝叶斯统计方法在空间流行病学中的广泛应用和OpenBUGS软件的不断更新[5],GeoBUGS的功能将会不断地增强,并将成为分析空间数据和绘制疾病地图的主流工具,掌握其绘制疾病地图的方法也会成为一门必备技术,这就要求我们在掌握GeoBUGS基本操作和贝叶斯统计基本理论的基础上反复摸索。

[1] 韦来生编著.贝叶斯统计.北京:高等教育出版社,2016.3.

[2] Carroll R,Lawson AB,Faes C,et al.Comparing INLA and OpenBUGS for hierarchical Poisson modeling in disease mapping.Spatial and Spatio-temporal Epidemiology,2015,14(15):5-54.

[3] Lyle W,Konigsberg,Frankenberg.Bayes in Biological Anthropology.American Journal of Physical Anthropology,2013,152(57):153-184.

[4] 刘桂芬,孟海英,张岩波.Bayes线性混合效应模型多中心临床试验研究.中国卫生统计,2005,22(4):200-203.

[5] 张继巍,高文龙,李娟生等.OpenBUGS软件介绍及应用.中国卫生统计,2017,34(1):170-172.

[6] Avril H,Daniel B.Bayesian disease mapping using product partition models.Statistics in Medicine,2008,27:3868-3893.

[7] Goicoa T,Ugarte MD,Etxeberria J,et al.Age-space-time CAR models in Bayesian disease mapping.Statistics in Medicine,2016,35:2391-2405.

[8] 徐丽,方亚.空间流行病学中的疾病制图常用方法.中国卫生统计,2015,32(2):338-341.

[9] GeoBUGS 3.2.3 user manual.

[10] OpenBUGS3.2.3 user manual.

[11] http://www.nhfpc.gov.cn/zwgkzt/tjnj/list.shtml.

[12] 谢世琴,柴微涛,等.ArcGIS制图表达在地图制图方面的应用.2014,2:11-13.

[13] 王劲松.流行病学地图软件EpiMap.医学信息,2000,13(7):374-375.

(责任编辑：刘 壮)