LASSO 与 Elastic Net：高维变量的正则化

上一章把 42 个 Bao 特征塞进一个普通逻辑回归，得到测试集 AUC 0.6966。这个数字并不差，但有一处令人不安。训练样本舞弊数只有 537，自变量却有 42 个，再加上常数项是 43 个待估系数。每个非舞弊样本对应的"信息量"被 42 个系数瓜分，普通逻辑回归会把训练集里的随机噪声当作信号学进去，在新数据上表现波动剧烈。会计文献里对此有专门的术语，叫"过拟合"。

惩罚回归是应对过拟合的标准工具。在原本的对数似然函数后面加一个对系数大小的惩罚项，强迫模型把那些"看起来有用、其实是噪声"的变量系数压向零。Tibshirani 在 1996 年提出 LASSO，用 $\sum |\beta_j|$ 形式的 L1 惩罚做变量选择；Hoerl 与 Kennard 早在 1970 年代提出 Ridge，用 $\sum \beta_j^2$ 形式的 L2 惩罚做整体收缩；Zou 与 Hastie 在 2005 年把两者按比例混合，提出 Elastic Net。本章把这三种惩罚都套到 Bao 的 42 个特征上，看哪一种最适合舞弊检测这种正样本极稀、变量高度相关的场景。

高维财务变量下的过拟合风险

第三章的逻辑回归用的是极大似然估计。最大化对数似然等价于最小化交叉熵损失 $-\sum_i [y_i \log \hat p_i + (1-y_i) \log (1-\hat p_i)]$。在低维数据上这个目标函数有唯一的极小值；在高维数据上它会沿着多个方向变得很平，模型可以在保持训练损失不变的前提下任意放大某些系数。会计变量之间存在天然的高相关，比如总资产 at 和总负债 lt 的相关系数通常在 0.95 以上，普通股权益 ceq 和留存收益 re 也高度相关。共线变量让损失函数沿着"一个系数升、另一个降"的方向几乎完全平坦，普通最大似然没有办法在这种平坦的方向上做出选择。

把这个问题翻译成具体数字。Bao 数据训练集 63,930 行，537 个舞弊样本，42 个特征。每个特征的标准差被先归一化到 1，意味着每个 $\beta_j$ 的"单位影响"是可比的。如果我们把 at 的系数设为 $+10$，lt 的系数设为 $-10$，由于这两个变量高度相关，它们的预测值之和与原本系数都为 0 时几乎一致，训练损失不会升高多少，但模型已经开始报告对总资产和总负债"非常敏感"。换到测试集 2009 至 2014 年的数据上，新公司的 at 与 lt 关系略有变化，模型的预测就会剧烈震荡。

惩罚回归的解决思路是在损失函数里加一项"系数不要乱动"的代价。形式化地，惩罚回归求解

其中 $\hat p_i = 1 / (1 + \exp(-x_i^\top \beta))$ 是逻辑回归的预测概率，$P(\beta)$ 是惩罚函数，$\lambda \ge 0$ 是惩罚强度。$\lambda = 0$ 退化为普通逻辑回归；$\lambda \to \infty$ 让所有非截距系数被压成零，模型退化为只用截距预测。中间的某个 $\lambda$ 在偏差与方差之间取得平衡，这就是交叉验证要找的最优值。

L1 在原点处不可微，这一点表面上是数学麻烦，实际上是它能做变量选择的根源。下一节会从几何上解释为什么 L1 会让某些系数恰好等于零。

L1 与 L2 惩罚的几何含义

先用一个二维数字例子建立直觉。假设只有两个特征 soft_assets 与 issue，无惩罚下的最大似然估计是 $\hat\beta = (0.6, 0.4)$。L1 惩罚要求 $|\beta_1| + |\beta_2| \le t$ 的约束在二维平面上是一个以原点为中心、四个顶点在坐标轴上的菱形；L2 惩罚要求 $\beta_1^2 + \beta_2^2 \le t$ 的约束是一个以原点为中心的圆。损失函数的等高线是绕 $\hat\beta$ 的椭圆。带惩罚的最优解是损失等高线第一次接触到约束边界的那一点。

椭圆与菱形相切，最容易切到菱形的尖角，也就是其中某个轴的顶点。一旦切到 $(\beta_1, 0)$ 这样的顶点，第二个系数就被精确压成零。L1 的菱形约束让损失函数的等高线倾向于在尖角处遇到约束，于是 LASSO 会把一部分变量的系数恰好设为零，自动完成变量选择。椭圆与圆相切的位置则是连续变化的，没有任何方向具有"角点优势"，Ridge 把所有系数都向零的方向收缩，但不会精确等于零。

把这个二维直觉推广到 42 维：LASSO 在 42 维超菱形上找等高线切点，会把多数系数压成零，最后只留下几个核心变量；Ridge 在 42 维超球面上找切点，所有 42 个系数都被收缩，但都保留下来。Elastic Net 用一个混合约束 $\alpha \sum|\beta_j| + (1-\alpha) \sum \beta_j^2 \le t$，几何形状介于菱形与圆之间，其顶点不那么尖，但仍然能让一部分系数为零。Zou 与 Hastie 的原文证明，当变量之间存在强相关时，Elastic Net 比 LASSO 更稳定，因为 LASSO 倾向于"在高度相关的一组变量里随机挑一个保留下来"，而 Elastic Net 会把它们一起保留或一起压低。

图 4·1　L1 与 L2 惩罚的几何对比。椭圆是损失函数的等高线，菱形与圆是惩罚约束的边界，最优解出现在两者首次相切处。L1 的菱形有四个角点，损失等高线最容易在角点处接触约束，使其中一个分量恰好为零，自动完成变量选择；L2 的圆面无角点，最优解处所有分量都被等比例缩小但都保留下来。完整 TikZ 结构图详见 PDF 全文。

停下来想一想。 如果所有特征都被先标准化到方差 1，惩罚的"单位"才是统一的。如果不标准化直接做 LASSO，那些天然量级很大的变量比如总资产 at，单位是百万美元，会被惩罚得过狠；量级小的变量比如比率类则几乎不受惩罚。所以惩罚回归的所有现成实现里都默认开启 standardize = TRUE。后面在 Bao 数据上做 LASSO 时，我们手动在训练集上拟合 z-score 标准化器，然后只用训练集的均值与标准差去 transform 验证集和测试集，避免把未来年份的均值信息泄漏到训练阶段。

Elastic Net 的混合策略

Elastic Net 的损失函数写为

其中 $\log L_i(\beta) = y_i \log \hat p_i + (1-y_i) \log (1-\hat p_i)$ 是单个样本的对数似然，$\alpha \in [0, 1]$ 是混合比例，$\lambda$ 仍然是整体惩罚强度。$\alpha = 1$ 退化为纯 LASSO，$\alpha = 0$ 退化为纯 Ridge。Bao 数据里有大量"原始值与衍生比率"的变量对，比如净利润 ni 与留存收益 / 总资产 reoa，某种程度上重复了同一信息，Elastic Net 的混合策略可以让这种成对的高相关变量平稳进入模型。

实践中的 $\alpha$ 选择有两条经验。第一条：如果研究目标是变量解读，倾向取 $\alpha = 1$，让 LASSO 给出一个稀疏的"被选中变量清单"。第二条：如果研究目标是预测稳健性，倾向取中间值 $\alpha = 0.5$，让模型在变量选择与系数平稳之间折中。本章对三个 $\alpha \in \{0, 0.5, 1\}$ 都做完整 CV，作为方法对照。

时间感知交叉验证

到这里读者可能会想：选 $\lambda$ 不是有现成的 cv.glmnet 吗？直接调它就好了。不行。cv.glmnet 的默认行为是把训练集随机切成 10 折，每一折当作验证集，剩下九折当作训练集。在静态横截面数据上这种做法标准且高效，在 Bao 这种 1991 至 2002 年跨 12 个会计年度的面板数据上会出大问题。

具体的问题是：随机折让模型在第 5 折训练时同时用到 1991 与 2002 年的样本，预测的"验证"对象包含 1995 年的样本。模型在做 1995 年预测时，已经"看到"了 2002 年发生的舞弊模式。换句话说，未来年份的样本被混入训练，模型的 CV 性能因此被高估。等到真正用 2009 至 2014 年测试时，模型从未见过该时段的样本，性能必然落差。这种泄漏在普通文献里通称"future information leakage"，在时间序列与因果推断里叫"穿越未来"。

正确的做法是 forward-chaining 时间感知 CV：第 1 折训练 1991 至 1995 年，验证 1996 年；第 2 折训练 1991 至 1996 年，验证 1997 年；以此类推到第 7 折训练 1991 至 2001 年，验证 2002 年。每一折的训练数据严格早于验证数据，模型不可能"看到未来"。这套切分方法的代价是不能用 cv.glmnet 默认 foldid 接口，因为它假设每行只属于一折，而 forward-chaining 中每一行可能在多折训练阶段重复出现，需要手写循环来实现。

在 Bao 数据上的实现

把上面的标准化、forward-chaining CV、三种惩罚拼到一起，就得到本章的核心代码。完整脚本在 code/04_penalized.R，这里展示骨架部分。

suppressPackageStartupMessages({
  library(tidyverse); library(glmnet); library(pROC); library(here)
})
set.seed(2026)

d <- read_csv(here::here("data", "bao2020_full.csv"),
              show_col_types = FALSE)

# 42 个 Bao 特征：28 原始 + 14 衍生比率
features <- c("act","ap","at","ceq","che","cogs","csho","dlc","dltis",
              "dltt","dp","ib","invt","ivao","ivst","lct","lt","ni",
              "ppegt","pstk","re","rect","sale","sstk","txp","txt",
              "xint","prcc_f","dch_wc","ch_rsst","dch_rec","dch_inv",
              "soft_assets","ch_cs","ch_cm","ch_roa","issue","bm",
              "dpi","reoa","EBIT","ch_fcf")

# Bao 时间切分 + 剔除 42 特征任一为 NA 的行
train <- d %>% filter(fyear >= 1991, fyear <= 2002) %>%
  drop_na(all_of(features))
test  <- d %>% filter(fyear >= 2009, fyear <= 2014) %>%
  drop_na(all_of(features))

# z-score 仅在训练集上拟合
mu <- colMeans(train[, features])
sg <- apply(train[, features], 2, sd); sg[sg == 0] <- 1
z_apply <- function(df) {
  m <- as.matrix(df[, features])
  sweep(sweep(m, 2, mu, "-"), 2, sg, "/")
}
X_train <- z_apply(train); y_train <- train$misstate
X_test  <- z_apply(test);  y_test  <- test$misstate

# 时间感知 CV：forward-chaining 1991..year-1 训练、year 验证
fold_years <- 1996:2002
lambda_grid <- exp(seq(log(1e-1), log(1e-6), length.out = 80))

run_time_cv <- function(alpha_val) {
  cv_auc <- matrix(NA_real_, nrow = length(fold_years),
                   ncol = length(lambda_grid))
  for (k in seq_along(fold_years)) {
    yr <- fold_years[k]
    tr_idx <- which(train$fyear < yr)
    va_idx <- which(train$fyear == yr)
    fit <- glmnet(X_train[tr_idx, ], y_train[tr_idx],
                  family = "binomial", alpha = alpha_val,
                  lambda = lambda_grid, standardize = FALSE)
    pp <- predict(fit, newx = X_train[va_idx, ], type = "response")
    for (j in seq_len(ncol(pp))) {
      if (length(unique(pp[, j])) < 2) next
      cv_auc[k, match(fit$lambda[j], lambda_grid)] <-
        as.numeric(auc(roc(y_train[va_idx], pp[, j], quiet = TRUE)))
    }
  }
  mean_auc <- colMeans(cv_auc, na.rm = TRUE)
  list(lambda = lambda_grid[which.max(mean_auc)],
       cv_auc = max(mean_auc, na.rm = TRUE))
}

cv_lasso <- run_time_cv(1.0)
cv_ridge <- run_time_cv(0.0)
cv_enet  <- run_time_cv(0.5)

表 4·1　时间感知 7 折 CV 选最优 λ

Ridge 的最优 $\lambda$ 比 LASSO 大约高三个数量级。这个差距来自 L1 与 L2 惩罚函数本身的尺度差异，不是异常现象。L2 是平方和，单个系数 $\beta_j = 0.1$ 贡献 0.01 的惩罚；L1 是绝对值之和，同样的系数贡献 0.1 的惩罚。要让两种惩罚施加相近的"压力"，L2 的 $\lambda$ 必须比 L1 的 $\lambda$ 大得多。两个 $\lambda$ 在数值上不可直接比较。

LASSO 选出的变量

LASSO 在最优 $\lambda$ 下保留了 28 个非零系数，剔除了 14 个。下表列出按系数绝对值排序的前十个变量。注意所有特征已先 z-score 标准化，所以系数的绝对值在不同变量之间可比。

表 4·2　LASSO 最优 λ 下系数绝对值前十的变量

最大的系数 soft_assets 是 Dechow 等人在 2011 年文献里反复强调的舞弊信号：软性资产，含应收账款、存货、商誉等可塑性强的科目，占总资产比例越高，公司账面操纵的空间就越大。第二位的 issue 也符合会计直觉，公司在外部融资压力下更可能美化报表来支持募资。第三位的 dlc 系数为负，意味着控制其它变量后，短期债务越高反而舞弊概率略低；这并不矛盾，可能是高短期债务的公司本身现金流压力更显眼，反而不容易藏住造假的迹象。dch_rec 是 Beneish 在 1999 年 M-Score 模型里使用的变量，本章 LASSO 也独立选中了它，与第二章规则模型的直觉一致。

被剔除的 14 个变量包括 at、lt、sale、ni、cogs、rect、ch_rsst、bm、EBIT 等"看起来很核心"的会计变量。这些变量自身的会计含义并未变弱，它们的信息已经被其它高度相关的变量捕获，比如 ceq、re、soft_assets、reoa 这一组留存收益与软性资产相关变量。LASSO 的稀疏选择是"信息冗余下的代表性挑选"，不能解读成"被剔除的变量与舞弊无关"。

Elastic Net 在相同 CV 下保留了 29 个非零系数，比 LASSO 多一个，整体上选出的变量与 LASSO 重合度极高。Ridge 不做变量选择，全部 42 个系数都非零。

性能评估与案例公司打分

最优 $\lambda$ 选定后，用全部训练数据 1991 至 2002 年重新拟合三个模型，再到测试集 2009 至 2014 年评估。下表列出三种惩罚的测试集性能。

表 4·3　测试集 2009–2014 性能；n = 27,628，阳性 107，1% 名额 k = 277

LASSO 的 AUC 在三种惩罚里最高，达到 0.6876。Elastic Net 紧随其后 0.6872，两者差距已小到 CV 的随机抖动量级；这与 LASSO 与 Elastic Net 选出的变量集合高度重合相符。Ridge 的 AUC 0.6599 落后约 3 个百分点，提示 Bao 数据中存在大量"几乎不携带信息"的弱相关变量，把它们的系数压到零比同时收缩所有系数更划算。Recall@1% 与 Precision@1% 这两个指标在三个模型间差距很小，因为测试集 1% 名额只有 277 个，舞弊总数 107，hits 在 6 至 7 之间波动。

回到 Bao 数据。 本章在测试集上挑出 AUC 最高的 LASSO 作为代表方法，给两家标志性舞弊公司打分。下表列出 Enron 2000 与 Tyco 2000 的预测概率。

表 4·4　LASSO 对标志性案例公司的舞弊概率打分

Enron 2000 被打到 0.7292，远高于全样本舞弊率 0.66%，也远高于第二章 Beneish M-Score 的"被 flag 即可疑"的二元判断。LASSO 抓住了 Enron 在 soft_assets、issue、prcc_f 三个核心变量上同时的极端表现，把它推到了测试集排序的前 1%。Tyco 2000 的概率只有 0.0925，依然高于全样本均值，但远不及 Enron 那么醒目。Tyco 的 2000 年报表在 soft_assets 和 issue 上的偏离没有 Enron 那么夸张，模型相应给出更保守的判断。

Python 实现

R 主线之外，本章配套一份 Python 实现，放在 code/04_penalized.py。Python 端使用 sklearn 的 LogisticRegression，通过 l1_ratio 参数在三种惩罚之间切换：1.0 对应 LASSO，0.0 对应 Ridge，0.5 对应 Elastic Net。求解器在纯 L1 与纯 L2 时使用 liblinear，在 Elastic Net 时使用 saga。其余流程与 R 端保持一致：z-score 标准化器仅在训练集上拟合，时间感知 CV 用 forward-chaining 切分，最优 $C$ 等价于 $1/\lambda$，通过 7 折 CV-AUC 选定。

import numpy as np, pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, ndcg_score

np.random.seed(2026)

# 标准化器仅在训练集上拟合
scaler = StandardScaler()
X_train = scaler.fit_transform(train[features].values)
X_test  = scaler.transform(test[features].values)

# l1_ratio: 1.0 = LASSO, 0.0 = Ridge, 0.5 = Elastic Net
def fit(C, l1_ratio):
    solver = "liblinear" if l1_ratio in (0.0, 1.0) else "saga"
    return LogisticRegression(C=C, l1_ratio=l1_ratio, solver=solver,
                              max_iter=4000, random_state=2026,
                              tol=1e-4).fit(X_train, y_train)

mod_lasso = fit(C_lasso, 1.0)
mod_ridge = fit(C_ridge, 0.0)
mod_enet  = fit(C_enet,  0.5)

Python 端三种惩罚在测试集上的 AUC 分别是：LASSO 0.6882、Ridge 0.6771、Elastic Net 0.6907。与 R 端的 0.6876 / 0.6599 / 0.6872 在小数点后两位上一致，差异在 0.005 以内。Python 端 Elastic Net 略胜 LASSO，R 端则反过来；这种排序的微小翻转源于 sklearn 的 saga 求解器与 glmnet 的坐标下降在 Elastic Net 子问题上的细微数值差异，不影响"惩罚回归整体在 0.68 至 0.69 区间"的结论。

本章累积对比表

表 4·5　方法对比表，截至第 4 章

LASSO 的 AUC 比第三章普通逻辑回归低了约 0.009。在 0.69 这个量级上 0.009 的差距并不大，可以视为 CV 选 $\lambda$ 的随机波动。更重要的是 LASSO 把 42 个变量压到 28 个，模型可解释性显著提升，部署时还能省下不必要的 14 个特征工程。下一章引入决策树与随机森林，把建模思路从"线性 + 惩罚"切到"非参数 + 集成"，看测试集 AUC 能不能突破 0.70 这道墙。

本章知识地图

表 4·6　第 4 章核心概念与常见误解