基于時空地理加權(quán)回歸模型的中國肺結(jié)核發(fā)病情況及影響因素研究

趙明揚，周乾宇，王榮榮，王宗熹，何雯倩，張文森，張恒榛，田卓旸，吳柯，王碧瑤，孫長青，*

肺結(jié)核是由結(jié)核分枝桿菌引起的具有強烈傳染性的慢性呼吸系統(tǒng)傳染病，是危害人體健康的重大傳染病之一［1］。肺結(jié)核患者在排菌期，即痰菌陽性時具有傳染性，1例未經(jīng)治療的排菌期活動性肺結(jié)核患者1年內(nèi)可感染15～20名接觸者。2019年全球結(jié)核報告顯示，全世界每年約有1 000萬人感染結(jié)核病，潛伏性結(jié)核感染人數(shù)占世界人口總數(shù)的1/4左右；中國新發(fā)結(jié)核病病例833 000例，發(fā)病率為58/10萬，位居世界第三，是結(jié)核病高負擔國家，肺結(jié)核防控形勢依舊嚴峻［2］。

隨著地理信息技術(shù)的發(fā)展和完善，大量基于空間計量學的空間統(tǒng)計分析方法被提出，地理信息系統(tǒng)和空間統(tǒng)計學在各領(lǐng)域得到廣泛應(yīng)用。英國學者Fotheringham提出的研究空間關(guān)系和空間相關(guān)關(guān)系的地理加權(quán)回歸（geographical weighted regression，GWR）模型可以直觀地探測空間關(guān)系的非平穩(wěn)性［3］，并在多學科領(lǐng)域得到廣泛應(yīng)用［4-6］。但GWR模型只將數(shù)據(jù)的空間特性納入模型，而忽略了時間特性對其的影響，在處理某些時間特性和空間特性均突出的數(shù)據(jù)資料時存在一定局限性。因此，香港大學黃波教授在GWR模型的基礎(chǔ)上提出了時空地理加權(quán)回歸（geographically and temporally weighted regression，GTWR）模型［7］。此模型能夠更好地對數(shù)據(jù)的時空分布及特性進行評估，有效解決了回歸模型的時空非平穩(wěn)性，一經(jīng)提出便得到廣泛關(guān)注，并在眾多學者的應(yīng)用過程中不斷完善［8-10］。

近年來，國內(nèi)外學者對肺結(jié)核的時空分布進行大量研究表明，肺結(jié)核具有較強的時間特性和空間特性［11-13］。但現(xiàn)有研究更多的是對肺結(jié)核發(fā)病的影響因素進行獨立的時間或空間回歸分析，研究結(jié)果存在局限性，例如使用差分整合滑動平均自回歸模型（autoregressive integrated moving average model，ARIMA）對肺結(jié)核數(shù)據(jù)進行時間序列分析時，忽略了其對應(yīng)的空間因素的影響［14-15］；在使用空間聚類和地理加權(quán)回歸時，肺結(jié)核數(shù)據(jù)所具有的時間特性便無法體現(xiàn)［16］。因此，將時間因素和空間因素同時進行回歸分析，采用GTWR模型對影響肺結(jié)核發(fā)病的影響因素進行擬合分析是十分必要的。

綜上所述，本文將使用GTWR模型探索中國肺結(jié)核分布的時間和空間異質(zhì)性，并分析肺結(jié)核發(fā)病情況與氣象和空氣質(zhì)量因素在時間和空間上的相關(guān)關(guān)系，為制訂相應(yīng)結(jié)核病防控措施提供科學參考。

1 材料與方法

1.1 數(shù)據(jù)來源 本研究使用的肺結(jié)核發(fā)病情況相關(guān)數(shù)據(jù)來源于公共衛(wèi)生科學數(shù)據(jù)中心（https://www.phsciencedata.cn/Share/）中的2016—2018年全國分地區(qū)（除香港、澳門、臺灣外的31個省份）肺結(jié)核分月統(tǒng)計數(shù)據(jù)，選取的指標為發(fā)病率（每10萬人中病例數(shù)）。氣象資料數(shù)據(jù)（包括月平均氣溫、濕度、風速等）來源于天氣后報網(wǎng)（http://www.tianqihoubao.com/）。空氣質(zhì)量指數(shù)數(shù)據(jù)（包括PM2.5、PM10、SO2、NO2、CO、O3等的月平均濃度）來源于空氣質(zhì)量指數(shù)歷史數(shù)據(jù)網(wǎng)站（https://www.aqistudy.cn/historydata/）。本研究所使用地理信息資料來源于國家基礎(chǔ)地理信息中心（http://www.ngcc.cn/），文中涉及的中國地圖均基于自然資源部標準地圖服務(wù)網(wǎng)站 （http://bzdt.ch.mnr.gov.cn/）下載的審圖號為GS（2019）1823號的標準地圖制作。

本研究使用數(shù)據(jù)均來源于公共數(shù)據(jù)庫，不適用倫理審查。

1.2 研究方法

1.2.1 多重共線性 多重共線性是指線性回歸模型中的解釋變量之間由于存在精確相關(guān)關(guān)系或高度相關(guān)關(guān)系而使模型估計失真或難以估計準確。因此，為保證回歸模型的合理性，應(yīng)在構(gòu)建模型前對備選自變量的共線性進行分析。方差膨脹因子（VIF）是常用的檢測自變量之間多重共線性的指標［17］。本研究將采用VIF對肺結(jié)核發(fā)病情況與影響因素之間的關(guān)系進行共線性檢驗，以避免由于影響因素之間的高度共線性而影響回歸分析結(jié)果，其計算公式如下：

其中，r為線性回歸中的決定系數(shù)，反映了回歸方程解釋因變量變化的百分比。VIF越大，說明解釋變量之間存在共線性的可能性越大，若VIF均在0～10，則影響因素之間不存在高度共線性，可直接進行回歸分析［18］。

1.2.2 空間自相關(guān) 使用空間計量學方法的前提是樣本數(shù)據(jù)之間存在空間異質(zhì)性，因此在構(gòu)建GWR和GTWR模型前需對自變量進行空間自相關(guān)分析。通常使用莫蘭指數(shù)（Moran'sⅠ）進行全局空間自相關(guān)分析，以確定所研究樣本點的某一屬性值與領(lǐng)域內(nèi)其他樣本點相同屬性值在空間上是否關(guān)聯(lián)。本研究將通過計算全局Moran'sⅠ以確定肺結(jié)核發(fā)病情況的空間自相關(guān)性，其計算公式如下：

其中，S0為所有樣本點之間空間權(quán)重的總和，zi為樣本點i的某一屬性值與其平均值的偏差。Moran'sⅠ的取值在-1～1，若指數(shù)為正值則表示樣本的某一屬性值在空間上呈現(xiàn)聚集狀態(tài)，且指數(shù)越趨近1則聚集程度越強；若指數(shù)為負值則表示樣本屬性值呈離散分布；指數(shù)為0則表示樣本屬性值呈隨機分布，無顯著特征。

1.2.3 回歸模型構(gòu)建 本研究分別構(gòu)建最小二乘法（OLS）模型、GWR模型和GTWR模型對肺結(jié)核發(fā)病情況進行實證分析，并比較模型優(yōu)度以確定GTWR模型是否為處理肺結(jié)核數(shù)據(jù)的最佳模型。

OLS模型是常用的傳統(tǒng)線性回歸模型，該模型僅對參數(shù)進行了平均或全局意義上的估計，但無法體現(xiàn)各參數(shù)在空間上的非平穩(wěn)性。模型計算公式為：

其中，Yi表示第i個樣本點的因變量，β0表示線性回歸方程的截距，βk表示第k個自變量的回歸系數(shù)，Xik表示第i個樣本點的第k個自變量，εi表示隨機誤差。

GWR模型是基于傳統(tǒng)線性回歸模型改進后的模型，其主要優(yōu)勢是能夠?qū)⒖臻g權(quán)重矩陣運用在線性回歸模型之中，可以更好的展現(xiàn)結(jié)果的空間結(jié)構(gòu)分異。模型計算公式為：

其中，Yi表示第i個樣本點的因變量，ui表示第i個樣本點的經(jīng)度坐標，vi表示第i個樣本點的緯度坐標，（ui，vi）表示第i個樣本點的空間經(jīng)緯度坐標，β0（ui，vi）表示第i個樣本點的常數(shù)項，βk（ui，vi）表示第k個自變量在第i個樣本點的回歸系數(shù)，Xik表示第i個樣本點的第k個自變量，εi表示隨機誤差。

GTWR模型是在GWR模型的基礎(chǔ)上將時間賦值到局部樣本點數(shù)據(jù)集上，求解局部樣本點i的參數(shù)，充分利用樣本數(shù)據(jù)的時間特性，提高參數(shù)估計的準確性。模型計算公式為：

其中，Yi表示第i個樣本點的因變量，ui表示第i個樣本點的經(jīng)度坐標，vi表示第i個樣本點的緯度坐標，ti表示第i個樣本點的時間坐標，（ui，vi，ti）表示第i個樣本點的時空維度坐標，β0（ui，vi，ti）表示第i個樣本點的常數(shù)項，βk（ui，vi，ti）表示第k個解釋變量在第i個樣本點的回歸系數(shù)，Xik表示第i個樣本點的第k個自變量，εi表示隨機誤差。

GWR與GTWR模型的參數(shù)方法如下：

其中，空間權(quán)重矩陣W是由空間帶寬、核函數(shù)、距離計算公式共同決定。根據(jù)既往文獻［7］，本研究將基于最小交叉驗證（CV）值、高斯（Gaussian）核函數(shù)和歐式距離（Euclidean distance）來共同構(gòu)建模型。模型優(yōu)度通過比較修正后的赤池信息量（AICc）與R2值來評估，R2值越大，AICc值越小，說明自變量對因變量的解釋度越強。

1.2.4 統(tǒng)計分析 使用均數(shù)、最小值、最大值、四分位數(shù)間距來描述GTWR模型的擬合系數(shù)?；贕TWR模型的擬合系數(shù)，分別繪制各個變量的核密度圖和時空分布圖。使用自然斷裂點法對相似度較高的數(shù)據(jù)進行分組，同時強行將“0”設(shè)置為區(qū)間值，以區(qū)分正系數(shù)和負系數(shù)，當擬合系數(shù)為正值時，表示自變量對因變量具有促進作用；當擬合系數(shù)為負值時，表示自變量對因變量具有抑制作用，且擬合系數(shù)的絕對值越大，作用程度越大。

本研究使用R 4.1.3軟件進行統(tǒng)計描述，使用Arc GIS 10.7軟件進行模型參數(shù)估計和模型構(gòu)建。

2 結(jié)果

2.1 肺結(jié)核發(fā)病情況的時空分布 2016—2018年全國肺結(jié)核發(fā)病率的空間分布見圖1：我國肺結(jié)核總發(fā)病率在逐年下降，且空間分布較為集中。肺結(jié)核發(fā)病率較高的地區(qū)主要集中在新疆、四川、西藏、青海、貴州、廣西等。其中，新疆的肺結(jié)核發(fā)病率連續(xù)3年處于最高水平；四川的肺結(jié)核發(fā)病率在2016年處于較高水平，但隨后2年發(fā)病率大幅下降；西藏、青海的肺結(jié)核發(fā)病率在2016年處于較低水平，但隨后2年發(fā)病率大幅增加。肺結(jié)核發(fā)病率較低的地區(qū)主要集中在寧夏、天津、上海、北京、海南等。

圖1 全國肺結(jié)核發(fā)病率空間分布圖Figure 1 Spatial distribution of pulmonary tuberculosis incidence in China during 2016—2018

2.2 模型對比結(jié)果 經(jīng)過多重共線性和空間自相關(guān)檢驗（全局Moran's I= 0.376），剔除氣溫這一具有強共線性的變量（VIFTmax=48.01，VIFTmin=48.34）后，分別建立OLS、GWR和GTWR模型，評估并比較模型優(yōu)度，結(jié)果如表1所示。GTWR模型的R2值與AdjustedR2值均比OLS和GWR模型要高，同時GTWR模型的AICc值均比OLS和GWR模型要小，表明GTWR模型能更好地解釋自變量對肺結(jié)核發(fā)病情況的影響，更能解釋具有時空特征的數(shù)據(jù)。

表1 OLS、GWR、GTWR模型比較的結(jié)果Table 1 Values of AICc，R2 and adjusted R2 of OLS，GWR and GTWR models

2.3 擬合系數(shù)的時空特性 使用均數(shù)、最小值、最大值、四分位數(shù)間距來描述GTWR模型的擬合系數(shù)，結(jié)果如表2所示。基于GTWR模型的擬合系數(shù)結(jié)果繪制每個變量的核密度圖，結(jié)果如圖2所示。擬合系數(shù)的核密度圖結(jié)果顯示：自變量風速呈現(xiàn)多峰分布，主峰約為-0.5；自變量濕度呈現(xiàn)多峰分布，主峰約為0.01，且各峰系數(shù)均大于0；自變量PM2.5呈現(xiàn)左偏峰分布，主峰約為-0.01；自變量PM10呈現(xiàn)單峰分布，主峰約為0.005；自變量SO2呈現(xiàn)多峰分布，主峰約為0；自變量NO2呈現(xiàn)多峰分布，主峰約為-0.03；自變量CO呈現(xiàn)右偏峰分布，主峰約為-0.5；自變量O3呈現(xiàn)多峰分布，主峰約為0。

圖2 各變量擬合系數(shù)的核密度分布圖Figure 2 Kernel density estimation of fitting coefficients of each variable

表2 GTWR模型的擬合系數(shù)Table 2 Fitting coefficients of meteorological and air quality factors included in GTWR model

基于GTWR模型的擬合系數(shù)，分別繪制各個變量的時空分布圖（圖3）。擬合系數(shù)的時空分布圖結(jié)果顯示：自變量風速呈現(xiàn)中部和東南部的擬合系數(shù)較低、東北和西部的擬合系數(shù)較高的分布格局，其中風速對青海、甘肅、湖南、江西等地區(qū)的肺結(jié)核發(fā)病率具有顯著的抑制作用，對西藏、新疆、黑龍江等地區(qū)的肺結(jié)核發(fā)病率具有顯著的促進作用；自變量濕度呈現(xiàn)西部和北部的擬合系數(shù)較低、東南和南部的擬合系數(shù)較高的分布格局，其中濕度對海南、廣東、廣西等地區(qū)的肺結(jié)核發(fā)病率具有顯著的促進作用；自變量PM2.5呈現(xiàn)東北部和中東部的擬合系數(shù)較低、西南部和西部的擬合系數(shù)較高的分布格局，其中PM2.5對廣西、云南、新疆、西藏等地區(qū)的肺結(jié)核發(fā)病率具有顯著的促進作用；自變量PM10呈現(xiàn)中部和西部的擬合系數(shù)較低、東南部和中東部的擬合系數(shù)較高的分布格局，其中PM10對浙江、上海、福建等地區(qū)的肺結(jié)核發(fā)病率具有顯著的促進作用；自變量SO2呈現(xiàn)中部和北部的擬合系數(shù)較低、西部和南部的擬合系數(shù)較高的分布格局，其中SO2對西藏、廣東、新疆等地區(qū)的肺結(jié)核發(fā)病率具有顯著的促進作用；自變量NO2呈現(xiàn)西北部和西南部的擬合系數(shù)較低，東北部、中部的擬合系數(shù)較高的分布格局，其中NO2對黑龍江、吉林等地區(qū)的肺結(jié)核發(fā)病率具有顯著的促進作用；自變量CO呈現(xiàn)西部和西南部的擬合系數(shù)較低、東南部和中西部的擬合系數(shù)較高的分布格局，其中CO對福建、青海、甘肅等地區(qū)的肺結(jié)核發(fā)病率具有顯著的促進作用；自變量O3呈現(xiàn)西北部和南部的擬合系數(shù)較低、東北部和西南部的擬合系數(shù)較高的分布格局，其中O3對新疆、浙江、福建等地區(qū)的肺結(jié)核發(fā)病率具有顯著的促進作用。

圖3 各變量擬合系數(shù)的時空分布圖Figure 3 Spatial and temporal distribution of fitting coefficients of each variable

3 討論

本研究構(gòu)建模型前預(yù)先進行了多重共線性檢驗，結(jié)果顯示，自變量氣溫具有強共線性。既往研究也表明氣溫與其他呼吸系統(tǒng)疾病危險因素具有交叉協(xié)同效應(yīng)。馬盼等［19］研究顯示，高溫、高濕共同作用下將對兒童的呼吸系統(tǒng)健康造成更嚴重的后果，而低溫、低濕共同作用下將增加老年人的呼吸系統(tǒng)疾病患病風險。一項基于醫(yī)院門診數(shù)據(jù)的時間序列分析結(jié)果顯示，氣溫和各種污染物濃度之間存在明顯的交互作用，尤其是低溫與高濃度污染物共同作用下，對呼吸系統(tǒng)門診就診人數(shù)的影響最大［20］。由于多重共線性數(shù)據(jù)可能會影響GTWR模型分析的精度，本研究未將氣溫這一具有強共線性的變量納入模型。

3.1 各變量擬合系數(shù)的核密度分布圖分析 各變量擬合系數(shù)的核密度圖結(jié)果顯示：自變量風速的主峰系數(shù)為負值，且各峰系數(shù)有正有負，表明自變量的增加對大多數(shù)地區(qū)的肺結(jié)核發(fā)病情況呈現(xiàn)顯著的保護作用，但在部分地區(qū)仍具有顯著的促進作用；自變量濕度、SO2的主峰系數(shù)為正值，且各峰系數(shù)均＞0，表明自變量的增加將顯著促進肺結(jié)核發(fā)病率的增加；自變量PM2.5、NO2、CO、O3呈現(xiàn)多峰分布，且各峰系數(shù)有正有負，表明自變量的增加對大多數(shù)地區(qū)的肺結(jié)核發(fā)病率的增加呈現(xiàn)顯著的促進作用；自變量PM10呈現(xiàn)單峰分布，主峰系數(shù)為正值，表明自變量的增加將顯著促進大多數(shù)地區(qū)肺結(jié)核發(fā)病率的增加。

3.2 各變量擬合系數(shù)的時空分布圖分析 GTWR模型的顯著優(yōu)勢是其擬合系數(shù)可以反映自變量對因變量影響的時空變異特征。本研究分別繪制了各變量擬合系數(shù)的時空分布圖，結(jié)果顯示：自變量風速對青海、甘肅、湖南、江西等城市的肺結(jié)核發(fā)病率的增加具有顯著的抑制作用，對西藏、新疆、黑龍江等地區(qū)的肺結(jié)核發(fā)病率的增加具有顯著的促進作用。既往研究表明，風速對肺結(jié)核發(fā)病存在促進作用［21］，較大的風速可降低氣溫，從而間接增加肺結(jié)核發(fā)病風險。針對本研究的雙向作用，可能存在的原因是：西藏、新疆屬于高海拔地區(qū)，東三省的平均氣溫本就偏低，這些地區(qū)更易受到風速的影響導(dǎo)致氣溫下降，從而使肺結(jié)核發(fā)病風險增加，因此這類地區(qū)應(yīng)在大風天氣注意保暖以預(yù)防疾病。而湖南、江西等地的平均溫度較高，溫度不易受風速影響，較大風速還會降低其空氣污染物濃度［22］，進而使肺結(jié)核發(fā)病風險降低。

自變量濕度對海南、廣東、廣西等地區(qū)的肺結(jié)核發(fā)病率增加具有顯著的促進作用。這與既往研究中較高濕度能增加肺結(jié)核發(fā)病風險的結(jié)論相符合［18，21］?？赡艽嬖诘脑蚴牵涸谝欢ǚ秶鷥?nèi)，較高的濕度將增加結(jié)核桿菌在空氣中停留存活的時間，從而增加了人群感染結(jié)核桿菌的風險。因此，針對這類地區(qū)應(yīng)注意對家具或生活用品的勤加晾曬，盡量保持生活環(huán)境的干燥。

自變量PM10對浙江、上海、福建等地區(qū)的肺結(jié)核發(fā)病率具有顯著的促進作用。這與既往研究中PM10暴露濃度增加將導(dǎo)致肺結(jié)核發(fā)病風險增加的結(jié)論相符［23］?？赡艽嬖诘脑蚴牵篜M10等可吸入顆粒物能夠伴隨呼吸進入人體，并在肺部沉積，可能會損傷肺泡及黏膜，導(dǎo)致肺部組織慢性纖維化等情況發(fā)生，降低肺部抵抗力；此外，PM10等空氣污染物還可以與空氣微生物結(jié)合［24］，明顯增加結(jié)核病菌侵入人體的概率。因此，針對這類地區(qū)應(yīng)注意在空氣污染程度較高時避免戶外活動。

既往研究表明空氣污染物（PM2.5、SO2、NO2、CO、O3）的濃度增加會顯著增加肺結(jié)核的發(fā)病風險［23，25-26］，本研究的大部分結(jié)果也能夠與該結(jié)論相互驗證。但上述幾類空氣污染物作為自變量也存在部分擬合系數(shù)為負數(shù)的現(xiàn)象，這部分結(jié)果與既往研究結(jié)論相悖。造成這一現(xiàn)象的原因可能是：由上述結(jié)論可知，擬合系數(shù)為負數(shù)的城市主要為江蘇、上海、廣東等沿海地區(qū)和西藏、新疆、青海等高海拔地區(qū)，此類地區(qū)的空氣質(zhì)量較好，空氣污染物濃度較低，因此空氣污染物對于此類地區(qū)肺結(jié)核發(fā)病情況的促進作用不顯著。但模型構(gòu)建過程中，在與其他地區(qū)、其他變量共同分析時，尤其是與空氣質(zhì)量較差地區(qū)做對比，且存在其他保護因素時，此類自變量在承擔部分結(jié)果解釋的功能后，最終導(dǎo)致過低濃度的空氣污染物呈現(xiàn)負擬合系數(shù)的現(xiàn)象。針對大部分城市的空氣污染物對肺結(jié)核發(fā)病人數(shù)的促進作用，應(yīng)采取積極措施改善空氣質(zhì)量，降低空氣污染物濃度；在空氣污染物濃度較高的時候盡量避免外出，在室內(nèi)即時關(guān)窗，以減少接觸空氣污染物，從而降低肺結(jié)核發(fā)病風險。

本研究存在一定的局限性。第一，為避免新型冠狀病毒肺炎造成的影響，尤其是疫情防控間接導(dǎo)致其他呼吸系統(tǒng)疾病發(fā)病率的減少［27］，本研究選取2016—2018年的數(shù)據(jù)進行分析，所反映的結(jié)果可能與當前疫情形勢下的實際情況有所偏差。第二，由于本研究數(shù)據(jù)存在較強的共線性，可能會導(dǎo)致模型結(jié)果存在一定偏差，研究組計劃后續(xù)使用更多數(shù)據(jù)集去驗證模型及研究結(jié)果。第三，既往研究顯示，氣象因素與空氣質(zhì)量因素不僅對當天肺結(jié)核發(fā)病情況有顯著影響外，還存在3～5 d的滯后效應(yīng)［23，28］，但本研究因數(shù)據(jù)局限并未能對其進行驗證。第四，影響肺結(jié)核發(fā)病情況的因素還包括生活方式、教育水平、經(jīng)濟狀況等，但由于相關(guān)數(shù)據(jù)缺乏可及性，本研究并未將其納入。因此，研究組計劃后續(xù)引入更新、更全面的數(shù)據(jù)集來提高模型的適用性和穩(wěn)健性。

本研究基于2016—2018年全國數(shù)據(jù)構(gòu)建的GTWR模型很好地展示了中國肺結(jié)核發(fā)病情況的時空分布，并詳細闡述了氣象因素和空氣質(zhì)量因素與肺結(jié)核發(fā)病情況的顯著相關(guān)關(guān)系及其時空特異性。在實際疾病預(yù)防工作中，應(yīng)結(jié)合不同地區(qū)的不同影響因素，針對性的制定疾病預(yù)防措施。例如：高海拔地區(qū)大風天氣注意保暖；濕度較高地區(qū)應(yīng)注意除濕；空氣質(zhì)量較差地區(qū)應(yīng)積極采取措施改善空氣質(zhì)量，并注意避免高空氣污染物濃度時的戶外活動。

作者貢獻：趙明揚負責提出概念、數(shù)據(jù)管理、形式分析、原稿寫作等工作；周乾宇負責方法學、監(jiān)督審查和編輯寫作等工作；王榮榮、王宗熹、何雯倩、張文森、張恒榛、田卓旸、吳柯、王碧瑤負責數(shù)據(jù)管理工作；孫長青負責資金支持、項目管理、監(jiān)督、審查和編輯寫作等工作。

本文無利益沖突。