GeoBUGS疾病制圖法在條件自回歸模型中的應(yīng)用*

蘭州大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)研究所(730000)

高文龍# 張繼巍# 拉扎提·木拉提 李學(xué)朝 秦天燕 李娟生△

教育部人文社科項(xiàng)目(項(xiàng)目號(hào)：15XJC910001)；中央高校基本科研業(yè)務(wù)專項(xiàng)資金(項(xiàng)目號(hào)：LZUjbky-2016-025)

#共同第一作者

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

GeoBUGS疾病制圖法在條件自回歸模型中的應(yīng)用*

蘭州大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)研究所(730000)

高文龍#張繼巍#拉扎提·木拉提 李學(xué)朝 秦天燕 李娟生△

貝葉斯統(tǒng)計(jì)起源于英國(guó)學(xué)者貝葉斯在1763年的一篇題為“機(jī)遇理論中一個(gè)問(wèn)題的解”的論文,他提出了著名的貝葉斯公式[1]。貝葉斯統(tǒng)計(jì)方法與經(jīng)典統(tǒng)計(jì)方法最根本的區(qū)別在于不僅利用總體信息和樣本信息進(jìn)行統(tǒng)計(jì)推斷,而且充分利用了參數(shù)的先驗(yàn)信息,它將每一個(gè)不確定的參數(shù)都看成一個(gè)隨機(jī)變量,通過(guò)給予先驗(yàn)分布,結(jié)合馬爾科夫鏈蒙特卡洛(markov chain monte carlo,MCMC)法進(jìn)行Gibbs抽樣,得出參數(shù)的后驗(yàn)分布,因此可以提高統(tǒng)計(jì)推斷的效果。其廣泛應(yīng)用于經(jīng)濟(jì)、金融、醫(yī)學(xué)、生物統(tǒng)計(jì)、自然科學(xué)和社會(huì)科學(xué)等各個(gè)領(lǐng)域[2-4]。隨著OpenBUGS軟件[5]的成功開發(fā),對(duì)GeoBUGS模塊的功能和界面做了相應(yīng)的調(diào)整和優(yōu)化,與其在WinBUGS中相比,增加了新的貝葉斯地圖模板和應(yīng)用案例,也為ArcView格式文件的導(dǎo)入提供了接口,促進(jìn)了其在疾病空間模型的構(gòu)造和空間地圖的繪制方面的發(fā)展[6-7]。

疾病制圖是空間流行病學(xué)研究的主要任務(wù),其目的在于將疾病危險(xiǎn)的空間變異可視化在地圖上,確定病例聚集地點(diǎn)和空間分布輪廓,揭示疾病空間分布聯(lián)系模式,為進(jìn)一步疾病病因和危險(xiǎn)因素的研究提供線索[8]。而貝葉斯統(tǒng)計(jì)分析軟件OpenBUGS中的模塊GeoBUGS就是用于空間數(shù)據(jù)的分析和空間地圖的繪制。但目前國(guó)內(nèi)有關(guān)GeoBUGS在空間地圖中的應(yīng)用報(bào)道的比較少,并且沒(méi)有一個(gè)詳細(xì)的介紹,使廣大讀者不能真正去掌握和應(yīng)用。本文主要通過(guò)實(shí)例詳細(xì)介紹如何實(shí)現(xiàn)GeoBUGS在疾病制圖中的應(yīng)用,希望能夠?qū)V大讀者起到拋磚引玉的作用。

GeoBUGS軟件介紹

1.GeoBUGS的功能:GeoBUGS是貝葉斯統(tǒng)計(jì)分析軟件OpenBUGS的一個(gè)附加模塊[9],專門用來(lái)分析空間數(shù)據(jù)并生成空間地圖。GeoBUGS通過(guò)Map下拉菜單為我們提供了一個(gè)窗口界面,我們可以通過(guò)鼠標(biāo)指令完成疾病地圖的繪制。

2.GeoBUGS包含的現(xiàn)成地圖有: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:用于導(dǎo)入用戶在ArcInfo軟件中自定義的多邊形地圖。

(4)Import Epimap:用于導(dǎo)入用戶在Epimap軟件中自定義的多邊形地圖。

(5)Import Splus:用于導(dǎo)入用戶在Splus中自定義的多邊形地圖。

(6)Export Splus:用于將GeoBUGS生成的空間地圖導(dǎo)出為Splus格式。

4.GeoBUGS操作步驟:

(1) 生成GeoBUGS格式的地圖:目前分析用的地圖普遍采用.shp格式,而GeoBUGS只能識(shí)別自帶的地圖,對(duì)于ArcInfo、EpiMap和Splus的地圖格式,需要進(jìn)行適當(dāng)?shù)霓D(zhuǎn)化才能識(shí)別。因此我們首先需生成一個(gè)GeoBUGS格式的目標(biāo)地圖,可以通過(guò)ArcInfo或R軟件來(lái)實(shí)現(xiàn)。

(2)生成空間鄰接矩陣:在OpenBUGS窗口界面,點(diǎn)擊“Map->Adjacency Tool”啟動(dòng) Adjacency Tool對(duì)話框,先在map標(biāo)簽中選擇自己將要分析的地圖,點(diǎn)擊adj map按鈕生成此地圖,此時(shí)adj matrix 按鈕和 show region按鈕被激活,然后點(diǎn)擊adj matrix 按鈕生成鄰接矩陣(注意:GeoBUGS在計(jì)算鄰接矩陣時(shí),存在0.1米的誤差[9];生成的矩陣在下一步的數(shù)據(jù)加載中將會(huì)用到);通過(guò)show region 按鈕和標(biāo)簽中的數(shù)字凸顯地圖中指定的區(qū)域。

(3)進(jìn)行Gibbs抽樣:具體過(guò)程參見OpenBUGS用戶指南[10]。

(4)生成疾病地圖:打開Map Tool對(duì)話框,如圖1所示:

①在Map標(biāo)簽的下拉菜單中選擇我們想要繪制的地圖;

②在Variable標(biāo)簽的空白框中輸入希望在地圖中呈現(xiàn)的模型參數(shù);

③根據(jù)Variable中變量的類型,在Quantity下拉菜單中選擇相應(yīng)的值。如果變量是data(例如SMR、期望值E或者協(xié)變量),則選擇Value;如果變量是stochastic quantity(如相對(duì)危險(xiǎn)度),我們可以選擇Quanntity菜單中的任意一個(gè);

a)當(dāng)為變量設(shè)置了summary monitor時(shí),我們只能選擇mean(summary);

圖1 Map Tool對(duì)話框

b)當(dāng)為變量設(shè)置了samples monitor,我們可以選擇Quantity中其他類型;(當(dāng)選擇Percentile時(shí),quantile標(biāo)簽會(huì)被激活,輸入適當(dāng)?shù)陌俜治粩?shù),就會(huì)繪制出變量后驗(yàn)分位數(shù)地圖;當(dāng)選擇prob greater或者prob less時(shí),threshold標(biāo)簽就會(huì)被激活,輸入相應(yīng)的值,將會(huì)繪制變量值大于等于或者小于等于指定值的后驗(yàn)概率地圖)

④設(shè)置分割點(diǎn):在GeoBUGS中有兩種地圖分割點(diǎn),分別為絕對(duì)值分割點(diǎn)(abs value)和百分位數(shù)分割點(diǎn)(percentile),對(duì)于絕對(duì)值分割點(diǎn),GeoBUGS選擇了一組基于變量絕對(duì)值的默認(rèn)間隔來(lái)繪圖(這些間隔一般為等距間隔);對(duì)于百分位數(shù)分割點(diǎn),GeoBUGS選擇變量先驗(yàn)分布的第10、第50和第90百分位數(shù)去繪制地圖。GeoBUGS也允許用戶自定義切割點(diǎn)的值。

⑤設(shè)置顏色:我們可以通過(guò)Palette下拉菜單來(lái)編輯地圖的顏色,也可以單擊地圖上注釋的小方塊來(lái)改變顏色。

⑥我們可以通過(guò) set cuts按鈕去更新目前選擇的地圖,也可以通過(guò)plot按鈕去生成一個(gè)新的地圖。

⑦地圖的輸出和保存:通過(guò)GeoBUGS工具生成的疾病地圖,通過(guò)“File<-Save As”將其保存為OpenBUGS能直接調(diào)用的.odc格式,以便我們重新編輯地圖的切割點(diǎn)和地圖顏色;另外,我們可以通過(guò)同時(shí)按住“Ctrl”和“Space”鍵選中地圖,將其復(fù)制、粘貼到Word或PowerPoint中。

GeoBUGS應(yīng)用實(shí)例

1.資料來(lái)源:資料來(lái)源于國(guó)家衛(wèi)計(jì)委網(wǎng)站《2009年中國(guó)衛(wèi)生統(tǒng)計(jì)年鑒》中女性乳腺癌的患病數(shù)據(jù)[11],利用乳腺癌的患病率(pi)和實(shí)查人數(shù)(ni)計(jì)算實(shí)際發(fā)病數(shù)(Yi=pi*ni),根據(jù)年齡分布計(jì)算其期望發(fā)病數(shù)Ei=∑(Nij*Pij),其中Nij和Pij分別是各年齡組的總?cè)藬?shù)和乳腺癌發(fā)病率;.shp格式的中國(guó)地圖來(lái)源于國(guó)家地理信息系統(tǒng)。目的是基于貝葉斯統(tǒng)計(jì)方法實(shí)現(xiàn)GeoBUGS在疾病制圖中的應(yīng)用。

2.研究方法與結(jié)果:首先利用R軟件將.shp格式的地圖轉(zhuǎn)換成GeoBUGS格式;然后依據(jù)貝葉斯統(tǒng)計(jì)推斷的基本原理,結(jié)合OpenBUGS軟件,進(jìn)行MCMC模擬,構(gòu)建貝葉斯條件自回歸模型(conditional autoregressive,CAR)如下:Log(mu[i])=Log(E[i])+alpha0+b[i]

式中alpha0反映的是各個(gè)區(qū)域間患病的相對(duì)基準(zhǔn)風(fēng)險(xiǎn);b[i]反應(yīng)的是與地域相關(guān)的潛在的患病風(fēng)險(xiǎn)因子。上述貝葉斯CAR模型的OpenBUGS代碼如下:

model { for (i in 1:N)

{O[i] ～ dpois(mu[i]) #觀察病例數(shù)服從泊松分布

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

RR[i] <- exp(alpha0 + b[i])} #地圖中第i個(gè)區(qū)域的相對(duì)危險(xiǎn)度

b[1:N] ～ car.normal(adj[],weights[],num[],tau) #設(shè)定b[i]通過(guò)car.normal先驗(yàn)分布來(lái)描述

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

alpha0 ～ dnorm(0,0.0001) #設(shè)定截距的無(wú)信息先驗(yàn)服從正態(tài)分布

tau ～ dgamma(0.5,0.0005) #設(shè)定參數(shù)精度的無(wú)信息先驗(yàn)服從伽馬分布

sigma <- sqrt(1 / tau) #通過(guò)方差求解標(biāo)準(zhǔn)差

b.mean <- sum(b[])} #通過(guò)sum求解隨機(jī)效應(yīng)b[i]的均數(shù)

其中:adj[]表示每個(gè)區(qū)域鄰接區(qū)域的編號(hào);weights[]表示各個(gè)區(qū)域間的權(quán)重因子;num[]表示每個(gè)區(qū)域相鄰區(qū)域的個(gè)數(shù);tau表示條件自回歸模型先驗(yàn)參數(shù)的精度;sumNumNeigh表示每個(gè)區(qū)域相鄰區(qū)域個(gè)數(shù)的合計(jì)。

本例采用兩條鏈,拋去前1000次迭代,以提高迭代的穩(wěn)定性和模型的收斂性。待模型收斂后得到參數(shù)后驗(yàn)分布的均數(shù)、標(biāo)準(zhǔn)差、中位數(shù)等信息如圖2和圖3。

最后利用GeoBUGS繪制的疾病地圖如圖4:

從圖中我們能直觀的看出,2008年女性乳腺癌發(fā)病RR≥6的有吉林、江蘇、浙江和貴州4省;4≤RR<5的只有廣西省;3≤RR<4的有山東、山西、湖北、青海、四川、云南和安徽7個(gè)省;而1≤RR<2的省(市)、自治區(qū)最多,總計(jì)11個(gè)。除此之外,我們能從圖中清楚的看出2008年全國(guó)各省(市)、自治區(qū)女性乳腺癌的空間分布輪廓和空間聚集現(xiàn)象。對(duì)那些發(fā)病相對(duì)危險(xiǎn)度高的地區(qū)我們應(yīng)該加強(qiáng)檢查力度,投入更多的衛(wèi)生資源,做到早發(fā)現(xiàn)、早診斷、早治療,降低疾病的發(fā)病風(fēng)險(xiǎn),減少疾病的發(fā)生。因此,依據(jù)地圖提供的信息,不僅有利于我們合理優(yōu)化配置衛(wèi)生資源,也有利于我們發(fā)現(xiàn)疾病的聚集傾向和分布輪廓,為進(jìn)一步探索疾病病因和危險(xiǎn)因素提供線索。

圖2 Gibbs抽樣結(jié)果

圖3 模型迭代結(jié)果

圖4 GeoBUGS輸出結(jié)果

除此之外,GeoBUGS也可以繪制出疾病期望發(fā)病數(shù)(Ei)、實(shí)際發(fā)病數(shù)(Oi)、標(biāo)準(zhǔn)化死亡比(SMR)和空間模型協(xié)變量等參數(shù)的地圖。當(dāng)然,GeoBUGS不僅適用于CAR模型,也實(shí)現(xiàn)了在Multivariate CAR模型、Gaussian kriging模型(Diggle等人)、Poisson-gamma convolution 模型(Best等人)和Shared comp-onent模型(Knorr-Held等人)等空間貝葉斯模型中的應(yīng)

用[9-10]。因此,我們可以根據(jù)自己的需要,繪制不同的疾病地圖。

討 論

本文著重介紹了貝葉斯空間分析工具GeoBUGS的功能及具體操作步驟,并通過(guò)實(shí)例分析,實(shí)現(xiàn)了GeoBUGS在疾病制圖中的應(yīng)用,使讀者對(duì)該工具有一個(gè)初步的認(rèn)識(shí),并能掌握其具體操作步驟和基本要求。相比于其他疾病制圖軟件(ArcGIS[12]、EpiMap[13]、ArcInfo),GeoBUGS不僅容易獲得,更重要的是其基于貝葉斯統(tǒng)計(jì)方法充分利用了先驗(yàn)信息[1],避免了頻率統(tǒng)計(jì)學(xué)方法中要求樣本之間相互獨(dú)立和無(wú)法利用先驗(yàn)信息的缺點(diǎn),使繪制的疾病地圖更加精確。

當(dāng)然,GeoBUGS工具也有自己的限制性缺點(diǎn),因?yàn)槠渲荒茏R(shí)別自帶的地圖格式,因而對(duì)很多我們需要分析的地圖(中國(guó)各省市地圖),必須通過(guò)其他軟件(ArcInfo、R、EpiMap)轉(zhuǎn)化成GeoBUGS格式才能使用。但隨著貝葉斯統(tǒng)計(jì)方法在空間流行病學(xué)中的廣泛應(yīng)用和OpenBUGS軟件的不斷更新[5],GeoBUGS的功能將會(huì)不斷地增強(qiáng),并將成為分析空間數(shù)據(jù)和繪制疾病地圖的主流工具,掌握其繪制疾病地圖的方法也會(huì)成為一門必備技術(shù),這就要求我們?cè)谡莆誈eoBUGS基本操作和貝葉斯統(tǒng)計(jì)基本理論的基礎(chǔ)上反復(fù)摸索。

[1] 韋來(lái)生編著.貝葉斯統(tǒng)計(jì).北京:高等教育出版社,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線性混合效應(yīng)模型多中心臨床試驗(yàn)研究.中國(guó)衛(wèi)生統(tǒng)計(jì),2005,22(4):200-203.

[5] 張繼巍,高文龍,李娟生等.OpenBUGS軟件介紹及應(yīng)用.中國(guó)衛(wèi)生統(tǒng)計(jì),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] 徐麗,方亞.空間流行病學(xué)中的疾病制圖常用方法.中國(guó)衛(wèi)生統(tǒng)計(jì),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制圖表達(dá)在地圖制圖方面的應(yīng)用.2014,2:11-13.

[13] 王勁松.流行病學(xué)地圖軟件EpiMap.醫(yī)學(xué)信息,2000,13(7):374-375.

(責(zé)任編輯：劉 壯)