一种新的Pearson相关系数稳健估计法

A New Robust Estimation Method of Pearson Correlation Coefficient

摘要: 传统Pearson相关系数计算公式具有不稳健性，离群值的存在会导致计算结果与实际不符。针对此问题，文章给出了一种稳健估计方法。在模拟样本量分别为20、50、100、200，污染率分别为1%、5%、10%情形下，比较传统相关系数值与稳健相关系数值，发现：稳健相关系数公式正确率均显著高于传统相关系数。在实例分析中进一步验证了稳健相关系数的可行性和有效性。文章研究结论可用于含离群值变量的相关系数稳健估计。

Abstract: The traditional Pearson correlation coefficient calculation formula is not robust, and the existence of outliers will cause the calculation results to be inconsistent with reality. To solve this problem, this paper presents a robust estimation method. When the simulated sample size is 20, 50, 100 and 200 respectively, the pollution rate is 1%, 5% and 10% respectively, it is found that the accuracy of the robust correlation coefficient formula is significantly higher than that of the traditional correlation coefficient. The feasibility and effectiveness of a robust correlation coefficient are further verified in the example analysis. The conclusions of this paper can be used for robust estimation of correlation coefficients with outlier variables.

1. 引言

相关分析就是通过对大量数字资料的观察，消除偶然因素的影响，探求现象之间相关关系的密切程度和表现形式。研究现象之间相关关系的理论方法称为相关分析法。在经济系统中，各个经济变量常常存在密切的关系，例如经济增长与财政收入、人均收入与消费支出等。在这些关系中，有一些是严格的函数关系，这类关系可以用数学表达式表示出来。例如在一定价格的条件下，商品销售额与销售量的依存关系。还有一些是非确定的关系，一个变量的产生会影响其他变量，使其产生变化。其变化具有随机的特性，但是仍然遵循一定的规律。对此，函数关系很容易解决，而那些非确定的相关关系才是我们关心的问题。因为在经济系统中，绝大多数经济变量之间的关系是非严格的、不确定的。

相关分析以现象之间是否相关、相关的方向和密切程度等为主要研究内容，它不区分自变量与因变量，也不关心各变量的构成形式。其主要分析方法有绘制相关图、计算相关系数和检验相关系数。但在做相关分析时，相关系数的值受离群值的影响较大，无此点时变量相关，有此点时变量可能就不相关了。在实际工作中，解决这一问题常用的方法是：采用及时复核数据或采用统计方法对数据中的离群值进行诊断而后再修正，但该两种方法在技术处理上均较为繁琐，因此，相关系数的稳健估计法被提出。在无离群值时，稳健相关系数与传统相关系数计算结果几乎无差异；在有离群值时，稳健相关系数几乎不受离群值干扰。从现有文献来看，国外学者对相关系数的稳健估计方法研究相对较多，有从理论上探讨相关系数的稳健检验和稳健估计[1]-[5]，也有学者将稳健相关系数应用于不同领域，如在二元正态性稳健图形检验中的应用[6]、在解决合流分析问题的初始数据中应用[7]、在对偶犹豫模糊软集与复杂直觉模糊集决策中的应用[8]-[10]等。而国内学者的研究相对较少，王志坚、王斌会(2014) [11]研究了自相关系数的稳健性。

2. 原理与方法

2.1. 传统Pearson相关系数的不稳健性

在所有相关分析中，最简单的是两个变量之间的线性相关，它只涉及两个变量。而且一变量数值发生变动，另一变量数值也随之发生大致均等变动，其各点的分布在平面图上近似地表现为一直线，这种相关关系就称为直线相关，也叫线性关系。

线性相关分析是用相关系数来表示两个变量间相互的线性关系，并判断其密切程度的统计方法。Pearson相关系数用来反映两个变量的线性相关关系。总体相关系数用$\rho$ 表示。其计算公式为：

$\rho =\frac{\mathrm{cov}\left(x,y\right)}{\sqrt{\mathrm{var}\left(x\right)\mathrm{var}\left(y\right)}}=\frac{{\sigma }_{xy}}{\sqrt{{\sigma }_x^2{\sigma }_y^2}}$ (1)

(1)式中，${\sigma }_x^2$ 为变量x的总体方差，${\sigma }_y^2$ 为变量y的总体方差，${\sigma }_{xy}$ 为变量x与变量y的总体协方差，相关系数$\rho$ 没有单位，在−1至+1范围内波动，其绝对值愈接近1，两个变量间的直线相关性愈密切；愈接近0，相关性愈不密切。

在实际中，我们通常要计算样本的线性相关系数，也即Person相关系数，其计算公式为：

$r=\frac{{\displaystyle \sum \left(x-\overline{x}\right)\left(y-\overline{y}\right)}}{\sqrt{{\displaystyle \sum {\left(x-\overline{x}\right)}^2}\cdot {\displaystyle \sum {\left(y-\overline{y}\right)}^2}}}$ (2)

(2)式中，$\overline{x}$ 为变量x样本均值，$\overline{y}$ 为变量y样本均值。

另外一方面，在实际问题研究中，研究者所观测到或搜集到的数据经常会遇到一些突发及非可控外界信息的干扰，导致数据序列常会出现一些偏离主体数据较大的值，即离群值。如：新冠肺炎的爆发、金融危机、政治动荡、新政策的颁布及观测仪器失灵等。许多统计学家经过长期观察发现：时间序列数据中常含有0.5%至10%的离群值。

而通过分析公式(1)，(2)不难发现，变量x方差${\sigma }_x^2$ 和均值$\overline{x}$ 、变量y的方差${\sigma }_y^2$ 和均值$\overline{y}$ 均不稳健，对离群值敏感，进而导致公式(1)，(2)具有不稳健性。因此，当变量存在离群值时，若用Pearson相关系数公式计算相关系数，往往会误导研究者对变量间相关程度的判断。事实上，即使单个离群值，都会导致相关系数值变得面目全非。

不失一般性，在样本量分别为20、50、100、200的情形下，且每种样本量下假设只有1个离群值，来观察离群值对相关系数的影响，数据均由R语言随机产生。

Figure 1. The influence of one outlier on the correlation coefficient at sample sizes of 20 and 50

图1. 在样本量20、50时，1个离群值对相关系数的影响

图1、图2显示，在不同样本量下，单个离群值对Person相关系数及回归线均有较大影响，但影响程度不一样，随着离群值在样本量中的比重减少，其影响也在减弱。我们定义：

$�差率=\frac{不含离群值相�系�-含离群值相�系�}{不含离群值相�系�}\times 100\%$

表1数据显示：第一，在不同样本量下，离群值对相关系数均有影响，第二，离群值的影响大小取决于其在样本数据中的比率，比率越高误差率越大。

2.2. Pearson相关系数稳健估计

接下来，构建Pearson相关系数的稳健估计公式。对一组样本观测值$\left(X_1,Y_1\right),\cdots ,\left(X_n,Y_n\right)$ 来说，首先来看变量X的观测值$X_1,\cdots ,X_n$ ，我们记：

Figure 2. The influence of one outlier on the correlation coefficient at sample sizes of 100 and 200

图2. 在样本量100、200时，1个离群值对相关系数的影响

Table 1. The influence of outliers on correlation coefficient under different sample size

表1. 不同样本量下离群值对相关系数的影响

样本量	不含离群值相关系数	含离群值相关系数	误差率(%)
20	0.9905558	0.7315845	26.14
50	0.9898884	0.8544762	13.68
100	0.9894738	0.9165406	7.37
200	0.9817959	0.9419404	4.06

$m=\left[\left(1-\beta \right)n\right]$

$W_i=\left|X_i-M_x\right|$

其中，$0\le \beta \le 0.5$ ，$\left[\left(1-\beta \right)n\right]$ 表示取$\left(1-\beta \right)n$ 值的整数部分，$M_x$ 取值为：

$M_x=\underset{\forall \left(i,j\right)}{median}\left(\frac{X_i+X_j}{2}\right)$

$W_{\left(1\right)}\le \cdots \le W_{\left(n\right)}$ 表示对$W_i$ 进行从小到大排序。记：

${\widehat{w}}_x=W_{\left(m\right)}$

$i_1$ 为满足不等式$\left(X_i-M_x\right)/{\widehat{w}}_x<-1$ 的$X_i$ 的个数，$i_2$ 为满足不等式$\left(X_i-M_x\right)/{\widehat{w}}_x>1$ 的$X_i$ 的个数。则得到以下两式：

$S_x={\displaystyle \sum _{i=i_1+1}^{n-i_2}{X}_{\left(i\right)}}$

${\widehat{\varphi }}_x=\frac{{\widehat{w}}_x\left(i_2-i_1\right)+S_x}{n-i_1-i_2}$

记$U_i=\left(X_i-{\widehat{\varphi }}_x\right)/{\widehat{w}}_x$ ，令：

$\psi \left(x\right)=\max \left[-1,\min \left(1,x\right)\right]$

同理，对变量Y的观测值$Y_1,\cdots ,Y_n$ ，我们记：

$m=\left[\left(1-\beta \right)n\right]$

$W_i=\left|Y_i-M_y\right|$

其中，$0\le \beta \le 0.5$ ，$\left[\left(1-\beta \right)n\right]$ 表示取$\left(1-\beta \right)n$ 值的整数部分，$M_y$ 取值为：

$M_y=\underset{\forall \left(i,j\right)}{median}\left(\frac{Y_i+Y_j}{2}\right)$

$W_{\left(1\right)}\le \cdots \le W_{\left(n\right)}$ 表示对$W_i$ 进行从小到大排序。记：

${\widehat{w}}_y=W_{\left(m\right)}$

$i_1$ 为满足不等式$\left(Y_i-M_y\right)/{\widehat{w}}_y<-1$ 的$Y_i$ 的个数，$i_2$ 为满足不等式$\left(Y_i-M_y\right)/{\widehat{w}}_y>1$ 的$Y_i$ 的个数。则得到以下两式：

$S_y={\displaystyle \sum _{i=i_1+1}^{n-i_2}{X}_{\left(i\right)}}$

${\widehat{\varphi }}_y=\frac{{\widehat{w}}_y\left(i_2-i_1\right)+S_y}{n-i_1-i_2}$

记$V_i=\left(Y_i-{\widehat{\varphi }}_y\right)/{\widehat{w}}_y$ 。令：

$\psi \left(y\right)=\max \left[-1,\min \left(1,y\right)\right]$

通过以上变换，我们得到稳健Pearson相关系数公式如下：

$r^R=\frac{{\displaystyle \sum \psi \left(U_i\right)\psi \left(V_i\right)}}{\sqrt{\psi {\left(U_i\right)}^2\psi {\left(V_i\right)}^2}}$ (3)

(3)式中$r^R$ 表示稳健Pearson相关系数r，下文中，我们称由(2)式计算的相关系数为传统相关系数，(3)式计算的相关系数为稳健相关系数。

3. 模拟研究

为验证文章给出的稳健相关系数的有效性，模拟样本量分别为20、50、100、200，污染率分别为1%、5%、10%，目的是从小样本到大样本，从低污染率到高污染率来观察稳健相关系数效果。由于在样本量为20、污染率为1%时，离群值不足1个，在此我们假设只有1个。同时，为了比较在相同样本量、相同污染率时，不同离群值大小对相关系数的影响，我们将在样本量为20、污染率为5%下的1个离群值设置成与样本量为20、污染率为1%下的离群值不同大小。比较结果如图3~6，传统相关系数与稳健相关系数的误差率比较如表2所示。

图3~图6、表2数据显示：第一，不管在哪种样本量、污染率下，稳健相关系数正确率均高于传统相关系数；第二，在高污染率下，稳健相关系数与传统相关系数正确率均有所下降，但传统相关系数下降速度远高于稳健相关系数；第三，离群值大小、位置对传统相关系数值与稳健相关系数均有影响，但对传统相关系数值的影响远大于稳健相关系数。以上三点说明稳健相关系数对离群值具有一定的耐抗性、计算结果具有可行性和有效性。

Figure 3. The sample size is 20, with contaminate rates of 1%, 5%, and 10%, respectively

图3. 样本量为20，污染率分别为1%、5%、10%

Figure 4. The sample size is 50, with contaminate rates of 1%, 5%, and 10%, respectively

图4. 样本量为50，污染率分别为1%、5%、10%

Figure 5. The sample size is 100, with contaminate rates of 1%, 5%, and 10%, respectively

图5. 样本量为100，污染率分别为1%、5%、10%

Figure 6. The sample size is 200, with contaminate rates of 1%, 5%, and 10%, respectively

图6. 样本量为200，污染率分别为1%、5%、10%

Table 2. Comparison of traditional correlation coefficient and robust correlation coefficient in simulation experiment

表2. 模拟实验中的传统相关系数与稳健相关系数比较

样本量	污染率	传统相关系数	稳健相关系数	传统相关系数 误差率(%)	稳健相关系数 误差率(%)
20	1%	0.82	0.96	——	——
	5%	0.73	0.91	26.14	8.15
	10%	0.62	0.86	37.23	13.36
50	1%	0.85	0.96	13.68	2.92
	5%	0.72	0.92	26.93	6.74
	10%	0.69	0.91	29.72	8.37
100	1%	0.92	0.97	7.37	1.49
	5%	0.79	0.95	20.37	3.73
	10%	0.81	0.95	18.11	3.99
200	1%	0.91	0.97	7.40	1.24
	5%	0.84	0.96	14.25	1.80
	10%	0.77	0.89	21.70	9.39

4. 实证分析

本研究选取R语言自带数据集做为实证分析样本数据，用于验证稳健相关系数的可行性和有效性。该数据来源于Daudin等人[12]给出的牛奶成分数据集milk，共包含8个变量、86个观测值。8个变量分别是：X1为密度、X2为脂肪含量、X3为蛋白质含量、X4为酪蛋白含量、X5为工厂测量的干酪干物质(简称：工厂测奶酪干)、X6为实验室测量的奶酪干物质(简称：实验室测奶酪干)、X7为乳干物质、X8为奶酪制品。除第一个测量单位外，其余均为克/升。选择该数据集作为样本数据主要有以下两点考虑：第一数据集里含有离群值，与本研究目的相一致；第二，Todorov等人(1994)和Atkinson (1994)用该数据集研究主成分分析的稳健性及研究正向搜索识别多个异常值的算法均取得了较好的研究成果。

接下来通过多变量散点图来观察8个变量间的相关关系，如图7所示。由于本研究重点关注离群值对相关系数的影响，因此在8个变量中，选取离群值较明显的散点图，通过观察比较图7，选取X3蛋白质含量与X6实验室测奶酪干作为一组来检验稳健相关系数的有效性。作为对比，选取一组含有好的高杠杆点数据，来检验稳健相关系数的可行性，因为通常在研究中我们认为一个好的高杠杆点对相关系数及回归模型的建立是几乎无影响，其中的变量X5工厂测奶酪干与X6实验室测奶酪符合要求。为比较稳健相关系数效果，本研究分别在含离群值、高杠杆点与不含离群值、高杠杆点情形下采用传统相关系数公式与稳健相关系数公式计算相关系数值，结果如图8与表3所示。

Figure 7. Multivariate scatter plot of milk dataset

图7. Milk数据集的多元散点图

Figure 8. Scatter plot of dried cheese and protein content measured in the laboratory (left) scatter plot of dried cheese measured in the factory (right)

图8. 实验室测奶酪干与蛋白质含量散点图(左) 工厂测奶酪干散点图(右)

分别在剔除离群值、高杠杆点情形下与包含离群值、高杠杆点情形下计算传统相关和稳健相关系数，结果如表3所示。

Table 3. Comparison between robust correlation coefficient and traditional correlation coefficient in empirical analysis

表3. 实证分析中的稳健相关系数与传统相关系数比较

	剔除离群值/高杠杆点相关系数		含离群值/高杠杆点相关系数
	传统相关系数	稳健相关系数	传统相关系数	稳健相关系数
离群值	0.96	0.95	0.64	0.89
高杠杆点	0.96	0.95	0.98	0.96

表3数据显示，在不含离群值、高杠杆点情形下稳健相关系数与传统相关系数几乎无差异，说明稳健相关系数具有可行性。而在含离群值情形下传统相关系数值改变很大(33.43%)，而稳健相关系数变化轻微(7.21%)。同时注意到，在一个好的高杠杆点情形下，传统相关系数和稳健相关系数在剔除或包含高杠杆点下几乎无差异，进一步证实了稳健相关系数的可行性和有效性。

5. 结束语

Pearson相关系数常用于度量不同变量之间的相关程度，研究者在用相关系数时，希望相关系数能反映大多数数据的变化趋势。但离群值的存在，往往会使得计算出来的相关系数值与实际不符，从而误导研究者对实际问题的判断。为解决这一问题，文章在前人研究的基础上，给出了一种Pearsn相关系数的稳健估计方法，模拟与实证分析均表明，在不含离群值，稳健估计方法与传统估计方法几乎无差别，在含离群值，稳健估计方法比传统估计方法具有耐抗性，能减弱异常值的影响，从而极大地提升了相关系数值的准确性。另外，研究发现，离群值大小和位置对传统相关系数值的影响远大于稳健相关系数。文章的研究结论可应用于实际问题中含离群值的两变量及多变量相关程度的判断。

基金项目

课题来源：广东省教学科学“十三五”规划2020年度课题《后疫时期高职网络课堂教学质量监控评价“6L”体系构建与实践》，课题编号：2020GXJK565。