2.14 高维统计与因果推断

前十三节的方法建立在一个共同的前提上：样本量明显大于参数个数，且数据生成机制满足模型假设。现代数据把这两条前提都推到边界。变量数常常达到数千甚至数百万（文本的词表、图像的像素、传感器阵列的通道），而样本量受限于采集成本；同时大量数据来自观测而非实验，变量之间的关联被共同原因与选择机制污染，回归系数不再等于因果效应。

本节讨论这两类问题。前一部分是高维统计：在 $p$ 接近或超过 $n$ 的情形下，最小二乘无解或极不稳定，需要引入稀疏性或低秩结构作为额外约束，正则化方法、随机矩阵理论与高维多重检验构成了这一领域的工具。后一部分是因果推断：从潜在结果框架与因果图两条路径给出因果效应的定义，说明识别所需的条件，介绍随机化实验之外的若干识别策略，并讨论敏感性分析与因果发现。两个部分处理的问题不同，却共享同一套概率语言，也都以本章前面的内容为工具。

2.14.1 高维数据与维数灾难

高维的含义

高维数据（high-dimensional data）通常指变量数 $p$ 与样本量 $n$ 同阶或更大。经典统计的渐近理论假设 $p$ 固定、$n\to \mathrm{\infty }$，高维情形下这一框架失效，需要 $p$ 随 $n$ 一起增长的渐近理论（$p/n\to \gamma$ 的比值渐近）。

高维带来的困难不止是计算量。维数灾难（curse of dimensionality）描述的是一组几何直觉在高维空间中的失效：

距离集中。在高维空间中，任意两点之间的距离趋于接近。若数据均匀分布在 $p$ 维单位立方体中，最远点距离与最近点距离之比随 $p$ 增大而趋于 $1$。基于距离的方法（最近邻、核方法、聚类）因此失去判别力：最近邻不再比任意其他点更近。

样本稀疏。要让 $p$ 维单位立方体中的样本覆盖整个空间，样本量需随 $p$ 指数增长。$p=10$ 时需要 ${10}^{10}$ 个样本才能达到每个维度上 $10$ 个点的密度，实际上不可行。空间中绝大部分区域没有观测，任何局部方法（局部回归、核密度估计）都缺乏支撑。

体积集中。高维球体的体积集中在表面附近的薄层中，多数点靠近边界而非中心。这一现象说明均匀分布在高维空间中的行为与低维直觉相反，基于体积的抽样与数值积分方法需要相应调整。

最小二乘的失效

高维情形对回归的直接影响是设计矩阵的秩不足。$p>n$ 时 ${\mathbf{X}}^T\mathbf{X}$ 是 $p\times p$ 的矩阵，秩最多为 $n$，因此不可逆，最小二乘解不唯一：存在无穷多组系数使残差平方和为零（可以完全拟合任何数据）。即使 $p<n$ 但两者接近，${\mathbf{X}}^T\mathbf{X}$ 也接近奇异，条件数很大，系数的方差极大。

这一局面有一个重要的推论：在高维情形下过拟合是默认状态。任何模型都可以在训练集上达到零误差，样本外表现完全由额外约束决定。约束的来源有两种：稀疏性假设（少数变量起作用）与低秩假设（数据位于低维子空间附近）。本节讨论的方法都建立在其中之一上。

import numpy as np

# 高维下的最小二乘：p > n 时可以完全拟合训练集，但样本外预测完全没有价值
rng = np.random.default_rng(16)
n, p = 30, 100
X = rng.normal(size=(n, p))
y = rng.normal(size=n)

beta_ls = np.linalg.lstsq(X, y, rcond=None)[0]
resid = y - X @ beta_ls
print(f"样本量 {n}，变量数 {p}")
print(f"训练残差平方和 = {(resid ** 2).sum():.6f}")     # 接近零，完全拟合

# 样本外数据上的表现
X_new = rng.normal(size=(50, p))
y_new = rng.normal(size=50)
pred = X_new @ beta_ls
print(f"样本外均方误差 = {((y_new - pred) ** 2).mean():.3f}")   # 远大于 1，比预测均值还差

# 各列相关矩阵的最大最小特征值，反映条件数
XtX = X.T @ X
eigs = np.linalg.eigvalsh(XtX)
print(f"X^T X 的秩 = {np.linalg.matrix_rank(XtX)}，小于 p = {p}")
print(f"条件数 = {np.sqrt(eigs[-1] / max(eigs[0], 1e-12)):.2e}")

2.14.2 稀疏性与正则化方法

稀疏假设

高维问题需要额外约束，最常用的约束是稀疏性（sparsity）：真正影响响应变量的只有少数几个特征，其余特征的系数为零或接近零。稀疏假设在文本分类（一篇文章的主题由少数词决定）、推荐系统（少数行为特征决定偏好）、传感器网络（少数位置反映关键状态）中都有实际依据。

稀疏性把估计问题的有效参数个数从 $p$ 降到 $s$（非零系数个数），使 $s<n$ 的经典理论重新适用。2.10 节的 Lasso 与压缩感知是这一思路的具体实现，本节补充它们的理论性质与算法。

Lasso 的选择性质

Lasso 估计满足

$$\hat{\mathit{\beta }}=\arg \underset{\mathit{\beta }}{min}\{\frac{1}{2n}\parallel \mathbf{y}-\mathbf{X}\mathit{\beta }{\parallel }^2+\lambda \parallel \mathit{\beta }{\parallel }_1\}$$

它的解满足最优性条件：对非零系数，梯度部分恰好抵消惩罚的次梯度；对零系数，梯度的绝对值不超过 $\lambda$。这一条件称为 KKT 条件（Karush-Kuhn-Tucker conditions），是坐标下降算法每步更新的依据。

理论结果需要设计矩阵满足一定条件，机制是防止真实非零变量与噪声变量之间出现较强的相关。不可表示条件（irrepresentable condition）要求非零变量与零变量的相关性不能太强，它是 Lasso 能够一致地识别支撑集（非零变量的位置）的关键前提。在这一条件下，选择合适的 $\lambda$ 可以使 Lasso 以高概率选出正确的变量集合，并达到接近已知支撑集时最优估计的收敛速度，这一性质称为 oracle 性质。

这些条件在相关系数较高的变量组上容易不成立，此时 Lasso 会在组内随机保留一个变量，选出的支撑集不稳定。弹性网通过加入 $L_2$ 惩罚缓解这一问题，代价是系数不再精确为零以外还引入了额外的收缩。

压缩感知

压缩感知（compressed sensing）处理的是从远少于参数个数的观测中恢复稀疏信号的欠定线性系统情形：给定 $\mathbf{y}=\mathbf{A}\mathbf{x}$，$\mathbf{A}$ 是 $m\times p$ 且 $m\ll p$，在 $\mathbf{x}$ 稀疏的假设下恢复它。这一问题在信号处理中称为采样与重建，在统计中就是高维稀疏回归。

恢复的可行性由受限等距性质（restricted isometry property，RIP）刻画：矩阵 $\mathbf{A}$ 在所有 $s$ 稀疏向量上近似保持长度，

$$(1-{\delta }_s)\parallel \mathbf{x}{\parallel }_2^2\le \parallel \mathbf{A}\mathbf{x}{\parallel }_2^2\le (1+{\delta }_s)\parallel \mathbf{x}{\parallel }_2^2$$

${\delta }_s$ 足够小时，$L_1$ 最小化（基追踪，basis pursuit）可以精确恢复 $s$ 稀疏信号，所需观测数约为 $m=O(s\log p)$。随机矩阵（高斯、伯努利、部分傅里叶）以高概率满足 RIP，这一结论是压缩感知实用性的理论基础。

算法上除了 $L_1$ 最小化，还有贪心类方法与迭代阈值类方法。正交匹配追踪（orthogonal matching pursuit，OMP）每次选出与当前残差最相关的列加入支撑集，再在支撑集上做最小二乘，迭代 $s$ 次。迭代硬阈值（iterative hard thresholding，IHT）在梯度步之后保留绝对值最大的 $s$ 个分量，直接控制稀疏度而非通过惩罚。近似消息传递（approximate message passing，AMP）利用矩阵的随机性给出闭式的迭代更新，在高维极限下有精确的性能预测，收敛速度通常快于前两者。

矩阵补全与低秩结构

另一类高维结构是低秩（low rank）：数据矩阵的秩远小于其维度。矩阵补全（matrix completion）在只观测到部分元素的情况下恢复低秩矩阵，典型的例子是用户对物品的评分矩阵（评分项稀疏，潜在维度少）。恢复的准则用核范数（矩阵奇异值之和）替代秩，作为秩的凸松弛：

$$\underset{\mathbf{M}}{min}\parallel \mathbf{M}{\parallel }_{\ast }\text{s.t.}{M}_{ij}=Y_{ij},(i,j)\in \mathrm{\Omega }$$

核范数是 $L_1$ 范数在奇异值上的版本，两者在稀疏与低秩之间形成完整的类比：稀疏对应 $L_1$，低秩对应核范数。矩阵分解（matrix factorization）用两个低秩因子 $\mathbf{M}=\mathbf{U}{\mathbf{V}}^T$ 直接参数化，用梯度方法优化，是推荐系统中的常用做法。

2.14.3 高维降维与协方差估计

主成分分析与稀疏主成分

主成分分析（principal component analysis，PCA）寻找方差最大的线性方向，计算归结为协方差矩阵的特征分解。线性代数部分已给出几何解释，此处补充高维下的两点变化。第一，$p>n$ 时样本协方差矩阵的秩最多为 $n-1$，非零特征值最多 $n-1$ 个，其余全为零，降维空间被样本量限制。第二，主成分方向是所有变量的线性组合，$p$ 很大时每个主成分都含有全部变量，解释性差。

稀疏主成分分析（sparse PCA）在载荷向量上加稀疏约束，使每个主成分只涉及少数变量。求解方式包括对载荷加 $L_1$ 惩罚后做迭代、或把问题写成带约束的优化用前向选择求解。稀疏载荷的价值在于解释：当主成分只涉及十来个变量时，可以给出这一方向的业务含义。代价是各个主成分方向不再正交（若同时施加稀疏约束），方差解释比例也随之下降。

其他降维方法包括因子分析（factor analysis，2.6 节）、典型相关分析与稀疏典型相关分析（2.6 节），它们分别针对协方差结构、两组变量关联与稀疏载荷的组合。

高维协方差估计

$p>n$ 时样本协方差矩阵 $\mathbf{S}$ 奇异，其逆不存在，依赖协方差逆的方法（判别分析、图模型、马氏距离）无法直接使用。即便 $p<n$，$p$ 较大时 $\mathbf{S}$ 的估计误差也很大，特征值被系统性压缩或膨胀（下面第 2.14.4 节的随机矩阵理论给出量化描述）。

应对方式是引入结构假设。对角线假设直接用对角矩阵或其缩放形式。因子结构用低秩加对角的形式（2.6 节的因子模型）。收缩估计把样本协方差向结构简单的目标（如对角矩阵）收缩：

$$\hat{\mathbf{\Sigma }}=(1-\delta )\mathbf{S}+\delta \mathbf{T}$$

$\delta$ 的选择以最小化期望误差为目标，Ledoit-Wolf 收缩给出闭式解。图形结构的假设见下。

图模型选择

高斯图模型（Gaussian graphical model）用精度矩阵 $\mathbf{\Theta }={\mathbf{\Sigma }}^{-1}$ 的非零元素表示变量之间的条件依赖（2.6 节）。估计精度矩阵的稀疏结构有两条主要路线。

图 Lasso（graphical Lasso）对精度矩阵刻画 $L_1$ 惩罚下的最大似然：

$$\hat{\mathbf{\Theta }}=\arg \underset{\mathbf{\Theta }\succ 0}{max}\{\log det\mathbf{\Theta }-\mathrm{tr}(\mathbf{S}\mathbf{\Theta })-\lambda \parallel \mathbf{\Theta }{\parallel }_1\}$$

求解用块坐标下降，每次更新一行一列。邻域选择（neighborhood selection）为每个变量做一次 Lasso 回归（用其余变量预测它），由回归系数的非零模式构造邻居集合，再通过对称化与一致性规则确定边。两种方法都利用稀疏假设把高维协方差估计分解为一组低维回归问题，计算复杂度显著低于直接优化。

import numpy as np

def graph_lasso_blockwise(S, lam, max_iter=100, tol=1e-5):
    """图 Lasso 的块坐标下降（简化版，用于展示更新逻辑）"""
    p = S.shape[0]
    W = S + lam * np.eye(p)
    for _ in range(max_iter):
        W_old = W.copy()
        for j in range(p):
            idx = np.arange(p) != j
            w11 = W[np.ix_(idx, idx)]
            s12 = S[idx, j]
            # 对子问题做一次近端梯度更新（此处用软阈值近似）
            beta = np.linalg.solve(w11, s12)
            beta = np.sign(beta) * np.maximum(np.abs(beta) - lam, 0)
            w12 = w11 @ beta
            W[np.ix_(idx, [j])] = w12.reshape(-1, 1)
            W[j, idx] = w12
        if np.max(np.abs(W - W_old)) < tol:
            break
    return W

rng = np.random.default_rng(28)
# 构造稀疏精度矩阵：三个变量构成一个团，其余独立
p = 6
Theta = np.eye(p)
Theta[0, 1] = Theta[1, 0] = -0.4
Theta[1, 2] = Theta[2, 1] = -0.4
Theta[3, 4] = Theta[4, 3] = -0.4
Sigma_true = np.linalg.inv(Theta)

n = 2000
X = rng.multivariate_normal(np.zeros(p), Sigma_true, size=n)
S = np.cov(X.T)
W = graph_lasso_blockwise(S, lam=0.05)
Theta_hat = np.linalg.inv(W)

print("估计的精度矩阵（绝对值 > 0.05 的位置视为边）：")
print(np.round(np.where(np.abs(Theta_hat) > 0.05, Theta_hat, 0), 3))
# 恢复出 0-1、1-2、3-4 三对条件依赖，与真实结构一致

2.14.4 随机矩阵理论

Marchenko-Pastur 定律

高维下样本协方差矩阵的特征值行为由随机矩阵理论描述。设 $\mathbf{X}$ 是 $n\times p$ 的矩阵，元素独立同分布于均值零、方差 ${\sigma }^2$ 的分布，且真实协方差为单位矩阵。当 $n,p\to \mathrm{\infty }$ 且 $p/n\to \gamma \in (0,1]$ 时，样本协方差矩阵 $\mathbf{S}=\frac{1}{n}{\mathbf{X}}^T\mathbf{X}$ 的特征值经验分布收敛到马尔琴科-帕斯图尔分布（Marchenko-Pastur law），其密度支撑在区间

$$[{\sigma }^2(1-\sqrt{\gamma }{)}^2,{\sigma }^2(1+\sqrt{\gamma }{)}^2]$$

上。

这一定律的含义是：即使真实协方差是单位矩阵（无任何相关性），样本协方差矩阵的特征值也不会都等于 ${\sigma }^2$，而是散布在一个区间内，区间宽度随 $\gamma =p/n$ 增大而变宽。最大样本特征值约为 ${\sigma }^2(1+\sqrt{\gamma }{)}^2$，远大于 $1$。

这一结论对 PCA 有直接影响。用样本特征值判断主成分个数时，噪声本身会贡献看似显著的大特征值。$\gamma =1$ 时最大噪声特征值可以达到 $4{\sigma }^2$，若阈值为 $1$，会把纯噪声方向当成信号。正确的做法是用 MP 分布的边界作为阈值，或改用基于随机矩阵理论的检验方法。

Tracy-Widom 分布与相变

最大特征值的涨落由 Tracy-Widom 分布（Tracy-Widom distribution）描述。把最大特征值标准化后，

$$\frac{{\lambda }_{max}-{\sigma }^2(1+\sqrt{\gamma }{)}^2}{{\sigma }^2{n}^{-2/3}(1+\sqrt{\gamma }{)}^{4/3}(1-\sqrt{\gamma }{)}^{-1/3}}\stackrel{d}{\to }{F}_1$$

其中 $F_1$ 是 Tracy-Widom 分布（与随机矩阵的普适性类有关，同一分布也出现在其他增长模型的最大值问题中）。有了这一分布，可以构造最大特征值是否超出噪声范围的检验。

信号是否存在还取决于其强度的下限。拜雷-埃斯金定律（Baik-Ben Arous-Péché law，BBP 相变）说明：若真实协方差有一个特征值 $\ell >{\sigma }^2(1+\sqrt{\gamma })$，对应的样本特征值会偏离噪声区间，其偏离程度不仅取决于 $\ell$，还伴随特征向量的相关。$\ell$ 低于这一阈值时，样本特征值被噪声区间吞没，PCA 无法察觉该方向；超过阈值后可以被检测，并被一致估计。这一相变现象说明高维下弱信号无法被检测，无论样本量多大，只要 $p/n$ 不变。

2.14.5 高维多重检验

高维 FDR 与本地 FDR

2.11 节给出 FDR 的概念与 BH 校正。高维情形下检验数量 $m$ 巨大（某些场景可达百万量级），p 值的分布形态成为关键信息。本地 FDR（local FDR）定义为给定 p 值（或给定检验统计量）时该假设为真的后验概率：

$$\text{lfdr}(t)=P(H_0\text{ 为真}\mid T=t)=\frac{{\pi }_0{f}_0(t)}{{\pi }_0{f}_0(t)+(1-{\pi }_0)f_1(t)}$$

其中 ${\pi }_0$ 是零假设为真的比例，$f_0$ 是零假设下统计量的密度，$f_1$ 是备择下的密度。本地 FDR 提供逐假设的判据，比整组的 FDR 阈值更精细。估计 ${\pi }_0$ 与 $f_1$ 需要额外假设（如 $f_0$ 已知、$f_1$ 光滑），实践中用经验零分布或 Ash 算法等去卷积方法估计。

高维下 FDR 控制的困难在于依赖结构：检验统计量之间往往正相关（相邻特征、相关指标），BH 的独立性假设不成立。应对方式是使用依赖结构下的修正（BY 方法）、重抽样校准（用置换分布估计零分布）或估计有效检验数（把相关导致的冗余折算掉）。

2.14.6 因果推断的问题与潜在结果框架

相关不等于因果

2.13 节的回归诊断已经提示过这一问题，这里给出正式的处理框架。两个直接原因使观测数据的关联偏离因果效应。混杂（confounding）：存在同时影响处理与结果的变量，使两者产生非因果的关联。选择（selection）：样本的处理状态与结果同时受某个变量影响，或被同一机制筛选。辛普森悖论（Simpson's paradox）是混杂导致结论反转的经典表现：分组数据中各组都显示正向关系，合并后关系反向，原因是分组变量的分布在不同处理组之间不同。

潜在结果

潜在结果（potential outcomes）框架（也称 Rubin 因果模型）为因果效应提供定义。对每个个体 $i$，定义两个潜在结果：$Y_i(1)$ 表示接受处理时的结果，$Y_i(0)$ 表示未接受处理时的结果。个体处理效应（individual treatment effect，ITE）为

$${\tau }_i=Y_i(1)-Y_i(0)$$

根本困难在于因果推断的基本问题：每个个体只能观测到其中一个潜在结果，另一个是反事实（counterfactual）。观测到的结果是 $Y_i=T_i{Y}_i(1)+(1-T_i)Y_i(0)$，其中 $T_i$ 是处理指示变量。

个体效应不可观测，因此推断的目标转向总体层面的量：平均处理效应（average treatment effect，ATE）

$${\tau }_{\text{ATE}}=E[Y(1)-Y(0)]$$

处理组的平均处理效应（ATT）与条件平均处理效应（conditional average treatment effect，CATE，$\tau (x)=E[Y(1)-Y(0)\mid \mathbf{X}=\mathbf{x}]$）也常用。CATE 描述效应的异质性，是精准干预与个性化决策的基础。

识别假设

从观测数据估计 ATE 需要三条假设：

可忽略性（ignorability，又称无混杂性、条件可交换性）：给定协变量 $\mathbf{X}$ 后，处理分配与潜在结果独立，

$$\{Y(1),Y(0)\}\perp T\mid \mathbf{X}$$

它的含义是所有同时影响处理与结果的变量都被观测并纳入了 $\mathbf{X}$。这一假设不能用数据检验，只能通过领域知识论证，这是因果推断中最需要谨慎的部分。

正性（positivity，又称重叠性）：对每个协变量取值，接受处理与不接受处理的概率都严格大于零，

$$0<P(T=1\mid \mathbf{X}=\mathbf{x})<1$$

违背时会出现某些类型个体只出现在一个处理组的情况，这些个体上的效应无法识别。

一致性（consistency，含 SUTVA）：个体的观测结果等于其处理状态对应的潜在结果，且一个个体的结果不受其他个体处理状态的影响（无干扰）。违背的例子包括疫苗的群体免疫效应、社交网络中的溢出效应。

三条假设成立时，ATE 可以用协变量调整后的期望表示：

$${\tau }_{\text{ATE}}=E_{\mathbf{X}}[E[Y\mid T=1,\mathbf{X}]-E[Y\mid T=0,\mathbf{X}]]$$

这就是后门调整的可估计形式，下节从图的角度重新表述。

随机化实验

随机化实验（randomized experiment）通过随机分配处理消除混杂：随机化保证 $T\perp \{Y(1),Y(0)\}$，可忽略性在无协变量条件下成立，正性在分配概率固定时成立。因此组间均值之差是 ATE 的无偏估计，不需要任何模型假设。这是实验设计在因果推断中地位最高的原因。

实践中的问题是随机化的完整执行与外部有效性：样本脱落、依从性不足、仅覆盖特定人群。前者破坏可忽略性（脱落往往与结果相关），后者限制结论的推广范围。

import numpy as np

# 辛普森悖论演示：分组内正效应，合并后变负
rng = np.random.default_rng(34)
group = rng.binomial(1, 0.5, size=4000)          # 0 为 A 类，1 为 B 类
# 处理分配依赖分组：A 类更可能被处理
T = rng.binomial(1, np.where(group == 0, 0.8, 0.2))
# 结果依赖分组与处理，组内效应为正
Y = 2.0 * group + 1.5 * T + rng.normal(scale=1.0, size=4000)

print(f"合并数据：处理组均值 {Y[T == 1].mean():.3f}，对照组均值 {Y[T == 0].mean():.3f}")
print(f"        合并差异 = {Y[T == 1].mean() - Y[T == 0].mean():+.3f}")

for g in [0, 1]:
    mask = group == g
    diff = Y[mask & (T == 1)].mean() - Y[mask & (T == 0)].mean()
    print(f"分组 {g} 内：处理组 {Y[mask & (T == 1)].mean():.3f}，"
          f"对照组 {Y[mask & (T == 0)].mean():.3f}，差异 {diff:+.3f}")
# 合并差异符号可能相反，分组内差异都接近 1.5

2.14.7 因果图与 do 算子

结构因果模型与有向无环图

2.8 节已引入结构因果模型与有向无环图的形式。因果图用箭头表示因果方向，与贝叶斯网络中箭头只表示条件依赖不同：$X\to Y$ 表示 $X$ 是 $Y$ 的原因，改变 $X$ 会改变 $Y$ 的分布。

d-分离（d-separation）给出图与条件独立性的对应：若一条路径被某个变量集合阻断，则该路径上的信息不传递。阻断规则识别三种结构的中间节点：链式 $X\to Z\to Y$ 与共因 $X\leftarrow Z\to Y$ 中的 $Z$ 被控制时路径阻断，对撞 $X\to Z\leftarrow Y$ 中的 $Z$ 被控制时路径反而打开（2.6 节的对撞偏倚）。

do 算子与 do-calculus

干预分布 $P(Y\mid do(X=x))$ 与观测条件分布 $P(Y\mid X=x)$ 一般不同。do 算子的计算在图上对应删除所有指向 $X$ 的边，把 $X$ 固定为 $x$，再在修改后的图上计算分布。do-calculus 给出三条从观测分布推导干预分布的规则，它们分别处理在干预下添加或删除观测条件、交换干预与观测顺序、删除与干预无关的变量的情形。三条规则配合图的拓扑结构，可以判断目标量是否可识别。

后门准则与前门准则

后门准则（back-door criterion）给出最常用的识别条件。变量集合 $\mathbf{Z}$ 满足后门准则，若它阻断了所有从 $X$ 到 $Y$ 的、以指向 $X$ 的箭头开头的路径，且不包含 $X$ 的后代。此时

$$P(Y\mid do(X=x))=\sum _{\mathbf{z}}P(Y\mid X=x,\mathbf{Z}=\mathbf{z})P(\mathbf{Z}=\mathbf{z})$$

后门调整（back-door adjustment）就是按混杂变量的分布加权平均条件关系。这一公式是分层分析与回归调整的因果依据：把 $\mathbf{Z}$ 控制在回归中，得到的系数对应调整后的效应。注意控制变量集合的选取由图结构决定，控制对撞变量或中介变量会引入偏倚甚至掩盖真实效应。

前门准则（front-door criterion）处理混杂变量不可观测的情形。当 $X$ 到 $Y$ 的效应完全通过一个可观测的中介 $M$ 时，且 $X$ 到 $M$ 与 $M$ 到 $Y$ 没有未观测混杂，则

$$P(Y\mid do(X=x))=\sum _{m}P(M=m\mid X=x)\sum _{x^{\prime }}P(Y\mid M=m,X=x^{\prime })P(X=x^{\prime })$$

它对两个可识别的部分分别估计后组合，绕过了不可观测的混杂。前门准则的实用价值在于它说明识别路径不唯一，图的局部结构决定了哪些量可以估计。

2.14.8 识别策略

工具变量

工具变量（instrumental variable，IV）适用于存在未观测混杂的场景。有效工具需满足：与处理相关（相关性）、只通过处理影响结果（排他性）、与潜在结果的条件分布独立（外生性）。估计用两阶段最小二乘（2.13 节），弱工具的诊断用第一阶段 F 统计量。IV 给出的效应对依从者（compliers，即处理状态会随工具变化的个体）而言是局部的，这一群体往往无法明确识别，解释效应时需要说明。

倾向得分方法

倾向得分（propensity score）是接受处理的概率：$e(\mathbf{x})=P(T=1\mid \mathbf{X}=\mathbf{x})$。它的核心性质是平衡性：在倾向得分相同的条件下，协变量在处理组与对照组的分布相同，可忽略性在单一维度上成立。

$$\{Y(1),Y(0)\}\perp T\mid e(\mathbf{X})$$

这一降维把高维协变量调整问题化为一维匹配或加权问题。倾向得分匹配（propensity score matching）为每个处理组个体寻找倾向得分接近的对照个体，常用最近邻匹配（取最近的若干个）、卡尺匹配（限制最大距离）、马氏距离匹配（在倾向得分之外加入协变量的距离）。匹配后需要检查平衡性：匹配样本中协变量的标准化差异是否足够小（经验规则是小于 $0.1$）。

逆概率加权（inverse probability weighting，IPW）用权重重构一个伪总体，使处理组与对照组在协变量上可比：

$${\hat{\tau }}_{\text{IPW}}=\frac{1}{n}\sum _{i=1}^n[\frac{T_i{Y}_i}{\hat{e}({\mathbf{x}}_i)}-\frac{(1-T_i)Y_i}{1-\hat{e}({\mathbf{x}}_i)}]$$

倾向得分接近 $0$ 或 $1$ 的个体会获得极大权重，使估计的方差失控。双重稳健估计（doubly robust estimation）把 IPW 与结果回归结合：两者只要有一个正确指定，估计就一致。这一性质使其成为实践中的默认选择，常见形式是增广逆概率加权（AIPW）。边际结构模型（marginal structural model）把 IPW 推广到时变处理的情形，用结构嵌套模型处理时变混杂。

import numpy as np
from scipy import optimize

def estimate_ate_ipw(X, T, Y, l2=1.0):
    """用逻辑回归估计倾向得分，再计算逆概率加权的 ATE"""
    n, p = X.shape
    Xc = np.column_stack([np.ones(n), X])

    def neg_loglik(beta):
        eta = Xc @ beta
        return -np.sum(T * eta - np.logaddexp(0, eta)) + l2 * beta @ beta / 2

    beta = optimize.minimize(neg_loglik, np.zeros(p + 1), method="BFGS").x
    e = 1 / (1 + np.exp(-Xc @ beta))
    e = np.clip(e, 1e-3, 1 - 1e-3)                 # 截断极端权重，控制方差
    return (T * Y / e).mean() - ((1 - T) * Y / (1 - e)).mean(), e

rng = np.random.default_rng(45)
n = 5000
# 混杂变量
x1 = rng.normal(size=n)
x2 = rng.normal(size=n)
# 处理分配依赖混杂
e_true = 1 / (1 + np.exp(-(0.8 * x1 - 0.5 * x2)))
T = rng.binomial(1, e_true)
# 结果同时依赖混杂与处理，真实 ATE = 2
Y = 1.0 * x1 + 0.5 * x2 + 2.0 * T + rng.normal(size=n)

naive = Y[T == 1].mean() - Y[T == 0].mean()
ate_ipw, e_hat = estimate_ate_ipw(np.column_stack([x1, x2]), T, Y)
print(f"朴素差异（有偏）  = {naive:+.3f}")
print(f"IPW 估计          = {ate_ipw:+.3f}（真值 2.0）")
print(f"倾向得分范围      = [{e_hat.min():.3f}, {e_hat.max():.3f}]")

准实验设计

部分设计利用数据生成过程中的制度性断点或时间结构来识别效应，不需要可忽略性假设。

断点回归（regression discontinuity）适用于处理由某个连续变量是否超过阈值决定的情形。阈值两侧的个体在协变量上近似可比，只有在阈值附近的一个小邻域内估计效应才可靠。模糊断点回归（fuzzy RD）处理阈值只改变处理概率而不完全决定处理状态的情形，用阈值作为工具变量。实施要点是带宽选择、断点附近的密度检验（排除操纵阈值）与协变量连续性检验。

双重差分（difference-in-differences，DID）比较处理组与对照组在处理前后的变化差异，识别假设是平行趋势：若无处理，两组的趋势相同。这一假设可以用处理前的数据部分检验。合成控制（synthetic control）在没有合适的对照组时，用多个未处理单元的加权组合构造一个反事实的对照，权重的选择使构造单元在处理前与处理单元的轨迹接近。

中介分析

中介分析（mediation analysis）把总效应分解为通过中介变量的间接效应与其余的直接效应。潜在结果框架下的定义需要引入嵌套的潜在结果 $Y(t,m)$：

$$\text{自然直接效应}=E[Y(1,M(0))-Y(0,M(0))],\text{自然间接效应}=E[Y(1,M(1))-Y(1,M(0))]$$

两者的和等于总效应。识别需要额外的假设：中介与结果之间没有未观测混杂，处理与中介之间没有未观测混杂。这些假设比可忽略性更强，敏感性分析尤其重要。

2.14.9 敏感性分析与因果发现

敏感性分析

可忽略性无法检验，因此需要评估结论对假设偏离的稳健程度。敏感性分析（sensitivity analysis）量化未测量混杂需要多强才能推翻结论。E 值（E-value）给出所需的效应量尺度：

$$\text{E 值}={\text{OR}}_{\text{obs}}+\sqrt{{\text{OR}}_{\text{obs}}({\text{OR}}_{\text{obs}}-1)}$$

它表示未测量混杂与处理和结果的关联强度至少需要达到这一水平，才能使观测到的效应完全由混杂解释。E 值远大于已知混杂的强度时，结论较为稳健；E 值接近 $1$ 时，任何微弱的混杂都能推翻结论。其他方法包括对混杂效应施加先验分布的贝叶斯敏感性分析、用负面对照结果做校准。

因果发现

因果发现（causal discovery）试图从数据推断图结构，前提假设比效应估计更强，输出也通常只能确定到马尔可夫等价类（同一等价类中的图蕴含相同的条件独立关系）。

PC 算法（Peter-Clark algorithm）基于条件独立检验：从完全图出发，逐个删除相互独立的变量对之间的边，再根据对撞结构的定向规则确定箭头方向。GES 算法（greedy equivalence search）在等价类空间中做贪心搜索，用评分（如 BIC）代替独立性检验，通常更稳定。两者都假设无未观测混杂（因果充分性）。

LiNGAM（linear non-Gaussian acyclic model）利用非高斯性识别方向：在线性、非高斯、无环的假设下，因果方向可以由独立性检验唯一确定，突破了马尔可夫等价类的限制。加性噪声模型（additive noise model）把方向识别建立在噪声独立性上：$Y=f(X)+\epsilon$ 中若 $\epsilon$ 与 $X$ 独立，则相反方向的同类表示不成立。因果表示学习（causal representation learning）把这一问题推广到从高维观测中学习潜在因果变量，是当前研究的前沿方向，与强化学习中的策略迁移问题相关。

2.14.10 本节小结与全章回顾

高维与因果：两条边界的处理方式

本节处理了经典方法的两条边界。高维情形下，$p$ 接近或超过 $n$ 使最小二乘无解且过拟合成为默认状态，解决方案是引入稀疏或低秩结构：Lasso 与压缩感知利用稀疏性，核范数最小化利用低秩性，两者的关系是 $L_1$ 与奇异值层面的对应。随机矩阵理论给出了高维噪声特征值的分布，Marchenko-Pastur 定律说明即使真实协方差为单位矩阵，样本特征值也会散布在一个区间内，这一结论直接修正了用特征值判断主成分个数的方法。高维协方差估计与图模型选择通过结构假设把高维问题分解为低维子问题。因果推断部分从潜在结果的不可观测性出发，给出可忽略性、正性与一致性三条识别假设，用因果图与 do 算子把识别问题转化为图上的路径判断，后门与前门准则给出两个可用的公式。识别策略按可利用的结构分为工具变量、倾向得分方法、断点回归、双重差分与合成控制，敏感性分析评估结论对未观测混杂的稳健程度。

第二章的内容到此结束。回顾全章的线索：2.1 到 2.3 建立了概率论的语言，用概率空间描述不确定性、用随机变量把事件编码为数值、用数字特征概括分布；2.4 与 2.5 按生成机制整理了常见分布族，为建模提供可直接使用的构件；2.6 处理多维情形，联合分布、条件分布与协方差结构是多元统计的基础；2.7 讨论大样本行为，收敛模式、极限定理与集中不等式给出了误差的量级；2.8 把独立同分布的假设放宽为依赖结构，随机过程与概率图模型提供了描述时序与网络的语言。统计推断部分从 2.9 的问题设定出发，2.10 与 2.11 分别给出估计与检验两套工具，2.12 用先验与后验重构推断流程并讨论计算方法，2.13 把它们组织成回归建模的框架，2.14 处理高维与因果这两条现代数据的边界。

这些工具与后续章节的分工是清楚的。最优化章节处理求解问题：极大似然、变分推断与正则化估计最终都要落到数值优化上，收敛性分析与算法选择是那里的内容。信息论章节处理度量问题：熵、KL 散度与互信息给出分布的差异度量，本节与 2.3 节中出现的 Fisher 信息与距离相关在那里得到统一处理。两章的内容都与本章交叉，共同构成机器学习的理论底座。

练习题

第 1 题 概念推导

设 $\mathbf{X}$ 为 $n\times p$ 的设计矩阵，$p>n$。说明为什么最小二乘解不唯一，并给出解集的一般形式；再说明加入 $L_2$ 惩罚后解为什么唯一，并写出解的表达式。

参考答案

解不唯一的原因：最小二乘的法方程为 ${\mathbf{X}}^T\mathbf{X}\mathit{\beta }={\mathbf{X}}^T\mathbf{y}$。$\mathbf{X}$ 的秩最多为 $min(n,p)=n$，因此 ${\mathbf{X}}^T\mathbf{X}$ 的秩最多为 $n$，而它是 $p\times p$ 的矩阵且 $p>n$，故不可逆，法方程的解不唯一。解集的一般形式为

$$\mathit{\beta }={\mathit{\beta }}_0+\mathbf{v},\mathbf{v}\in \mathrm{Null}(\mathbf{X})$$

其中 ${\mathit{\beta }}_0$ 是任意一个特解（例如最小范数解 ${\mathbf{X}}^{+}\mathbf{y}$，${\mathbf{X}}^{+}$ 是 Moore-Penrose 伪逆）。零空间维数为 $p-\mathrm{rank}(\mathbf{X})\ge p-n>0$，说明有无穷多组系数给出相同的拟合值与相同的残差平方和。

加入 $L_2$ 惩罚后唯一：目标函数变成

$$\parallel \mathbf{y}-\mathbf{X}\mathit{\beta }{\parallel }^2+\lambda \parallel \mathit{\beta }{\parallel }^2$$

对 $\mathit{\beta }$ 求导得 $({\mathbf{X}}^T\mathbf{X}+\lambda \mathbf{I})\mathit{\beta }={\mathbf{X}}^T\mathbf{y}$。${\mathbf{X}}^T\mathbf{X}$ 半正定，$\lambda \mathbf{I}$ 正定，故 ${\mathbf{X}}^T\mathbf{X}+\lambda \mathbf{I}$ 正定、可逆，解唯一：

$${\hat{\mathit{\beta }}}_{\text{ridge}}={({\mathbf{X}}^T\mathbf{X}+\lambda \mathbf{I})}^{-1}{\mathbf{X}}^T\mathbf{y}$$

除了保证唯一性，惩罚项还使解在零空间方向上被压缩到零，即所有解中范数最小的那一个（$\lambda \to 0$ 时趋近最小范数解）。从贝叶斯角度看，$L_2$ 惩罚对应系数的高斯先验，先验在零空间方向上确定了原本无法确定的部分。

第 2 题 计算推理

某平台有 $n=200$ 条用户记录，$p=800$ 个特征。请判断普通最小二乘是否可以用于建模；若使用 Lasso，说明非零系数的个数上限与样本量的关系；并估计随机矩阵理论下最大噪声特征值相对于真实方差的比例（取 $p/n=4$ 的近似值，说明计算依据）。

参考答案

最小二乘的可行性：$p=800>n=200$，设计矩阵秩不足以使 ${\mathbf{X}}^T\mathbf{X}$ 可逆，最小二乘解不唯一。即使强行取最小范数解（伪逆），模型也能完全拟合训练数据，样本外表现没有保证，因此不适合直接使用。

Lasso 的非零系数上限：Lasso 的解由 KKT 条件确定，在一般位置下非零系数的个数不超过样本量 $n$。原因是每个非零系数对应一个等式约束，而可用的约束个数受设计矩阵秩的限制。因此 $n=200$ 时最多选出约 $200$ 个变量（实际通常远少于此）。这说明 Lasso 天然完成了降维：最终模型的有效复杂度不超过样本量。

最大噪声特征值：随机矩阵理论中，真实协方差为单位矩阵（方差 ${\sigma }^2$）时，样本协方差的最大特征值约为 ${\sigma }^2(1+\sqrt{\gamma }{)}^2$，其中 $\gamma =p/n$。取 $\gamma =4$（$p$ 远大于 $n$ 的极端情形），

$$\frac{{\lambda }_{max}}{{\sigma }^2}\approx (1+\sqrt{4}{)}^2=9$$

即纯噪声也能产生达到真实方差 $9$ 倍的特征值。若用固定的阈值（例如 $2{\sigma }^2$）判断是否存在主成分，会选出大量纯噪声方向。正确的做法是根据 $\gamma$ 计算 MP 分布的上界作为阈值，或使用假设检验。这一例子说明高维下特征值的大小本身不是信号的证据，必须与噪声的预期分布比较。

第 3 题 代码验证

构造一个含混杂的观测数据，比较三种估计方式：未调整的组间均值差、后门调整（在回归中加入混杂变量）、逆概率加权。用一个已知真值的数据生成过程验证三种估计的偏差；再演示当混杂变量未被观测时，后门调整失效的情形。

参考答案

import numpy as np
from scipy import optimize

def estimate_ate_ipw(X, T, Y, l2=1.0):
    n, p = X.shape
    Xc = np.column_stack([np.ones(n), X])

    def neg_loglik(beta):
        eta = Xc @ beta
        return -np.sum(T * eta - np.logaddexp(0, eta)) + l2 * beta @ beta / 2

    beta = optimize.minimize(neg_loglik, np.zeros(p + 1), method="BFGS").x
    e = np.clip(1 / (1 + np.exp(-Xc @ beta)), 1e-3, 1 - 1e-3)
    return (T * Y / e).mean() - ((1 - T) * Y / (1 - e)).mean()

rng = np.random.default_rng(52)
n = 20_000
u = rng.normal(size=n)                 # 混杂变量（可观测）
w = rng.normal(size=n)                 # 未观测混杂
T = rng.binomial(1, 1 / (1 + np.exp(-(1.2 * u + 1.0 * w))))
Y = 1.5 * u + 1.2 * w + 2.0 * T + rng.normal(size=n)

# 情形一：只观测到 u，w 未观测
naive = Y[T == 1].mean() - Y[T == 0].mean()
Xc = np.column_stack([np.ones(n), u, T])
coef = np.linalg.lstsq(Xc, Y, rcond=None)[0]
print(f"未调整差异            = {naive:+.3f}")
print(f"回归调整（控制 u）     = {coef[2]:+.3f}")
print(f"IPW（控制 u）          = {estimate_ate_ipw(u.reshape(-1, 1), T, Y):+.3f}")
print(f"真值                  = +2.000")
# 三者都偏离真值，因为 w 未观测，可忽略性不成立

# 情形二：u 与 w 都观测
X2 = np.column_stack([u, w])
Xc2 = np.column_stack([np.ones(n), u, w, T])
coef2 = np.linalg.lstsq(Xc2, Y, rcond=None)[0]
print(f"\n控制 u 与 w 后的回归系数 = {coef2[3]:+.3f}")
print(f"控制 u 与 w 后的 IPW     = {estimate_ate_ipw(X2, T, Y):+.3f}")

两种情形对比说明可忽略性的作用：混杂变量全部观测到时，回归调整与 IPW 都接近真值 $2.0$；存在未观测混杂时，两种方法都给出有偏估计。增加样本量不能消除这一偏差，这是敏感性分析存在的原因。实践中还应检查倾向得分的取值范围，出现接近 $0$ 或 $1$ 的极端值时需要截断或改用匹配方法，因为少数极端权重会主导整个估计。

第 4 题 综合应用

某平台想评估新版推荐算法对用户停留时长的影响，但无法做随机实验。已知数据包含用户属性、历史行为与新算法是否启用（由平台的灰度策略决定，与用户的历史活跃度相关）。请设计一个估计方案：列出需要假设的条件、选择识别策略、说明必须检查的诊断项，并给出结论报告时应包含的内容。

参考答案

需要假设的条件：核心假设是可忽略性：给定用户属性与历史行为后，算法分配与潜在停留时长独立，即不存在未观测的同时影响分配与结果的因素。灰度策略按用户分组推进时，若分组依据是活跃度、地域、设备等可观测属性，这一假设有可能成立；若分组还依据运营对用户偏好的主观判断，则该判断若未被记录则假设可能违背。另外需要正性（每种用户类型在两组中都出现）与一致性（无跨用户干扰，要注意社交分享可能带来的溢出效应）。

识别策略：以用户属性与历史行为为协变量，采用双重稳健估计（结果回归加逆概率加权的 AIPW 形式），兼顾模型错配风险。若灰度策略基于某个明确的阈值（例如活跃度分数超过某个值才启用），改用它做断点回归，识别条件更弱，只在阈值附近取样本。若灰度是分批推进的，可以用双重差分比较早用组与晚用组在推进前后的变化差异，需要检验平行趋势。

必须检查的诊断项：其一，倾向得分的重叠情况，画出两组的倾向得分分布，检查是否有大量个体落在极端区域，必要时截断权重或改用匹配。其二，协变量平衡性，计算加权或匹配后各协变量的标准化差异，应小于 $0.1$。其三，若用双重差分，检验处理前的趋势是否平行。其四，若用断点回归，检验阈值附近的协变量连续性与密度（排除操纵阈值）。其五，做敏感性分析，计算 E 值并与已知混杂的强度对比。

报告内容：给出效应的点估计与区间；说明识别假设的内容与论证依据；列出平衡性与重叠性的诊断结果；报告敏感性分析的结论，说明需要多强的未观测混杂才能推翻结论；说明效应针对的群体（若用 IV 则是依从者，若用断点回归则是阈值附近的个体），避免把局部效应当作全平台效应；最后给出结论对假设偏离的稳健性评估，而非单一的显著性判断。

常见错误

错误 1 · 在高维下用 $R^2$ 或训练误差评价模型

原因：$p$ 接近或超过 $n$ 时模型可以完全拟合训练数据，$R^2$ 接近 $1$、训练误差接近零，这些数值不包含任何关于样本外表现的信息。$p>n$ 时甚至可以构造出对随机标签也完全拟合的模型。

解决：评价一律使用交叉验证或独立的留出集，报告样本外误差。特征筛选必须放在交叉验证的每一折内部进行，否则选择过程本身会泄漏信息，使误差估计偏乐观。

错误 2 · 用样本特征值的大小判断主成分个数

原因：高维下噪声本身会产生较大的样本特征值。Marchenko-Pastur 定律说明即使真实协方差为单位矩阵，最大样本特征值也可达 ${\sigma }^2(1+\sqrt{p/n}{)}^2$，$p/n=1$ 时约为 $4{\sigma }^2$，$p/n=4$ 时约为 $9{\sigma }^2$。用固定的经验阈值会把噪声方向当作信号。

解决：用噪声分布的上界作为阈值（按 $\gamma =p/n$ 计算 MP 边界），或用基于 Tracy-Widom 分布的检验。更稳妥的做法是结合平行分析（用随机数据生成同样规模的矩阵并比较特征值）与领域知识判断，不单纯依赖数值大小。

错误 3 · 把关联直接解释为因果而不检查识别条件

原因：混杂存在时回归系数会系统偏离因果效应，偏离方向与大小取决于混杂结构，样本量增加不能消除偏差。对撞变量被当作混杂控制进去时，反而会引入偏倚，甚至使原本存在的效应消失或反转。

解决：给出因果结论前明确写出识别假设，并用图结构判断应控制哪些变量：控制共因与链式结构的中间变量，避免控制对撞变量与中介变量。无法确定时报告多种控制方案下的结果，说明结论对变量选择的敏感性，并做敏感性分析。

错误 4 · 把倾向得分的极端权重直接代入估计

原因：倾向得分接近 $0$ 或 $1$ 的个体会获得极大权重，少数几个观测可能主导整个估计。若某个用户类型几乎只出现在一组中（违背正性），权重会趋向无穷，估计失去意义。

解决：先检查倾向得分的分布与重叠性，对超出范围的值截断（如限制在 $[0.01,0.99]$）或用稳定权重（乘上处理概率的边际比例）。重叠性严重不足时放弃 IPW，改用匹配限制在共同支撑区域内估计，或明确说明效应只在重叠区域可识别。