基于广义可加模型的时空数据可视化

王艳,柯蓉

（上海对外经贸大学统计与信息学院,上海 201620）

气候变化、空气污染和健康危害是全球极为关注的研究领域,公共卫生科学家观察到气候变化对各种健康结果的影响,研究者通过广义线性模型或广义可加模型进行分析,使用参数、非参数样条估计大气参数和空气污染的健康效应.可加模型有助于非线性变量的建模,在制定环境空气质量的质量标准方面,国家大气监测和评估计划发挥了重要作用.

广义可加模型适用于响应变量与解释变量之间为非线性或非单调关系的数据,允许解释变量(协变量)和响应之间的关系由平滑曲线或样条来描述.广义可加模型(GAMs)是回归模型,其中平滑样条用于代替协变量的线性系数[1-2],通常是连续预测值和连续正态分布结果之间的参数和非参数形式.GAM 最广泛的应用是在地理环境领域[3],空气污染的短期和长期影响一直被认为是不利健康结果的风险,许多使用GAM和不同分布样条的时间序列模型,提供了环境空气质量和健康风险之间的时间相互关系,以评估空气污染对呼吸系统死亡率的影响.O3、NO2与德黑兰的呼吸系统死亡有关,可通过采取可测量的方法减少环境空气污染[4-5].Pearce 等[6]运用广义可加模型框架对气象-污染物关系进行了评估,表明澳大利亚墨尔本的污染物体积分数和气象参数之间存在关系.Wood 等[7-8]讨论了分布式滞后模型,并试图解释滞后在健康相关数据集中的使用,通常在空气质量相关模型中,滞后值是基于短期或对数期效应的目标来选择的.Li 等[9]采用带有泊松分布族对数连接函数的GAM模型进行分析,得出2014—2015 年合肥市空气污染与儿童上呼吸道感染(URTI)的关系.Huo 等[10]运用广义可加模型(GAM)研究分析了1985 年至2015 年中国东南部长潭水库中气候变暖和养分含量增加对沉积色素记录的藻类群落的影响,表明气候变暖和人为养分负荷增加迫使该水库浮游植物动态变化,并促使藻类生物量增加.Li 等[11]通过GAMM和GAM克服使用自然样条预测的可加效应中的变化和随机效应,研究显示由于PM10 和NO2增加10 μg/m3,每日心血管死亡率上升0.31 个百分点,空气污染的空间效应也随着北京市的地理位置而变化.Zhang 等[13]描述了新冠疫情期间空气污染物体积分数的趋势和COVID-19 的发病率,并应用广义可加性模型来评估中国235 个城市中短期暴露于空气污染与新冠肺炎每日确诊病例之间的关联.

本文利用GAM在环境数据方面的各种应用的模型优势,特别是在空气质量污染和健康领域,通过R 软件分别应用于两个具有时间、空间变量以及存在交互效应的环境数据中,通过可视化分析进一步来解释模型.重点是利用R 软件实现GAM在时空数据统计分析中的应用,以图形可视化的形式清晰地表明环境因素与预测结果的关系.

1 数据源与研究方法

1.1 数据来源

实例1 数据源于环境计算网站,来自美国国家海洋和大气局的研究基地莫纳罗亚山天文台的传感器数据,研究选取了从1958 年3 月1 日到2013 年3 月1 日二氧化碳体积分数变化检测数据,55 年间二氧化碳体积分数波动范围从315.71×10-6增加到397.23×10-6.最近一年的检测数据显示大气中二氧化碳体积分数迅速超过了415×10-6,达到452×10-6的新高.这是人类历史上地球大气中二氧化碳体积分数前所未有的峰值.人类已经生活在空气中二氧化碳体积分数超过450×10-6的时代,“350×10-6安全线”一去不复返.澳大利亚国立大学威尔·斯特芬教授在《自然》杂志发表文章指出目前全球气候的变化,或许已经不可逆转,需要面对整个地球不稳定的风险.本文通过GAM模型来拟合时间序列的数据,以尝试区分年度内和年度间差异,发现变化趋势及规律.

实例2 数据使用R 软件中sp 包的meuse 数据集,关于荷兰默兹河沿岸重金属土壤污染的地理空间数据,在默兹河泛滥平原表层土壤中测量到的4 种重金属镉、铜、铅、锌,以及一些协变量,重金属分布由河流携带的污染沉积物控制,主要在靠近河岸和低海拔地区.数据集包含具有x 和y 坐标的数据框、土壤中重金属含量以及其他空间协变量.本研究主要通过GAM 建立具有多个变量相互作用的模型,并使用这些交互作用来对地理空间数据进行拟合,对复杂的表面使用三维图形可视化显示出来.

1.2 研究方法

1.2.1 广义线性模型 在GLM中因变量的分布可以是非正态的,也可以是非连续的,与一般的线性模型相比,响应变量可以遵循泊松分布.预测变量的线性组合预测因变量,以通过连接函数与因变量“相关联”.在一般的线性模型中因变量的值被假定为遵循渐近分布,预测值的线性组合没有进行转换.GLM是一般线性模型的扩展,在环境研究的应用中应用广泛.线性模型也是可加模型的一个特例,通过使用线性参数化来解释关系.但是在某些情况下,线性的假设非常不准确,可能导致潜在的不正确推断.而标准GLM依赖于数据服从指数族分布的假设.因此相应的密度可以写成

式（1）中:θi是未知参数,a,b,c 是依赖于特定分布的固定函数,且满足E(Yi|xi)=µi=b′(θ),Var(Yi|xi)=在标准GLM中期望值可表示为为线性预测,然而将预测限制为协变量的线性组合通常是不够的,因此引入半参数预测[2].

1.2.2 广义可加模型 广义可加模型在处理空气污染研究、相关的复杂非线性方面特别有效.模型可写成以下形式[1]

式(2)中:yi是 第i 个时间点污染体积分数,β0是响应变量的总体平均值,s j(xi,j)是第j 个协变量第i个值的平滑函数,n 是协变量的总数,εi是具有Var(εi)=σ2的残差,假设服从正态分布.平滑函数通过模型选择和惩罚回归样条的自动平滑参数选择来表示,优化了拟合并最小化模型中的维数.相互作用项也可以建模为薄板回归样条或张量积平滑,平滑参数的选择通过限制最大似然(REML),置信区间估计使用无条件贝叶斯方法.广义可加模型非常适合处理复杂的变量相互作用,方法的微小改进都可能提高所有气象数据的空间分辨率.

1.2.3 建立模型 对数据1 从1958—2013 年CO2体积分数测量日期进行预处理后,选取部分数据集进行建模,首先选时间变量作为平滑项,模型表示为

由于数据存在年内的波动,再通过增加一个变量平滑项,同时使用循环三次样条回归的方法,R 中函数的样条回归中参数节点k 设置为12 个节点数,模型可表示为

周期性光滑项由基本函数组成,样条曲线的端点被约束为相等.

数据2 是荷兰默兹河沿岸重金属土壤污染的地理空间数据,一般情况下可加模型由一个或多个平滑组成,每个平滑包含一个变量.但是在处理空间数据时,将相互作用表示为光滑表面,是空间数据建模最合适的方法之一,选择使用复杂的平面而不是单条平滑线来表示.在线性模型中通过添加两个变量相乘的项,可能导致结果高于或低于单独两个值之和所预测的结果.GAM模型中变量与结果之间的关系会在平滑范围内变化,变量值之间的交互作用也不同.或者将相互作用与其他项(线性或非线性)混合在一起,包括离散的分类项以及交互项和线性项等灵活的结合方式.运用R 中mgcv 软件包[2]建立具有相互作用的4 个模型,并使这些交互作用对地理空间数据进行拟合,以三维图形可视化实现.进一步分析出平滑变量和分类变量之间的相互作用,并为空间和时间等不同类型变量进行建模.

具体R 实现步骤:

(1)从sp 包导入土壤重金属污染数据预处理后,首先利用x 和y 坐标位置的交互作用预测土壤中镉浓度建立模型1.可加模型比线性模型的参数估计更为复杂,线性模型中的每个变量只有一个系数或参数.GAM中的平滑实际是由许多基本函数构成,每个平滑是多个基本函数的总和,每个基本函数乘以一个系数,每个系数都是模型中的一个参数.最后得出模型1 的因变量和自变量之间的单个非线性关系具有29 个基础参数和一个截距项.

古代封建专制政权家天下，皇家也是基因论的忠实拥趸。皇帝号称真龙天子，代替上天行事。皇子阿哥之流自然顺理成章地“善游”。甭管他们是否天天只知拈弓弹雀、提笼架鸟、不务正业，只要生在皇家，便是龙子龙种，将来是一定要子承父业做皇帝的。打着君权神授的幌子，蒙蔽百姓，妄想家天下，历朝历代的皇帝概莫能外。秦始皇是第一个从皇帝名号上动此脑筋的，所以他自称“始皇帝”，妄想子子孙孙二世三世乃至万世地传承下去，可惜只传到二世，秦王朝就寿终正寝了。

(2)通过mgcv 软件包中的gam 函数为模型添加其他预测变量,建立模型2 预测土壤中的镉浓度,除x,y 曲面外,还包括对海拔(elev)和离河岸距离(dist)的影响进行平滑处理.数据中的因子变量Landuse,指定土壤采样地点的土地利用或覆盖类型.

(3)通过建立模型3 同时考虑分类连续变量的交互作用,加入dist 变量来预测铜含量,即每个土地利用水平Landuse 具有不同的平滑度,拟合具有dist 和Landuse 变量之间的因子平滑关系的模型,得到一个整体的交互作用.平滑公式将由R 函数自动调整,尤其当类别非常多或某些类别中的数据点很少时,可以很好地控制不是关注的主要变量类别影响.

(4)除了二维平滑交互和分类连续交互两种方式外,张量平滑可以对不同尺度变量(如空间和时间)的交互进行建模.很多情况下只有一个平滑参数没有意义,数据中具有x 和y 以及高度的变量,以米为单位.理论上希望得到变化的尺度在x 和y 的方向上相似,但是随着海拔的升高,每米的变化可能会大不相同,其中小的差异可能意味着非常不同的环境.同样对距离和om 变量(以g/kg 表示的有机物)也没有意义,它们具有不同的度量单位.张量平滑形式适合具有不同尺度或单位的变量交互,具有两个平滑参数,为每个平滑指定不同数量的基础函数k 值.另一个优势是可将相互作用与单个变量效应分开.因此建立模型5 预测镉污染随位置和海拔的变化,使用平滑项对每个单变量分别建模,然后为每个交互项使用张量平滑.同时单变量平滑是可加的,交互作用是在这基础上的附加作用,这种分离效果使得复杂的模型更容易理解.

2 结果与分析

2.1 时间序列GAM 模型

首先将日期转换为连续时间变量,选取子数据集(2000 年到2010 年)对年际趋势变化进行描述性统计,如图1 所示,从2000 年1 月1 日到2010 年12 月1 日二氧化碳的体积分数含量变化呈小范围周期性波动及整体向上趋势,接着采用可加模型来拟合数据,并对拟合情况进行诊断.如图2 所示.

图1 年际趋势Fig.1 Inter-annual trend

图2 中的时间变量平滑项实际上变为一般的线性项,是由于惩罚线性样条曲线的作用,其中模型的有效自由度(edf)为1,自由度大于1 则表示模型的非线性明显.mgcv 中的默认设置是薄板回归样条,平滑项是一些基础函数的总和.通过模型的诊断图,会发现残差图出现了上升和下降的趋势,存在某种依赖性结构，可能与年内波动有关,再通过改变平滑的参数方式引入循环三次样条回归方法,得到新的拟合平滑项，如图3 所示.

图2 单一时间平滑及诊断图Fig.2 Single time smoothing and diagnostic graph

图3 两个平滑项及诊断图Fig.3 Two smoothing items and a diagnostic diagram

通过设置参数bs=“cc”,k=12,因为三次回归样条具有一定数量的结,一年12 个月则将k 设为12.从两个平滑项拟合来看,可以看到模型月份变量平滑项很好地解释了图2 中残差图的上升和下降的波动趋势,可以看出s(month)平滑项分解出时间序列组成部分的波动效应.诊断结果显示模型的偏差解释从89.5%增加到接近100%,s(month)和s(time)的edf 分别为9.367、8.847,非线性拟合效果提升.图4 显示了季节性因素与整个数据长期趋势相对应.对于平滑项没有推论出的单个系数(正负值、效应大小等),因此可从图形中解释平滑项的效果或根据二氧化碳预测值推断变化趋势,在年内呈周期性波动、年际间整体上升.

图4 长期趋势及预测Fig.4 Long-term trends and projections

2.2 空间数据GAM 模型

模型1 使用x 和y 位置坐标的相互作用预测土壤中镉浓度,将x,y 组合的效果合并在一个平滑项,模型的偏差解释为66.7%.而在线性模型中,变量及其组合项是分开的.如图5 所示,轮廓图体现交互作用,轴表示预测变量x 和y 的值.内部是预测值的地形图,等高线表示相等的预测值.虚线表示预测的不确定性,如果预测是一个较高或较低的标准误差,轮廓线将产生移动.将plot 函数中的scheme 参数设置为1 获得三维透视图,设置为2 会生成一个热力图,浅色代表镉浓度较大值,深色代表镉浓度较小值.

图5 模型1 可视化Fig.5 Model 1 visualization

图6 模型1 预测Fig.6 Model 1 prediction

图7 是改变旋转角度、缩放后的透视图,通过se 参数显示预测的置信区间,即利用标准误差与平均预测值的差值绘制高低预测曲面.从图可直观的看出随着位置y 的增大,曲面的高度在不断增加,镉浓度在增大;区域的中心位置镉浓度最低.

模型2 预测土壤中镉含量,除了x,y 交互变量还加入海拔高度(elev)和距河流距离(dist)变量进行平滑处理,将模型2 的交互作用项分别绘制为图8 中的等高线图、透视图、热力图,以及图7 第三幅的预测置信曲面图,三种类别可视化图表明随着海拔高度和距离增加镉含量在不断减小,模型2 的偏差解释达到了84.4%,较单一平滑项模型1 有所提高.

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

图8 模型2 可视化Fig.8 Model 2 visualization

模型3 通过土地使用类别变量拟合具有单独平滑项的模型,只根据距离(dist)预测铜含量,模型的偏差解释为58.3%.在模型中为分类变量的每个值拟合不同的平滑,另一种类别连续的交互为”因子平滑”,将bs 参数设置为fs,指定两个变量作为平滑的一部分,观察不同连续类别交互类型之间的差异.可以看出当使用bs 参数在GAM上调用因子平滑拟合的函数时,默认情况下将绘制一张具有多个平滑度的图,使用vis.gam 函数可视化因子平滑,产生类似阶梯的透视图(图9 所示),显示了不同土地利用方式的污染分布.

图9 分类交互透视图Fig.9 Categorizes interactive perspectives

模型4 使用张量平滑对不同尺度变量的交互进行建模,预测镉污染随位置和海拔的变化.即使变量尺度不同也可将x,y,elev 各自的作用与相互作用分离,x 和y 在相同尺度上交互、elev 是单独的平滑项、三个不同尺度的相互作用是一个单独项.使用具有平滑和张量交互作用的拟合模型,以分离出变量的独立影响,并评估这些平滑项的重要性,由模型结果得出三项都具有显著性,模型的解释达到84.7%,可视化结果如图10 所示.通过对4 个模型的比较可知,通过增加多变量平滑项的GAM模型2 以及使用张量平滑项的模型4 拟合效果最好,预测模型的偏差解释都达到了80%以上.

图10 张量平滑Fig.10 Tensor smooths

3 小结

广义可加模型一直被用作衡量空气污染短期和长期影响的重要统计工具.为了分析短期效应,带惩罚样条的GAM被认为是研究环境、气候和健康联系的最佳方法.R 中的mgcv 包为各种数据集提供了GAM的实现,包括各种样条.在拟合GAM模型时,自由度也是一个重要的考虑因素,自由度取决于平滑参数,如果平滑参数都设置为零,那么模型的自由度就是模型系数的维数.在预测变量中存在很大噪声的情况下,GAM提供了一种灵活的解决方法,在预测值和自变量之间存在非线性关系的情况下,也表现出了最佳拟合.将GAM应用于数据集需要注意小的改变可进一步改进模型,使用残差图可以容易地识别出可能的异常值,并且为了更好地拟合模型,可以进一步消除这些异常值.由于大多数环境数据是非正态的,GAM及可视化交互提供了比传统线性模型更有效的分析方法.