偏相关

在概率论和统计学中，偏相关（Partial correlation）度量了两个随机变量在移除了其他控制变量（Confounding variables）的影响后的相关程度。

在确定两个感兴趣变量之间的数值关系时，如果存在与这两个变量都相关的混杂变量，直接计算它们的相关系数可能会产生误导性结果（伪相关）。通过计算偏相关系数控制混杂变量，可以避免这种误导性信息。这正是多元线性回归中包含其他右侧变量的动机；虽然多元回归给出了效应值的无偏结果，但它并不直接给出一个数值来度量两个感兴趣变量之间关系的强度。

例如，给定关于个人消费、收入和财富的经济数据，如果希望分析消费与收入之间的关系，而在计算相关系数时未能控制“财富”这一变量，可能会得到误导性的结果。因为收入可能与财富相关，而财富又与消费相关；直接测量的消费与收入之间的相关性可能会被这些其他相关性所“污染”。使用偏相关可以避免这个问题。

与普通的相关系数一样，偏相关系数的取值范围也是从 -1 到 1。

- 值 -1 表示在控制某些变量后，存在完全的负相关（即一个变量值越高，另一个变量值越低）。
- 值 1 表示完全的正线性关系。
- 值 0 表示没有线性关系。

如果随机变量的联合分布是多元正态分布、其他椭圆分布、多元超几何分布、多元负超几何分布、多项分布或狄利克雷分布，则偏相关与条件相关一致，但一般情况下两者并不相同。

## 正式定义

形式上，给定一组 $n$ 个控制变量 $\mathbf{Z} = \{Z_1, Z_2, \dots, Z_n\}$，变量 $X$ 和 $Y$ 之间的偏相关系数 $\rho_{XY\cdot\mathbf{Z}}$ 定义为：$X$ 关于 $\mathbf{Z}$ 的线性回归的残差 $e_X$ 与 $Y$ 关于 $\mathbf{Z}$ 的线性回归的残差 $e_Y$ 之间的相关系数。

一阶偏相关（即 $n=1$ 时）是相关系数与可移除相关系数之积的差，除以可移除相关系数的相异系数（coefficient of alienation）之积。

## 计算方法

### 使用线性回归

计算样本偏相关的一种简单方法是求解两个相关的线性回归问题，并计算残差之间的相关性。

设 $X$ 和 $Y$ 为取实数值的随机变量，$\mathbf{Z}$ 为 $n$ 维向量值随机变量。设 $x_i, y_i, \mathbf{z}_i$ 表示从实随机变量 $X, Y, \mathbf{Z}$ 的联合分布中抽取的 $N$ 个独立同分布（i.i.d.）观测值中的第 $i$ 个（其中 $\mathbf{z}_i$ 增加了一个 1 以允许回归中的常数项）。

求解线性回归问题相当于找到 $(n+1)$ 维回归系数向量 $\mathbf{w}_X^*$ 和 $\mathbf{w}_Y^*$，使得：

$$\mathbf{w}_X^* = \arg\min_{\mathbf{w}} \left\{ \sum_{i=1}^N (x_i - \langle\mathbf{w}, \mathbf{z}_i \rangle)^2 \right\}$$

$$\mathbf{w}_Y^* = \arg\min_{\mathbf{w}} \left\{ \sum_{i=1}^N (y_i - \langle\mathbf{w}, \mathbf{z}_i \rangle)^2 \right\}$$

其中 $N$ 是观测数量，$\langle\mathbf{w}, \mathbf{z}_i \rangle$ 是向量 $\mathbf{w}$ 和 $\mathbf{z}_i$ 之间的标量积。

残差为：

$$e_{X,i} = x_i - \langle\mathbf{w}_X^*,\mathbf{z}_i \rangle$$

$$e_{Y,i} = y_i - \langle\mathbf{w}_Y^*,\mathbf{z}_i \rangle$$

样本偏相关系数即为这些新生成的残差值之间的样本相关系数：

$$\begin{align} \hat{\rho}_{XY\cdot\mathbf{Z}}&=\frac{N\sum_{i=1}^N e_{X,i}e_{Y,i}-\sum_{i=1}^N e_{X,i}\sum_{i=1}^N e_{Y,i}} {\sqrt{N\sum_{i=1}^N e_{X,i}^2-\left(\sum_{i=1}^N e_{X,i}\right)^2}~\sqrt{N\sum_{i=1}^N e_{Y,i}^2-\left(\sum_{i=1}^N e_{Y,i}\right)^2}}\\ &=\frac{N\sum_{i=1}^N e_{X,i}e_{Y,i}} {\sqrt{N\sum_{i=1}^N e_{X,i}^2}~\sqrt{N\sum_{i=1}^N e_{Y,i}^2}}. \end{align}$$

在第一个表达式中，减号后的三项都等于 0，因为它们包含普通最小二乘法回归的残差和（必定为 0）。

#### 示例

考虑三个变量 $X, Y, Z$ 的以下数据：

X｜Y｜Z

2｜1｜0

4｜2｜0

15｜3｜1

20｜4｜1

计算变量 $X$ 和 $Y$ 之间的皮尔逊相关系数约为 0.970；而使用上述公式计算 $X$ 和 $Y$ 之间的偏相关系数（控制 $Z$）则为 0.919。以下是使用 R语言 进行的计算：

```R
> x <- c(2, 4, 15, 20)
> y <- c(1, 2,  3,  4)
> z <- c(0, 0,  1,  1)

# 将 x 对 z 进行回归并计算残差
> res_x <- lm(x ~ z)$residuals

# 将 y 对 z 进行回归并计算残差
> res_y <- lm(y ~ z)$residuals

# 计算残差的相关性
> cor(res_x, res_y)
# [1] 0.919145

# 显示这与 x 和 y 之间的直接相关性不同
> cor(x, y)
# [1] 0.9695016
```

### 使用递推公式

求解线性回归问题在计算上可能非常昂贵。实际上，$n$ 阶偏相关（即 $|\mathbf{Z}| = n$）可以很容易地从三个 $(n-1)$ 阶偏相关计算得出。零阶偏相关 $\rho_{XY\cdot\varnothing}$ 定义为普通的相关系数 $\rho_{XY}$。

对于任意 $Z_0 \in \mathbf{Z}$，有：

$$\rho_{XY\cdot \mathbf{Z} } = <pre> \frac{\rho_{XY\cdot\mathbf{Z}\setminus\{Z_0\}} - \rho_{XZ_0\cdot\mathbf{Z}\setminus\{Z_0\}}\rho_{Z_0Y\cdot\mathbf{Z}\setminus\{Z_0\}}} {\sqrt{1-\rho_{XZ_0\cdot\mathbf{Z}\setminus\{Z_0\}}^2} \sqrt{1-\rho_{Z_0Y\cdot\mathbf{Z}\setminus\{Z_0\}}^2}}$$

如果单纯地将此计算实现为递归算法，会产生指数级的时间复杂度。然而，由于该计算具有重叠子问题性质，使用动态规划或简单地缓存递归调用的结果可以产生 $\mathcal{O}(n^3)$ 的复杂度。

特别地，当 $Z$ 是单个变量时，公式简化为：

$$\rho_{XY\cdot Z}= <pre> \frac{\rho_{XY}-\rho_{XZ}\rho_{ZY}} {\sqrt{1-\rho_{XZ}^2}\sqrt{1-\rho_{ZY}^2}}$$

### 使用矩阵求逆

偏相关也可以用联合精度矩阵（Precision matrix）来表示。考虑一组随机变量 $\mathbf{V} = \{X_1,\dots, X_n\}$。我们想要计算在给定其他所有变量 $\mathbf{V} \setminus \{X_i,X_j\}$ 的情况下，变量 $X_i$ 和 $X_j$ 之间的偏相关。

假设（联合/完整）协方差矩阵 $\Sigma = (\sigma_{ij})$ 是正定矩阵，因此是可逆矩阵。定义精度矩阵为 $\Omega = (p_{ij}) = \Sigma^{-1}$，则：

$$\rho_{X_i X_j\cdot \mathbf{V} \setminus \{X_i,X_j\}} = - \frac{p_{ij}}{\sqrt{p_{ii}p_{jj}}}$$

计算此式需要求协方差矩阵 $\Sigma$ 的逆 $\Sigma^{-1}$，时间复杂度为 $\mathcal{O}(n^3)$。值得注意的是，只需一次矩阵求逆即可获得 $\mathbf{V}$ 中所有变量对之间的偏相关系数。

## 解释

### 几何解释

<image>

设三个变量 $X, Y, Z$（其中 $Z$ 是“控制”或“额外变量”）是从 $n$ 个变量 $\mathbf{V}$ 的联合概率分布中选出的。进一步设 $\mathbf{v}_i, 1 \le i \le N$ 是从 $\mathbf{V}$ 的联合概率分布中抽取的 $N$ 个 $n$ 维独立同分布观测值。

考虑由观测值形成的 $N$ 维向量 $\mathbf{x}$、$\mathbf{y}$ 和 $\mathbf{z}$。

可以证明，来自 $X$ 对 $\mathbf{Z}$ 的线性回归的残差 $e_{X,i}$，如果视为一个 $N$ 维向量 $\mathbf{e}_X$（在附图中记为 $\mathbf{r}_X$），与 $\mathbf{Z}$ 生成的向量 $\mathbf{z}$ 的标量积为零。这意味着残差向量位于垂直于 $\mathbf{z}$ 的 $(N-1)$ 维超平面 $S_{\mathbf{z}}$ 上。

同理也适用于 $Y$ 的残差向量 $\mathbf{e}_Y$。所需的偏相关就是向量 $\mathbf{x}$ 和 $\mathbf{y}$ 在垂直于 $\mathbf{z}$ 的超平面上的投影 $\mathbf{e}_X$ 和 $\mathbf{e}_Y$ 之间夹角 $\phi$ 的余弦。

### 作为条件独立性测试

假设所有涉及的变量都是多元高斯分布的，偏相关 $\rho_{XY\cdot\mathbf{Z}}$ 为零当且仅当 $X$ 在给定 $\mathbf{Z}$ 的条件下与 $Y$ 条件独立。这一性质在一般情况下不成立。

为了检验样本偏相关 $\hat{\rho}_{XY\cdot\mathbf{Z}}$ 是否意味着总体偏相关不为 0，可以使用偏相关的费希尔变换（Fisher's z-transform）：

$$z(\hat{\rho}_{XY\cdot\mathbf{Z}}) = \frac{1}{2} \ln\left(\frac{1+\hat{\rho}_{XY\cdot\mathbf{Z}}}{1-\hat{\rho}_{XY\cdot\mathbf{Z}}}\right)$$

零假设是 $H_0: \rho_{XY\cdot\mathbf{Z}} = 0$，针对双尾备择假设 $H_A: \rho_{XY\cdot\mathbf{Z}} \neq 0$。如果下式成立，则拒绝 $H_0$：

$$\sqrt{N - |\mathbf{Z}| - 3}\cdot |z(\hat{\rho}_{XY\cdot\mathbf{Z}})| > \Phi^{-1}(1-\alpha/2)$$

其中 $\Phi$ 是具有零均值和单位标准差的高斯分布的累积分布函数，$\alpha$ 是显著性水平，$N$ 是样本大小。

## 半偏相关

半偏相关（Semipartial correlation，或 part correlation）统计量与偏相关统计量类似；两者都比较两个变量在控制某些因素后的变异。然而，计算半偏相关时，第三个变量仅针对 $X$ 或 $Y$ 中的一个保持常数，而不是同时针对两者。偏相关比较的是一个变量的唯一变异与另一个变量的唯一变异，而半偏相关比较的是一个变量的唯一变异（移除了与 $Z$ 相关的变异）与另一个变量的未过滤变异。

半偏相关可以被视为更具有实际相关性，因为它相对于因变量（响应变量）的总变异进行了缩放。相反，它在理论上的用处较小，因为它对自变量的唯一贡献的作用不那么精确。

$X$ 与 $Y$ 的半偏相关的绝对值总是小于或等于 $X$ 与 $Y$ 的偏相关。原因如下：假设从 $X$ 中移除了 $X$ 与 $Z$ 的相关性，得到残差向量 $e_x$。在计算半偏相关时，$Y$ 仍然包含唯一方差和由于其与 $Z$ 关联而产生的方差。但是 $e_x$ 由于与 $Z$ 不相关，只能解释 $Y$ 方差中唯一的一部分，而不能解释与 $Z$ 相关的那部分。相比之下，对于偏相关，只有 $e_y$（$Y$ 方差中与 $Z$ 无关的部分）需要被解释，因此原本 $e_x$ 无法解释的那类方差在 $e_y$ 中已经不存在了。

## 时间序列分析中的应用

在时间序列分析中，时间序列的偏自相关函数（Partial autocorrelation function, PACF，有时也称为“偏相关函数”）定义为滞后 $h$ 时的偏相关：

$$\varphi(h)= \rho_{X_0X_h\,\cdot\, \{X_1,\,\dots\,,X_{h-1} \}}$$

该函数用于确定自回归（AR）过程的适当滞后长度。

## 带有收缩估计的偏相关

当样本量小于变量数量（即高维数据设置）时，估计偏相关具有挑战性。在这种情况下，样本协方差 $\hat{\Sigma}$ 是病态的，求其逆矩阵 $\hat{\Omega}$ 会出现问题。

收缩估计（Shrinkage estimation）方法可以改善 $\hat{\Sigma}$ 或 $\hat{\Omega}$ 的估计，从而产生更可靠的偏相关估计。一个例子是 Ledoit-Wolf 收缩估计量：

$$\hat{\Sigma}^{[\lambda]} = \lambda T + (1 - \lambda) \Sigma$$

其中 $\hat{\Sigma}$ 是样本协方差矩阵，$T$ 是目标矩阵（例如对角矩阵），收缩强度 $\lambda\in (0,1)$。

Ledoit-Wolf 收缩下的偏相关估计为：

$$\hat{P}_{ij}^{[\lambda]} = \frac{\hat{\Omega}_{ij}^{[\lambda]}}{\sqrt{\hat{\Omega}_{ii}^{[\lambda]} \hat{\Omega}_{jj}^{[\lambda]} }} <pre>$$

其中 $\hat{\Omega}_{ij}^{[\lambda]}$ 是 $\hat{\Sigma}_{ij}^{[\lambda]}$ 的逆矩阵元素。这种方法应用于金融和基因组学等多个领域。

## 参见

- 线性回归
- 条件独立
- 复相关 (Multiple correlation)
- 偏自相关函数 (Partial autocorrelation function)