非线性效应的多元meta分析及其应用*

陈金建 杨 军 陈仁超 刘起勇 欧春泉

【提 要】 目的 介绍多元meta随机效应模型评价多个城市暴露-反应非线性关系的理论，并通过实例加以详细阐述。 方法 多元meta随机效应模型广泛应用于处理多重相关参数的复杂效应研究，本文基于该方法评估多个地区逐日气温对人群非意外死亡人数的非线性影响。建模过程主要基于二阶段策略，首先采用分布滞后非线性模型(distributed lag non-linear model,DLNM)构建气温-死亡-滞后效应的二维交叉基，分别估计各个城市气温与死亡的非线性暴露-反应关系，并提取各城市回归系数与协方差矩阵；然后，应用多元meta随机效应模型对各城市估计的参数进行合并，得到气温对人群非意外死亡的综合效应。结果 我国16个城市的气温对人群非意外死亡人数存在滞后影响，总体效应关系呈现U形，提示低温与高温对人群健康影响较大，且低温效应大于高温效应。结论 多元随机效应meta模型可拟合不同研究间多个复杂的相关参数，非常适合于暴露-反应呈现非线性的效应综合研究。

气候变化背景下，环境危险因素的健康影响愈来愈受到各级政府、公众与科学家关注。以往大量研究表明气温与健康结局(如：发病率与死亡率)呈现V、U或J型的非线性的暴露-反应关系[1-3]。即当气温高于或低于某个临界值时，健康结局发生率会随着温度的升高或降低而上升。由于不同地区间气候、人群适应性与社会经济结构存在较大差距，难以将单个城市的研究结论推广至其他地区。因此有必要采用meta分析的技术将不同城市或地区的单独研究结果进行合并，从而得到更有说服力的证据。单变量的meta分析模型仅针对温度两端的效应(低温与高温)进行合并，无法获得多个地区气温与健康的综合反应关系曲线，同时对温度滞后效应的拟合不够充分，未能客观反映气温对死亡的非线性效应。2012年，Gasparrini[3-4]将传统多元meta分析的技术引入气温健康效应研究中，重新阐述该理论实现过程，并制作了”mvmeta”的R软件包，为实现多个地区暴露-反应非线性效应的合并提供便利。

因此，本研究将详细介绍多元随机效应meta分析的理论过程，并以中国16个城市为实例应用该模型进行分析。

原理与方法

暴露-反应关系的合并通常采用标准的二阶段建模分析策略：第一阶段得到单个研究城市的效应曲线，本文运用分布滞后非线性模型(distributed lag non-linear model,DLNM)拟合各城市非线性的气温与人群死亡结局的暴露-反应关系[5]；然后第二阶段使用扩展的多元meta分析方法综合多个城市的非线性暴露-反应关系。

1.第一阶段模型

(1)广义相加模型基本结构[5-7]：

(1)

其中g为log函数；μt≡E(Y)，Y可为多种指数分布族(exponential family)，如正态分布、二项分布、Gamma分布、Poisson分布等；在评估环境危险因素健康效应的时间序列研究中，因变量Yt(t=1,2,…,n)通常是某一地区人群中某种阳性事件的逐日累计人数(如每日死亡人数、每日患病人数等)，与当地人口基数相比，其为小概率事件，且时常出现过度离散现象(方差大于均数)，故通常采用广义Poisson分布；自变量xj通常是同期的逐日空气污染物浓度、温度、相对湿度等环境因子，uk表示其他混杂因素，βj、γk为相应的系数；fj表示自变量xj的各种基函数，用于拟合非线性的自变量，常用的参数类型的基函数有正交函数、多项式、线性阈值函数、自然三次样条和B样条等，非参数类型的基函数有局部回归函数和惩罚样条函数。

(2)分布滞后非线性模型

DLNM模型以广义线性模型和广义相加模型等传统模型的思想为基础，利用交叉基(cross-basis)过程，可同时灵活拟合暴露与健康非线性及其滞后效应，因此非常适合于温度健康效应研究。具体地，定义长度为N的向量X=[x1,…,xt,…,xN]T,L是最大滞后天数，l=[0,…,l,…,L]T，设定φ1(.),φ2(.),…,φvx(.)是拟合暴露-反应维度的基函数；ψ1(.),ψ2(.),…,ψvl(.)是拟合暴露滞后效应维度的基函数；vx,vl是暴露/滞后基函数的维度。可得如下公式[5]：

(2)

β=Mη

V(β)=MV(η)MT

其中：M=1(L+1)C⊗Ivx是简化矩阵；β=(β1,β2,…,βvx)是简化后的参数

2.第二阶段模型

(1)多元meta模型

(3)

其中Si+Ψ=Σi。公式(3)包含暴露在研究内与研究间独立的效应变异，前者可由抽样误差Nk(θi,Si)表示，服从k维正态分布；后者来自Nk(θ,Ψ)总体,Ψ为研究间待估计的协方差矩阵。当存在有p个可用来预测结局的协变量xi=[x1i,x2i,…,xpi]T时，模型(3)可进一步扩展为多元meta回归模型：

(4)

其中：

Xi(k×kp)=I(k)⊗xiT=

另外，若假定各个研究与总体的真实效应之间差异的随机效应为ui=Nk(0,Ψ)，则模型(4)可用另一种形式表示：

(5)

(2)系数的估计

模型(3)系数的估计主要基于最大似然和限制性最大似然法。首先，模型(3)的边际对数似然函数为：

(6)

其中：Σi(ξ)简记为Σi，n为观测值总数，若无缺失时其等于km。

(7)

但若Ψ未知时，另一种基于n-q个线性误差对比的估计法更为合理，q为β中固定效应系数个数，采用限制性最大似然法得到：

(8)

(3)假设检验与模型预测

随机效应meta分析关注点在于不同研究间是否存在显著的异质性，其原假设H0为:Ψ=0，备择假设H1:Ψ≠0。将Cochran Q异质性检验扩展到多元情况得到：

(9)

进一步地，可得到以下两个评价异质性的具体指标：

(10)

其中H2为异质性大小；I2为协变量不能解释部分占比。

实例分析

从中国疾病预防控制中心获得我国16个省会城市2007-2013年逐日非意外死亡人数数据，同时期逐日平均温度、相对湿度与大气气压资料来自于中国气象数据共享网。16个省会城市包括哈尔滨、长春、沈阳、北京、天津、银川、济南、郑州、上海、南京、合肥、成都、武汉、重庆、长沙、广州，各城市气象与死亡数据统计描述见我们前期发表的文章[8-9]。首先采用分布滞后非线性模型拟合各城市温度-死亡之间的非线性暴露-反应关系，其中温度和滞后维度均采用自由度为5的三次样条函数，同时采用三次样条函数控制逐日湿度、逐日大气气压、时间变量的混杂，这些变量自由度分别为3、3、8/年；以上自由度采用Q-AIC最小化准则进行筛选得到；温度最大滞后天数设定为21天[4]。得到各城市单独效应后，采用多元meta随机效应模型进行效应合并。

表1 各个城市第一阶段模型的系数估计

*：b1～b5是由于拟合气温变量时自由度为5，待估计的参数为5个

表1展示了第一阶段中估计得到的气温的模型系数。在得到各个城市系数估计及相应的协方差矩阵后，可利用第二阶段多元meta随机效应进行效应合并。

图1 16个省会城市气温对非意外死亡的效应曲线

图1提示16个城市温度-死亡效应曲线呈U型，最适(死亡风险最低)温度为17℃，低温与高温均使人群死亡风险上升。相对于最适温度，第1百分位数气温(低温)与第99百分位数气温(高温)的相对危险度分别为2.039(95%CI:1.744～2.384)与1.390(95%CI:1.171～1.650)。进一步，以秦岭-淮河南北分界线为界，得到南方低温与高温的相对危险度分别为2.309(95%CI:1.773～2.995)与1.175(95%CI:1.030～1.341)，北方分别为(95%CI:1.521～2.277)与1.566(95%CI:1.268～1.2.099)。

讨 论

本文介绍的多元meta随机效应模型具有综合效应关系的功能。它是传统的标准两阶段设计的扩展，可对信息量更丰富的多元关联参数进行综合估计。Sauerbrei和Royston[10]提出另一种通过对暴露的效应进行网格化，以期构建总体非线性关系的meta分析方法，称为meta平滑(meta-smoothing)。该方法应用前提是假设各网格的效应为独立的，但通常情况该假设条件难以满足，因此在暴露效应值的标准误和置信区间的估计可能存在偏倚。多元meta随机效应模型通过估计构建协方差矩阵来允许各城市的效应估计之间存在相关性。可针对模型拟合与比较、异质性和预测分别进行统计检验，并且在参数估计方面有更好的统计性能，比如可通过估计研究间协方差矩阵结构来提高精度。但也有一定的不足[3]，需要估计的模型参数繁多，包括k维固定效应参数向量，p个个体的参数个数为kp，协方差参数个数为k(k+1)/2；另外，要求所有个体研究取相同的暴露-反应函数，可能造成部分信息的缺失。

实例应用中探究了中国16个城市的每日气温和非意外死因之间的关系，得出气温-死亡综合的暴露-反应关系曲线呈U形，提示低温与高温会在一定程度上增大人群死亡率，并且低温比高温的影响更大，这与以往研究相一致[1，11-12]。

综上所述，多元meta分析可非常适合应用于多个地区暴露与健康呈现复杂的非线性关系的环境流行病学研究。另外，该方法还可应用于多站点临床试验、多中心队列、多个癌症登记研究或多城市时间序列研究等；其研究目的还包括合并生存曲线[13-14]、纵向剖面[15]，ROC曲线[16]和热浪效应[17]等不同研究问题，具有非常良好的应用前景。