張輝國,梁韻婷,胡錫健

布魯氏菌病(簡稱布病)是一種常見的人獸共患傳染病,既影響身體健康,又影響畜牧業(yè)發(fā)展與社會穩(wěn)定。目前,人間布魯氏菌病在我國甲乙類法定報告?zhèn)魅静≈邪l(fā)病率及發(fā)病人數(shù)已升至第8位,是一個不可忽視的公共衛(wèi)生問題。

在我國,布病多高發(fā)于東北、西北及部分華北地區(qū),并呈現(xiàn)出明顯的空間聚集性。布病的空間分布格局離不開氣候、自然環(huán)境及畜牧業(yè)等因素。一方面,這些高發(fā)地區(qū)由于適宜的氣候、自然環(huán)境等條件形成了發(fā)達的畜牧業(yè)并進一步影響了布病發(fā)病率;另一方面,氣候、自然環(huán)境等條件又能直接影響布氏菌的存活。目前,學術(shù)界圍繞布病空間流行病學分析已經(jīng)開展了廣泛的研究,從研究方法看主要涉及全局與局部空間自相關(guān)分析[1-3]、時空加權(quán)泊松回歸模型[4]、地理加權(quán)回歸模型[5]、貝葉斯時空模型[6],而考慮空間相關(guān)性并建立空間回歸模型的定量研究較為缺乏。因此,本文提出了貝葉斯LASSO矩陣指數(shù)空間模型(BL-MESS),用于探索2020年中國大陸31個省、直轄市、自治區(qū)布病發(fā)病率與氣候因素、自然環(huán)境及畜牧業(yè)因素的關(guān)系,目的在于從關(guān)系復雜的眾多變量中篩選出關(guān)鍵影響因素,為精準疫情防控策略的制定和實施提供依據(jù)。

1 資料與方法

1.1 資料來源與預處理 2020年全國各地區(qū)布病發(fā)病率(1/10萬)來自2021年《中國衛(wèi)生健康統(tǒng)計年鑒》,各地區(qū)豬、牛、山羊與綿羊年末存欄量(萬頭)來自2021年《中國農(nóng)村統(tǒng)計年鑒》,牧草地面積(萬公頃)來自各地區(qū)第三次國土調(diào)查主要數(shù)據(jù)公報,省面積(萬公頃)來自中華人民共和國行政區(qū)劃統(tǒng)計表,平均氣溫(℃)、平均風速(m/s)、年降水量(mm)、平均相對濕度(%)、年日照時數(shù)(h)、平均氣壓(hpa)來自于2020年中國氣象要素年度空間插值數(shù)據(jù)集(https://www.resdc.cn/DOI/DOI.aspx?DOIID=96),平均海拔高度來自NASA SRTM3 V4.1 90 m精度和ASTER-GDEM V2 30 m精度DEM數(shù)據(jù)。以省、自治區(qū)、直轄市為地理單元,在Arcgis軟件中分別對平均氣溫等氣候因素及平均海拔高度進行匯總并計算出各地理單元的平均值。對平均海拔高度做平方處理后對除布病發(fā)病率以外的其他變量進行標準化處理以消除各量綱影響,對布病發(fā)病率除以100再進行反正弦平方根變換以滿足誤差分布的正態(tài)性與方差齊性。

1.2 空間自相關(guān)分析 全局Moran’sI指數(shù)是衡量事物空間要素屬性值在全區(qū)域內(nèi)的聚合或離散程度的指標,值域位于[-1,1]之間,當取值大于0時,意味著較大的觀測值通常被較大的值所包圍,反之亦然,表明存在空間正自相關(guān)性;當取值小于0時,意味著較小的觀測值通常被較大的值所包圍,反之亦然,表明存在空間負自相關(guān)性;當取值接近于0時,意味著沒有明顯的空間分布規(guī)律,表明不存在空間自相關(guān)性。絕對值越大,意味著空間自相關(guān)程度越高。利用Geoda軟件計算全局Moran’sI指數(shù)并設顯著性水平為0.05,通過蒙特卡羅試驗對Moran’sI指數(shù)進行顯著性檢驗,置換次數(shù)為999,當Z>1.96時具有統(tǒng)計學意義。

1.3 相關(guān)性分析 Pearson相關(guān)系數(shù)是衡量兩變量間線性相關(guān)程度的指標,值域位于[-1,1],當取值為正時表明兩變量呈正相關(guān),取值為負時表明兩變量呈負相關(guān),而取值為0時表明兩變量間無線性關(guān)系,同時絕對值越大,相關(guān)程度越高。此外,條件數(shù)是衡量多重共線性嚴重程度的一個重要指標,當取值小于100時表明共線性程度較小,取值大于1 000時表明存在嚴重的共線性,取值為100-1 000時則存在中等程度的共線性。在R軟件中應用Pearson相關(guān)性分析了解風險因素間的相關(guān)性并可視化,同時計算自變量相關(guān)矩陣的條件數(shù)。

1.4 研究方法 考慮到2020年全國布病數(shù)據(jù)的空間相關(guān)性及模型預測因子的稀疏性假設,本研究構(gòu)建貝葉斯LASSO矩陣指數(shù)空間模型(BL-MESS),同時實現(xiàn)系數(shù)估計與變量選擇,從一系列潛在發(fā)病風險因子中篩選出關(guān)鍵影響因素。

空間自回歸模型[7]為:

Y=αWY+Xβ+ε

其中α為空間自回歸參數(shù),取值范圍為(-1,1),正值表明存在空間正相關(guān)性,負值表明存在空間負相關(guān)性,取零表明無空間相關(guān)性,將該式移項并用exp(ρW)代替In-αW即為MESS模型,進一步可得兩空間參數(shù)存在對應關(guān)系α=1-eρ。

根據(jù)LeSage和Pace的設定[8],MESS全模型的貝葉斯分層表示如下

exp(ρW)Y=Xβ+ε,ε～Nn(0n,σ2In)

π(β)～Np(c,T)

π(ρ)～N(0,ζ)

為了實現(xiàn)變量收縮,在MESS模型的基礎上結(jié)合了貝葉斯LASSO[9],用拉普拉斯密度的正態(tài)尺度混合表示代替貝葉斯MESS中系數(shù)的多元正態(tài)分布。此外,將隨機搜索變量選擇先驗作為空間系數(shù)ρ的先驗,該處理的目的是在空間相關(guān)性不顯著時實現(xiàn)空間系數(shù)的變量選擇,使其退化為經(jīng)典線性回歸模型,同時該處理能使空間系數(shù)限制在一個更合理的范圍內(nèi)。綜上,BL-MESS全模型的分層表示如下:

exp(ρW)Y=Xβ+ε,ε～Nn(0n,σ2In)

π(ρ)～N(0,ζ)

γ～Ber(0.5)

根據(jù)以上模型設定,分別推導兩類模型各參數(shù)的全條件后驗分布,由于ρ的全條件后驗分布是非標準的分布函數(shù)形式,因此對參數(shù)ρ采用Metropois-Hastings抽樣,而對其余參數(shù)采用Gibbs抽樣。根據(jù)以上方法在R軟件中生成足夠多的樣本點,將前面一定范圍內(nèi)的樣本點作為預熱,取后面若干值的平均值作為參數(shù)的估計并通過80%置信區(qū)間指導變量選擇。

2 結(jié) 果

2.1 空間自相關(guān)分析 使用Queen鄰接計算空間權(quán)重矩陣,由于海南省在地理上不與任何省份鄰接,空間權(quán)重矩陣對應行全為0元素,為了避免由此造成的計算問題,令海南與廣東、廣西相鄰?？臻g自相關(guān)分析結(jié)果表明2020年布病發(fā)病率全局Moran’sI指數(shù)為0.529,呈正空間自相關(guān)性,具有明顯空間聚集性,Z值為5.549,P值為0.001,具有統(tǒng)計學意義。因此,模型的建立需要考慮空間相關(guān)性。

2.2 風險因素的相關(guān)性分析 由圖1可得,氣候因素內(nèi)、自然環(huán)境因素內(nèi)、畜牧業(yè)因素內(nèi)及各因素間均存在多組相關(guān)性較高的變量,此外由于空間滯后項為對應變量的線性組合,因此各變量及其空間滯后項不可避免地成為了相關(guān)變量。通過計算可得多組變量的Pearson相關(guān)系數(shù)大于0.9且解釋變量相關(guān)矩陣的條件數(shù)遠大于1 000,表明存在嚴重的多重共線性。而多重共線性的存在會導致參數(shù)估計值標準差較大、回歸方程不穩(wěn)定、參數(shù)估計值的正負號與實際不符等問題。

2.3 空間回歸模型分析 表1為兩模型的擬合結(jié)果,從系數(shù)估計結(jié)果初步可得BL-MESS模型篩選了6個變量,除截距項外分別為平均氣溫、平均海拔高度、牧草地面積占比、山羊與綿羊年末存欄量,而其余變量由于置信區(qū)間含零,表明對布病發(fā)病影響不顯著從而被排除。由顯著變量組成的相關(guān)矩陣所計算的條件數(shù)僅為61.8,小于100說明多重共線性程度小。將該6個變量作為自變量重新擬合MESS模型可得,各系數(shù)均顯著且擬合優(yōu)度為0.909,表示布病發(fā)病率有91%左右的變異程度可由該6個變量的變異程度解釋。擬合結(jié)果顯示,MESS模型的擬合優(yōu)度略高于BL-MESS模型,然而由SAR模型與MESS模型空間系數(shù)的對應關(guān)系可得,MESS模型呈現(xiàn)出空間負自相關(guān)性,與實際數(shù)據(jù)相悖。此外,由圖2可得該模型中不少變量系數(shù)具有較長的置信區(qū)間,表明系數(shù)估計精度低、可靠性差。因此,綜合各項指標可得,BL-MESS模型具有更高的可信度。

表1 MESS和BL-MESS模型估計結(jié)果

圖2 MESS與BL-MESS模型下系數(shù)的后驗平均值和相應的80%置信區(qū)間

由于模型包含空間滯后項,模型的回歸系數(shù)并不能直接表征各變量對布病發(fā)病率的影響程度。借鑒LeSage和Pace提出的方法[11],將模型的總體效應分解為直接效應與間接效應。表2結(jié)果顯示,氣候因素中只有平均氣溫的平均直接效應和平均總體效應顯著為正。自然環(huán)境因素來看,平均海拔高度各項效應均顯著為負,而牧草地面積占比的直接效應、總體效應顯著為正。從畜牧業(yè)因素來看,僅有山羊年末存欄量的各項效應及綿羊年末存欄量的直接效應顯著為正,其余因素均不顯著。

表2 BL-MESS模型解釋變量的效應分解

3 討 論

過去國內(nèi)外針對布病的空間分析更多集中于空間聚集性分析,并以此了解疾病在空間上的分布特點與規(guī)律,而近年來結(jié)合空間信息建立空間回歸模型受到越來越多學者的青睞。本研究利用MESS模型在理論與計算上的良好優(yōu)勢并融合貝葉斯LASSO,以中國大陸31個省份為例建立BL-MESS模型,從空間關(guān)系角度探索布病發(fā)病風險的關(guān)鍵影響因素。

本文通過空間回歸模型分析篩選出5個宏觀因素,分別為平均氣溫、平均海拔高度、牧草地面積占比、山羊年末存欄量、綿羊年末存欄量。結(jié)果表明平均氣溫對布病發(fā)病有明顯促進作用,這種現(xiàn)象主要受多種因素影響:首先,溫度是牧草生長發(fā)育的重要因素,適宜的溫度能提高牧草結(jié)實率和產(chǎn)量從而有利于養(yǎng)殖牲畜,為布氏菌提供天然宿主。其次,氣溫對綿羊等牲畜的繁殖有明顯影響,夏季高溫使得母羊妊娠受到影響,易引起胚胎死亡,研究表明流產(chǎn)史是影響布病的重要因素,此外病畜流產(chǎn)物會污染草場與水源,促進布魯氏菌病的傳播,而處理流產(chǎn)母羊也是人感染布病的風險因素。另外,平均海拔高度升高能抑制布病發(fā)病,這也解釋了為何青海和西藏牧草地面積占比較高,而布病發(fā)病率卻相對較低,這在其他研究中也能得到類似的結(jié)論[12]。牧草地面積占比也是布病發(fā)病風險的危險因素,在其他條件不變的情況下,牧草地面積占比越高,越能保證牲畜有足夠的采食量,從而越有利于布病的繁殖與傳播。而山羊、綿羊年末存欄量各項系數(shù)的顯著性則提示了動物防護的重要性,許多國家經(jīng)驗表明控制畜間布病能有效降低人間布病的發(fā)病率。

根據(jù)文獻對比,BL-MESS模型篩選所得的顯著影響因素與已有研究結(jié)果一致,而牛、豬年末存欄量的不相關(guān)性也得到了印證[13],這進一步驗證了該方法能有效識別重要風險因素及無關(guān)變量,且該法相對于MESS模型具有顯著降低參數(shù)標準差、提高估計精度的優(yōu)勢。此外,空間溢出效應反映了鄰近地區(qū)影響本地區(qū)發(fā)病率的途徑,由結(jié)果分析可得低海拔的地區(qū)及山羊養(yǎng)殖規(guī)模較大的省份有布病外溢的風險,從而為控制布病疫情擴散提供了新方向和理論依據(jù)。綜上所述,由于布病深受氣候因素與畜牧業(yè)因素影響,密切關(guān)注氣候變化與動物防護應是今后布病防護的重要內(nèi)容。

已有研究表明氣象因素對布病發(fā)生驅(qū)動作用明顯[14-15],然而除了平均氣溫本研究暫未發(fā)現(xiàn)其它氣候因素與布病發(fā)病率的關(guān)聯(lián)。由于現(xiàn)有數(shù)據(jù)來源限制,目前較難獲得地級行政區(qū)及以下的布病發(fā)病率數(shù)據(jù),這可能使得在全國尺度分析氣候因素對布病發(fā)病率造成一定局限性。

利益沖突：無