基于数据挖掘技术的PM2.5污染与居民死亡人数的暴露-反应关系

张书源, 程全国, 邢红彬

(沈阳大学 a. 环境学院, b. 大学生创新创业指导中心, 辽宁 沈阳 110044)

大气污染是影响人类健康的主要环境危害因素之一。在公认的大气污染物中,可吸入颗粒物PM10、PM2.5与人群健康效应终点的流行病学联系最为密切[1-2]。近年来,对可吸入颗粒物的定量健康危害评价已成为国际社会关注的热点之一[3]。

大量的研究结果表明,空气可吸入颗粒物污染与居民每日死亡总人数、呼吸系统与心血管等疾病的超额死亡数显著相关[4-5],即使在空气颗粒物污染浓度低于标准的情况下,随着颗粒物浓度的升高,呼吸系统和心脑血管疾病的发病率和死亡率亦呈现上升趋势[6-8]。

数据挖掘是指从大量的数据中通过算法与技术搜索隐藏于其中的信息的过程。数据挖掘中经常用到的算法和技术有统计分析、决策树、贝叶斯网络、粗糙集、人工神经网络等。随着生态文明建设和生态环境保护进入了数据驱动的新时代,数据挖掘技术在生态环境保护领域中应用的广度和深度都不断加强。

目前国内外研究人员运用流行病学和毒理学的研究方法探索大气污染对死亡率的影响,并建立暴露反应关系模型[9-10],对于大气污染与各死因类型的相关分析研究相对较少。本文以辽宁省本溪市2019年的死因数据、大气污染和气象监测数据作为数据来源,通过数据挖掘技术分析PM2.5污染与各死因类型的显著性与相关性,并以此为基础建立PM2.5污染与每日死亡人数的暴露-反应关系模型,为促进数据挖掘技术在生态环境保护领域中的应用、制定有效的环境法规、保护易感人群提供科学依据。

1 数据与方法

1.1 数据收集与处理

分析所用数据包括数2019年1月1日—12月31日本溪市的死因数据、大气污染和气象数据。死因数据来源于本溪市疾控中心,死因按照国际疾病分类ICD-10进行分类;大气污染数据来自中国空气质量在线监测平台(https:∥www.aqistudy.cn/historydata/)的监测数据,包括PM2.5、PM10、SO2、CO、NO2、O3的每日平均质量浓度(μg·m-3);气象数据来自国家气象科学数据中心(http:∥data.cma.cn/)本溪站的监测数据,包括每日平均温度(℃)、相对湿度(%)、平均气压(kPa)、平均风速(m·s-1)。

数据处理按照数据挖掘技术中统计分析的一般步骤,即数据清洗、缺失值填充、z-score数据标准化方法的顺序对原始数据进行处理。由于死因数据中的ICD-10死因分类类别繁杂,为了避免特征维度的扩张与方便统计分析软件识别各类死因,根据ICD-10分类标准对本溪市2019年1月1日—12月31日的死因进行分类,并将死因进行数值化编码,如表1所示。

表1 死因分类数值化编码Table 1 Numerical coding for the classification of causes of death

1.2 统计方法

大气污染流行病学研究中常用时间序列的广义相加泊松(Poisson)回归模型评价大气污染短期暴露对人群的健康影响,由于对于总人口来说,每日死亡人数是小概率事件,作为一种时间序列模型,其实际分布近似泊松分布。时间序列的广义相加泊松回归模型是对传统广义对数线性模型的进一步拓展,模型中除拟合普通的线性项(如大气污染物浓度)外,还可将一些与因变量之间存在复杂非线性关系的变量(如时间序列资料中的长期趋势、季节和其他一些与时间长期变异有关的混杂因子)以不同函数加和的形式拟合模型[11]。广义相加Poisson回归模型的一般形式为

lnE(Y)=βX。

式中:E(Y)为日死亡人数预测值;X为大气污染物浓度及气象因素数据等一些其他变量构成的矩阵,包括截距列;β为模型自变量的系数向量。

通过PM2.5污染对各死因的差异显著性分析,判断PM2.5污染对不同类型的死因是否具有显著性影响,在此基础上将死亡人数的自然对数与大气污染物指标、气象指标进行相关性分析,并考虑滞后日效应,选择自然对数与PM2.5和其他指标存在显著关联的变量进行线性回归分析,从而建立PM2.5污染与居民每日死亡人数的暴露-反应关系。

2 研究结果

2.1 描述性统计分析

本溪市2019年1月1日—12月31日的PM2.5日均质量浓度在各死因分类(编码1～20)中的分布情况如表2所示。

表2 PM2.5日均质量浓度在各死因分类中的分布情况Table 2 Distribution of PM2.5 daily average mass concentration in the classification of death causes

从表2中可以看出,2019年本溪市居民死因中的肿瘤,内分泌、营养和代谢疾病,循环系统疾病,呼吸系统疾病及消化系统疾病为主要死因,因循环系统疾病导致死亡的人数最多,因皮肤和皮下组织疾病导致死亡的人数最少,PM2.5质量浓度在死因分类各组中的均值在35～50 μg·m-3之间波动。

2.2 差异显著性分析

显著性P>0.05表示差异性不显著;0.01≤P≤0.05表示差异性显著;P<0.01表示差异性极显著。本文在研究本溪市PM2.5污染对每日死亡人数的影响时,首先进行基于死因分组的大气污染物指标均值差异显著性分析,判断PM2.5等大气污染物对不同死因是否具有差异显著性。

利用SPSS软件,将数值编码后的死因分类(编码1～20)作为因子,将大气污染物数据中PM2.5、PM10、SO2、CO、NO2、O3的每日平均质量浓度作为因变量进行差异显著性分析,结果如表3所示。从表3可以看出,PM2.5、PM10、SO2、CO、NO2污染的显著性大于0.05,O3污染的显著性小于0.05且接近0.01。分析结果说明只有O3污染对不同类型的死因具有显著影响,PM2.5及其他大气污染物对不同类型的死因影响不显著,可以认为PM2.5对不同类型的死因影响相同。因此本文在分析本溪市PM2.5污染对每日死亡人数的影响时,应对全部死因类型的死亡人数进行分析。

表3 大气污染物在各死因组间的均值差异显著性分析Table 3 Significant analysis of the difference in the mean value of air pollutants among the death cause groups

差异显著性分析除了对研究的总体均值的差异进行显著性检验以外,还需要对2个独立样本所属总体的总体方差的差异进行显著性检验,统计学上称为方差齐性检验。

利用SPSS软件,基于死因分组的大气污染物方差齐性检验结果如表4所示,从表4中可以看出,除O3外,其他大气污染物的显著性均大于0.05,可以认为除O3外的其他大气污染物在各组死因中的方差相等,O3在各组死因中的方差差别较大,说明PM2.5在各组死因中数据波动一致,进一步证明PM2.5对不同类型的死因影响不显著。

表4 大气污染物在各死因组间的方差齐性检验Table 4 Test for the homogeneity of variance of air pollutants among the death cause groups

2.3 相关性分析

相关性分析是指对2个或多个具备相关性的变量元素进行分析,从而衡量2个变量因素的相关密切程度。相关系数是研究变量之间线性相关程度的量,较为常用的是皮尔逊相关系数,系数为正值表示变量间为正相关关系,负值则表示为负相关关系,系数绝对值越大,表示关系紧密程度越高,本文相关性分析结果中的相关系数均为皮尔逊相关系数。

为探究大气污染物指标、气象指标与每日全因死亡人数(D)的相关性,利用SPSS软件对本溪市2019年1月1日—12月31日的大气污染物指标中的PM2.5、PM10、SO2、CO、NO2、O3每日平均质量浓度以及气象指标中的平均温度、相对湿度、平均气压、平均风速,与每日全因死亡人数作为因子,进行双变量相关性分析,并考虑滞后日效应。滞留日是在时间序列资料分析中常用的概念,其含义是用今日的健康指标与前期若干日的气象或大气污染质量浓度值进行分析,从而研究前几日的气象情况或大气污染对以后健康问题的影响。本文使用滞后0(即当天)～7d的大气污染物数据、气象数据与每日全因死亡人数进行相关性分析,分析结果如表5所示。

表5 滞后0～7 d日全因死亡人数(D)与各指标相关性分析Table 5 Analysis of the correlation between the number of daily all-cause deaths (D) and various indicators after a lag of 0～7 days

表5中,相关系数右上角的星号(*)代表显著性水平,即假设变量间具有线性相关的可能性,有星号则表示可能有关系,否则表示无关系。通过表5可知,滞后0、3、6、7 d,日全因死亡人数与PM2.5的日平均质量浓度显著性均小于0.01,表明假设结果极显著,具有相关性,日全因死亡人数与PM2.5污染具有明显的统计学意义。同时,PM10、SO2、CO、NO2日平均质量浓度及平均温度、平均气压也与日全因死亡人数具有统计学意义。

为进一步探究日全因死亡人数与PM2.5污染的内在关联,需要继续进行偏相关分析。偏相关分析是指当2个变量同时与第3个变量相关时,将第3个变量的影响剔除,只分析另外2个变量之间相关程度的过程。通过偏相关分析,控制气象因素(平均温度、相对湿度、平均气压、平均风速)的影响,只分析日全因死亡人数与大气污染物指标的内在相关性,利用SPSS软件分别对当天、滞后3、6、7 d的数据进行偏相关分析,分析结果如表6。通过表6可以看出,PM2.5污染的相关性由大到小为:当天、滞后6 d、滞后3 d、滞后7 d,显著性由大到小为:滞后7 d、滞后3 d、滞后6 d、当天。当天PM2.5污染的显著性相对其他组最小,且小于0.01,相关系数相对其他组最大。说明当天的PM2.5数据与日全因死亡人数具有最大线性相关的可能性,并且相关程度也最高。因此本文将采用死亡日当天的大气污染物数据、气象数据与日死亡人数进行回归分析。

表6 日全因死亡人数(D)与各指标偏相关分析Table 6 Analysis of partial correlation between daily all-cause deaths (D) and various indicators

利用SPSS软件将死因数据中日死亡人数的自然对数lnD的计算结果与大气污染物指标及气象指标进行相关分析,结果见表7。通过表7可以看出,lnD与大气污染物指标中的PM2.5、PM10、SO2、NO2的日平均质量浓度以及气象指标中的平均温度、平均气压P0的显著性小于0.01,具有显著的线性相关可能性与明显的统计学意义;lnD与CO日平均质量浓度、平均湿度的显著性大于0.01,线性相关可能性软弱;lnD与O3日平均质量浓度、平均风速无线性相关性。因此,PM2.5、PM10、SO2、NO2日平均质量浓度、日平均温度、日平均气压P0均可作为预测日死亡人数模型的主要变量。

表7 ln D与各指标相关分析Table 7 Correlation analysis of ln D and each index

2.4 主成分分析

主成分分析是一种将多个指标化为少数几个不相关的综合指标(即主成分)的统计分析方法,该方法适于分析多指标的大量数据间的关系和趋势,也是数学上用来降维的一种方法。主成分载荷是主成分分析中原始指标与主成分之间的相关系数,主成分载荷的正、负决定原始指标对主成分影响的方向是正面还是负面的,主成分载荷的绝对值大小表示原始指标对主成分影响的强度,绝对值越大表示原始指标对主成分越具有代表性。

为使模型的自变量降维,从而达到简化模型的作用,利用SPSS软件对PM2.5、PM10、SO2、NO2的日平均质量浓度、日平均温度、日平均气压P06个评价指标进行主成分分析,提取出了2个不相关的综合指标作为主成分进行分析,主成分载荷矩阵如表8所示。从表8中可以看出,在第1主成分中的最大载荷为ρPM2.5,为0.934,在第2成分中的最大载荷为P0,为0.955,说明第1主成分中ρPM2.5对该主成分最具有代表性,第2主成分中P0对该主成分最具有代表性。因此,为使模型自变量降维,选取PM2.5的日平均质量浓度ρPM2.5和日平均气压P0与lnD进行回归分析。

表8 评价指标主成分载荷矩阵Table 8 Evaluation index principal component load matrix

2.5 模型拟合

利用SPSS软件的回归分析功能,将本溪市2019年1月1日—12月31日的死亡数据中日死亡人数的自然对数lnD作为因变量,将评价指标中PM2.5日平均质量浓度ρPM2.5、日平均气压P0作为自变量进行回归分析,拟合模型系数及95%置信区间如表9所示。拟合模型调整后的R2为0.992,大于0.8,说明拟合模型能够解释因变量99.2%的变异,拟合精度高。表9显示评价指标ρPM2.5回归系数为0.001 1(0.000 3～0.001 9),评价指标P0的回归系数为0.027 8(0.027 4～0.028 2),2个回归系数估计值的95%置信限均分别在0的同侧,且2个回归系数P值均小于0.05,表明2个回归系数均具有显著性与明显的统计学意义。根据分析结果,以泊松回归模型为基本模型,拟合本溪市PM2.5污染与死亡人数的暴露-反应关系模型表达式为

表9 系数及95%置信区间结果Table 9 Coefficient and 95% confidence interval results

lnD=0.001 1ρPM2.5+0.027 8P0。

3 结 论

本文研究发现,本溪市PM2.5污染对各类死因的影响不显著,O3对各类死因的影响显著;PM2.5日平均质量浓度、日平均气压与当日居民全因死亡人数的自然对数之间存在显著关联。当PM2.5日平均质量浓度上升10 μg·m-3时,带入拟合的日死亡人数泊松回归模型中计算可得本溪市日死亡人数上升1.1%,而我国其他城市的同类研究结果在0.36%～0.85%之间[12-14]。研究结果表明,本溪市居民死亡率受PM2.5污染的影响程度相比国内其他城市更为严重。