多重稳健的高维缺失数据插补研究

熊 巍,王 娟,潘 晗,田茂再

(1.对外经济贸易大学 统计学院,北京 100029;2.北京大学 数学科学学院,北京 100871;3.中国人民大学 统计学院,北京 100872)

一、引言

在实际应用中,由于客观条件的限制以及数据收集渠道、数据结构的差异化,某些数据通常被标记为未知、空白等,这种数据即为缺失数据。缺失数据问题普遍存在于抽样调查、社会科学、教育学、流行病学等许多领域。新冠病毒感染疫情爆发以来,由于不同地区统计口径、时间口径以及数据收集手段的差异,某些地区的病例详情在一段时间内会出现缺失,给研究者们开展实际分析带来重重困难。若不对缺失数据加以处理,统计分析工作中的有效样本量将会缩小,导致统计效率下降及分析结果偏误。另外,随着科技的发展以及可获取数据维度的增高,高维缺失数据已成为常态,特别是高维随机设计阵中缺失项通常具有厚尾特征,使得利用常规处理方法得到的结果往往不稳健,难以有效代表总体。面对高维特征,如何利用不完全数据对总体做出科学准确的推断,已成为数据科学时代统计研究的难点和热点问题。据此,本文围绕高维缺失数据,综合使用加法模型和增强的逆概率加权估计并融入协变量平衡倾向性评分,提出一类新的稳健有效的估计量,不仅发掘未充分利用的潜在信息,提升数据分析效率,还能避免高维非参数问题中的维数灾难,为有效处理高维缺失数据提供新的研究视角和可靠的统计工具。

本文重点研究随机缺失机制下高维数据的稳健插补问题。自Rubin提出缺失机制的概念以来[1],国内外学者对缺失数据问题的研究日益增多,并涌现出大量成果探讨缺失数据的处理方法,最常用的是插补法。研究者们可以利用最可能的值去替代缺失值,如人工填写、均值插补等;也可以借助统计模型或机器学习方法实现缺失值的填充,如回归插补、基于极大似然估计、聚类和关联规则的插补等。而在众多的缺失数据插补研究中,最具有代表性的两类插补方法为基于倾向性评分的插补方法以及逆概率加权法(IPW)。倾向性评分由Rosenbaum和Rubin提出,用来度量实验设计中不同治疗组间的处理效应,其基本原理是将多个协变量的影响用一个倾向评分值加以表示,并根据评分值进行不同组间的匹配(PSM)或加权(PSW),使观察数据达到“接近随机分配数据”的效果,实现有效降维[2]。倾向性评分模型自提出以来被广泛应用于各个领域,如抽样调查中对总体参数的推断问题,金勇进和刘晓宇考虑到权数对模型推断的影响,提出将权数引入倾向性评分模型及预测模型,构造了一种双重稳健估计方法,有效提升了处理效应估计结果的准确率[3]。逆概率加权法由Horvitz和Thompson提出,要求先估计缺失概率,再利用缺失概率加权确定最终估计结果,但该方法仅当缺失概率已知或模型被正确参数化时表现良好[4]。为此,Robins等提出了IPW估计量的增强版本,并表明该类估计量具有双稳健性,即缺失概率结构和均值回归函数的参数结构只要有一个设定正确,IPW估计量就具有相合性[5]。这一优良性质激发了学者们借鉴逆概率加权思想对随机缺失机制下的模型选择与估计问题进行更深入的研究[6-7]。

随着可获取数据维度的不断增高,数据的缺失比例也在逐步攀升,若每一特征都以一定概率缺失,缺失数据的规模将相当可观。然而,高维数据的复杂性伴随着其稀疏特性使得传统插补方法几乎失效,如Garcia等采用的基于EM算法的插补方法在维度增大时,计算量大幅度地膨胀,估计量的有效性极大削减[8]。为提高高维缺失数据插补的有效性同时减缓运算压力,近年来降维技术广泛地应用于缺失值插补中。祝丽萍和邵伟利用多元函数提出一种半参数降维方法,将高维数据转化为一元数据,并利用核回归实现了缺失值插补,在避免维数灾难的同时,保证了估计量的无偏性,但该方法在数据转换过程中会损失原始数据信息[9];熊巍等针对高维成分数据结合Lasso分位回归及修正的EM算法提出了一种稳健的近似零值插补方法,但该方法需对原始数据进行一定的预处理[10]。Wang等利用增强的IPW并结合充分降维方法,将高维协变量映射到低维空间,实现了响应变量的均值插补及双稳健特性,然而一旦模型被错误设定,该方法的估计效率会大打折扣[11]。可见上述方法或破坏原始数据结构,或存在模型误设的风险,不宜直接应用于高维缺失数据,特别在实际应用中,鲜少有先验信息表明响应变量的条件期望是线性函数或属于一个有限维参数空间。邰凌楠等结合倾向评分匹配和模型平均,提出了逆概率多重加权估计,并在理论性质及实际应用中证实了倾向匹配方法可以显著增强估计量的稳健性及有效性[12]。然而如上方法需要准确的倾向性评分估计,轻微的倾向评分模型的误设都会导致处理效应的严重偏差[13-14]。这也凸显了倾向性评分的矛盾性:其本质为降低协变量的维数,但其估计需要高维协变量的合理建模。虽然已有各种方法对PSM和PSW进行改进和完善[7,15],但尚缺乏稳健性。而Imai和Ratkovic提出的协变量平衡倾向性评分(CBPS)方法[16]通过选择使协变量平衡最大化的参数值能够消除倾向评分模型误设带来的影响,当PSM和PSW表现不好时,CBPS却能显著提升性能。目前,CBPS方法已应用于缺失插补领域,然而相关方法仍基于线性模型[17]。综上,已有高维缺失插补文献或基于线性模型,或仅考虑IPW估计和CBPS中的一种方法,不够稳健,难以处理高维随机设计阵中具有厚尾特征的随机缺失项。于是本文综合利用IPW估计和CBPS方法的优势,并通过引入一类更加灵活的非参数加法模型[18],以期对高维缺失问题的理论建立带来有效改善,实现双重降维及多重稳健性。

加法模型(AM)如式(1)所示:

(1)

其中,α为截距项,{mk(·),k=1,2,…,d}为未知的光滑函数集,ε为误差项,并满足E(ε)=0,var(ε)=σ2<∞。为保证式(1)的可识别性,假定每一特征满足E(mk(Xk))=0。由式(1),加法模型包含了线性模型,又显著提高了线性模型的灵活度,允许协变量以一种全新的解析模式进入线性模型,具有如下优势:一是能够捕捉响应变量与协变量之间的非线性关系,提高建模灵活度;二是在因变量分布不确定或不符合正态假设下依然适用;三是保持了非参数模型数据驱动的优势,同时避免了缺失值插补下协变量较多引发的维数灾难;四是有效获取各协变量的边际效应,在构建的插补估计量中能够充分利用协变量信息,提高插补精度。可见,应用加法模型对高维缺失数据进行插补是一种合理的选择。

本文的主要贡献及创新性如下:首先,已有缺失数据插补的文献多基于线性模型,少数对高维缺失数据插补的稳健性进行研究。特别是文献中利用非参数加法模型对高维缺失数据进行插补的理论研究尚属空白。为此,创新性地将增强的IPW方法与加法模型AM融合,应用协变量平衡倾向评分法CBPS估计缺失概率,提出一种适用于高维缺失数据的可加协变量平衡倾向评分插补方法(CBPS-AM)。该方法不仅具有多重稳健性,起到双重降维的作用,还能实现建模的灵活性。其次,借鉴广义矩估计方法(GMM)和Backfitting算法给出了CBPS估计算法[19]。该算法简洁有效,能够提高数据使用效率与插补精度。最后,基于广义矩估计理论并综合应用Slutsky定理和加法模型的性质,证实了满足一定条件下,所提CBPS-AM估计量具有相合性和渐近正态性。模拟研究与实证数据分析表明,CBPS-AM方法适用于多种场景下的缺失数据插补问题,可为高维缺失数据的稳健插补提供合理的研究思路和有效的理论框架,也能为极端突发事件的预测分析提供一定的启示。

二、CBPS-AM方法与渐近理论

假设有一个样本容量为n的观测数据集{(xij,yi):i=1,2,…,n,j=1,2,…,d},X=(X1,X2,…,Xd)为d维协变量,令xi=(xi1,xi2,…,xid)′,Y为存在部分缺失的一维响应变量。本文仅考虑最常用的随机缺失机制,即二元指示变量T的取值仅与X有关:

P(T=1|Y,X)=P(T=1|X)=π(X)

(2)

其中,π(·)为选择概率(缺失概率)函数,π(x)=P(T=1|X=x)。当响应变量Y存在随机缺失时,现有样本的分布并不能代表总体的真实分布,因此仅利用完全观测到的数据对总体进行推断会引起偏误。为避免数据缺失带来的统计推断不准确、插补失效等问题,本节从逆概率加权估计出发,逐步提出适用于高维缺失数据的CBPS-AM插补方法并研究其理论性质。

(一)逆概率加权估计IPW及协变量平衡倾向评分法CBPS

当响应变量出现随机缺失时,Horvitz和Thompson提出IPW法估计总体均值μ[4]:

(3)

(4)

此外,在观测性研究中,式(2)中的π(X)也称为倾向性评分(PS)值,通常未知。为进行估计,研究者们经常假定其具有某种参数结构,如Logistic模型:

(5)

其中,β为未知的d维参数向量。若再假定πβ(·)关于向量β二阶连续可导,易得估计方程:

(6)

(7)

(二)可加协变量平衡倾向得分估计方法CBPS-AM及算法

已有高维缺失插补方法或基于线性模型,或仅考虑IPW估计和CBPS中的一种方法,不够灵活稳健,难以处理高维随机设计阵中具有厚尾特征的随机缺失项。于是本文综合利用IPW估计和CBPS方法的优势,并融入加法模型,创新性地构建CBPS-AM方法,有效避免维数灾难的同时,实现高维缺失数据的稳健插补。本节给出CBPS-AM方法以及相应算法:

第1步,利用CBPS估计选择概率π(xi)。相应估计方程为:

(8)

对于式(8),利用广义矩估计方法(GMM),得到如下广义矩估计量:

(9)

第3步,得到CBPS-AM插补估计值。基于以上两步,总体均值的CBPS-AM估计为:

(10)

于是,CBPS-AM算法汇总如下:

输入:一组观测数据(yi,xij),i=1,2,…,n,j=1,2,…,d,其中Y为响应变量并满足随机缺失机制,Ti为指示Yi是否缺失的二元指示变量,以及任意小常数γ。

(3)重复;

(4)forl=1 toK,j=1 tod,执行计算

综上,CBPS-AM方法具有如下优势:

其一,该方法以增强的IPW为基础,具有双稳健性,只要缺失概率与响应变量回归模型有一结构设定正确,估计量就具有相合性。当上述两个模型均被错误设定时,IPW方法表现非常差,而本文提出的CBPS-AM克服了这一问题,且在实际数据关系难以准确设定时仍能得到较好的估计结果(见模拟例子)。

其二,应用CBPS方法估计缺失概率,可以有效改进缺失概率结构误设所带来的偏误,具有本文所谓的“稳健性”;并且CBPS中通过将高维协变量汇总为一个评分值,能够有效处理高维协变量情形,此为本文所谓的“第一重降维作用”。

其三,使用加法模型模拟均值回归函数灵活高效,并能避免高维非参数的维数灾难问题,此为“第二重降维作用”。估计过程中假定各分量间相互独立,估计效率可以大大提升。可见随机缺失机制下,利用本文的CBPS-AM方法进行插补能够达到多重稳健、双重降维的效果。

(三)CBPS-AM渐近理论

下面给出所提估计量CBPS-AM的渐近性质,对此给出如下所需要的条件和假设以及记号。令v(x)≡var(Y|X=x),hj表示核函数中的带宽,j=1,2,…,d。

(A1)选择概率(倾向得分)满足0

(A2)πβ(x)关于β二阶可导,且有下界,即πβ(x)≥c0(c0为一正常数);

(A3)设计矩阵X的概率密度为f(x),0

(A4)v(x)为连续函数且严格正,并具有二阶连续导数;

(A5)对于每一个可加函数mj(x)均二阶连续可微,j=1,2,…,d;

(A6)核函数K(·)是定义在[-1,1]区间上的对称密度函数,且满足李普希茨(Lipschitz)条件;

(A7)当n→∞时,有n1/5hj→δj(δj为正常数),j=1,2,…,d。

基于以上条件,可以得到如下定理1和定理2。

定理2在满足条件(A1)～(A7)时,本文所提的CBPS-AM插补估计量具有渐近正态性:

三、模拟研究

例1(低维情形,d=4) 响应变量Y分别来自如下三种不同模型:

模型1:Y=X1+X2+X3+X4+ε,存在模型误设;

模型2:Y=sin(X1)+sin(X2)+sin(X3)+sin(X4)+ε,存在模型误设;

模型3:Y=sin(X1)+sin(X2)+sin(X3)+sin(X4)+ε,模型设定正确。

log(π(X)/(1-π(X)))=-0.2(X1+X2+X3+X4)-0.6

表1 不同插补方法在例1下的模拟结果(BIAS×0.01(RMSE))

第三,各插补方法结果在异方差情形与同方差情形下相似,除CBPS-AM与AM方法,其余五种方法均表现较差。一方面说明,与增强的逆概率加权方法相比,传统的IPW方法不具有稳健性,估计量的相合性难以保证,而基于线性模型的插补方法由于缺乏灵活性,也可能导致较高误差;另一方面表明,缺失概率估计中,相比于Logistic回归模型,更宜选择CBPS方法。

综上,本文提出的CBPS-AM方法在模型和缺失概率结构存在误设与否、均值回归函数线性或非线性、误差方差是同方差还是异方差情形下,与现有插补方法相比都具有更好的估计结果,表明了CBPS-AM方法在缺失插补中的优越性及多重稳健性。此外,随样本容量的增加,CBPS-AM方法的均方误差也逐渐减小,与定理结论保持一致。

例2(高维情形,d=100) 本例考虑高维情形,响应变量Y来自如下模型:

Y=g1(X1)+g2(X2)+g3(X3)+g4(X4)+1.5g1(X5)+1.5g2(X6)+1.5g3(X7)+1.5g4(X8)+

log(π(X)/(1-π(X)))=0.2(X1+X2…+X5-X6-X7…-X10+2X11+2X12)+0.6

考虑n=50,200,1 000。令t=0,1,当t=0时,协变量之间相互独立;当t=1时,协变量之间中度相关,相关系数为0.5,插补结果如表2所示。

综上,本文提出的CBPS-AM方法在高维缺失数据中具有良好的适应性和稳健性,不仅能够有效避免高维“维数灾难”,还适用于自变量高度相关、误差项厚尾及偏态分布的情形。

表2 例2不同插补方法模拟结果(BIAS×0.01(RMSE))

四、实证研究

(一)HIV临床实验数据的实证研究

艾滋病由人类免疫缺陷病毒(简称“HIV”)引起,已成为21世纪威胁人类健康最大的传染病之一。HIV会破坏人体免疫系统,使其丧失抵抗各种疾病的能力,进而严重危害生命。人类免疫系统的CD4细胞在抵御HIV的入侵中起着重要作用,当CD4被HIV感染裂解时,其数量会急剧减少,HIV将迅速增加,终至艾滋病发作。于是需要有效的治疗手段尽量减少人体内HIV的数量,同时激活更多的CD4,以提高人体免疫能力。为探究不同治疗方案对患者的疗效是否存在显著差异,本文将CBPS-AM方法应用于HIV数据集以建立科学的综合评价以及针对疗效的合理预测。

HIV数据集包含2 139个HIV感染患者[20],记录了每个患者的人口统计学特征,如年龄、性别、种族、有无病史及给予抗转录病毒治疗后的一些生理指标在不同阶段的表现。本文主要研究不同治疗方案对患者的治疗效果是否存在显著差异,响应变量为:接受抗转录病毒治疗后96±5周的CD4数量(单位为:每立方毫米)。CD4 96±5细胞数量可以衡量艾滋病的发病风险,其水平越低患病风险越大。由于跟踪病患会出现死亡及中途退出,CD4 96±5存在缺失,总体缺失率高达37.26%。本文还考虑了6个连续协变量,分别为:年龄、体重、CD4 0±5、CD4 20±5、CD8 0±5以及CD8 20±5。试验小组随机地将这些病患分成四组,对应四种不同日用药,见表3。治疗现状(CD4 96±5)与CD4和CD8计数在0±5周、20±5周记录有关。

表3 病患分组及治疗方案

图1 四种治疗方案下96±5周CD4计数水平估计

缺失的96±5周的CD4计数受协变量的影响,如若跟踪患者不能接受有效的治疗,感染HIV低基线CD4计数的患者可能更容易退出。本文假定跟踪受访是随机缺失的,即只与六个协变量有关,而与缺失的96±5周CD4计数无关,这在实际应用中也是较为合理的。

为对比分析,利用观测数据分别计算如上所述七个估计量,四种不同治疗方案下各估计结果在图1中给出,横坐标表示四种不同治疗方案,纵坐标代表经过治疗后的CD4水平。图1表明,七种插补方法得到的估计结果相差不大。对比各治疗方案疗效可以发现在艾滋病的治疗中,组合治疗均优于单一疗法,单一疗法中新药DID疗效更佳;和新药DID搭配服用的治疗方案3在所有方案中疗效最佳。

表4 各插补方法在HIV临床试验数据的预测能力

表4给出了200次模拟下预测误差指标的平均值,进一步表明:(1)直接删除缺失值预测性能最差,因其损耗原始数据信息,这也阐明了缺失值插补的必要性;(2)与模拟研究结果类似,基于CBPS和线性模型的插补方法预测性能相对较差,而基于加法模型的插补方法不论在预测的准确性和稳定性上都表现较优,尤其本文提出的CBPS-AM方法具有最小的预测误差。因此,CBPS-AM方法能最大程度恢复数据原貌,更好地提升数据的预测性能。

(二)中国新冠病毒感染疫情数据的实证研究

2020年初,一场突如其来的新冠病毒感染疫情在全球大范围蔓延。面对严峻形势,中国政府积极采取防控措施切断传染源,对重点人群进行核酸检测,做到“应检尽检”,启用“小汤山”模式对确诊患者实施医疗救治,使得疫情得到有效控制,疫情态势逐渐向好。防控过程中,为有效把握疫情发展态势,降低疑似病例流动率,及时掌握疫情动态变化是相关部门的重点关注问题。然而这类极端事件更易出现数据缺失,若利用缺失数据进行统计预测可能高估疾病流行的严重程度,引发公众恐慌,对遏制疫情产生极大的挑战。于是本文将CBPS-AM方法应用于新冠病毒感染疫情数据,以期给出对于疗效和疫情动态的合理预测。

本文采用的数据时间跨度为2020年2月6日至2020年4月27日,变量包括省份、市/区、省累计确诊数、治愈数、死亡数、疑似病例数以及市/区确诊、治愈、死亡和疑似病例数,并人为产生100个来自标准正态分布的噪声变量以满足高维设定。本例主要研究缺失情形下,如何对日市/区的疫情状况进行准确预测和插补。为此,在加法模型的建模中,本文将响应变量设定为市/区确诊病例数。为保证算法正常运行,本文对省份及市/区数据进行如下编码:以武汉为中心,将湖北省设定为1,根据武汉到其他省会城市的距离,将其他各省份编码为2,3,4,…,各市区对应编码为1.01,1.02,…,2.01,2.02,…。为验证本文的CBPS-AM方法的插补效果,随机抽取2月13日、2月25日、3月11日、3月17日和4月26日五天的疫情数据,变量描述统计如表5所示。

表5 各变量描述统计结果

表6 各插补方法在疫情数据中的预测能力比较(BIAS(SSE))

表6展示了七种插补方法的偏差与平均误差平方和。从偏差来看,CBPS估计量及CBPS1估计量分别在2月13日和2月25日插补偏差最小,表明了CBPS方法在选择概率估计中的稳健性;而本文提出的CBPS-AM方法在余下三个日期中均有最小的插补偏差,体现了加法模型应用的合理性。另外,基于线性回归的插补估计量不论在偏差还是误差平方和方面均表现较差,也说明了线性模型在实际应用中的局限性。综上,本文提出的CBPS-AM方法不仅能实现稳健准确的缺失数据插补,在极端突发事件的预测性能上也具有明显的优势。

五、结论

本文构建了一种适用于高维缺失数据的插补方法CBPS-AM,并证明了该方法的相合性及渐近正态性。与现有高维缺失数据插补方法相比,CBPS-AM方法结合了逆概率加权与协变量平衡倾向评分方法的优势,有效避免了模型误设和缺失概率结构误设对插补结果造成的影响,保证了算法的多重稳健性;同时在随机缺失机制下创新性地应用加法模型,不仅实现了建模的灵活性,也与CBPS方法一起实现了双重降维。模拟研究和实证分析表明,本文提出的CBPS-AM方法具有广泛适用性,适用于多种情形和场景下的均值插补问题,包括:(1)线性或非线性均值回归函数结构;(2)模型存在误设与否(包括缺失概率结构);(3)误差同方差或异方差,以及误差项来自正态分布或偏态、厚尾分布;(4)协变量高度相关;(5)医药卫生等领域的高维及非高维缺失数据的插补问题。

充分认识到数据时代数据缺失问题在统计工作中的重要地位,本文得到如下启示:第一,缺失值插补已成为统计分析工作的重要环节,提出的CBPS-AM方法可以有效提升数据质量,为后续统计分析工作提供强有力的支持;第二,CBPS-AM方法为高维缺失数据的稳健插补提供了一套合理的研究框架;第三,未来研究中,一方面可以借助于定理2中的渐近方差进一步探究CBPS-AM估计量的统计推断问题(区间估计和假设检验),另一方面可以将CBPS-AM理论框架拓展至其他非参数模型,如变系数模型、部分线性可加模型等。为更好地进行高维降维,可以将本文提出的方法与充分降维、投影寻踪方法等进行结合,也可以将该方法拓展到分位数回归插补中,实现更稳健的插补,或通过并行运算实现大数据的稳健插补。