王艷,柯蓉

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

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

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

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

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

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

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

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

1.2 研究方法

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

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

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

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

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

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

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

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

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

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

古代封建專制政權(quán)家天下，皇家也是基因論的忠實擁躉?；实厶柗Q真龍?zhí)熳?，代替上天行事。皇子阿哥之流自然順理成章地“善游”。甭管他們是否天天只知拈弓彈雀、提籠架鳥、不務(wù)正業(yè)，只要生在皇家，便是龍子龍種，將來是一定要子承父業(yè)做皇帝的。打著君權(quán)神授的幌子，蒙蔽百姓，妄想家天下，歷朝歷代的皇帝概莫能外。秦始皇是第一個從皇帝名號上動此腦筋的，所以他自稱“始皇帝”，妄想子子孫孫二世三世乃至萬世地傳承下去，可惜只傳到二世，秦王朝就壽終正寢了。

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

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

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

2 結(jié)果與分析

2.1 時間序列GAM 模型

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

圖1 年際趨勢Fig.1 Inter-annual trend

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

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

圖3 兩個平滑項及診斷圖Fig.3 Two smoothing items and a diagnostic diagram

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

圖4 長期趨勢及預(yù)測Fig.4 Long-term trends and projections

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

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

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

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

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

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

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

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

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

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

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

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

3 小結(jié)

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