基于MRCD估计的多元线性回归模型的稳健估计

颜海波, 邓 罡, 姜云卢*

(1.暨南大学公共管理学院, 广东广州510632; 2.暨南大学经济学院, 广东广州510632)

回归分析是指在统计理论指导下，对统计数据进行数学处理，构建回归模型，借以反映变量间相互依赖关系和进行外推建模的一种建模技术。当回归分析中的自变量和因变量均有2个或2个以上，且自变量和因变量之间是线性关系时，则将这种模型称为多元线性回归模型(multivariate linear model)。自提出以来，多元线性回归模型已被广泛应用于医学、生物、气象、图像识别、交通运输、物流、经济管理等方面的影响因素分析和预测研究[1-5]。

由于存储设备容量越来越大，价格越来越低，同时云存储技术也在快速发展，各种数据都能得以保留。随着各行业的快速发展，产生的数据越来越多，包含了各种结构，形式也更加多样。因此，需要处理的数据更加复杂且与传统数据有很大的差别，不仅数据出现异常值的概率大大增加，数据的维数也在不断增加，经常出现高维的情况。在多元线性回归模型的估计中，最常用的方法是最小二乘估计。由于最小二乘估计十分依赖于对样本均值向量和样本协方差矩阵进行运算，而样本均值向量和样本协方差矩阵对重尾分布和异常值非常敏感，因此，最小二乘估计十分容易受到异常值的影响，使得估计出现误差，甚至得到错误的结论。另一方面，在高维数据特别是维数大于样本量的数据中，样本协方差矩阵不再可逆，传统的估计方法也不再适用。

最早把稳健估计方法应用于多元线性回归模型的是Koenker 和Portnoy[20]，他们把以凸函数作为影响函数ρ的M估计应用于每个响应向量的坐标。随后Bilodeau和 Duchesne[21]推广了Davies[22]的S估计，Roelant等[23]推广了Croux等[24]的GS估计，Ben等[25]推广了τ估计。Rousseeuw[13]于1984年提出了MVE(minimum volumn ellipsoid)和MCD(minimum covariance determinant)方法，该方法的崩溃点值最高达到50%，但由于计算过于复杂且当时计算机技术并未达到要求，未能获得广泛应用。Gao[26]于2017年提出了基于深度函数的多元线性回归模型，并指出其在稳健线性回归中具有潜在的应用。为了解决现有方法无法在大样本下快速计算的缺点，Rousseeuw和Van Driessen[27]于1999年提出改良后的快速MCD方法(fast-MCD)，也使得MCD方法真正用于稳健统计中，并得到广泛认可。随后Rousseeuw等[28]提出了使用MCD估计来得到多元线性回归的参数估计，基于MCD估计得到的参数估计比普通最小二乘法有更好的稳健性，且与LTS等方法相比可以更好地处理杠杆点和异常点。Agulló等[15]推广了Rousseeuw的MCD估计。然而MCD估计也有局限性，其只有在p5p，随着p的增大，MCD估计的精度会不断降低，其稳健性也随之降低。因此，当两组变量的维数之和大于样本个数时，基于MCD估计方法的稳健多元线性回归模型就不再适用。此外，随着样本维数的增加，当自变量或者因变量个数大于样本量时，其样本协方差矩阵将不再可逆，经典的最小二乘估计方法也随之失效，无法得到系数估计。Boudt等[29]于2020年提出了高维稳健的协方差矩阵稳健估计方法MRCD(minimum regularized covariance determinant)，使得当变量个数大于样本量时也能够获得出比其他方法更稳健的估计。姜云卢等[30]将MRCD估计方法和主成分聚类方法相结合，并证明其比传统方法能更好地抵御异常值，且在高维情况下更为有效。因此，本文选择使用MRCD方法来估计总体均值向量和协方差矩阵，再进行多元线性回归模型的参数估计。本文提出的方法比传统最小二乘估计方法更加稳健，且在高维情况下克服了MCD估计的多元线性回归方法失效的缺点。

本文结构如下：第1章介绍多元线性回归模型的传统最小二乘估计及其在含异常值数据下的局限性，介绍基于MCD估计的多元线性回归模型估计及其在高维情况下的局限性，提出基于MRCD估计的多元线性回归模型，同时通过数值模拟与传统方法和基于MCD估计方法的多元线性回归进行比较；第2章使用实证数据对基于MRCD估计的多元线性回归模型进行进一步分析；第3章基于理论、数值模拟和实证分析的结果提出结论与建议。

1 多元线性回归模型和基于MCD估计的多元线性回归模型的局限性

1.1 传统的多元线性回归模型的局限性

从上述传统多元线性回归模型的最小二乘估计过程中可以看出，模型估计的关键步骤是估计总体均值和总体协方差矩阵，样本均值和样本协方差矩阵在正态分布的总体下是总体均值和协方差矩阵最优的估计，传统的多元线性回归模型的估计就是在正态总体的假定下进行的。然而，在重尾分布的情况下，样本均值向量和样本协方差矩阵对样本中的异常值非常敏感，并不稳健，这就使得传统方法得到的回归系数估计值精度下降，回归系数出现偏差，甚至会得到错误的结论。

1.2 基于MCD估计方法的稳健多元线性回归估计的局限性

① 从样本数据中随机抽取含h个样本的数据子集，重复步骤②～④。

② 计算选定数据子集的均值向量和协方差矩阵：

Si(H)=(h-1)-1(Xi-mi(H))T(Xi-mi(H))。

式中：H表示选定的数据子集；Xi表示第i次迭代的子集数据；mi(H)是第i个子集的均值向量；Si(H)是第i个子集的协方差矩阵。

③ 通过mi(H)和Si(H)计算所有n个数据样本到数据中心miH的马氏距离：

④ 选出所有n个样本中离样本中心mi(H)的马氏距离Di(X)最小的h个样本。

⑤ 当detSi(H)收敛时停止迭代，记SM(H)为最后一次迭代的协方差矩阵。

⑥ 计算均值向量和协方差矩阵的MCD估计：

mMCD=mM(HMCD),

SMCD=cαSM(HMCD)。

通过以上步骤即可得到MCD估计方法的稳健均值和协方差矩阵，并利用其进行进一步的多元线性回归模型估计。基于MCD估计方法的多元线性回归估计受离群值的影响较小，能得到更有效的系数。但与此同时，MCD估计方法也有一些局限性，其只有在p5p的情况下，MCD估计有较高的准确性，而随着样本维数p的增大，MCD估计的精度会不断降低。因此，当自变量或者 因变量的维数大于样本个数时，基于MCD估计方法的多元线性回归模型估计就不再适用。

1.3 基于MRCD估计方法的高维稳健多元线性回归估计方法

1.3.1 MRCD估计方法简介

MRCD估计方法由Boudt等[29]于2020年提出，是一种高维稳健的多变量总体均值向量和协方差矩阵估计方法，其在最小协方差矩阵行列式的子集基础上进行正则化，从而得到估计值，是在MCD估计的基础上进行改进的。与MCD估计方法类似，MRCD估计方法仍然具有稳健性，能够有效地应对异常值和偏态数据，同时，MRCD估计方法还克服了MCD估计方法在p>n的情形下失效的缺点，可以用于任何维度数据的均值向量和协方差矩阵的估计。

在实际运用中，含异常值的数据和高维数据越来越频繁地出现，在医学统计、生物统计、经济管理、图像识别等领域，对高维数据的处理和分析更为普遍，依赖于样本均值和样本协方差矩阵的传统方法不可避免地失效。使用MRCD估计方法可以获得高维数据的总体均值向量和协方差矩阵的估计，使应对含异常值的数据和高维数据成为可能。

1.3.2 基于MRCD估计的多元线性回归模型

与MCD估计方法类似，MRCD估计方法也使用了马氏距离来寻找样本数据的子集，并以此数据子集来估计总体均值向量和协方差矩阵。该方法先进行分位数标准化，再寻找出正则化协方差矩阵行列式最小的数据子集，之后对寻找出的数据子集进行反标准化，从而得到总体均值向量和协方差矩阵的稳健估计。寻找正则化协方差矩阵行列式最小的数据子集使用的是C-step算法，即先随机选择一个包含h1个样本观测的子集，计算均值向量和协方差矩阵，并对h个样本观测计算马氏距离，选择距离最小的h2个观测作为第二个子集，这样进行迭代，直到满足收敛条件。挑选MRCD数据子集的C-step算法[13]步骤如下：

① 初始化6个均值向量mi和6个稳健无病态的协方差矩阵Si(i=1,…,6)；

③ 计算每个子集Hi对应的正则化系数ρi，其中0≤ρi<1，ρi定义为满足ρiI+(1-ρi)cαSW(Hi)的状态良好的矩阵的最小值；

⑤ 对于每个初始化的子集Hi，如果ρi≤ρ，则重复步骤②～④，直至ρiI+(1-ρi)cαSW(Hi)收敛，将最后得到的子集定义为Hm；

⑥ 从6个候选子集中选择使得ρiI+(1-ρi)cαSW(Hi)具有最小行列式的子集，记为HMRCD。

进一步，即可得到整个MRCD估计方法，其简单计算步骤如下：

① 随机选择包含h个样本的数据子集，并进行分位数标准化，

式中：vX为每一个变量中位数向量；DX是每一个变量的绝对中位差(median absolute deviation)估计组成的对角阵。

② 确定目标矩阵T和惩罚项系数，并得到正则化协方差矩阵公式K(H)：

T=cJp+(1-c)Ip,

SU(H)=(h-1)-1(Ui-mi(H))T(Ui-mi(H)),

KH=ρT+1-ρcαSUH。

(1)

④ 使用C-step算法挑选MRCD数据子集

当T、Q和Λ确定时，也可以写为

⑤ 根据选出的子集HMRCD计算MRCD估计的均值向量和协方差矩阵：

mMRCD=vX+DXmU(HMRCD),

通过以上步骤即可得到MRCD估计方法的高维稳健均值向量和稳健协方差矩阵估计，并可以进行进一步的多元线性回归模型参数估计。基于MRCD估计方法的高维稳健多元线性回归模型估计不仅受离群值的影响较小，还克服了在自变量和因变量的个数之和大于样本个数时，基于MCD估计方法的稳健多元回归模型就不再适用的问题。

1.4 数值模拟

假设自变量的维度为p，因变量的维度为q，使用R语言分别生成n=100和n=200，以及p=q=10、p=q=25、p=q=30、p=q=50、p=q=60、p=q=75和p=q=120的服从混合分布Xp,n～(1-α)Np(0,1)+α(γNp(0,1)+10p)和εq,n～(1-β)Nq(0,1)+β(ηNq(0,1)+10q)的数据[32]。式中：10p是各元素均为10的n×p维矩阵；10q是各元素均为10的n×q维矩阵；α是自变量污染比例；β是因变量的污染比例；γ决定自变量异常点的协方差矩阵；η决定因变量异常点的协方差矩阵。数据中的自变量异常值来源于γNp(0,1)+10p的部分，因变量异常值来源于ηNq(0,1)+10q的部分。α和β分别取0、0.1和0.2，γ和η分别取0和1。对同一组数据，分别进行传统多元线性回归模型估计、基于MCD估计方法的多元线性回归模型估计和基于MRCD估计方法的多元线性回归模型估计，每次实验重复100次，并求其均值进行比较。其中，MRCD估计方法使用Boudt等[29]提供的R语言程序包(rrcov)。

在进行结果比较时，每种估计方法的表现使用重复实验的平均SMSE(standard mean squared error)来衡量。SMSE的定义为

由表1、表2可以看出，在低维无异常值的情况下，3种方法的多元线性回归模型估计结果都只有很小的SMSE。在自变量和因变量都存在异常值的情况下，传统方法的多元线性回归模型估计有更大的SMSE，而基于MCD方法和基于MRCD方法的多元线性回归模型估计都取得了更好的效果，即基于MCD方法和基于MRCD方法的多元线性回归模型估计在低维情形下都具有稳健性，在n≥5(p+q)的情形下，基于MRCD方法的多元线性回归模型估计的稳健性要稍弱于基于MCD方法的多元线性回归模型估计。当数据的维度不断上升时，特别是当n<2(p+q)时，基于MCD方法的估计效果开始变差，稳健性也比基于MRCD方法的多元线性回归模型估计差，而在n≤p+q的时候，基于MCD方法的稳健协方差矩阵估计失效，无法进行多元线性回归模型估计，而在np或n>q时，传统方法多元线性回归模型估计失效，而基于MRCD估计方法的估计仍然保持了较小的SMSE。

表1 n=100时3种方法多元线性回归模型估计的数值模拟结果

表2 n=200时3种方法多元线性回归的数值模拟结果

同时，在自变量和因变量的异常值情况相同，且在相同的p和q下，随着n的增大，基于3种方法的多元线性回归系数SMSE均有所下降，但是当p和q较大时，传统方法多元线性回归模型估计和基于MCD方法的估计SMSE仍然很大，说明仍然有很大的误差。而基于MRCD方法的多元线性回归模型估计仍然保持了较小的SMSE。

另一方面，在p、q和n相同的情况下，自变量和因变量同时包含异常值的情况下传统方法和基于MCD方法的多元线性回归模型估计的SMSE均比仅有自变量包含异常值的情况下的SMSE要大，而基于MRCD方法的多元线性回归模型估计在自变量和因变量都包含异常值的情况下仍然保持了较小的SMSE。

这三方面的结果说明，基于MRCD估计方法的高维稳健多元线性回归模型估计在有异常值的高维数据下，特别是数据维数大于样本数的情况下仍然是稳健的。

2 实证分析

本文选择UCI machine learning repository的concrete slump test data set(CST，混凝土坍落度的数据集)[33]，初始数据集有78个数据样本，包含7个自变量和3个因变量，几年后数据发布者又获得了25个新的样本，数据集一共包含103个样本。高性能混凝土是一种复杂的材料，其坍塌度受到了含水量和混凝土成分含量的影响，因此自变量描述了混凝土的成分，包含水泥、粉煤灰和水等成分指标；因变量描述了混凝土的性质，包含凹陷和抗压程度等3个指标。

首先对处理后的CST数据集进行异常值检测。对每一个样本观测值，利用MRCD方法估计的均值向量和协方差矩阵计算稳健马氏距离

图1 CST数据集的稳健马氏距离异常值检测

由图1可以看出，CST数据集中存在部分异常点，如果直接建立多元线性回归模型将得到错误的估计系数，带来错误的结论。

图2 CST数据集的稳健异常值

由图2可以看出，基于MRCD方法的高维稳健多元线性回归模型估计识别出了6个坏的杠杆点，具有用于研究稳健回归的价值。

进一步，为了考察基于MRCD方法的高维稳健多元线性回归模型的预测效果，本文使用5折交叉验证方法来比较经典估计、MCD估计和MRCD估计在多元线性回归模型中的预测效果。比较的标准为RMSPE(root mean squared prediction error),计算公式为

将样本平均分为 5份，每次将其中一份作为测试集，用其他样本进行训练再对测试集样本进行预测，得到预测值。如此重复5次，即可将RMSPE的平均值作为测试结果。其中，每次测试的测试集分别记为T1,T2,…,T5。预测的结果如表3所示。

由表3可以看出，在该数据集中，经典估计受到了异常值的影响，导致RMSPE较大。而MRCD方法抵御了异常值的影响，取得了较小的RMSPE。这说明在数据集有异常值和坏的杠杆点时，基于MRCD方法的高维稳健多元线性回归模型表现最好，符合数值模拟中自变量和因变量都有异常值时基于MRCD方法的高维稳健多元线性回归模型估计效果更好的结果。因此可以认为，本文提出的方法在实际应用中是有效的。

表3 交叉验证预测结果

3 结语

随着经济社会和计算机技术的高速发展，各行各业每天都会产生大量结构复杂且形式多样的数据，包含异常值和高维的数据也越来越频繁地出现在实际应用中。这些数据蕴含着大量有意义的信息，因此挖掘高维数据的潜在信息有着非常重要的实际意义。但在处理高维数据时，现有方法往往导致传统分析方法的失效，多元线性回归模型估计的关键在于估计多变量总体的均值向量及协方差矩阵，但基于样本协方差矩阵的最小二乘法多元线性回归估计对于模型异常值极为敏感，且会在自变量或因变量的维数大于样本量时失效。基于传统方法在这两方面的局限性，本文提出了基于MRCD方法的多元线性回归模型。数值模拟的结果表明，当数据存在异常值，特别是自变量和因变量均存在异常值时，基于MRCD方法的多元线性回归模型估计在高维情况下是稳健的；实证分析的结果表明，基于MRCD方法的多元线性回归模型在实际情况下能更好地抵御异常值，且能得到更好的预测效果。在医学统计、生物统计、经济管理等领域中，常常遇到回归模型中有多个自变量和多个因变量的情况，而这些领域中又经常会遇到维数大于样本量的数据，可以考虑使用基于MRCD方法的高维稳健多元线性回归模型。