因子正交旋转（方差最大化）与因子得分计算

小象实战讲义 · 多元统计分析

在上一节中，我们学习了因子分析的基本模型和因子载荷矩阵的求解方法。然而，直接求解得到的因子载荷矩阵往往难以解释其经济含义。本节将介绍两种关键技术：因子正交旋转 和 因子得分计算。通过正交旋转，我们可以使因子载荷矩阵的结构更加简化，从而清晰地揭示每个公共因子与原始变量之间的关系。通过因子得分计算，我们可以将每个样本在不可观测的公共因子上的“得分”估算出来，为后续的分类、评价等应用提供基础。

💡 核心导读

• 因子旋转的必要性：原始因子载荷矩阵通常结构复杂，难以解释。通过正交旋转（特别是方差最大化旋转），可以使每个变量仅在一个公共因子上有较大载荷，从而简化结构、明确因子含义。
• 方差最大化旋转 (Varimax)：一种最常用的正交旋转方法，通过最大化因子载荷平方的方差，使载荷向 0 和 1 两极分化，实现因子结构的简化。
• 因子得分的估算：公共因子是不可观测的，但我们可以通过汤姆生 (Thomson) 回归法，将公共因子表达为原始变量的线性组合，从而估算每个样本在各个公共因子上的得分。
• 主成分分析与因子分析的核心区别：本节最后将系统梳理这两种降维方法在模型假设、求解目标、结果唯一性等方面的本质差异。

一、因子旋转：为什么需要旋转？

在因子分析中，我们最终的目标是解释公共因子 $F_1, F_2, \ldots, F_m$ 的实际意义。回顾因子模型： $$X = A F + \varepsilon$$ 其中 $A$ 是 $p \times m$ 的因子载荷矩阵。直接通过主因子法或主成分法求解得到的 $A$，其元素 $a_{ij}$ 的绝对值可能分布较为均匀，导致我们难以判断某个公共因子 $F_j$ 主要与哪些原始变量 $X_i$ 相关。

因子旋转 的目的，就是通过一个正交变换，使旋转后的因子载荷矩阵具有 “简单结构”：

• 每个变量 $X_i$ 仅在一个（或少数几个）公共因子上有较大的载荷（绝对值接近1）。
• 在其他公共因子上的载荷较小（绝对值接近0）。
• 这样，每个公共因子 $F_j$ 的经济意义就主要由那些在其上载荷较大的变量来定义，解释起来更加清晰。

1.1 因子旋转的数学原理

因子载荷矩阵 $A$ 不是唯一的。对于任意 $m \times m$ 的正交矩阵 $\Gamma$（满足 $\Gamma \Gamma^T = I_m$），令： $$A^* = A \Gamma, \quad F^* = \Gamma^T F$$ 则原模型可改写为： $$X = (A \Gamma) (\Gamma^T F) + \varepsilon = A^* F^* + \varepsilon$$ 容易验证，新的公共因子 $F^*$ 仍然满足 $E(F^*) = 0$，$Cov(F^*) = I_m$，且与特殊因子 $\varepsilon$ 不相关。因此，$A^*$ 和 $F^*$ 构成了原模型的另一组等价解。

因子旋转，就是寻找一个合适的正交矩阵 $\Gamma$，使得旋转后的载荷矩阵 $A^*$ 具有更理想的简单结构。

二、方差最大化旋转 (Varimax)

在众多正交旋转方法中，方差最大化旋转 (Varimax) 是最常用的一种。其核心思想是：让旋转后每个公共因子上的载荷平方的 方差最大化。方差越大，说明载荷值越向 0 和 1 两极分化，结构越简单。

2.1 Varimax 的目标函数

设旋转后的载荷矩阵为 $A^* = (a_{ij}^*)$，其维度为 $p \times m$。为了消除不同变量 $X_i$ 的共性方差 $h_i^2$ 大小不同的影响，我们先对载荷进行标准化： $$d_{ij} = \frac{a_{ij}^*}{h_i}, \quad i=1,\ldots,p; \ j=1,\ldots,m$$ 其中 $h_i^2 = \sum_{j=1}^m (a_{ij}^*)^2$ 是变量 $X_i$ 的共性方差（旋转前后保持不变）。

对于第 $j$ 个公共因子，定义其载荷平方的相对方差为： $$V_j = \frac{1}{p} \sum_{i=1}^p \left( d_{ij}^2 - \bar{d}_j \right)^2, \quad \text{其中} \quad \bar{d}_j = \frac{1}{p} \sum_{i=1}^p d_{ij}^2$$

Varimax 旋转的目标是选择正交矩阵 $\Gamma$，使得所有 $m$ 个公共因子的方差之和最大： $$V = \sum_{j=1}^m V_j \quad \rightarrow \quad \max_{\Gamma}$$

2.2 二维旋转的几何解释与迭代算法

当 $m=2$ 时，旋转矩阵 $\Gamma$ 可以表示为平面旋转矩阵： $$\Gamma = \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix}$$ 旋转角度 $\theta$ 可以通过求解 $V$ 对 $\theta$ 的导数并令其为零得到解析解。

对于 $m > 2$ 的一般情况，Varimax 采用 逐对旋转 的迭代算法：

1. 将所有 $m$ 个因子两两配对，共有 $C_m^2$ 对。
2. 对每一对因子 $(l, k)$，构造一个仅在 $(l, l), (l, k), (k, l), (k, k)$ 位置非零的 $m \times m$ 正交矩阵 $\Gamma_{lk}$（即一个二维旋转矩阵嵌入到 $m$ 维单位矩阵中）。
3. 对当前载荷矩阵 $A$ 右乘 $\Gamma_{lk}$，并计算旋转后 $V$ 的增量。选择使 $V$ 增加最大的旋转角度 $\theta$ 进行更新。
4. 完成所有配对的一次旋转称为 一轮循环。
5. 重复进行多轮循环，直到 $V$ 的增量小于预设的容差（例如 $10^{-6}$），此时认为算法收敛。

由于 $|a_{ij}^*| \le 1$，$V$ 是一个有上界的单调递增序列，因此算法必然收敛。

2.3 Python 实战：验证 Varimax 旋转

下面我们使用 Python 实现一个简化的 Varimax 旋转算法（此处实现的是未规范化 raw 版本，未做 2.1 节所述的 Kaiser 行归一化），并验证旋转前后因子方差贡献保持不变（正交旋转不改变总解释方差）。

"""
多元统计分析 6.2 因子正交旋转与因子得分计算
验证方差最大化正交旋转算法 (Varimax) 与旋转前后方差贡献保持
"""
import numpy as np

def varimax(Phi, gamma=1.0, q=500, tol=1e-6):
    """
    简化版 Varimax 正交旋转算法 (未规范化 raw 版本)
    
    参数:
        Phi: 原始因子载荷矩阵, 形状 (p, k)
        gamma: Crawford-Ferguson 族参数, 1.0=Varimax, 0.0=Quartimax
        q: 最大迭代次数
        tol: 收敛容差
    
    返回:
        Phi_rot: 旋转后的因子载荷矩阵
        Rot_mat: 旋转矩阵 (正交矩阵)
    """
    p, k = Phi.shape
    # 初始化旋转矩阵为单位矩阵
    R = np.eye(k)
    V_old = 0.0
    for i in range(q):
        # 当前旋转后的载荷矩阵
        Lambda = Phi @ R
        # 计算梯度方向 (简化版 Kaiser 公式)
        u, s, vh = np.linalg.svd(Phi.T @ (Lambda**3 - (gamma / p) * Lambda @ np.diag(np.sum(Lambda**2, axis=0))))
        # 更新旋转矩阵
        R = u @ vh
        Lambda_new = Phi @ R
        V_new = 0.25 * np.sum(np.sum(Lambda_new**4, axis=0) - (gamma / p) * np.sum(Lambda_new**2, axis=0) ** 2)
        # 检查收敛
        if abs(V_new - V_old) < tol:
            break
        V_old = V_new
    else:
        print(f"警告: 迭代 {q} 次后仍未收敛, 建议增大 q")
    # 返回旋转后的载荷矩阵和旋转矩阵
    return Phi @ R, R

# 示例：一个 5x2 的因子载荷矩阵 (5个变量，2个公共因子)
Lambda_raw = np.array([
    [0.85, 0.20],  # 变量1
    [0.88, 0.22],  # 变量2
    [0.82, 0.18],  # 变量3
    [0.30, 0.80],  # 变量4
    [0.25, 0.78]   # 变量5
])

print("原始因子载荷矩阵 Lambda_raw:")
print(np.round(Lambda_raw, 4))
print()

# 执行 Varimax 旋转
Lambda_rot, Rot_mat = varimax(Lambda_raw)

print("旋转后的因子载荷矩阵 Lambda_rot:")
print(np.round(Lambda_rot, 4))
print()

print("旋转矩阵 R (满足 R * R^T = I):")
print(np.round(Rot_mat, 4))
print()

# 计算旋转前后各因子的方差贡献 (载荷平方和)
var_raw = np.sum(Lambda_raw**2, axis=0)
var_rot = np.sum(Lambda_rot**2, axis=0)

print("旋转前因子方差贡献:", np.round(var_raw, 4).tolist())
print("Varimax 旋转后因子方差贡献:", np.round(var_rot, 4).tolist())
print("旋转前后总方差贡献对比:", round(float(np.sum(var_raw)), 4), "vs", round(float(np.sum(var_rot)), 4))
print("旋转矩阵正交性误差 (R * R^T 偏离度):", round(float(np.linalg.norm(Rot_mat @ Rot_mat.T - np.eye(2))), 8))

运行结果：

原始因子载荷矩阵 Lambda_raw:
[[0.85 0.2 ]
 [0.88 0.22]
 [0.82 0.18]
 [0.3  0.8 ]
 [0.25 0.78]]

旋转后的因子载荷矩阵 Lambda_rot:
[[0.8409 0.2355]
 [0.87   0.2567]
 [0.8117 0.2142]
 [0.2662 0.8119]
 [0.2171 0.7898]]

旋转矩阵 R (满足 R * R^T = I):
[[ 0.9991  0.0419]
 [-0.0419  0.9991]]

旋转前因子方差贡献: [2.3218, 1.3692]
Varimax 旋转后因子方差贡献: [2.2409, 1.4501]
旋转前后总方差贡献对比: 3.691 vs 3.691
旋转矩阵正交性误差 (R * R^T 偏离度): 0.0

结果分析：

1. 结构简化：旋转后，前三个变量在第一个因子上载荷较大（0.810.87），在第二个因子上载荷较小（0.210.26）；后两个变量在第二个因子上载荷较大（0.790.81），在第一个因子上载荷较小（0.220.27）。这体现了“简单结构”。
2. 方差贡献保持：旋转前后两个因子的方差贡献之和均为 3.691，验证了正交旋转不改变总解释方差，只是在因子间重新分配。
3. 正交性：旋转矩阵 $R$ 满足 $R R^T = I$，误差为 0（在浮点精度内）。

在实际应用中，我们可以直接使用 factor_analyzer 库中的 Rotator 类进行 Varimax 旋转，其算法更加稳健。

三、因子得分计算

因子分析的另一重要应用是计算每个样本在公共因子上的 得分。例如，在教育评估中，我们通过学生各科成绩（可观测变量）来估算其阅读能力、逻辑能力等潜在因子（公共因子）的得分。

3.1 问题的挑战

在主成分分析中，主成分 $Z$ 是原始变量 $X$ 的线性组合：$Z = W^T X$。因此，给定 $X$ 的观测值，可以直接计算 $Z$。

但在因子分析中，模型是 $X = A F + \varepsilon$，公共因子 $F$ 是“因”，原始变量 $X$ 是“果”。当 $m < p$ 时，载荷矩阵 $A$ 不可逆，无法直接由 $X$ 解出 $F$。

3.2 汤姆生 (Thomson) 回归法

汤姆生法的核心思想是：将每个公共因子 $F_j$ 对 $p$ 个原始变量 $X_1, \ldots, X_p$ 做线性回归。假设 $X$ 和 $F$ 都已标准化（均值为0，方差为1），则回归模型没有截距项： $$F_j = b_{j1} X_1 + b_{j2} X_2 + \cdots + b_{jp} X_p + e_j, \quad j=1,\ldots,m$$ 其中 $e_j$ 为回归残差。

记回归系数矩阵 $B = (b_{ji})_{m \times p}$，则上述 $m$ 个方程可写为： $$F = B X + e$$ 其中 $e$ 为回归残差向量；取条件期望（忽略残差）时，得分估计量为 $\hat{F} = B X = A^T R^{-1} X$。

3.3 回归系数的求解

根据因子载荷的统计意义，$a_{ij} = Cov(X_i, F_j) = Corr(X_i, F_j)$（标准化后）。将回归模型代入： $$\begin{aligned} a_{ij} &= Cov(X_i, F_j) \\ &= Cov\left(X_i, \sum_{k=1}^p b_{jk} X_k\right) \\ &= \sum_{k=1}^p b_{jk} Cov(X_i, X_k) \\ &= \sum_{k=1}^p b_{jk} r_{ik} \end{aligned}$$ 其中 $r_{ik}$ 是 $X_i$ 与 $X_k$ 的相关系数。

将上述关系写为矩阵形式： $$A = R B^T \quad \Rightarrow \quad B^T = R^{-1} A$$ 其中 $R$ 是 $X$ 的 $p \times p$ 相关矩阵。

因此，回归系数矩阵为： $$B = A^T R^{-1}$$

3.4 因子得分估算公式

最终，因子得分的估算公式为： $$\hat{F} = B X = A^T R^{-1} X$$ 对于一个新的样本，只要将其标准化后的观测值 $x$ 代入上式，即可得到其在 $m$ 个公共因子上的得分向量 $\hat{f}$。

3.5 Python 实战：因子得分计算

假设我们已经通过因子分析得到了载荷矩阵 $A$ 和相关矩阵 $R$，下面演示如何计算因子得分。

"""
因子得分计算示例 (汤姆生回归法)
"""
import numpy as np
from scipy import linalg

# 设定随机种子确保可复现
np.random.seed(42)

# 假设我们有 p=5 个变量，m=2 个公共因子
p, m = 5, 2

# 1. 生成一个随机的因子载荷矩阵 A (已旋转后的简化结构)
A = np.array([
    [0.9, 0.0],
    [0.8, 0.1],
    [0.7, 0.2],
    [0.1, 0.8],
    [0.0, 0.9]
])

# 2. 生成特殊因子方差 (个性方差)
D_epsilon = np.diag([0.19, 0.35, 0.47, 0.35, 0.19])  # 个性方差矩阵

# 3. 根据因子模型，生成总体相关矩阵 R = A A^T + D_epsilon
R = A @ A.T + D_epsilon

print("因子载荷矩阵 A:")
print(np.round(A, 4))
print("\n总体相关矩阵 R:")
print(np.round(R, 4))

# 4. 生成 n=100 个样本的标准化数据 X (均值为0，方差为1)
n = 100
# 生成多元正态数据，协方差矩阵为 R
mean = np.zeros(p)
X = np.random.multivariate_normal(mean, R, size=n)
# 标准化 (样本中心化，但这里总体均值为0，只需除以标准差)
X = (X - X.mean(axis=0)) / X.std(axis=0, ddof=1)

print(f"\n生成 {n} 个样本的标准化数据 X，形状: {X.shape}")

# 5. 计算汤姆生回归法的因子得分
# 公式: F_hat = X * (R^{-1} A)^T = X * B^T
R_inv = linalg.inv(R)
B = A.T @ R_inv  # 回归系数矩阵，形状 (m, p)
F_hat = X @ B.T   # 因子得分矩阵，形状 (n, m)

print("\n前5个样本的因子得分:")
print(np.round(F_hat[:5], 4))

# 6. 验证因子得分的性质
# a) 因子得分均值为0 (近似)
print(f"\n因子得分均值 (应接近0): {np.round(F_hat.mean(axis=0), 6)}")
# b) 正交因子模型下，Thomson 回归因子得分近似不相关（本例样本相关阵接近对角）；
#    注意其理论协方差为 A^T R^{-1} A（收缩估计），对角元小于1，并不等于单位矩阵
F_corr = np.corrcoef(F_hat, rowvar=False)
print("因子得分之间的样本相关矩阵:")
print(np.round(F_corr, 4))
# c) 因子得分与原始变量的相关系数应接近因子载荷 A
# 计算因子得分与原始变量的相关系数矩阵 (p x m)
corr_XF = np.zeros((p, m))
for i in range(p):
    for j in range(m):
        corr_XF[i, j] = np.corrcoef(X[:, i], F_hat[:, j])[0, 1]
print("\n原始变量与因子得分的样本相关系数矩阵 (应接近 A):")
print(np.round(corr_XF, 4))
print("\n与真实载荷矩阵 A 的差异 (绝对值平均):", np.round(np.mean(np.abs(corr_XF - A)), 4))

代码说明：

1. 我们首先生成一个具有简单结构的因子载荷矩阵 $A$。
2. 根据因子模型 $R = A A^T + D_\varepsilon$ 生成总体相关矩阵 $R$。
3. 从该多元正态分布中生成 $n=100$ 个样本。
4. 使用汤姆生公式 $\hat{F} = X (A^T R^{-1})^T$ 计算因子得分。
5. 验证因子得分的性质：均值为0、因子间近似不相关（Thomson 回归得分的协方差为 $A^T R^{-1} A$，不必等于单位阵，本例因载荷结构特殊而接近对角占优）、与原始变量的相关系数接近真实载荷 $A$。

四、主成分分析与因子分析的核心区别

学习至此，我们有必要系统梳理主成分分析 (PCA) 与因子分析 (FA) 的本质区别，这也是面试和实际应用中经常被问到的问题。

核心思想差异：

• PCA 是 数据的视角：从空间上转换观看数据的角度，寻找数据变异最大的方向。
• FA 是 变量的视角：从显在变量去“提炼”背后的潜在因子，解释变量间的相关关系。

📝 动手练一练

1. 因子旋转的意义
   假设你对一组社会经济指标进行了因子分析，得到了两个公共因子，但旋转前的载荷矩阵显示大多数变量在两个因子上都有中等大小的载荷。请问：

   • 这种情况下直接解释因子的经济含义会遇到什么困难？
   • 进行 Varimax 旋转后，你期望看到载荷矩阵发生怎样的变化？
   • 旋转会改变每个变量被公共因子解释的方差（共性方差）吗？为什么？

   参考答案：
   直接解释困难在于无法清晰判断每个因子主要代表哪些变量。Varimax 旋转后，期望看到载荷矩阵出现“简单结构”：每个变量仅在一个因子上有高载荷，在另一个因子上载荷接近0。旋转不会改变共性方差 $h_i^2$，因为正交旋转不改变载荷的平方和 $\sum_j a_{ij}^2$，而 $h_i^2$ 正是这个平方和。

2. 因子得分计算
   在因子模型中，已知： $$A = \begin{pmatrix} 0.8 & 0.2 \\ 0.7 & 0.3 \\ 0.1 & 0.9 \\ 0.3 & 0.8 \end{pmatrix}, \quad R = \begin{pmatrix} 1.0 & 0.6 & 0.3 & 0.4 \\ 0.6 & 1.0 & 0.4 & 0.5 \\ 0.3 & 0.4 & 1.0 & 0.7 \\ 0.4 & 0.5 & 0.7 & 1.0 \end{pmatrix}$$ 现有一样本，其标准化后的观测值为 $x = [1.2, -0.5, 0.8, 0.1]^T$。请使用汤姆生法计算该样本在两个公共因子上的得分。

   参考答案：

   import numpy as np
   A = np.array([[0.8,0.2],[0.7,0.3],[0.1,0.9],[0.3,0.8]])
   R = np.array([[1.0,0.6,0.3,0.4],
                 [0.6,1.0,0.4,0.5],
                 [0.3,0.4,1.0,0.7],
                 [0.4,0.5,0.7,1.0]])
   x = np.array([1.2, -0.5, 0.8, 0.1])
   R_inv = np.linalg.inv(R)
   B = A.T @ R_inv  # 回归系数矩阵
   f_hat = B @ x    # 因子得分
   print("因子得分:", np.round(f_hat, 4))

   计算结果约为 $[0.30, 0.53]$。

本章小结

本节深入探讨了因子分析中两个至关重要的应用环节：因子正交旋转 与 因子得分计算。

要点回顾：

1. 因子旋转 是为了获得具有“简单结构”的载荷矩阵，使每个公共因子的经济意义更加明确。Varimax 旋转通过最大化载荷平方的方差，实现载荷向 0 和 1 两极分化。
2. 因子得分 是将不可观测的公共因子表达为原始变量的线性组合，从而估算每个样本在潜在因子上的“得分”。汤姆生回归法是最常用的方法，公式为 $\hat{F} = A^T R^{-1} X$。
3. 主成分分析与因子分析 虽然都是降维技术，但在模型假设、求解目标、结果唯一性等方面存在本质区别。PCA 侧重数据变异的最大化，FA 侧重变量间相关结构的解释。

行动清单：

1. 动手实践：使用 factor_analyzer 库的 Rotator 类对你自己的数据集进行 Varimax 旋转，观察旋转前后载荷矩阵的结构变化。
2. 代码复现：根据汤姆生公式，编写一个函数 compute_factor_scores(X, A)，输入标准化数据 $X$ 和载荷矩阵 $A$，返回因子得分矩阵。
3. 对比分析：对同一份数据分别进行 PCA 和 FA，比较两者的成分/因子载荷矩阵、得分以及解释的方差，深入理解两者的异同。

通过本节学习，你已掌握了因子分析从模型求解到结果解释、再到实际应用（得分计算）的完整流程。下一节我们将结合真实数据集，使用 Python 完整演示因子分析的全过程。

— 小象教研组