基于廣義可加模型的時(shí)空數(shù)據(jù)可視化

王艷,柯蓉

（上海對(duì)外經(jīng)貿(mào)大學(xué)統(tǒng)計(jì)與信息學(xué)院,上海 201620）

氣候變化、空氣污染和健康危害是全球極為關(guān)注的研究領(lǐng)域,公共衛(wèi)生科學(xué)家觀察到氣候變化對(duì)各種健康結(jié)果的影響,研究者通過廣義線性模型或廣義可加模型進(jìn)行分析,使用參數(shù)、非參數(shù)樣條估計(jì)大氣參數(shù)和空氣污染的健康效應(yīng).可加模型有助于非線性變量的建模,在制定環(huán)境空氣質(zhì)量的質(zhì)量標(biāo)準(zhǔn)方面,國(guó)家大氣監(jiān)測(cè)和評(píng)估計(jì)劃發(fā)揮了重要作用.

廣義可加模型適用于響應(yīng)變量與解釋變量之間為非線性或非單調(diào)關(guān)系的數(shù)據(jù),允許解釋變量(協(xié)變量)和響應(yīng)之間的關(guān)系由平滑曲線或樣條來描述.廣義可加模型(GAMs)是回歸模型,其中平滑樣條用于代替協(xié)變量的線性系數(shù)[1-2],通常是連續(xù)預(yù)測(cè)值和連續(xù)正態(tài)分布結(jié)果之間的參數(shù)和非參數(shù)形式.GAM 最廣泛的應(yīng)用是在地理環(huán)境領(lǐng)域[3],空氣污染的短期和長(zhǎng)期影響一直被認(rèn)為是不利健康結(jié)果的風(fēng)險(xiǎn),許多使用GAM和不同分布樣條的時(shí)間序列模型,提供了環(huán)境空氣質(zhì)量和健康風(fēng)險(xiǎn)之間的時(shí)間相互關(guān)系,以評(píng)估空氣污染對(duì)呼吸系統(tǒng)死亡率的影響.O3、NO2與德黑蘭的呼吸系統(tǒng)死亡有關(guān),可通過采取可測(cè)量的方法減少環(huán)境空氣污染[4-5].Pearce 等[6]運(yùn)用廣義可加模型框架對(duì)氣象-污染物關(guān)系進(jìn)行了評(píng)估,表明澳大利亞墨爾本的污染物體積分?jǐn)?shù)和氣象參數(shù)之間存在關(guān)系.Wood 等[7-8]討論了分布式滯后模型,并試圖解釋滯后在健康相關(guān)數(shù)據(jù)集中的使用,通常在空氣質(zhì)量相關(guān)模型中,滯后值是基于短期或?qū)?shù)期效應(yīng)的目標(biāo)來選擇的.Li 等[9]采用帶有泊松分布族對(duì)數(shù)連接函數(shù)的GAM模型進(jìn)行分析,得出2014—2015 年合肥市空氣污染與兒童上呼吸道感染(URTI)的關(guān)系.Huo 等[10]運(yùn)用廣義可加模型(GAM)研究分析了1985 年至2015 年中國(guó)東南部長(zhǎng)潭水庫(kù)中氣候變暖和養(yǎng)分含量增加對(duì)沉積色素記錄的藻類群落的影響,表明氣候變暖和人為養(yǎng)分負(fù)荷增加迫使該水庫(kù)浮游植物動(dòng)態(tài)變化,并促使藻類生物量增加.Li 等[11]通過GAMM和GAM克服使用自然樣條預(yù)測(cè)的可加效應(yīng)中的變化和隨機(jī)效應(yīng),研究顯示由于PM10 和NO2增加10 μg/m3,每日心血管死亡率上升0.31 個(gè)百分點(diǎn),空氣污染的空間效應(yīng)也隨著北京市的地理位置而變化.Zhang 等[13]描述了新冠疫情期間空氣污染物體積分?jǐn)?shù)的趨勢(shì)和COVID-19 的發(fā)病率,并應(yīng)用廣義可加性模型來評(píng)估中國(guó)235 個(gè)城市中短期暴露于空氣污染與新冠肺炎每日確診病例之間的關(guān)聯(lián).

本文利用GAM在環(huán)境數(shù)據(jù)方面的各種應(yīng)用的模型優(yōu)勢(shì),特別是在空氣質(zhì)量污染和健康領(lǐng)域,通過R 軟件分別應(yīng)用于兩個(gè)具有時(shí)間、空間變量以及存在交互效應(yīng)的環(huán)境數(shù)據(jù)中,通過可視化分析進(jìn)一步來解釋模型.重點(diǎn)是利用R 軟件實(shí)現(xiàn)GAM在時(shí)空數(shù)據(jù)統(tǒng)計(jì)分析中的應(yīng)用,以圖形可視化的形式清晰地表明環(huán)境因素與預(yù)測(cè)結(jié)果的關(guān)系.

1 數(shù)據(jù)源與研究方法

1.1 數(shù)據(jù)來源

實(shí)例1 數(shù)據(jù)源于環(huán)境計(jì)算網(wǎng)站,來自美國(guó)國(guó)家海洋和大氣局的研究基地莫納羅亞山天文臺(tái)的傳感器數(shù)據(jù),研究選取了從1958 年3 月1 日到2013 年3 月1 日二氧化碳體積分?jǐn)?shù)變化檢測(cè)數(shù)據(jù),55 年間二氧化碳體積分?jǐn)?shù)波動(dòng)范圍從315.71×10-6增加到397.23×10-6.最近一年的檢測(cè)數(shù)據(jù)顯示大氣中二氧化碳體積分?jǐn)?shù)迅速超過了415×10-6,達(dá)到452×10-6的新高.這是人類歷史上地球大氣中二氧化碳體積分?jǐn)?shù)前所未有的峰值.人類已經(jīng)生活在空氣中二氧化碳體積分?jǐn)?shù)超過450×10-6的時(shí)代,“350×10-6安全線”一去不復(fù)返.澳大利亞國(guó)立大學(xué)威爾·斯特芬教授在《自然》雜志發(fā)表文章指出目前全球氣候的變化,或許已經(jīng)不可逆轉(zhuǎn),需要面對(duì)整個(gè)地球不穩(wěn)定的風(fēng)險(xiǎn).本文通過GAM模型來擬合時(shí)間序列的數(shù)據(jù),以嘗試區(qū)分年度內(nèi)和年度間差異,發(fā)現(xiàn)變化趨勢(shì)及規(guī)律.

實(shí)例2 數(shù)據(jù)使用R 軟件中sp 包的meuse 數(shù)據(jù)集,關(guān)于荷蘭默茲河沿岸重金屬土壤污染的地理空間數(shù)據(jù),在默茲河泛濫平原表層土壤中測(cè)量到的4 種重金屬鎘、銅、鉛、鋅,以及一些協(xié)變量,重金屬分布由河流攜帶的污染沉積物控制,主要在靠近河岸和低海拔地區(qū).數(shù)據(jù)集包含具有x 和y 坐標(biāo)的數(shù)據(jù)框、土壤中重金屬含量以及其他空間協(xié)變量.本研究主要通過GAM 建立具有多個(gè)變量相互作用的模型,并使用這些交互作用來對(duì)地理空間數(shù)據(jù)進(jìn)行擬合,對(duì)復(fù)雜的表面使用三維圖形可視化顯示出來.

1.2 研究方法

1.2.1 廣義線性模型 在GLM中因變量的分布可以是非正態(tài)的,也可以是非連續(xù)的,與一般的線性模型相比,響應(yīng)變量可以遵循泊松分布.預(yù)測(cè)變量的線性組合預(yù)測(cè)因變量,以通過連接函數(shù)與因變量“相關(guān)聯(lián)”.在一般的線性模型中因變量的值被假定為遵循漸近分布,預(yù)測(cè)值的線性組合沒有進(jìn)行轉(zhuǎn)換.GLM是一般線性模型的擴(kuò)展,在環(huán)境研究的應(yīng)用中應(yīng)用廣泛.線性模型也是可加模型的一個(gè)特例,通過使用線性參數(shù)化來解釋關(guān)系.但是在某些情況下,線性的假設(shè)非常不準(zhǔn)確,可能導(dǎo)致潛在的不正確推斷.而標(biāo)準(zhǔn)GLM依賴于數(shù)據(jù)服從指數(shù)族分布的假設(shè).因此相應(yīng)的密度可以寫成

式（1）中:θi是未知參數(shù),a,b,c 是依賴于特定分布的固定函數(shù),且滿足E(Yi|xi)=μi=b′(θ),Var(Yi|xi)=在標(biāo)準(zhǔn)GLM中期望值可表示為為線性預(yù)測(cè),然而將預(yù)測(cè)限制為協(xié)變量的線性組合通常是不夠的,因此引入半?yún)?shù)預(yù)測(cè)[2].

1.2.2 廣義可加模型 廣義可加模型在處理空氣污染研究、相關(guān)的復(fù)雜非線性方面特別有效.模型可寫成以下形式[1]

式(2)中:yi是 第i 個(gè)時(shí)間點(diǎn)污染體積分?jǐn)?shù),β0是響應(yīng)變量的總體平均值,s j(xi,j)是第j 個(gè)協(xié)變量第i個(gè)值的平滑函數(shù),n 是協(xié)變量的總數(shù),εi是具有Var(εi)=σ2的殘差,假設(shè)服從正態(tài)分布.平滑函數(shù)通過模型選擇和懲罰回歸樣條的自動(dòng)平滑參數(shù)選擇來表示,優(yōu)化了擬合并最小化模型中的維數(shù).相互作用項(xiàng)也可以建模為薄板回歸樣條或張量積平滑,平滑參數(shù)的選擇通過限制最大似然(REML),置信區(qū)間估計(jì)使用無條件貝葉斯方法.廣義可加模型非常適合處理復(fù)雜的變量相互作用,方法的微小改進(jìn)都可能提高所有氣象數(shù)據(jù)的空間分辨率.

1.2.3 建立模型 對(duì)數(shù)據(jù)1 從1958—2013 年CO2體積分?jǐn)?shù)測(cè)量日期進(jìn)行預(yù)處理后,選取部分?jǐn)?shù)據(jù)集進(jìn)行建模,首先選時(shí)間變量作為平滑項(xiàng),模型表示為

由于數(shù)據(jù)存在年內(nèi)的波動(dòng),再通過增加一個(gè)變量平滑項(xiàng),同時(shí)使用循環(huán)三次樣條回歸的方法,R 中函數(shù)的樣條回歸中參數(shù)節(jié)點(diǎn)k 設(shè)置為12 個(gè)節(jié)點(diǎn)數(shù),模型可表示為

周期性光滑項(xiàng)由基本函數(shù)組成,樣條曲線的端點(diǎn)被約束為相等.

數(shù)據(jù)2 是荷蘭默茲河沿岸重金屬土壤污染的地理空間數(shù)據(jù),一般情況下可加模型由一個(gè)或多個(gè)平滑組成,每個(gè)平滑包含一個(gè)變量.但是在處理空間數(shù)據(jù)時(shí),將相互作用表示為光滑表面,是空間數(shù)據(jù)建模最合適的方法之一,選擇使用復(fù)雜的平面而不是單條平滑線來表示.在線性模型中通過添加兩個(gè)變量相乘的項(xiàng),可能導(dǎo)致結(jié)果高于或低于單獨(dú)兩個(gè)值之和所預(yù)測(cè)的結(jié)果.GAM模型中變量與結(jié)果之間的關(guān)系會(huì)在平滑范圍內(nèi)變化,變量值之間的交互作用也不同.或者將相互作用與其他項(xiàng)(線性或非線性)混合在一起,包括離散的分類項(xiàng)以及交互項(xiàng)和線性項(xiàng)等靈活的結(jié)合方式.運(yùn)用R 中mgcv 軟件包[2]建立具有相互作用的4 個(gè)模型,并使這些交互作用對(duì)地理空間數(shù)據(jù)進(jìn)行擬合,以三維圖形可視化實(shí)現(xiàn).進(jìn)一步分析出平滑變量和分類變量之間的相互作用,并為空間和時(shí)間等不同類型變量進(jìn)行建模.

具體R 實(shí)現(xiàn)步驟:

(1)從sp 包導(dǎo)入土壤重金屬污染數(shù)據(jù)預(yù)處理后,首先利用x 和y 坐標(biāo)位置的交互作用預(yù)測(cè)土壤中鎘濃度建立模型1.可加模型比線性模型的參數(shù)估計(jì)更為復(fù)雜,線性模型中的每個(gè)變量只有一個(gè)系數(shù)或參數(shù).GAM中的平滑實(shí)際是由許多基本函數(shù)構(gòu)成,每個(gè)平滑是多個(gè)基本函數(shù)的總和,每個(gè)基本函數(shù)乘以一個(gè)系數(shù),每個(gè)系數(shù)都是模型中的一個(gè)參數(shù).最后得出模型1 的因變量和自變量之間的單個(gè)非線性關(guān)系具有29 個(gè)基礎(chǔ)參數(shù)和一個(gè)截距項(xiàng).

古代封建專制政權(quán)家天下，皇家也是基因論的忠實(shí)擁躉?；实厶?hào)稱真龍?zhí)熳樱嫔咸煨惺隆；首影⒏缰髯匀豁樌沓烧碌亍吧朴巍?。甭管他們是否天天只知拈弓彈雀、提籠架鳥、不務(wù)正業(yè)，只要生在皇家，便是龍子龍種，將來是一定要子承父業(yè)做皇帝的。打著君權(quán)神授的幌子，蒙蔽百姓，妄想家天下，歷朝歷代的皇帝概莫能外。秦始皇是第一個(gè)從皇帝名號(hào)上動(dòng)此腦筋的，所以他自稱“始皇帝”，妄想子子孫孫二世三世乃至萬世地傳承下去，可惜只傳到二世，秦王朝就壽終正寢了。

(2)通過mgcv 軟件包中的gam 函數(shù)為模型添加其他預(yù)測(cè)變量,建立模型2 預(yù)測(cè)土壤中的鎘濃度,除x,y 曲面外,還包括對(duì)海拔(elev)和離河岸距離(dist)的影響進(jìn)行平滑處理.數(shù)據(jù)中的因子變量Landuse,指定土壤采樣地點(diǎn)的土地利用或覆蓋類型.

(3)通過建立模型3 同時(shí)考慮分類連續(xù)變量的交互作用,加入dist 變量來預(yù)測(cè)銅含量,即每個(gè)土地利用水平Landuse 具有不同的平滑度,擬合具有dist 和Landuse 變量之間的因子平滑關(guān)系的模型,得到一個(gè)整體的交互作用.平滑公式將由R 函數(shù)自動(dòng)調(diào)整,尤其當(dāng)類別非常多或某些類別中的數(shù)據(jù)點(diǎn)很少時(shí),可以很好地控制不是關(guān)注的主要變量類別影響.

(4)除了二維平滑交互和分類連續(xù)交互兩種方式外,張量平滑可以對(duì)不同尺度變量(如空間和時(shí)間)的交互進(jìn)行建模.很多情況下只有一個(gè)平滑參數(shù)沒有意義,數(shù)據(jù)中具有x 和y 以及高度的變量,以米為單位.理論上希望得到變化的尺度在x 和y 的方向上相似,但是隨著海拔的升高,每米的變化可能會(huì)大不相同,其中小的差異可能意味著非常不同的環(huán)境.同樣對(duì)距離和om 變量(以g/kg 表示的有機(jī)物)也沒有意義,它們具有不同的度量單位.張量平滑形式適合具有不同尺度或單位的變量交互,具有兩個(gè)平滑參數(shù),為每個(gè)平滑指定不同數(shù)量的基礎(chǔ)函數(shù)k 值.另一個(gè)優(yōu)勢(shì)是可將相互作用與單個(gè)變量效應(yīng)分開.因此建立模型5 預(yù)測(cè)鎘污染隨位置和海拔的變化,使用平滑項(xiàng)對(duì)每個(gè)單變量分別建模,然后為每個(gè)交互項(xiàng)使用張量平滑.同時(shí)單變量平滑是可加的,交互作用是在這基礎(chǔ)上的附加作用,這種分離效果使得復(fù)雜的模型更容易理解.

2 結(jié)果與分析

2.1 時(shí)間序列GAM 模型

首先將日期轉(zhuǎn)換為連續(xù)時(shí)間變量,選取子數(shù)據(jù)集(2000 年到2010 年)對(duì)年際趨勢(shì)變化進(jìn)行描述性統(tǒng)計(jì),如圖1 所示,從2000 年1 月1 日到2010 年12 月1 日二氧化碳的體積分?jǐn)?shù)含量變化呈小范圍周期性波動(dòng)及整體向上趨勢(shì),接著采用可加模型來擬合數(shù)據(jù),并對(duì)擬合情況進(jìn)行診斷.如圖2 所示.

圖1 年際趨勢(shì)Fig.1 Inter-annual trend

圖2 中的時(shí)間變量平滑項(xiàng)實(shí)際上變?yōu)橐话愕木€性項(xiàng),是由于懲罰線性樣條曲線的作用,其中模型的有效自由度(edf)為1,自由度大于1 則表示模型的非線性明顯.mgcv 中的默認(rèn)設(shè)置是薄板回歸樣條,平滑項(xiàng)是一些基礎(chǔ)函數(shù)的總和.通過模型的診斷圖,會(huì)發(fā)現(xiàn)殘差圖出現(xiàn)了上升和下降的趨勢(shì),存在某種依賴性結(jié)構(gòu)，可能與年內(nèi)波動(dòng)有關(guān),再通過改變平滑的參數(shù)方式引入循環(huán)三次樣條回歸方法,得到新的擬合平滑項(xiàng)，如圖3 所示.

圖2 單一時(shí)間平滑及診斷圖Fig.2 Single time smoothing and diagnostic graph

圖3 兩個(gè)平滑項(xiàng)及診斷圖Fig.3 Two smoothing items and a diagnostic diagram

通過設(shè)置參數(shù)bs=“cc”,k=12,因?yàn)槿位貧w樣條具有一定數(shù)量的結(jié),一年12 個(gè)月則將k 設(shè)為12.從兩個(gè)平滑項(xiàng)擬合來看,可以看到模型月份變量平滑項(xiàng)很好地解釋了圖2 中殘差圖的上升和下降的波動(dòng)趨勢(shì),可以看出s(month)平滑項(xiàng)分解出時(shí)間序列組成部分的波動(dòng)效應(yīng).診斷結(jié)果顯示模型的偏差解釋從89.5%增加到接近100%,s(month)和s(time)的edf 分別為9.367、8.847,非線性擬合效果提升.圖4 顯示了季節(jié)性因素與整個(gè)數(shù)據(jù)長(zhǎng)期趨勢(shì)相對(duì)應(yīng).對(duì)于平滑項(xiàng)沒有推論出的單個(gè)系數(shù)(正負(fù)值、效應(yīng)大小等),因此可從圖形中解釋平滑項(xiàng)的效果或根據(jù)二氧化碳預(yù)測(cè)值推斷變化趨勢(shì),在年內(nèi)呈周期性波動(dòng)、年際間整體上升.

圖4 長(zhǎng)期趨勢(shì)及預(yù)測(cè)Fig.4 Long-term trends and projections

2.2 空間數(shù)據(jù)GAM 模型

模型1 使用x 和y 位置坐標(biāo)的相互作用預(yù)測(cè)土壤中鎘濃度,將x,y 組合的效果合并在一個(gè)平滑項(xiàng),模型的偏差解釋為66.7%.而在線性模型中,變量及其組合項(xiàng)是分開的.如圖5 所示,輪廓圖體現(xiàn)交互作用,軸表示預(yù)測(cè)變量x 和y 的值.內(nèi)部是預(yù)測(cè)值的地形圖,等高線表示相等的預(yù)測(cè)值.虛線表示預(yù)測(cè)的不確定性,如果預(yù)測(cè)是一個(gè)較高或較低的標(biāo)準(zhǔn)誤差,輪廓線將產(chǎn)生移動(dòng).將plot 函數(shù)中的scheme 參數(shù)設(shè)置為1 獲得三維透視圖,設(shè)置為2 會(huì)生成一個(gè)熱力圖,淺色代表鎘濃度較大值,深色代表鎘濃度較小值.

圖5 模型1 可視化Fig.5 Model 1 visualization

圖6 模型1 預(yù)測(cè)Fig.6 Model 1 prediction

圖7 是改變旋轉(zhuǎn)角度、縮放后的透視圖,通過se 參數(shù)顯示預(yù)測(cè)的置信區(qū)間,即利用標(biāo)準(zhǔn)誤差與平均預(yù)測(cè)值的差值繪制高低預(yù)測(cè)曲面.從圖可直觀的看出隨著位置y 的增大,曲面的高度在不斷增加,鎘濃度在增大;區(qū)域的中心位置鎘濃度最低.

模型2 預(yù)測(cè)土壤中鎘含量,除了x,y 交互變量還加入海拔高度(elev)和距河流距離(dist)變量進(jìn)行平滑處理,將模型2 的交互作用項(xiàng)分別繪制為圖8 中的等高線圖、透視圖、熱力圖,以及圖7 第三幅的預(yù)測(cè)置信曲面圖,三種類別可視化圖表明隨著海拔高度和距離增加鎘含量在不斷減小,模型2 的偏差解釋達(dá)到了84.4%,較單一平滑項(xiàng)模型1 有所提高.

圖7 置信曲面圖Fig.7 Confidence surface

圖8 模型2 可視化Fig.8 Model 2 visualization

模型3 通過土地使用類別變量擬合具有單獨(dú)平滑項(xiàng)的模型,只根據(jù)距離(dist)預(yù)測(cè)銅含量,模型的偏差解釋為58.3%.在模型中為分類變量的每個(gè)值擬合不同的平滑,另一種類別連續(xù)的交互為”因子平滑”,將bs 參數(shù)設(shè)置為fs,指定兩個(gè)變量作為平滑的一部分,觀察不同連續(xù)類別交互類型之間的差異.可以看出當(dāng)使用bs 參數(shù)在GAM上調(diào)用因子平滑擬合的函數(shù)時(shí),默認(rèn)情況下將繪制一張具有多個(gè)平滑度的圖,使用vis.gam 函數(shù)可視化因子平滑,產(chǎn)生類似階梯的透視圖(圖9 所示),顯示了不同土地利用方式的污染分布.

圖9 分類交互透視圖Fig.9 Categorizes interactive perspectives

模型4 使用張量平滑對(duì)不同尺度變量的交互進(jìn)行建模,預(yù)測(cè)鎘污染隨位置和海拔的變化.即使變量尺度不同也可將x,y,elev 各自的作用與相互作用分離,x 和y 在相同尺度上交互、elev 是單獨(dú)的平滑項(xiàng)、三個(gè)不同尺度的相互作用是一個(gè)單獨(dú)項(xiàng).使用具有平滑和張量交互作用的擬合模型,以分離出變量的獨(dú)立影響,并評(píng)估這些平滑項(xiàng)的重要性,由模型結(jié)果得出三項(xiàng)都具有顯著性,模型的解釋達(dá)到84.7%,可視化結(jié)果如圖10 所示.通過對(duì)4 個(gè)模型的比較可知,通過增加多變量平滑項(xiàng)的GAM模型2 以及使用張量平滑項(xiàng)的模型4 擬合效果最好,預(yù)測(cè)模型的偏差解釋都達(dá)到了80%以上.

圖10 張量平滑F(xiàn)ig.10 Tensor smooths

3 小結(jié)

廣義可加模型一直被用作衡量空氣污染短期和長(zhǎng)期影響的重要統(tǒng)計(jì)工具.為了分析短期效應(yīng),帶懲罰樣條的GAM被認(rèn)為是研究環(huán)境、氣候和健康聯(lián)系的最佳方法.R 中的mgcv 包為各種數(shù)據(jù)集提供了GAM的實(shí)現(xiàn),包括各種樣條.在擬合GAM模型時(shí),自由度也是一個(gè)重要的考慮因素,自由度取決于平滑參數(shù),如果平滑參數(shù)都設(shè)置為零,那么模型的自由度就是模型系數(shù)的維數(shù).在預(yù)測(cè)變量中存在很大噪聲的情況下,GAM提供了一種靈活的解決方法,在預(yù)測(cè)值和自變量之間存在非線性關(guān)系的情況下,也表現(xiàn)出了最佳擬合.將GAM應(yīng)用于數(shù)據(jù)集需要注意小的改變可進(jìn)一步改進(jìn)模型,使用殘差圖可以容易地識(shí)別出可能的異常值,并且為了更好地?cái)M合模型,可以進(jìn)一步消除這些異常值.由于大多數(shù)環(huán)境數(shù)據(jù)是非正態(tài)的,GAM及可視化交互提供了比傳統(tǒng)線性模型更有效的分析方法.