新疆大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院(830046)

畢圣賢 別思羽 張輝國 胡錫健△

【提 要】 目的 分析2014-2018年全國31個省、直轄市、自治區(qū)全年日照時(shí)長、平均氣溫與牛、豬、羊年底頭數(shù)對布魯氏菌病(簡稱布病)發(fā)病數(shù)影響的時(shí)空變化特征，為今后制定科學(xué)防治措施提供理論依據(jù)。方法 建立時(shí)空加權(quán)泊松回歸模型，對系數(shù)函數(shù)進(jìn)行置信區(qū)間檢驗(yàn)與可視化分析，揭示變量間相互影響的時(shí)空非平穩(wěn)性，并與地理加權(quán)泊松回歸模型擬合結(jié)果相比較。結(jié)果 全國布病發(fā)病數(shù)時(shí)空分布特征明顯，且與五個宏觀因素密切相關(guān)；時(shí)空加權(quán)泊松回歸模型能有效吸收時(shí)間信息與地理信息，回歸關(guān)系與系數(shù)函數(shù)具有時(shí)空非平穩(wěn)性，擬合優(yōu)度檢驗(yàn)(R2，AICc，RSS)結(jié)果優(yōu)于地理加權(quán)泊松回歸模型。結(jié)論 全國布病防治工作取得一定進(jìn)展，但今后仍需重點(diǎn)關(guān)注北方疫情熱點(diǎn)區(qū)域，并做好南方疫情冷點(diǎn)區(qū)域的監(jiān)測工作，避免疫情向南方蔓延。

布魯氏菌病(brucellosis，簡稱布病)是最常見的人畜共患慢性傳染病之一，通過接觸感染的牛羊等動物與未消毒的奶制品傳播，其病死率低，但易復(fù)發(fā)[1-3]。2000年以來，我國人間布魯氏菌病疫情進(jìn)入了快速擴(kuò)張階段[4]，成為公共衛(wèi)生安全的重大威脅。

目前對布病的流行病學(xué)研究集中在時(shí)空相關(guān)性方面[5-8]，通過時(shí)空回歸模型研究布病發(fā)病與其感染人群所處環(huán)境間時(shí)空異質(zhì)性研究尚未見報(bào)道。Huang等基于傳統(tǒng)地理加權(quán)回歸理論提出了時(shí)空回歸模型理論[9]，并進(jìn)行了實(shí)證分析；張亞慧等與孫舒曼等分別利用地理加權(quán)泊松回歸模型與時(shí)空加權(quán)回歸模型分析了全國艾滋病發(fā)病數(shù)分布特征[10-11]；趙媛通過線性回歸方程、spearman相關(guān)分析與空間自相關(guān)分析研究了寧夏布魯氏菌病的分布特征與影響因素[12]。本文在上述研究基礎(chǔ)上，建立時(shí)空加權(quán)泊松回歸模型，分析2014-2018年全國31個省、直轄市、自治區(qū)布病發(fā)病數(shù)與全年日照時(shí)長、平均氣溫、與牛豬羊年底頭數(shù)的異質(zhì)性關(guān)系，揭示變量間相互影響的時(shí)空非平穩(wěn)性，并與地理加權(quán)回歸模型擬合結(jié)果相比較。

資料來源與方法

1.資料來源與預(yù)處理

2014-2018年全國各地區(qū)布病發(fā)病率(每10萬人)來自2015-2017年《中國衛(wèi)生和計(jì)劃生育統(tǒng)計(jì)年鑒》與2018-2019年《中國衛(wèi)生健康統(tǒng)計(jì)年鑒》，各地區(qū)年底人口數(shù)(萬人)、降雨量(毫米)、全年日照時(shí)長(小時(shí))、平均氣溫(攝氏度)、平均相對濕度(%)，牛、豬、羊年底頭數(shù)(萬頭)來自國家統(tǒng)計(jì)局官網(wǎng)《中國統(tǒng)計(jì)年鑒》，布病發(fā)病人數(shù)由發(fā)病率與各地區(qū)人口數(shù)相乘得到。31個省、直轄市、自治區(qū)空間位置的平面坐標(biāo)通過WGS_1984_Albers投影坐標(biāo)系統(tǒng)轉(zhuǎn)換得到。對降雨量、全年日照時(shí)長等七個宏觀因素進(jìn)行標(biāo)準(zhǔn)化處理，消除單位影響，建立線性回歸方程，并通過逐步回歸的方法篩選自變量，顯著性水平取0.05，最終保留下的五個自變量是全年日照時(shí)長、平均氣溫與牛、豬、羊年底頭數(shù)。

2.研究方法

廣義線性泊松回歸模型用來研究疾病發(fā)病數(shù)與各影響因素間的關(guān)系。本文以五個宏觀因素為解釋變量，2014-2018年各省、直轄市、自治區(qū)布病發(fā)病數(shù)為響應(yīng)變量，分別建立時(shí)空加權(quán)泊松回歸模型(GTWPR)與地理加權(quán)泊松回歸模型(GWPR)，通過研究不同模型系數(shù)函數(shù)間的關(guān)系，揭示回歸關(guān)系與回歸系數(shù)的時(shí)空非平穩(wěn)性。

GTWPR模型[9]選取對數(shù)函數(shù)作為連接函數(shù)，其線性預(yù)測部分表達(dá)形式為：

ηij=ln(uij)

其中uij為第i地區(qū)第j年布病發(fā)病數(shù)，Xij1、Xij2、Xij3、Xij4、Xij5分別表示第i地區(qū)第j年的全年日照時(shí)長、平均氣溫與牛、豬、羊年底頭數(shù)指標(biāo)，βp(ui，vi，tj)，p=1，2，3，4，5表征五個宏觀因素對布病發(fā)病數(shù)的影響程度，(ui，vi，tj)反映系數(shù)函數(shù)的時(shí)空位置。

GWPR模型[13]線性預(yù)測部分表達(dá)式為：

ηi=ln(ui)

ui為第i地區(qū)布病發(fā)病數(shù)，系數(shù)函數(shù)只包含地理信息。兩類加權(quán)回歸模型核光滑方法均選用Gauss核函數(shù)，光滑參數(shù)通過AICc準(zhǔn)則確定，非平穩(wěn)性檢驗(yàn)通過bootstrap置信區(qū)間算法完成。

GTWPR模型與GWPR模型系數(shù)函數(shù)的估計(jì)與bootstrap算法檢驗(yàn)由R軟件完成；2014-2018年全國布魯氏菌病時(shí)空分布圖通過excel 2013繪制；其他作圖部分通過argis 10.2軟件實(shí)現(xiàn)。

結(jié) 果

1.布病發(fā)病數(shù)的時(shí)空分布

2014-2018年全國布病報(bào)告病數(shù)共239228例，總體呈下降趨勢。其中山西、陜西呈遞減趨勢；遼寧、黑龍江、河南、甘肅、寧夏、新疆在2015年達(dá)到峰值后穩(wěn)步下降；河北、吉林報(bào)告病數(shù)逐年遞減后在2018年出現(xiàn)小幅增長；內(nèi)蒙古與山東布病疫情在2016年分別達(dá)到峰谷與峰頂后最終恢復(fù)到2014年?duì)顟B(tài)；南方各地區(qū)報(bào)告病例數(shù)則一直維持在較低水平?？傮w來看，全國布病疫情時(shí)空差異明顯，呈北重南輕分布。以甘肅、寧夏、陜西、河南、山東等地區(qū)構(gòu)成的疫情緩沖帶，將全國劃分成新疆、內(nèi)蒙古與山西、黑龍江與吉林等三個北方疫情熱點(diǎn)區(qū)域，與南方疫情冷點(diǎn)區(qū)域，結(jié)果見圖1。

圖1 2014-2018年全國布病時(shí)空分布圖

2.系數(shù)的估計(jì)

兩類時(shí)空加權(quán)回歸模型系數(shù)函數(shù)估計(jì)范圍見表1。結(jié)果顯示，GTWPR模型五個系數(shù)函數(shù)估計(jì)范圍均表現(xiàn)從負(fù)到正，且四分位間距較大，提示在該模型中影響布病發(fā)病的五個宏觀因素間存在較強(qiáng)的時(shí)空非平穩(wěn)性；在GWPR模型中，β1、β4與β5系數(shù)函數(shù)估計(jì)結(jié)果表現(xiàn)正值，β2與β3系數(shù)函數(shù)表現(xiàn)負(fù)值，提示在該模型中全年日照時(shí)長、豬年底頭數(shù)與羊年底頭數(shù)對布病發(fā)病起正相關(guān)作用，平均氣溫與牛年底頭數(shù)對布病發(fā)病起負(fù)相關(guān)作用。同時(shí)GWPR模型系數(shù)函數(shù)的取值范圍都在GTWPR模型的系數(shù)函數(shù)取值空間內(nèi)，提示地理加權(quán)模型估計(jì)結(jié)果能在時(shí)空加權(quán)模型中得到解釋。

表1 兩種模型系數(shù)估計(jì)范圍

3.非平穩(wěn)性檢驗(yàn)

檢驗(yàn)回歸關(guān)系的非平穩(wěn)性，實(shí)質(zhì)就是檢驗(yàn)?zāi)Ｐ蛥?shù)是否為零，即檢驗(yàn)假設(shè)：H0：βk(u，v)=0或βk(u，v，t)=0對某個參數(shù)成立。本文通過比較模型所有系數(shù)函數(shù)置信區(qū)間是否含零進(jìn)行總體顯著性檢驗(yàn)，GTWPR模型與GWPR模型系數(shù)函數(shù)置信區(qū)間通過bootstrap方法得到，置信度為95%，結(jié)果見表2。結(jié)果顯示，兩類加權(quán)回歸模型都通過總體顯著性檢驗(yàn)，存在時(shí)空非平穩(wěn)性。同時(shí)，兩類模型的β1、β4與β5系數(shù)函數(shù)置信區(qū)間表現(xiàn)正值，β2與β3系數(shù)函數(shù)置信區(qū)間表現(xiàn)負(fù)值，提示五個宏觀因素對布病發(fā)病存在不同的正負(fù)相關(guān)作用。

檢驗(yàn)時(shí)空加權(quán)回歸模型中回歸參數(shù)的非平穩(wěn)性，即檢驗(yàn)回歸參數(shù)在不同時(shí)間位置與空間位置上是否有相同取值。推廣Alireza等[14]通過比較置信區(qū)間寬度檢驗(yàn)系數(shù)函數(shù)空間非平穩(wěn)性的方法，檢驗(yàn)GTWPR模型與GWPR模型系數(shù)函數(shù)的非平穩(wěn)性，結(jié)果見表3。檢驗(yàn)結(jié)果顯示，GTWPR模型系數(shù)函數(shù)非平穩(wěn)性檢驗(yàn)指標(biāo)優(yōu)于GWPR模型。

表2 回歸關(guān)系的非平穩(wěn)性檢驗(yàn)

表3 系數(shù)函數(shù)的非平穩(wěn)性檢驗(yàn)

4.擬合優(yōu)度檢驗(yàn)

擬合優(yōu)度檢驗(yàn)是衡量模型擬合優(yōu)良性的重要準(zhǔn)則[15]。分別對GTWPR模型與GWPR模型進(jìn)行擬合優(yōu)度檢驗(yàn)，結(jié)果見表4。結(jié)果顯示，GTWPR模型AICc信息準(zhǔn)則、殘差平方和RSS與均方誤差MSE小于GWPR模型，可決定系數(shù)大于GWPR模型，表明GTWPR模型能更好擬合布病發(fā)病數(shù)與五個宏觀因素的回歸關(guān)系。

表4 擬合優(yōu)度檢驗(yàn)

5.系數(shù)的可視化分析

對GTWPR模型五個系數(shù)函數(shù)的估計(jì)均值進(jìn)行可視化分析。結(jié)果顯示，五個宏觀因素在不同空間位置上對布病發(fā)病的影響存在顯著差異。全年日照時(shí)長與布病發(fā)病呈正相關(guān)關(guān)系的是甘肅、青海等17個地區(qū)，呈負(fù)相關(guān)關(guān)系的是新疆、西藏、內(nèi)蒙古等6個地區(qū)，東北、華北與江西、湖南等地區(qū)相關(guān)關(guān)系不明顯，提示多數(shù)地區(qū)日照時(shí)長對布病發(fā)病有正相關(guān)作用。平均氣溫與布病發(fā)病呈正相關(guān)關(guān)系的是西藏、陜西等5個地區(qū)，寧夏、湖北、湖南、上海相關(guān)關(guān)系不顯著，布病報(bào)告病例數(shù)最多的內(nèi)蒙古、新疆、山西、河北等地區(qū)，均表現(xiàn)強(qiáng)負(fù)相關(guān)關(guān)系，提示較低的氣溫是影響布病發(fā)病的重要因素。牛年底頭數(shù)與豬年底頭數(shù)系數(shù)函數(shù)從西向東分別呈遞增與遞減趨勢，西藏與西北地區(qū)在牛年底頭數(shù)上的負(fù)相關(guān)作用最明顯，西北地區(qū)在豬年底頭數(shù)上的正相關(guān)關(guān)系最明顯。羊年底頭數(shù)方面，與布病發(fā)病呈強(qiáng)負(fù)相關(guān)關(guān)系的是西藏、云南、貴州、廣西、廣東、海南等疫情冷點(diǎn)區(qū)域，呈強(qiáng)正相關(guān)關(guān)系的是甘肅、青海兩省，東南沿海地區(qū)按與疫情熱點(diǎn)區(qū)域的距離正相關(guān)作用不斷增強(qiáng)。畜牧業(yè)的三項(xiàng)宏觀因素系數(shù)函數(shù)表現(xiàn)結(jié)果提示，西藏、云南、廣東等南方地區(qū)受羊年底頭數(shù)影響最大，處于疫情冷點(diǎn)；新疆、內(nèi)蒙古、河北、山西等北方疫情熱點(diǎn)地區(qū)單項(xiàng)畜牧業(yè)因素系數(shù)函數(shù)特征不顯著；東南沿海地區(qū)受牛、羊年底頭數(shù)的影響，疫情防控形勢嚴(yán)峻，但豬年底頭數(shù)對東南地區(qū)布病疫情影響不大，可能與該地區(qū)成熟的圈養(yǎng)模式有關(guān)，減少了家畜感染布魯氏菌病的機(jī)會[5]。

討 論

時(shí)空加權(quán)泊松回歸模型將時(shí)間信息與空間信息引入到布魯氏菌病發(fā)病的回歸關(guān)系研究中，探究布病發(fā)病的時(shí)空分布特征與影響因素，能充分反映回歸關(guān)系的時(shí)空異質(zhì)性。逐步回歸方法篩選出五個宏觀因素，即全年日照時(shí)長、平均氣溫與牛、豬、羊年底頭數(shù)對布病發(fā)病數(shù)均有顯著影響。

時(shí)空分布圖顯示，全國各地區(qū)布病發(fā)病數(shù)具有時(shí)空聚集性，在時(shí)間上呈遞減趨勢，在空間上呈北重南輕、中間過渡分布態(tài)勢，表明隨著2014年加強(qiáng)布病疫苗接種工作以來，全國布病疫情時(shí)空擴(kuò)散趨勢已得到較好控制，但北方尤其是內(nèi)蒙古地區(qū)仍是高發(fā)區(qū)域，存在疫情擴(kuò)散風(fēng)險(xiǎn)。系數(shù)估計(jì)結(jié)果、非平穩(wěn)性檢驗(yàn)與擬合優(yōu)度檢驗(yàn)表明，GTWPR模型能更好擬合布病發(fā)病數(shù)與五個宏觀因素的回歸關(guān)系，且GWPR模型系數(shù)估計(jì)結(jié)果能在GTWPR模型中得到合理解釋，進(jìn)一步說明了GTWPR模型在解釋布病發(fā)病異質(zhì)性特征方面的優(yōu)越性。

異質(zhì)性分析結(jié)果表明，西南地區(qū)受羊年底頭數(shù)較少影響，處于疫情冷點(diǎn)；東南沿海地區(qū)受牛、羊年底頭數(shù)共同影響，雖然近幾年疫情防控壓力增加，但因其牛、羊基數(shù)較少，仍處于疫情冷點(diǎn)；北方尤其是新疆、內(nèi)蒙古等疫情熱點(diǎn)區(qū)域，受五個宏觀因素綜合影響，人間布病的感染居高不下。其原因可能是(1)畜間傳播是導(dǎo)致布病發(fā)病的主要因素，北方尤其是內(nèi)蒙古、新疆等地畜牧業(yè)發(fā)達(dá)，人畜同院、散養(yǎng)放養(yǎng)的傳統(tǒng)飼養(yǎng)模式導(dǎo)致布病疫情多發(fā)；(2)北方地區(qū)充足的日照與較適宜的溫度有利于布魯氏菌的存活繁殖，也為發(fā)展畜牧業(yè)提供良好條件，從而增加傳染源與傳播的宿主密度，為布病疫情的擴(kuò)散帶去了便利[4]；(3)人間布病的感染不僅受畜牧業(yè)的影響，還受環(huán)境因素、經(jīng)濟(jì)因素、社會因素與相關(guān)行業(yè)從業(yè)人員衛(wèi)生防護(hù)意識、身體免疫力等的綜合影響，東南沿海地區(qū)經(jīng)濟(jì)發(fā)達(dá)，從事牛羊養(yǎng)殖與屠宰職業(yè)人群較多，增加了布病防控壓力。五個宏觀因素的整體分析表明，畜間發(fā)病是引起人間布病的必要不充分條件，這與楊祖榮等[16]和吳晨等[7]通過分布滯后線性模型與空間回歸模型研究布病發(fā)病影響因素的結(jié)論一致。

圖2 GTWPR模型系數(shù)函數(shù)分布

布病防控是一項(xiàng)長期工程，需聯(lián)系五項(xiàng)宏觀因素并結(jié)合各地區(qū)實(shí)際采取有效措施。要密切監(jiān)測北方牧區(qū)畜間疫情，提高監(jiān)測水平與疫苗接種率，對病畜要及時(shí)捕殺、焚燒或填埋，并嚴(yán)格做好圈舍的消毒工作；要增強(qiáng)東南沿海地區(qū)牛羊暴露人群的健康教育，提升防護(hù)意識，同時(shí)加強(qiáng)西南地區(qū)跨地區(qū)畜牧產(chǎn)品檢疫與管理，切實(shí)做好布病防控工作[17]。

綜上所述，2014-2018年全國布病發(fā)病呈遞減趨勢，時(shí)空分布特征明顯，且受氣候因素與畜牧業(yè)因素等的多重影響。布病防控工作雖取得一定進(jìn)展，但仍需重點(diǎn)關(guān)注新疆、內(nèi)蒙古等疫情熱點(diǎn)區(qū)域的發(fā)病情況，并做好南方疫情冷點(diǎn)區(qū)域的監(jiān)測工作，避免疫情向南方蔓延。