第 13 章 模型选择、交叉验证与预测评估

实际建模很少只有一个候选方案。同一个响应变量可以采用不同变换，可以加入不同解释变量，也可以用线性项、非线性项或交互项描述均值结构。候选模型越灵活，通常越容易降低训练误差；但模型还会拟合样本中的偶然波动，因此训练误差不能直接代表新数据上的预测误差。

模型选择的关键在于先说明评价目标，再使用与该目标一致的证据。本章讨论三类证据：信息准则衡量似然拟合与复杂度之间的权衡，交叉验证估计指定损失下的样本外表现，独立测试集则用于评价整个建模过程。第 12 章说明了模型怎样被拟合；本章进一步说明，拟合完成后应怎样比较模型，而不把训练拟合误认为泛化能力。

13.1 学习目标

完成本章后，读者应当能够：

  1. 区分训练集、验证集和测试集各自承担的作用；
  2. 解释 AIC、AICc 和 BIC 的构成、适用条件及局限；
  3. 使用同一组数据分折公平比较多个候选模型；
  4. 手写 K 折交叉验证，并正确汇总逐观测预测损失；
  5. 利用杠杆值计算线性模型的留一交叉验证误差；
  6. 根据预测任务选择损失函数和评价指标；
  7. 识别预处理、变量筛选和调参过程中的数据泄漏；
  8. 为时间、空间或分组数据设计符合预测场景的验证方案。

13.2 模型选择的目标与数据角色

“选择模型”至少可能指向三种不同目标。

目标 主要问题 合适的证据
解释与推断 哪些参数及其不确定性有实质含义 研究设计、模型假设、估计稳定性与敏感性分析
结构识别 哪个候选模型最接近设定的数据生成结构 信息准则及其成立条件、模型诊断
预测 哪种方法对未来或未见观测损失更小 与应用场景一致的验证集或测试集表现

这三个目标可以重合，也可能给出不同选择。预测效果相近时，简洁模型通常更容易解释和部署；一个具有明确经济含义的变量，即使只能小幅改善预测，也可能值得保留。反过来，交叉验证选中的变量组合并不自动具有因果含义。变量能否作因果解释主要由研究设计和识别假设决定，而不是由预测分数决定。

模型比较还应明确比较的单位。若最终任务是预测新住户，就应让验证观测对应未参与拟合的新住户；若任务是预测同一企业的未来季度，验证方案必须保留时间顺序。随机拆分观测是否合理，取决于未来预测任务与数据依赖结构，而不是一个固定的软件选项。

目标确定以后，还要区分数据在建模流程中的角色。训练、验证和测试并不是三种不同来源的数据，而是对数据使用方式的约束。

数据角色 用途 可以影响哪些决定
训练集 估计给定模型的参数 系数、预处理量和其他拟合结果
验证集 比较候选模型或选择调节参数 模型形式、变量集合、阈值和超参数
测试集 评价已经完成的建模过程 原则上不再反过来修改模型

K 折交叉验证让不同观测轮流充当验证数据。模型和调节参数确定后，通常在全部训练数据上重新拟合，最后在独立测试集上评价一次。如果查看测试结果后继续改变变量、模型或阈值，这批数据就已经参与了选择，应当称为验证数据；新的无偏评价需要另一批未使用的数据。

13.3 贯穿案例与训练误差的乐观偏差

下面构造一组教学用租金数据。条件均值包含面积、距离的曲线关系、收入和房屋状态；noise1 与noise2 是无关变量。模拟数据使真实结构已知，便于检查各种选择方法的行为，但案例中的具体系数不应解释为真实市场规律。

simulate_rent <- function(n) {
  income <- rlnorm(n, meanlog = log(8), sdlog = 0.5)
  area <- runif(n, min = 30, max = 120)
  distance <- runif(n, min = 1, max = 25)
  new_building <- rbinom(n, size = 1, prob = 0.35)
  noise1 <- rnorm(n)
  noise2 <- runif(n)

  mean_rent <- 4 + 0.075 * area - 0.22 * distance +
    0.009 * distance^2 + 0.40 * log(income) +
    0.70 * new_building

  data.frame(
    rent = mean_rent + rnorm(n, sd = 1.5),
    income = income,
    area = area,
    distance = distance,
    new_building = new_building,
    noise1 = noise1,
    noise2 = noise2
  )
}

set.seed(2026)
rent_train <- simulate_rent(450)
rent_test <- simulate_rent(3000)

candidate_formulas <- list(
  simple = rent ~ area + distance,
  main = rent ~ area + distance + log(income) + new_building,
  curved = rent ~ area + distance + I(distance^2) +
    log(income) + new_building,
  crowded = rent ~ area * distance + I(distance^2) +
    log(income) + new_building + noise1 + noise2
)

candidate_description <- data.frame(
  model = names(candidate_formulas),
  description = c(
    "仅含面积和距离的线性项",
    "加入收入和房屋状态",
    "进一步加入距离二次项",
    "再加入交互项和两个无关变量"
  )
)

knitr::kable(
  candidate_description,
  row.names = FALSE,
  col.names = c("模型", "均值结构")
)
模型 均值结构
simple 仅含面积和距离的线性项
main 加入收入和房屋状态
curved 进一步加入距离二次项
crowded 再加入交互项和两个无关变量

rent_test 在这里已经生成，但在模型选择完成前不会使用。3000 个测试观测远多于实际项目通常能够保留的数量，这是模拟教学中的安排：较大的独立样本可以减少最终性能估计的偶然波动。

这个案例也可以用来说明训练误差为什么不足以评价预测能力。设训练数据为\(\mathcal D=\{(\mathbf x_i,y_i):i=1,\ldots,n\}\)，其中 \(\mathbf x_i\) 是第 \(i\) 个观测的解释变量向量；损失函数写作 \(L\{y,f(\mathbf x)\}\)。给定训练数据拟合出的预测函数记作 \(\hat f_{\mathcal D}\)。训练风险为

\[ \hat R_{\mathrm{train}} =\frac{1}{n}\sum_{i=1}^n L\{y_i,\hat f_{\mathcal D}(\mathbf x_i)\}, \]

而我们真正希望了解的是同一目标总体中新观测 \((\mathbf X^\ast,Y^\ast)\) 的预测风险

\[ R(\hat f_{\mathcal D}) =\operatorname{E}\left[L\{Y^\ast,\hat f_{\mathcal D}(\mathbf X^\ast)\} \mid \mathcal D\right]. \]

训练数据既参与估计又参与评价，使训练风险具有乐观偏差。对于嵌套的普通最小二乘模型，增加变量不会增大残差平方和；即使新增变量与响应没有真实关系，它也可能吸收一部分样本噪声。

candidate_fits <- lapply(
  candidate_formulas,
  lm,
  data = rent_train
)

training_table <- data.frame(
  model = names(candidate_fits),
  coefficients = vapply(candidate_fits, function(fit) {
    length(coef(fit))
  }, integer(1)),
  training_RMSE = vapply(candidate_fits, function(fit) {
    sqrt(mean(residuals(fit)^2))
  }, numeric(1))
)

knitr::kable(
  training_table,
  row.names = FALSE,
  digits = 4,
  col.names = c("模型", "回归系数个数", "训练 RMSE")
)
模型 回归系数个数 训练 RMSE
simple 3 1.437
main 5 1.426
curved 6 1.319
crowded 9 1.317

训练 RMSE 随模型扩展而下降，因此它无法决定何时停止增加变量。模型选择方法要估计或校正的，正是这种重复使用训练数据产生的乐观程度。

13.4 信息准则

13.4.1 AIC、AICc 与 BIC

设模型的最大对数似然为 \(\ell(\widehat{\boldsymbol\theta})\)，估计参数总数为 \(k\)，用于拟合的观测数为\(n\)。Akaike 信息准则为 (Akaike 1974)

\[ \operatorname{AIC} =-2\ell(\widehat{\boldsymbol\theta})+2k. \]

第一项奖励样本内拟合，第二项校正因估计参数而产生的乐观偏差。AIC 的理论目标与样本外对数预测损失有关；它并不是把“拟合优度”和“简单性”任意相加。AIC 只具有相对意义，同一组候选模型中数值较小者得到更多支持，AIC 本身不能检验模型是否正确。

当 \(n\) 与 \(k\) 相比不够大时，常使用小样本修正

\[ \operatorname{AICc} =\operatorname{AIC} +\frac{2k(k+1)}{n-k-1}, \qquad n>k+1. \]

样本量增大时，修正项趋近于 0。不同资料对 AICc 的适用模型和参数计数可能采用略有差异的约定，使用软件结果时应核对其定义。

贝叶斯信息准则（Bayesian information criterion，BIC，也称 Schwarz 准则）为(Schwarz 1978)

\[ \operatorname{BIC} =-2\ell(\widehat{\boldsymbol\theta})+k\log n. \]

在通常的数据规模下，BIC 对新增参数的惩罚比 AIC 更强。BIC 的大样本理论更接近从一组固定候选模型中识别简洁结构，并依赖候选集合、参数维数和模型正则性等条件。AIC 偏重预测信息损失，BIC 偏重模型识别；两者目标不同，出现不同选择并不构成计算错误。BIC 与模型证据的关系将在第19 章讨论。

13.4.2 参数计数与可比性

信息准则中的 \(k\) 是似然中被估计的参数总数，不一定等于回归系数个数。例如，lm() 的高斯似然还要估计误差方差，因此 logLik.lm() 返回的自由度通常是回归系数个数再加 1。直接调用 AIC() 和BIC() 可以减少手工计数不一致，但仍应理解软件采用的似然和常数项约定。

信息准则比较还需要满足以下条件：

  1. 候选模型使用相同的响应变量和同一批观测；
  2. 对数似然来自可以比较的概率模型，并采用一致的常数项约定；
  3. 参数估计已经可靠收敛，数值优化失败时的信息准则没有解释基础；
  4. 比较的是预先定义的候选集合，而不是把大量尝试中最小的数值当作没有选择偏差的证据。

尤其要注意，lm(y ~ x) 与 lm(log(y) ~ x) 的 AIC 一般不能直接比较，因为二者的响应尺度和似然不同；缺失值使候选模型使用不同样本时，也不宜直接比较输出。

log_likelihood <- vapply(candidate_fits, function(fit) {
  as.numeric(logLik(fit))
}, numeric(1))

parameter_count <- vapply(candidate_fits, function(fit) {
  as.numeric(attr(logLik(fit), "df"))
}, numeric(1))

sample_size <- vapply(candidate_fits, nobs, numeric(1))
aic_value <- vapply(candidate_fits, AIC, numeric(1))
bic_value <- vapply(candidate_fits, BIC, numeric(1))
aicc_value <- aic_value +
  2 * parameter_count * (parameter_count + 1) /
  (sample_size - parameter_count - 1)

information_table <- data.frame(
  model = names(candidate_fits),
  k = parameter_count,
  logLik = log_likelihood,
  AIC = aic_value,
  delta_AIC = aic_value - min(aic_value),
  AICc = aicc_value,
  BIC = bic_value,
  delta_BIC = bic_value - min(bic_value)
)

knitr::kable(
  information_table,
  row.names = FALSE,
  digits = 2,
  col.names = c(
    "模型", "$k$", "logLik", "AIC", "$\\Delta$AIC",
    "AICc", "BIC", "$\\Delta$BIC"
  ),
  escape = FALSE
)
模型 \(k\) logLik AIC \(\Delta\)AIC AICc BIC \(\Delta\)BIC
simple 4 -801.6 1611 70.99 1611 1628 58.67
main 6 -798.2 1608 68.18 1609 1633 64.07
curved 7 -763.1 1540 0.00 1541 1569 0.00
crowded 10 -762.3 1545 4.41 1545 1586 16.74

报告差值比报告孤立的 AIC 或 BIC 更有意义。这里的 curved 模型获得最小准则值；crowded 的训练误差略低，但新增交互项和无关变量带来的似然改善不足以抵消复杂度惩罚。这个结论仍只是在给定候选集合中的相对比较，不能代替残差诊断或研究设计。

13.4.3 调整 \(R^2\) 与自适应搜索

在线性模型中，调整决定系数也会对新增参数作校正。若模型含 \(p\) 个回归系数（包括截距），则

\[ R_{\mathrm{adj}}^2 =1-\frac{\operatorname{RSS}/(n-p)} {\operatorname{TSS}/(n-1)}. \]

普通 \(R^2\) 在增加变量后不会下降，调整 \(R^2\) 则可能下降。它适合在相同响应和样本上的线性模型之间描述经自由度校正的样本内拟合，但不是指定损失下的样本外误差，也不能推广为所有模型的统一比较尺度。

信息准则可以评价一个已经拟合的候选模型，却不会消除模型搜索本身的选择偏差。逐步选择每次只考察局部的加入或删除操作，结果可能依赖起始模型和允许搜索的范围，也不保证找到所有候选组合中的最优值。如果先搜索大量变量组合，再把选中模型的常规标准误、\(p\) 值和置信区间当作模型预先确定时的结果，通常会低估选择带来的不确定性。用于解释和推断时，应报告搜索范围与选择规则，并采用独立验证、样本拆分、选择后推断或其他与研究目的相符的方法评估稳定性。

13.5 交叉验证：从一次划分到重复分折

13.5.1 验证集方法与 K 折估计量

最直接的样本外评价是把数据随机分成训练集和验证集：只在训练集拟合，在验证集计算损失。这样保留了评估数据的独立角色，但结果可能对某一次划分十分敏感，而且减少了实际用于拟合的数据量。

K 折交叉验证在数据利用率和计算量之间作出折中。把观测索引划分为互不重叠的集合\(\mathcal I_1,\ldots,\mathcal I_K\)。对第 \(k\) 折，用其余观测拟合 \(\hat f^{(-k)}\)，再对\(i\in\mathcal I_k\) 产生预测\(\hat y_i^{(-k)}\)。交叉验证风险为

\[ \hat R_{\mathrm{CV}} =\frac{1}{n}\sum_{k=1}^K\sum_{i\in\mathcal I_k} L\{y_i,\hat y_i^{(-k)}\}. \]

这个写法按观测汇总损失。若各折大小不完全相同，先计算每折平均损失再简单平均，会给予小折更大的单观测权重。保存每个观测的折外预测，再统一计算损失，可以避免这一问题。

常见的 \(K\) 为 5 或 10。较大的 \(K\) 让每次拟合使用更多数据，却需要更多次拟合，也会改变误差估计的方差。不存在对所有数据都最优的固定折数；应结合样本量、拟合成本和数据结构选择。交叉验证方法的系统讨论可参见 Arlot and Celisse (2010)。

13.5.2 共同分折与折外预测

公平比较要求同一次交叉验证中的所有候选模型使用完全相同的分折。否则，模型差异与随机划分差异混在一起，尤其在模型性能接近时可能改变排序。下面的函数把分折编号作为显式输入，不在函数内部重新抽样。

cv_lm <- function(formula, data, fold_id) {
  n <- nrow(data)
  if (length(fold_id) != n || anyNA(fold_id)) {
    stop("fold_id 必须为每个观测提供一个非缺失的分折编号。")
  }

  fold_levels <- sort(unique(fold_id))
  if (length(fold_levels) < 2L) {
    stop("交叉验证至少需要两个分折。")
  }

  full_frame <- model.frame(formula, data = data, na.action = na.fail)
  y <- model.response(full_frame)
  if (!is.numeric(y) || !is.null(dim(y))) {
    stop("这个教学函数只处理数值型单响应线性模型。")
  }

  prediction <- rep(NA_real_, n)
  for (fold in fold_levels) {
    in_validation <- fold_id == fold
    fit <- lm(
      formula,
      data = data[!in_validation, , drop = FALSE],
      na.action = na.fail
    )
    prediction[in_validation] <- predict(
      fit,
      newdata = data[in_validation, , drop = FALSE]
    )
  }

  if (any(!is.finite(prediction))) {
    stop("折外预测含有非有限值；请检查模型矩阵和各训练折。")
  }

  error <- y - prediction
  list(
    prediction = prediction,
    squared_error = error^2,
    absolute_error = abs(error),
    MSE = mean(error^2),
    RMSE = sqrt(mean(error^2)),
    MAE = mean(abs(error)),
    fold_MSE = tapply(error^2, fold_id, mean)
  )
}

这个函数适合完整、独立同分布观测下的教学示例。若某一训练折的分类变量水平不全，或者模型矩阵在某折秩亏，predict() 可能产生警告或不可估预测；这不是应当隐藏的软件麻烦，而是分折方案或模型设定需要检查的信号。

set.seed(1301)
K <- 5L
fold_id <- sample(rep(seq_len(K), length.out = nrow(rent_train)))

cv_results <- lapply(
  candidate_formulas,
  cv_lm,
  data = rent_train,
  fold_id = fold_id
)

cv_table <- data.frame(
  model = names(cv_results),
  CV_MSE = vapply(cv_results, function(x) x$MSE, numeric(1)),
  CV_RMSE = vapply(cv_results, function(x) x$RMSE, numeric(1)),
  CV_MAE = vapply(cv_results, function(x) x$MAE, numeric(1))
)

knitr::kable(
  cv_table,
  row.names = FALSE,
  digits = 4,
  col.names = c("模型", "CV MSE", "CV RMSE", "CV MAE")
)
模型 CV MSE CV RMSE CV MAE
simple 2.087 1.445 1.153
main 2.077 1.441 1.147
curved 1.765 1.329 1.063
crowded 1.785 1.336 1.064

这一次分折中，curved 的交叉验证误差最低。crowded 在训练集上拟合更紧，却在折外预测中略差，符合无关变量增加估计方差的预期。不同指标若给出不同排序，应回到应用中的损失含义，而不是事后选择最有利于某个模型的指标。

13.5.3 重复交叉验证与比较稳定性

一次随机分折可能偶然有利于某个模型。重复交叉验证使用多组分折，可以观察模型排序对划分的敏感程度。每一次重复仍须让所有模型共享同一分折，这样模型间的误差差异是配对比较。

repeated_cv_lm <- function(formulas, data, K = 5L,
                           repeats = 30L, seed = 1L) {
  set.seed(seed)
  result <- matrix(
    NA_real_,
    nrow = repeats,
    ncol = length(formulas),
    dimnames = list(NULL, names(formulas))
  )

  for (r in seq_len(repeats)) {
    fold_id <- sample(rep(seq_len(K), length.out = nrow(data)))
    for (j in seq_along(formulas)) {
      result[r, j] <- cv_lm(
        formulas[[j]], data, fold_id
      )$RMSE
    }
  }
  result
}

cv_repeated <- repeated_cv_lm(
  candidate_formulas,
  rent_train,
  K = 5,
  repeats = 30,
  seed = 1302
)

reference_model <- names(which.min(colMeans(cv_repeated)))
paired_difference <- sweep(
  cv_repeated,
  1,
  cv_repeated[, reference_model],
  FUN = "-"
)
selected_index <- apply(cv_repeated, 1, which.min)
repeat_summary <- data.frame(
  model = colnames(cv_repeated),
  mean_RMSE = colMeans(cv_repeated),
  partition_SD = apply(cv_repeated, 2, sd),
  mean_paired_difference = colMeans(paired_difference),
  paired_difference_SD = apply(paired_difference, 2, sd),
  times_selected = tabulate(
    selected_index,
    nbins = ncol(cv_repeated)
  )
)

knitr::kable(
  repeat_summary,
  row.names = FALSE,
  digits = 4,
  col.names = c(
    "模型", "平均 CV RMSE", "随机划分标准差",
    "相对参考模型的平均差", "配对差标准差", "成为最优的次数"
  )
)
模型 平均 CV RMSE 随机划分标准差 相对参考模型的平均差 配对差标准差 成为最优的次数
simple 1.450 0.0055 0.1085 0.0060 0
main 1.447 0.0064 0.1059 0.0046 0
curved 1.341 0.0058 0.0000 0.0000 30
crowded 1.349 0.0075 0.0075 0.0043 0 

boxplot(
  cv_repeated,
  ylab = "五折 CV RMSE",
  xlab = "候选模型",
  col = "grey90",
  border = "grey35"
)
points(
  seq_len(ncol(cv_repeated)),
  colMeans(cv_repeated),
  pch = 19,
  col = "firebrick"
)
候选模型在30次五折交叉验证中的RMSE

图13.1: 候选模型在30次五折交叉验证中的RMSE

表中的参考模型（curved）是平均 CV RMSE 最小的候选模型；相对差值为“当前模型减参考模型”，正值表示当前模型折外误差更大。由于各模型共享同一次划分，配对差比两个互不相关的误差标准差更直接地反映模型差异。“随机划分标准差”和“配对差标准差”都只描述固定数据集上更换划分造成的波动，不是总体预测风险估计的标准误。各次重复仍使用同一批观测，不能当作独立重复样本。模型间差异若与配对波动处于同一量级，更稳妥的结论是预测表现接近，再结合简洁性、成本和解释目的选择。

13.5.4 线性模型的留一交叉验证快捷公式

留一交叉验证（leave-one-out cross-validation，LOOCV）取 \(K=n\)，每次删除一个观测。直接实现需要拟合 \(n\) 次模型。对于普通最小二乘，可以利用第 12 章的杠杆值避免重复拟合。

设设计矩阵满列秩，且各 \(h_{ii}<1\)。记全样本残差为 \(\hat e_i=y_i-\hat y_i\)，帽子矩阵的第 \(i\) 个对角元为 \(h_{ii}\)。删除第 \(i\) 个观测后对它产生的预测残差满足

\[ \hat e_{i,(-i)}=y_i-\hat y_{i,(-i)} =\frac{\hat e_i}{1-h_{ii}}. \]

因此平方损失下的留一误差为

\[ \operatorname{CV}_{\mathrm{LOO}} =\frac{1}{n}\sum_{i=1}^n \left(\frac{\hat e_i}{1-h_{ii}}\right)^2. \]

分母解释了高杠杆观测为什么会显著影响留一误差：当 \(h_{ii}\) 接近 1 时，删除该观测后对它的预测可能发生很大变化。下面用逐个重拟合核对快捷公式。

loo_fit <- candidate_fits$curved
loo_shortcut_mse <- mean(
  (residuals(loo_fit) / (1 - hatvalues(loo_fit)))^2
)

loo_prediction <- vapply(
  seq_len(nrow(rent_train)),
  function(i) {
    fit_minus_i <- lm(
      candidate_formulas$curved,
      data = rent_train[-i, , drop = FALSE]
    )
    predict(fit_minus_i, newdata = rent_train[i, , drop = FALSE])
  },
  numeric(1)
)

loo_brute_force_mse <- mean(
  (rent_train$rent - loo_prediction)^2
)

stopifnot(abs(loo_shortcut_mse - loo_brute_force_mse) < 1e-10)

knitr::kable(
  data.frame(
    method = c("杠杆值快捷公式", "逐个删除并重拟合"),
    LOOCV_MSE = c(loo_shortcut_mse, loo_brute_force_mse)
  ),
  row.names = FALSE,
  digits = 8,
  col.names = c("计算方法", "LOOCV MSE")
)
计算方法 LOOCV MSE
杠杆值快捷公式 1.787
逐个删除并重拟合 1.787

两种结果在数值容差内相同，而第一种方法只需一次完整拟合。这个例子体现了计算统计中的一个基本原则：先利用模型的代数结构，再决定是否进行昂贵的重复计算。该公式针对未惩罚普通最小二乘；不能未经推导直接套到 GLM、Lasso 或其他拟合方法上。

13.6 预测损失与评价指标

模型选择依赖损失函数。不同损失强调不同类型的错误，评价指标应在查看模型结果前依据应用目标确定。

13.6.1 连续响应

对验证观测 \((y_i,\hat y_i)\)，均方误差、均方根误差和平均绝对误差分别为

\[ \operatorname{MSE} =\frac{1}{n}\sum_{i=1}^n(y_i-\hat y_i)^2, \qquad \operatorname{RMSE}=\sqrt{\operatorname{MSE}}, \]

\[ \operatorname{MAE} =\frac{1}{n}\sum_{i=1}^n|y_i-\hat y_i|. \]

MSE 对大误差给予平方惩罚，RMSE 与响应变量同量纲，MAE 对极端误差相对不敏感。需要注意，先合并全部折外平方误差再开方，与先计算每折 RMSE 再平均并不相同。本章代码采用前者，使每个观测获得相同权重。

若高估和低估的成本不同，可以使用不对称绝对损失或分位数损失；若不同观测代表的业务规模不同，可以采用预先规定的加权损失。指标应反映决策代价，而不只是沿用软件默认值。

13.6.2 二元结果：概率预测与分类决定

设 \(y_i\in\{0,1\}\)，模型给出的事件概率为 \(\hat p_i\)。Brier 分数 (Brier 1950) 和平均对数损失分别为

\[ \operatorname{Brier} =\frac{1}{n}\sum_{i=1}^n(y_i-\hat p_i)^2, \]

\[ \operatorname{LogLoss} =-\frac{1}{n}\sum_{i=1}^n \left\{y_i\log\hat p_i+(1-y_i)\log(1-\hat p_i)\right\}. \]

两者直接评价概率预测，并且都属于适当评分规则：在期望意义下，诚实报告真实概率最有利(Gneiting and Raftery 2007)。对数损失对极度自信却错误的预测惩罚更重；程序实现时要避免直接对 0 取对数。

准确率、灵敏度和特异度评价的是给定阈值下的分类决定。准确率依赖阈值和类别比例，在低发生率事件中，“全部预测为不发生”也可能得到很高准确率。AUC 衡量阳性观测排序高于阴性观测的能力，但不直接反映概率校准，也不包含具体决策成本。完整报告通常同时考虑概率误差、校准、区分能力和实际阈值下的代价。

下面在独立测试样本上比较只有截距的基准模型与包含风险变量的 Logistic 模型。这里的测试数据只承担评价作用，模型形式已经预先给定。

binary_metrics <- function(y, probability, threshold = 0.5) {
  if (length(y) != length(probability) || length(y) == 0L) {
    stop("y 与 probability 必须等长且不能为空。")
  }
  if (anyNA(y) || !all(y %in% c(0, 1))) {
    stop("y 必须只包含非缺失的 0 和 1。")
  }
  if (any(!is.finite(probability)) ||
      any(probability < 0 | probability > 1)) {
    stop("预测概率必须是 0 与 1 之间的有限数。")
  }
  if (length(threshold) != 1L || !is.finite(threshold) ||
      threshold < 0 || threshold > 1) {
    stop("threshold 必须是 0 与 1 之间的有限标量。")
  }

  eps <- .Machine$double.eps
  probability_safe <- pmin(pmax(probability, eps), 1 - eps)
  classification <- as.integer(probability >= threshold)

  n1 <- sum(y == 1)
  n0 <- sum(y == 0)
  if (n1 == 0L || n0 == 0L) {
    auc <- NA_real_
  } else {
    rank_sum <- sum(rank(probability, ties.method = "average")[y == 1])
    auc <- (rank_sum - n1 * (n1 + 1) / 2) / (n1 * n0)
  }

  sensitivity <- if (n1 == 0L) {
    NA_real_
  } else {
    mean(classification[y == 1] == 1)
  }
  specificity <- if (n0 == 0L) {
    NA_real_
  } else {
    mean(classification[y == 0] == 0)
  }

  c(
    accuracy = mean(classification == y),
    sensitivity = sensitivity,
    specificity = specificity,
    Brier = mean((y - probability)^2),
    log_loss = -mean(
      y * log(probability_safe) +
        (1 - y) * log1p(-probability_safe)
    ),
    AUC = auc
  )
}

set.seed(1310)
n_credit <- 1600
credit_data <- data.frame(
  debt_ratio = rnorm(n_credit),
  income_stability = rnorm(n_credit),
  late_payments = rpois(n_credit, lambda = 1.2)
)
credit_probability <- with(
  credit_data,
  plogis(
    -2.7 + 0.9 * debt_ratio - 0.6 * income_stability +
      0.35 * late_payments
  )
)
credit_data$default <- rbinom(n_credit, 1, credit_probability)

credit_train_id <- sample(seq_len(n_credit), 1100)
credit_train <- credit_data[credit_train_id, ]
credit_test <- credit_data[-credit_train_id, ]

credit_base <- glm(
  default ~ 1,
  family = binomial(),
  data = credit_train
)
credit_full <- glm(
  default ~ debt_ratio + income_stability + late_payments,
  family = binomial(),
  data = credit_train
)

probability_base <- predict(
  credit_base, newdata = credit_test, type = "response"
)
probability_full <- predict(
  credit_full, newdata = credit_test, type = "response"
)

classification_table <- rbind(
  intercept_only = binary_metrics(
    credit_test$default, probability_base
  ),
  risk_variables = binary_metrics(
    credit_test$default, probability_full
  )
)

knitr::kable(
  classification_table,
  digits = 4,
  col.names = c(
    "准确率", "灵敏度", "特异度", "Brier",
    "对数损失", "AUC"
  )
)
准确率 灵敏度 特异度 Brier 对数损失 AUC
intercept_only 0.852 0.0000 1.0000 0.1275 0.4256 0.5000
risk_variables 0.858 0.0541 0.9977 0.1133 0.3655 0.7636

事件发生率较低时，截距模型在 0.5 阈值下也能取得较高准确率，但它不能区分个体风险。包含风险变量的模型主要通过更低的 Brier 分数和对数损失、更高的 AUC 显示改进。若任务是决定哪些申请需要人工审核，阈值还应根据漏判与误报成本确定，不能由 0.5 的习惯值代替。

13.7 交叉验证的边界：泄漏、调参与相关数据

数据泄漏是指拟合过程使用了在真实预测时不可获得的信息。泄漏会使验证误差过于乐观，而且往往不会触发程序报错。交叉验证中的基本原则是：任何从数据中估计的步骤，都必须只在当前训练折中估计，再原样应用到验证折。

操作 错误做法 正确做法
标准化 用全样本均值和标准差后再分折 每折只用训练数据计算中心和尺度
缺失值处理 先用全样本估计填补值 在训练折估计填补规则并应用到验证折
变量筛选 先用全样本相关性筛变量 把筛选步骤放进每个训练折
降维或特征提取 先在全样本做 PCA 每折只在训练数据估计载荷
调节参数选择 用测试集反复选择 \(\lambda\) 在训练数据内部交叉验证，测试集最后使用

像 log(income) 这样不需要从样本估计参数的逐点变换，可以直接应用；中心化、标准化、分位点截断和主成分方向则依赖样本，必须在折内学习。第 14 章选择惩罚参数时，同一原则尤其重要。

13.7.1 调参与嵌套交叉验证

如果用交叉验证选择了模型或调节参数，再用同一交叉验证误差报告最终性能，这个数值也参与了选择，通常带有乐观偏差。数据充足时可以保留独立测试集。没有测试集且需要评价整个调参过程时，可以使用嵌套交叉验证：外层分折估计最终预测性能，内层分折只负责选择模型或超参数。

外层第 \(k\) 折的正确流程是：

  1. 暂时封存外层验证折；
  2. 只在外层训练数据中进行内层交叉验证并选择超参数；
  3. 用选定设置在全部外层训练数据上重新拟合；
  4. 对外层验证折预测；
  5. 汇总所有外层折的预测损失。

嵌套交叉验证估计的是一整套建模程序的表现，而不是某个预先固定模型的表现。它计算成本更高，却避免把调参收益重复计入性能评价。

13.7.2 时间、分组与空间数据

普通随机 K 折交叉验证隐含观测可以交换、训练样本与未来目标来自相同机制等条件。下列数据需要按结构设计分折。

数据结构 主要风险 合适的验证思路
时间序列或面板的未来预测 随机分折把未来信息放入训练集 按时间排序，采用扩展窗口或滚动窗口
同一企业、医院或住户的重复观测 同一单位同时进入训练与验证造成泄漏 以单位为整体分组分折
空间数据 邻近位置高度相似，使随机验证过于乐观 留出空间区块或目标地区
类别极不平衡 某些折几乎没有少数类 在不破坏分组结构的前提下分层分折

验证方案应模拟模型部署后真正遇到的信息边界。例如，用截至某年末的数据预测下一年，训练折就不能包含预测时点以后才会公布的变量或修订值。即使代码完全正确，错误的时间边界仍会产生过于乐观的评价。

13.8 从候选模型到最终报告

13.8.1 在独立测试集上完成案例

租金案例先用信息准则和重复交叉验证比较候选模型，没有查看 rent_test。这里按照平均交叉验证 RMSE选择模型，在全部训练数据上重新拟合，然后只对测试集评价一次。

selected_model <- names(which.min(colMeans(cv_repeated)))
final_fit <- lm(
  candidate_formulas[[selected_model]],
  data = rent_train
)
final_prediction <- predict(final_fit, newdata = rent_test)

final_test_table <- data.frame(
  selected_model = selected_model,
  training_observations = nrow(rent_train),
  test_observations = nrow(rent_test),
  test_RMSE = sqrt(mean((rent_test$rent - final_prediction)^2)),
  test_MAE = mean(abs(rent_test$rent - final_prediction))
)

knitr::kable(
  final_test_table,
  row.names = FALSE,
  digits = 4,
  col.names = c(
    "选定模型", "训练样本量", "测试样本量",
    "测试 RMSE", "测试 MAE"
  )
)
选定模型 训练样本量 测试样本量 测试 RMSE 测试 MAE
curved 450 3000 1.544 1.228

测试误差是对完整选择流程的一次外部检查。实际报告还应说明测试样本如何获得、是否代表目标总体、预测发生在什么时间，以及测试结果的不确定性。一个精确计算出的测试 RMSE 也只对应特定总体、时间范围和损失函数。

13.8.2 报告模型比较的证据

一份可核验的模型比较，需要让读者从报告中恢复下列关键决定。

报告部分 需要说明的问题
预测任务 对谁、在什么时点预测；届时可以取得哪些变量；损失函数为何符合决策目标
候选集合 比较了哪些结构和超参数；候选集合是在什么时候确定的
数据与计算管道 训练、验证和测试数据如何划分；分折如何处理时间、分组或类别比例；预处理、筛选和调参是否都在训练折内完成
比较结果 折外指标及其实际差异；随机划分敏感性；拟合警告与收敛状态
最终评价 测试集是否只使用一次；结论适用于什么总体和时间范围；还有哪些抽样波动与已知局限

信息准则、交叉验证和测试误差提供的是不同层面的证据。它们可以相互印证，却不应被混成一个没有明确目标的“综合评分”。

13.9 本章小结

训练误差因重复使用拟合数据而偏于乐观。AIC 通过参数惩罚估计与对数预测损失有关的相对信息，AICc提供小样本修正，BIC 则采用随样本量增长的复杂度惩罚。信息准则只能比较使用相同数据和可比似然的候选模型，数值较小也不表示模型已经通过诊断。

K 折交叉验证为每个观测产生折外预测，再按预先规定的损失函数评价。多个模型必须共享相同分折，所有数据驱动的预处理和调参步骤必须在训练折内完成。重复交叉验证可以检查结果对随机划分的敏感性；线性模型还可利用残差和杠杆值一次计算留一误差。

预测指标必须服从任务。连续响应的 MSE、RMSE 和 MAE 强调不同误差；二元结果的 Brier 分数与对数损失评价概率预测，准确率等阈值指标评价具体分类决定，AUC 评价排序。时间、分组和空间相关数据不能机械使用普通随机分折。最后保留的测试集评价的是整个建模流程，而不是又一次调参机会。

13.10 思考题

  1. 解释、结构识别和预测三个目标为什么可能选择不同模型。各举一个统计应用说明选择标准。
  2. 为什么嵌套线性模型的训练残差平方和不会随变量增加而上升？这为什么不能证明新增变量提高了预测能力？
  3. AIC 与 BIC 的惩罚项有何差别？当二者选择不同模型时，为什么不能只把数值更小的那个准则称为“更正确”？
  4. 两个回归模型因缺失值而使用不同观测。为什么不应直接比较其 AIC？请提出一种公平的处理方案。
  5. 说明为什么多个候选模型应使用同一组交叉验证分折。若分别随机分折，模型差异中会混入什么？
  6. 推导普通最小二乘的留一残差公式\(\hat e_{i,(-i)}=\hat e_i/(1-h_{ii})\)，并解释高杠杆值对留一预测误差的影响。
  7. 某违约数据的事件率为 2%。一个模型把所有人都预测为不违约，准确率为 98%。为什么这个结果不足以说明模型有用？还应报告哪些指标？
  8. 在预测下一季度销售额时，随机 K 折交叉验证可能造成哪些信息泄漏？请设计一个保留时间顺序的验证方案。
  9. 为什么重复交叉验证结果的标准差不能直接解释为总体预测风险估计的标准误？
  10. 比较独立测试集和嵌套交叉验证的用途。什么情况下更适合使用后者？

13.11 上机实验（Lab）

  1. Lab 1：信息准则与无关变量。 在租金训练数据中继续加入 5、20 和 50 个独立噪声变量，比较训练 RMSE、AIC、AICc、BIC 和五折交叉验证误差。记录各方法从什么时候开始拒绝更复杂的模型。

  2. Lab 2：折数与随机划分。 分别使用 \(K=2,5,10\) 和留一法比较四个租金模型。对前三种设置各重复 50 次，报告计算时间、平均误差和分折间波动，并解释差异。

  3. Lab 3：制造并识别数据泄漏。 构造包含缺失值的数据，比较“全样本填补后交叉验证”与“每个训练折内估计填补值”得到的误差。说明哪一步使用了验证信息，以及偏差方向是否每次都相同。

  4. Lab 4：分类阈值与成本。 把违约训练数据进一步分成拟合集和验证集，令漏判一个违约客户的成本为误报成本的 8 倍。用验证集选择阈值，再到原测试集评价一次，并与 0.5 阈值的结果比较。

  5. 拓展 Lab：滚动时间验证。 生成带趋势、季节性和结构变化的月度序列，比较随机五折验证与扩展窗口验证的预测误差。说明哪一种方案更接近真实的下一期预测任务。

13.12 延伸阅读

AIC 的原始思想见 Akaike (1974)，BIC 见 Schwarz (1978)。Stone (1977) 讨论了交叉验证与 AIC 的渐近联系，Arlot and Celisse (2010) 系统综述了交叉验证的理论与实践。概率预测的评分规则及其解释可参见 Gneiting and Raftery (2007)。

参考文献

Akaike, Hirotugu. 1974. “A New Look at the Statistical Model Identification.” IEEE Transactions on Automatic Control 19 (6): 716–23. https://doi.org/10.1109/TAC.1974.1100705.
Arlot, Sylvain, and Alain Celisse. 2010. “A Survey of Cross-Validation Procedures for Model Selection.” Statistics Surveys 4: 40–79. https://doi.org/10.1214/09-SS054.
Brier, Glenn W. 1950. “Verification of Forecasts Expressed in Terms of Probability.” Monthly Weather Review 78 (1): 1–3.
Gneiting, Tilmann, and Adrian E. Raftery. 2007. “Strictly Proper Scoring Rules, Prediction, and Estimation.” Journal of the American Statistical Association 102 (477): 359–78. https://doi.org/10.1198/016214506000001437.
Schwarz, Gideon. 1978. “Estimating the Dimension of a Model.” The Annals of Statistics 6 (2): 461–64. https://doi.org/10.1214/aos/1176344136.
Stone, Mervyn. 1977. “An Asymptotic Equivalence of Choice of Model by Cross-Validation and Akaike’s Criterion.” Journal of the Royal Statistical Society: Series B (Methodological) 39 (1): 44–47.