第 18 章 贝叶斯回归模型

回归模型是经济统计中最重要的建模工具之一。第 15 章建立了先验、似然和后验的基本框架，第 16 章介绍了低维后验近似，第 17 章进一步讨论了复杂后验的抽样方法。本章把这些思想和计算工具汇合到回归模型中。与经典回归相比，贝叶斯回归的核心变化是：回归系数不再只用点估计和标准误描述，而是用后验分布描述。这样可以更自然地回答下面这类问题：

  1. 某项消费券政策使人均消费增加的概率有多大？
  2. 政策带动的消费额超过财政补贴面额的概率有多大？
  3. 对一个新家庭或一家新企业，预测值的不确定性有多大？
  4. 当样本量不大或变量较多时，怎样把合理的先验信息转化为稳定估计？

18.1 学习目标

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

  1. 从似然和先验推导高斯线性回归的后验均值与精度；
  2. 根据变量尺度设置系数先验，并用先验预测检查其实际含义；
  3. 用后验样本计算系数、条件效应、政策概率和预测区间；
  4. 在误差方差未知时实现多链 Gibbs 抽样并检查计算质量；
  5. 写出贝叶斯逻辑回归的梯度和曲率，实施 MAP 与 Laplace 近似；
  6. 区分条件均值、未来观测、关联解释和因果解释。

18.2 高斯线性回归：模型、先验与后验

18.2.1 似然与系数先验

考虑线性模型

\[ \mathbf Y\mid\boldsymbol\beta,\sigma^2 \sim N_n(\mathbf X\boldsymbol\beta,\sigma^2\mathbf I_n). \]

其中 \(\mathbf Y\) 是 \(n\) 维随机响应向量，其观测实现记为\(\mathbf y=(y_1,\ldots,y_n)^\top\)；\(\mathbf X\) 是 \(n\times p\)设计矩阵，\(\boldsymbol\beta=(\beta_0,\ldots,\beta_{p-1})^\top\) 是回归系数，\(\sigma^2\) 是误差方差，\(\mathbf I_n\) 是 \(n\) 阶单位矩阵。若 \(\mathbf X\) 第一列为 1，则对应系数是截距项。

为了先看清主线，本节暂把 \(\sigma^2\) 视为已知，并给 \(\boldsymbol\beta\) 指定多元正态先验：

\[ \boldsymbol\beta\sim N_p(\mathbf b_0,\mathbf V_0). \]

其中 \(\mathbf b_0\) 是先验均值向量，\(\mathbf V_0\) 是先验协方差矩阵。这个共轭结果并不需要记忆。把高斯似然与先验的对数密度相加，并略去与 \(\boldsymbol\beta\) 无关的项，得到

\[ \log p(\boldsymbol\beta\mid\mathbf y) = -\frac{1}{2}\boldsymbol\beta^\top \left(\frac{\mathbf X^\top\mathbf X}{\sigma^2}+\mathbf V_0^{-1}\right) \boldsymbol\beta +\boldsymbol\beta^\top \left(\frac{\mathbf X^\top\mathbf y}{\sigma^2}+\mathbf V_0^{-1}\mathbf b_0\right) +C. \]

二次项给出后验精度，线性项给出“精度乘均值”。配方后，后验仍为多元正态分布：

\[ \boldsymbol\beta\mid \mathbf y \sim N_p(\mathbf b_n,\mathbf V_n), \]

其中

\[ \mathbf V_n = \left( \frac{1}{\sigma^2}\mathbf X^\top\mathbf X + \mathbf V_0^{-1} \right)^{-1}, \]

\[ \mathbf b_n = \mathbf V_n \left( \frac{1}{\sigma^2}\mathbf X^\top\mathbf y + \mathbf V_0^{-1}\mathbf b_0 \right). \]

18.2.2 后验精度、收缩与正则化

这里 \(\mathbf V_n\) 是后验协方差矩阵，\(\mathbf b_n\) 是后验均值向量。注意后验精度矩阵\(\mathbf V_n^{-1}\) 可以写成

\[ \mathbf V_n^{-1} = \frac{1}{\sigma^2}\mathbf X^\top\mathbf X+\mathbf V_0^{-1}. \]

这说明后验信息来自两部分：数据提供的精度 \(\mathbf X^\top\mathbf X/\sigma^2\)，以及先验提供的精度\(\mathbf V_0^{-1}\)。如果先验方差很大，\(\mathbf V_0^{-1}\) 很小，后验会更接近普通最小二乘；如果先验方差较小，系数会更明显地向先验均值收缩。在已知 \(\sigma\) 的模型中，\(\mathbf V_n\) 只由设计矩阵、误差尺度和先验决定，观测响应 \(\mathbf y\) 通过 \(\mathbf b_n\) 改变后验位置。

这种收缩与 Ridge 回归有直接联系。若 \(\mathbf b_0=\mathbf 0\) 且 \(\mathbf V_0=\tau^2\mathbf I_p\)，后验众数等价于最小化

\[ \frac{1}{2\sigma^2} (\mathbf y-\mathbf X\boldsymbol\beta)^\top (\mathbf y-\mathbf X\boldsymbol\beta) + \frac{1}{2\tau^2}\boldsymbol\beta^\top\boldsymbol\beta. \]

第一项衡量拟合误差，第二项惩罚系数大小。\(\tau^2\) 越小，惩罚越强。若 Ridge 实现不惩罚截距，对应的贝叶斯模型也应给截距单独设置更宽的先验；不能在没有说明尺度与截距处理的情况下，把两种方法的调节参数逐字对应。

Lasso 也有相应的 MAP 解释。若斜率系数相互独立且

\[ p(\beta_j\mid b)=\frac{1}{2b}\exp\left(-\frac{|\beta_j|}{b}\right), \qquad j=1,\ldots,p-1, \]

则负对数先验贡献为 \(b^{-1}\sum_{j=1}^{p-1}\lvert\beta_j\rvert\)，与高斯负对数似然相加后得到一范数惩罚目标。只有把似然、样本量、误差方差和软件目标函数的缩放统一后，\(b\) 才能与第 14 章的\(\lambda\) 对应。还要注意，出现精确零的是后验众数；连续双指数先验下的完整后验并不会在 0 点产生正概率质量。

18.3 消费券案例：从先验到政策结论

18.3.1 数据、模型与变量尺度

下面使用一个模拟案例。某城市想评估“夜间商圈消费券”是否真的带来增量消费。消费券面额为 15 元，部分居民随机收到消费券。研究者记录一周内夜间商圈消费额，并同时记录收入指数、到商圈距离和当周雨天数。数据是模拟数据，不对应真实调查。

政策评价包含两个不同层次的问题：消费券效应是否为正，以及平均带动消费额是否超过 15 元面额。第二个问题不能只靠回归系数的点估计回答，需要计算效应超过政策阈值的后验概率。

set.seed(1801)
n <- 240

coupon <- rbinom(n, size = 1, prob = 0.5)
income_index <- rnorm(n)
distance <- runif(n, min = 0.3, max = 8)
rainy_days <- rbinom(n, size = 4, prob = 0.35)

baseline <- 62 + 9 * income_index - 1.8 * distance - 3.2 * rainy_days
low_income <- as.integer(income_index < -0.5)
extra_effect <- 15 + 8 * low_income
night_spend <- baseline + coupon * extra_effect + rnorm(n, sd = 18)

coupon_dat <- data.frame(
  night_spend = night_spend,
  coupon = coupon,
  income_index = income_index,
  distance = distance,
  rainy_days = rainy_days,
  low_income = low_income
)

coupon_dat$income_c <- coupon_dat$income_index -
  mean(coupon_dat$income_index)
coupon_dat$distance_c <- coupon_dat$distance - mean(coupon_dat$distance)
coupon_dat$rainy_c <- coupon_dat$rainy_days - mean(coupon_dat$rainy_days)

head(coupon_dat)
#>   night_spend coupon income_index distance rainy_days low_income income_c
#> 1       48.62      0     -0.26605    1.961          0          0  -0.2433
#> 2       30.01      0      0.03238    1.559          0          0   0.0551
#> 3       61.52      0     -0.98459    1.801          1          1  -0.9619
#> 4       73.67      1     -2.25775    5.541          1          1  -2.2350
#> 5       51.91      0     -1.13338    4.265          1          1  -1.1107
#> 6       48.78      1      1.37992    6.570          2          0   1.4026
#>   distance_c rainy_c
#> 1    -2.4984 -1.4292
#> 2    -2.8999 -1.4292
#> 3    -2.6587 -0.4292
#> 4     1.0817 -0.4292
#> 5    -0.1939 -0.4292
#> 6     2.1101  0.5708

设回归模型为

\[ Y_i = \beta_0 +\beta_1\operatorname{coupon}_i +\beta_2\operatorname{income}_i +\beta_3\operatorname{distance}_i +\beta_4\operatorname{rain}_i +\varepsilon_i, \]

\[ \varepsilon_i\overset{\mathrm{iid}}{\sim}N(0,\sigma^2), \qquad i=1,\ldots,n. \]

其中 \(Y_i\) 是第 \(i\) 个居民一周夜间商圈的随机消费额，观测实现记为 \(y_i\)；\(\operatorname{coupon}_i\) 表示是否收到消费券，其余解释变量在进入模型前减去样本均值。中心化不改变相应斜率，却使截距表示一个具有平均协变量、未收到消费券的居民的预期消费，更便于设置先验。系数 \(\beta_1\) 表示在当前无交互项模型中，控制这些变量后收到消费券的平均消费差异。

为了先核对解析公式，本节暂把模拟机制中的 \(\sigma=18\) 当作已知。真实数据不会告诉我们这个数，也不应把同一数据估计出的残差标准差悄悄当成已知；后文将把 \(\sigma^2\) 纳入后验并进行多链抽样。

18.3.2 先验设置与先验预测

“系数服从 \(N(0,10^2)\)”是否宽松，取决于自变量的一单位究竟代表 1 元、1 万元还是一个标准差。本例在原有计量单位上保留斜率，但先中心化连续协变量，再为每个系数分别设置尺度。先验标准差分别表示：截距约 20 元、消费券效应约 20 元、收入指数每单位效应约 15 元、距离每公里效应约 4 元、每个雨天效应约 6 元。

先验预测把这些分散的系数判断放回消费额尺度。下面从观测到的协变量设计中随机取居民画像，但不使用其消费额，再生成先验预测消费。极端预测过多，通常意味着系数先验组合起来过宽；这比逐个查看先验标准差更容易发现问题。

X <- model.matrix(night_spend ~ coupon + income_c + distance_c + rainy_c,
                  data = coupon_dat)
y <- coupon_dat$night_spend
p <- ncol(X)

ols_fit <- lm(night_spend ~ coupon + income_c + distance_c + rainy_c,
              data = coupon_dat)
sigma_known <- 18

b0 <- c(60, 0, 0, 0, 0)
prior_sd <- c(20, 20, 15, 4, 6)
V0 <- diag(prior_sd^2)

rmvnorm_chol <- function(n, mean, Sigma) {
  z <- matrix(rnorm(n * length(mean)), nrow = n)
  sweep(z %*% chol(Sigma), 2, mean, FUN = "+")
}

set.seed(1800)
S_prior <- 4000
beta_prior <- rmvnorm_chol(S_prior, b0, V0)
profile_id <- sample(seq_len(nrow(X)), S_prior, replace = TRUE)
mu_prior <- rowSums(beta_prior * X[profile_id, , drop = FALSE])
y_prior <- rnorm(S_prior, mean = mu_prior, sd = sigma_known)

prior_predictive_summary <- data.frame(
  median = median(y_prior),
  q025 = unname(quantile(y_prior, 0.025)),
  q975 = unname(quantile(y_prior, 0.975)),
  P_below_0 = mean(y_prior < 0),
  P_above_200 = mean(y_prior > 200)
)

knitr::kable(
  prior_predictive_summary, row.names = FALSE, digits = 3,
  caption = "给定观测协变量设计的先验预测检查"
)
表18.1: 给定观测协变量设计的先验预测检查
median q025 q975 P_below_0 P_above_200
60.64 -9.787 130.7 0.048 0

先验预测出现少量负值，反映高斯响应模型允许整个实数轴。若研究对象中零消费很多，或负值具有不可忽略的先验概率，应考虑两部模型、截断分布或对数尺度模型，而不是只靠收紧所有先验来掩盖支持集问题。

18.3.3 已知方差下的解析后验

posterior_precision <- crossprod(X) / sigma_known^2 + solve(V0)
Vn <- chol2inv(chol(posterior_precision))
bn <- solve(
  posterior_precision,
  crossprod(X, y) / sigma_known^2 + solve(V0, b0)
)
post_sd <- sqrt(diag(Vn))

coupon_summary <- data.frame(
  term = colnames(X),
  mean = as.vector(bn),
  sd = post_sd,
  q025 = as.vector(bn) + qnorm(0.025) * post_sd,
  q975 = as.vector(bn) + qnorm(0.975) * post_sd
)

knitr::kable(coupon_summary, row.names = FALSE, digits = 3)
term mean sd q025 q975
(Intercept) 50.333 1.650 47.100 53.566
coupon 18.796 2.315 14.259 23.334
income_c 7.264 1.104 5.100 9.428
distance_c -1.279 0.504 -2.268 -0.291
rainy_c -4.121 1.156 -6.386 -1.856

表中的 coupon 行是本案例最关心的政策系数。后验均值可以像普通回归系数一样解释，但可信区间的含义是：在给定模型、数据和先验后，系数位于该区间内的后验概率约为 95%。这与经典置信区间的解释不同，后者是关于重复抽样程序长期覆盖率的陈述。

18.3.4 从后验样本回答政策问题

贝叶斯回归的一个优点是：许多政策问题可以直接转化为后验概率。若消费券面额为 15 元，则政策部门可以关注

\[ P(\beta_1>0\mid \mathbf y) \]

和

\[ P(\beta_1>15\mid \mathbf y). \]

第一个概率衡量消费券带来正向消费增量的证据，第二个概率衡量平均消费带动额超过券面额的证据。因为\(\boldsymbol\beta\mid\mathbf y\) 是多元正态分布，可以直接模拟后验样本。

set.seed(1802)
S <- 6000
beta_draw <- rmvnorm_chol(S, as.vector(bn), Vn)
colnames(beta_draw) <- colnames(X)

coupon_effect <- beta_draw[, "coupon"]
face_value <- 15
coupon_index <- match("coupon", colnames(X))
coupon_mean_exact <- as.numeric(bn[coupon_index])
coupon_sd_exact <- sqrt(Vn[coupon_index, coupon_index])

policy_prob <- data.frame(
  quantity = c("效应为正的后验概率",
               "效应超过券面额的后验概率",
               "效应的后验均值",
               "效应的后验标准差"),
  simulation = c(mean(coupon_effect > 0),
                 mean(coupon_effect > face_value),
                 mean(coupon_effect),
                 sd(coupon_effect)),
  analytic = c(
    1 - pnorm(0, coupon_mean_exact, coupon_sd_exact),
    1 - pnorm(face_value, coupon_mean_exact, coupon_sd_exact),
    coupon_mean_exact,
    coupon_sd_exact
  )
)

knitr::kable(policy_prob, row.names = FALSE, digits = 3)
quantity simulation analytic
效应为正的后验概率 1.000 1.000
效应超过券面额的后验概率 0.951 0.949
效应的后验均值 18.794 18.796
效应的后验标准差 2.291 2.315

本例中 \(\beta_1\) 的边际后验是正态分布，所以两项概率也可由正态分布函数精确计算。模拟值与解析值的差异只是有限次后验模拟造成的 Monte Carlo 误差。这里仍保留抽样计算，因为同一组样本随后可以处理交互效应、预测和其他没有简洁公式的非线性量；有解析答案时将它作为程序核验基准，是比只看代码能否运行更可靠的习惯。

这些概率比“显著或不显著”更贴近政策讨论。如果 \(P(\beta_1>0\mid\mathbf y)\) 很高，但\(P(\beta_1>15\mid\mathbf y)\) 不高，说明消费券可能确实刺激消费，但平均带动额未必足以覆盖券面额。即使第二个概率很高，也不能直接推出政策具有正净收益：核销率、商户利润、财政筹资成本和消费跨期替代都不在这个回归系数中。

后验样本也可以用来比较不同居民画像下的平均消费。下面比较两个居民：一个住得较近、收入指数较高；另一个住得较远、收入指数较低且雨天较多。我们分别计算收到和未收到消费券时的预测均值差异。

profile_dat <- data.frame(
  coupon = c(0, 1, 0, 1),
  income_index = c(0.8, 0.8, -0.8, -0.8),
  distance = c(1.2, 1.2, 6.5, 6.5),
  rainy_days = c(1, 1, 3, 3)
)

profile_dat$income_c <- profile_dat$income_index -
  mean(coupon_dat$income_index)
profile_dat$distance_c <- profile_dat$distance - mean(coupon_dat$distance)
profile_dat$rainy_c <- profile_dat$rainy_days - mean(coupon_dat$rainy_days)

X_profile <- model.matrix(~ coupon + income_c + distance_c + rainy_c,
                          data = profile_dat)
mu_profile <- beta_draw %*% t(X_profile)

lift_near <- mu_profile[, 2] - mu_profile[, 1]
lift_far <- mu_profile[, 4] - mu_profile[, 3]

profile_summary <- data.frame(
  profile = c("近距离、高收入", "远距离、低收入"),
  mean_lift = c(mean(lift_near), mean(lift_far)),
  q025 = c(quantile(lift_near, 0.025), quantile(lift_far, 0.025)),
  q975 = c(quantile(lift_near, 0.975), quantile(lift_far, 0.975))
)

knitr::kable(profile_summary, row.names = FALSE, digits = 3)
profile mean_lift q025 q975
近距离、高收入 18.79 14.28 23.24
远距离、低收入 18.79 14.28 23.24

在当前线性模型中，消费券效应被设为相同的 \(\beta_1\)，所以两个画像的平均增量相同。如果研究者怀疑低收入或远距离居民的消费券反应不同，应加入交互项，例如 coupon:income_index 或 coupon:distance。两个结果相同来自模型施加的共同效应假设，并不是数据证明了人群之间没有差异。后验解释始终以回归方程所表达的结构为条件。

18.3.5 用交互项表达效应异质性

模拟机制中，收入指数低于 \(-0.5\) 的居民获得了更大的消费增量。主效应模型只能估计跨居民平均后的消费券系数。令 \(L_i=\mathbb I\{\operatorname{income}_i<-0.5\}\)，加入\(\operatorname{coupon}_iL_i\) 后，非低收入组效应为 \(\beta_1\)，低收入组效应为\(\beta_1+\beta_{1L}\)。加入交互项后，\(\beta_1\) 的含义随之变为参照组中的消费券效应。

X_het <- model.matrix(
  night_spend ~ coupon * low_income + income_c + distance_c + rainy_c,
  data = coupon_dat
)

b0_het <- setNames(rep(0, ncol(X_het)), colnames(X_het))
b0_het["(Intercept)"] <- 60
prior_sd_het <- setNames(rep(15, ncol(X_het)), colnames(X_het))
prior_sd_het[c("(Intercept)", "coupon", "income_c",
               "distance_c", "rainy_c", "coupon:low_income")] <-
  c(20, 20, 15, 4, 6, 12)
V0_het <- diag(prior_sd_het^2)

precision_het <- crossprod(X_het) / sigma_known^2 + solve(V0_het)
Vn_het <- chol2inv(chol(precision_het))
bn_het <- solve(
  precision_het,
  crossprod(X_het, y) / sigma_known^2 + solve(V0_het, b0_het)
)

set.seed(1806)
beta_het_draw <- rmvnorm_chol(8000, as.vector(bn_het), Vn_het)
colnames(beta_het_draw) <- colnames(X_het)

effect_nonlow <- beta_het_draw[, "coupon"]
effect_low <- effect_nonlow + beta_het_draw[, "coupon:low_income"]

heterogeneity_summary <- data.frame(
  group = c("非低收入组", "低收入组"),
  mean = c(mean(effect_nonlow), mean(effect_low)),
  q025 = c(quantile(effect_nonlow, 0.025),
           quantile(effect_low, 0.025)),
  q975 = c(quantile(effect_nonlow, 0.975),
           quantile(effect_low, 0.975)),
  P_above_face_value = c(
    mean(effect_nonlow > face_value),
    mean(effect_low > face_value)
  )
)

knitr::kable(
  heterogeneity_summary, row.names = FALSE, digits = 3,
  caption = "交互模型中的分组消费券效应"
)
表18.2: 交互模型中的分组消费券效应
group mean q025 q975 P_above_face_value
非低收入组 16.49 11.10 22.00 0.705
低收入组 23.37 15.88 30.81 0.986

这里的消费券是随机分配的，所以在随机化、无干扰和执行一致等条件下，效应差异可以作因果解释。观察性数据中的回归系数没有这种自动保证；改变先验或使用 MCMC 都不能补回缺失的识别条件。

18.3.6 先验敏感性分析

贝叶斯回归需要报告先验。一个实用做法是比较不同先验下关键结论是否稳定。下面比较三种先验：

  1. weakly_regularizing：消费券效应先验标准差较大，表示弱正则化先验；
  2. skeptical：消费券效应先验均值为 0，标准差较小，表示对政策效果较保守；
  3. positive_prior：消费券效应先验均值为 10，表示先前试点或专家判断认为可能存在正效应。

“弱正则化”比“中性”更准确：任何有信息的先验都会改变模型，只是影响程度不同。

fit_bayes_lm <- function(b0, prior_sd, X, y, sigma) {
  V0 <- diag(prior_sd^2)
  precision <- crossprod(X) / sigma^2 + solve(V0)
  Vn <- chol2inv(chol(precision))
  bn <- solve(
    precision,
    crossprod(X, y) / sigma^2 + solve(V0, b0)
  )
  list(mean = as.vector(bn), cov = Vn)
}

prior_list <- list(
  weakly_regularizing = list(b0 = c(60, 0, 0, 0, 0),
                             sd = c(20, 20, 15, 4, 6)),
  skeptical = list(b0 = c(60, 0, 0, 0, 0),
                   sd = c(20, 6, 15, 4, 6)),
  positive_prior = list(b0 = c(60, 10, 0, 0, 0),
                        sd = c(20, 8, 15, 4, 6))
)

sens <- data.frame()
for (nm in names(prior_list)) {
  fit_s <- fit_bayes_lm(prior_list[[nm]]$b0,
                        prior_list[[nm]]$sd,
                        X, y, sigma_known)
  idx <- which(colnames(X) == "coupon")
  m <- fit_s$mean[idx]
  se <- sqrt(fit_s$cov[idx, idx])
  sens <- rbind(
    sens,
    data.frame(
      prior = nm,
      coupon_mean = m,
      coupon_sd = se,
      prob_gt_0 = 1 - pnorm(0, mean = m, sd = se),
      prob_gt_face_value = 1 - pnorm(face_value, mean = m, sd = se)
    )
  )
}

knitr::kable(sens, row.names = FALSE, digits = 3)
prior coupon_mean coupon_sd prob_gt_0 prob_gt_face_value
weakly_regularizing 18.80 2.315 1 0.949
skeptical 16.55 2.173 1 0.763
positive_prior 18.34 2.238 1 0.932

如果不同先验下政策概率差异很小，说明数据对该结论提供了较强信息。如果差异很大，报告时就应说明结论依赖哪些先验判断。对本科阶段而言，先验敏感性分析比追求某个“唯一正确”的先验更重要。

18.4 未知方差：从条件后验到预测

18.4.1 条件后验与 Gibbs 更新

前面为了得到简洁公式，把 \(\sigma^2\) 当作已知。在更完整的贝叶斯线性回归中，误差方差也可以作为未知参数。设

\[ \boldsymbol\beta\sim N_p(\mathbf b_0,\mathbf V_0), \]

\[ \sigma^2\sim \operatorname{IG}(a_0,d_0), \]

其中 \(\operatorname{IG}(a_0,d_0)\) 表示逆伽马分布，\(a_0\) 是形状参数，\(d_0\) 是尺度参数。这里采用约定：若 \(\sigma^2\sim\operatorname{IG}(a,d)\)，则 \(1/\sigma^2\sim\operatorname{Gamma}(a,d)\)，其中 Gamma 的第二个参数为 rate。

这里的 \(\boldsymbol\beta\) 先验与 \(\sigma^2\) 独立。另一种教材中常见的共轭设定是\(\boldsymbol\beta\mid\sigma^2\sim N_p(\mathbf b_0,\sigma^2\mathbf V_0)\)；它会得到不同的条件后验。推导和编程时必须先确认采用哪一种先验，不能把两套更新公式混用。

在独立正态先验下，可以用 Gibbs 抽样交替更新：

  1. 给定当前 \(\sigma^2\)，从\(\boldsymbol\beta\mid\sigma^2,\mathbf y\) 中抽样；
  2. 给定当前 \(\boldsymbol\beta\)，从\(\sigma^2\mid\boldsymbol\beta,\mathbf y\) 中抽样；
  3. 重复以上步骤，用烧入期减弱初值影响，再对多条正式样本链进行诊断和汇总。

给定 \(\sigma^2\) 时，\(\boldsymbol\beta\) 的条件后验仍是正态分布：

\[ \boldsymbol\beta\mid\sigma^2,\mathbf y \sim N_p(\mathbf b_\sigma,\mathbf V_\sigma), \]

其中

\[ \mathbf V_\sigma = \left( \frac{1}{\sigma^2}\mathbf X^\top\mathbf X+\mathbf V_0^{-1} \right)^{-1}, \]

\[ \mathbf b_\sigma = \mathbf V_\sigma \left( \frac{1}{\sigma^2}\mathbf X^\top\mathbf y+\mathbf V_0^{-1}\mathbf b_0 \right). \]

给定 \(\boldsymbol\beta\) 时，令残差平方和

\[ \operatorname{RSS}(\boldsymbol\beta) = (\mathbf y-\mathbf X\boldsymbol\beta)^\top (\mathbf y-\mathbf X\boldsymbol\beta), \]

则

\[ \sigma^2\mid\boldsymbol\beta,\mathbf y \sim \operatorname{IG} \left( a_0+\frac{n}{2}, d_0+\frac{\operatorname{RSS}(\boldsymbol\beta)}{2} \right). \]

18.4.2 多链实现与计算诊断

取 \(a_0=3\)，并令 \(d_0=(a_0-1)20^2\)，使先验均值\(\operatorname{E}(\sigma^2)=20^2\)。这里的 20 元来自对消费额波动量级的外部判断，不是用同一数据拟合后再当作先验。逆伽马先验在这里主要为了得到可直接抽样的条件后验，并不意味着它对所有尺度参数都是默认选择；形状和尺度取值看似很小，也可能在接近零处产生很强约束。实际建模仍应检查 \(\sigma\) 尺度上的先验和先验预测。下面把一次 Gibbs 运行封装成函数，并从分散初值运行四条链。

a0 <- 3
d0 <- (a0 - 1) * 20^2
V0_inv <- solve(V0)
XtX <- crossprod(X)
Xty <- crossprod(X, y)

run_lm_gibbs <- function(beta0, sigma20, n_iter, seed) {
  set.seed(seed)
  beta_cur <- beta0
  sigma2_cur <- sigma20
  beta_save <- matrix(NA_real_, n_iter, p)
  sigma2_save <- numeric(n_iter)
  colnames(beta_save) <- colnames(X)

  for (s in seq_len(n_iter)) {
    precision <- XtX / sigma2_cur + V0_inv
    R <- chol(precision)
    rhs <- Xty / sigma2_cur + V0_inv %*% b0
    mean_beta <- backsolve(
      R, forwardsolve(t(R), rhs)
    )
    beta_cur <- as.numeric(
      mean_beta + backsolve(R, rnorm(p))
    )

    resid <- y - drop(X %*% beta_cur)
    shape <- a0 + length(y) / 2
    rate <- d0 + sum(resid^2) / 2
    sigma2_cur <- 1 / rgamma(1, shape = shape, rate = rate)

    beta_save[s, ] <- beta_cur
    sigma2_save[s] <- sigma2_cur
  }

  list(beta = beta_save, sigma2 = sigma2_save)
}
S_gibbs <- 6000
burnin_gibbs <- 1000
beta_ols <- coef(ols_fit)
beta_initials <- rbind(
  b0,
  beta_ols,
  beta_ols + c(20, -20, 5, 2, 3),
  beta_ols + c(-20, 20, -5, -2, -3)
)
sigma2_initials <- c(10^2, 18^2, 30^2, 50^2)

gibbs_fits <- lapply(seq_len(4), function(j) {
  run_lm_gibbs(
    beta0 = beta_initials[j, ],
    sigma20 = sigma2_initials[j],
    n_iter = S_gibbs,
    seed = 1810 + j
  )
})

keep_gibbs <- (burnin_gibbs + 1L):S_gibbs
beta_post_chains <- lapply(
  gibbs_fits, function(f) f$beta[keep_gibbs, , drop = FALSE]
)
sigma_post_chains <- lapply(
  gibbs_fits, function(f) sqrt(f$sigma2[keep_gibbs])
)

后续诊断沿用第 17 章定义的 ESS 与秩标准化 \(\widehat R\)。为使本章代码能够独立运行，下一个隐藏代码块保留了相同算法的紧凑实现，正文不再重复展开。正式项目应使用经过充分测试的软件计算秩标准化 \(\widehat R\)、bulk ESS 和 tail ESS。

coupon_chains <- do.call(
  cbind, lapply(beta_post_chains, function(x) x[, "coupon"])
)
sigma_chains <- do.call(cbind, sigma_post_chains)
policy_indicator_chains <- coupon_chains > face_value

diagnostic_series <- list(
  "$\\beta_1$" = coupon_chains,
  "$\\sigma$" = sigma_chains,
  "$\\mathbb I\\{\\beta_1>15\\}$" = policy_indicator_chains
)

gibbs_diagnostics <- do.call(rbind, lapply(
  names(diagnostic_series),
  function(nm) {
    chains <- diagnostic_series[[nm]]
    total_ess <- sum(apply(chains, 2, ess_demo))
    data.frame(
      quantity = nm,
      rank_folded_rhat = rank_folded_rhat_demo(chains),
      total_ess = total_ess,
      mcse_mean = sd(as.numeric(chains)) / sqrt(total_ess)
    )
  }
))

knitr::kable(
  gibbs_diagnostics, row.names = FALSE, digits = 4,
  col.names = c(
    "后验量", "秩折叠 $\\widehat R$", "总 ESS", "均值 MCSE"
  ),
  caption = "消费券 Gibbs 抽样的多链诊断",
  escape = FALSE
)
表18.3: 消费券 Gibbs 抽样的多链诊断
后验量 秩折叠 \(\widehat R\) 总 ESS 均值 MCSE
\(\beta_1\) 1 18568 0.0189
\(\sigma\) 1 17145 0.0070
\(\mathbb I\{\beta_1>15\}\) 1 19285 0.0019 
old_par <- par(no.readonly = TRUE)
par(mfrow = c(2, 1), mar = c(3.5, 4, 2, 1), family = "sans")
chain_col <- c("#0072B2", "#D55E00", "#009E73", "#CC79A7")
matplot(coupon_chains[1:1500, ], type = "l", lty = 1,
        col = chain_col, xlab = "迭代次数", ylab = "消费券系数")
matplot(sigma_chains[1:1500, ], type = "l", lty = 1,
        col = chain_col, xlab = "迭代次数", ylab = expression(sigma))
消费券系数与误差标准差的四条正式样本链。

图18.1: 消费券系数与误差标准差的四条正式样本链。 

invisible(par(old_par))

诊断通过后再合并各链。政策概率使用指标序列的 ESS 和 MCSE，而不是借用消费券系数的 ESS；这是因为 Monte Carlo 精度取决于最终计算的函数。

beta_gibbs <- do.call(rbind, beta_post_chains)
sigma_gibbs <- unlist(sigma_post_chains, use.names = FALSE)

gibbs_summary <- data.frame(
  quantity = c("消费券系数后验均值", "消费券系数2.5%分位数",
               "消费券系数97.5%分位数", "效应超过15的后验概率",
               "误差标准差后验均值"),
  value = c(
    mean(beta_gibbs[, "coupon"]),
    quantile(beta_gibbs[, "coupon"], 0.025),
    quantile(beta_gibbs[, "coupon"], 0.975),
    mean(beta_gibbs[, "coupon"] > face_value),
    mean(sigma_gibbs)
  )
)

knitr::kable(gibbs_summary, row.names = FALSE, digits = 3)
quantity value
消费券系数后验均值 18.722
消费券系数2.5%分位数 13.624
消费券系数97.5%分位数 23.762
效应超过15的后验概率 0.927
误差标准差后验均值 20.065

这里 \(\sigma\) 的后验均值略高于模拟时的 18，不只是“抽样误差”。Gibbs 模型仍然省略了收入组别与消费券的交互，未解释的异质性会进入残差并抬高误差尺度。增加迭代次数可以减小 MCSE，却不能修复模型遗漏。

18.4.3 后验预测与基本模型检查

给定未来协变量 \(\widetilde{\mathbf x}\) 后，\(\widetilde{\mathbf x}^\top\boldsymbol\beta\) 是条件均值；未来消费额还包含新的误差\(\widetilde\varepsilon\sim N(0,\sigma^2)\)。下面同时计算两种区间。

\[ \widetilde y\mid\widetilde{\mathbf x},\boldsymbol\beta,\sigma^2 \sim N(\widetilde{\mathbf x}^\top\boldsymbol\beta,\sigma^2). \]

mu_gibbs <- beta_gibbs %*% t(X_profile)
set.seed(1820)
y_future <- matrix(
  rnorm(
    length(mu_gibbs), mean = as.vector(mu_gibbs),
    sd = rep(sigma_gibbs, times = ncol(mu_gibbs))
  ),
  nrow = nrow(mu_gibbs)
)

profile_labels <- c(
  "近距离高收入：无券", "近距离高收入：有券",
  "远距离低收入：无券", "远距离低收入：有券"
)

predictive_comparison <- data.frame(
  profile = profile_labels,
  mean_response = colMeans(mu_gibbs),
  mean_q025 = apply(mu_gibbs, 2, quantile, 0.025),
  mean_q975 = apply(mu_gibbs, 2, quantile, 0.975),
  predictive_q025 = apply(y_future, 2, quantile, 0.025),
  predictive_q975 = apply(y_future, 2, quantile, 0.975)
)

knitr::kable(
  predictive_comparison, row.names = FALSE, digits = 2,
  caption = "条件均值可信区间与未来消费额预测区间"
)
表18.4: 条件均值可信区间与未来消费额预测区间
profile mean_response mean_q025 mean_q975 predictive_q025 predictive_q975
近距离高收入：无券 62.26 56.85 67.77 22.29 102.49
近距离高收入：有券 80.98 75.29 86.80 41.05 120.89
远距离低收入：无券 35.77 29.48 42.20 -3.34 75.56
远距离低收入：有券 54.50 48.48 60.55 14.57 94.71

预测区间明显更宽，因为它同时包含参数不确定性和个体消费波动。固定 \(\sigma\) 的后验不仅低估方差不确定性，也容易让读者误把均值的不确定性当成未来观测的不确定性。

后验预测还能检查模型是否再现了数据中与研究问题有关的结构。下面从合并后的正式样本中等距选取1000 组参数，在原设计矩阵上生成复制数据 \(\mathbf y^{\mathrm{rep}}\)。除了均值和标准差，还检查低收入组与非低收入组的消费券效应差；这个统计量直接针对主效应模型可能遗漏的异质性。

pp_draw_id <- unique(round(seq(
  1, nrow(beta_gibbs), length.out = 1000
)))
mu_observed_pp <- beta_gibbs[pp_draw_id, , drop = FALSE] %*% t(X)
sigma_observed_pp <- sigma_gibbs[pp_draw_id]

set.seed(1821)
y_rep_observed <- matrix(
  rnorm(
    length(mu_observed_pp),
    mean = as.vector(mu_observed_pp),
    sd = rep(sigma_observed_pp, times = ncol(mu_observed_pp))
  ),
  nrow = length(pp_draw_id)
)

coupon_heterogeneity <- function(z) {
  effect_low <- mean(z[coupon == 1 & low_income == 1]) -
    mean(z[coupon == 0 & low_income == 1])
  effect_nonlow <- mean(z[coupon == 1 & low_income == 0]) -
    mean(z[coupon == 0 & low_income == 0])
  effect_low - effect_nonlow
}

pp_statistics <- cbind(
  mean = rowMeans(y_rep_observed),
  sd = apply(y_rep_observed, 1, sd),
  coupon_heterogeneity = apply(
    y_rep_observed, 1, coupon_heterogeneity
  )
)
observed_statistics <- c(
  mean = mean(y),
  sd = sd(y),
  coupon_heterogeneity = coupon_heterogeneity(y)
)
upper_tail <- colMeans(sweep(
  pp_statistics, 2, observed_statistics, FUN = ">="
))

posterior_predictive_check <- data.frame(
  statistic = names(observed_statistics),
  observed = observed_statistics,
  replicated_q025 = apply(pp_statistics, 2, quantile, 0.025),
  replicated_median = apply(pp_statistics, 2, median),
  replicated_q975 = apply(pp_statistics, 2, quantile, 0.975),
  upper_tail_probability = upper_tail
)

knitr::kable(
  posterior_predictive_check, row.names = FALSE, digits = 3,
  caption = "消费券模型的基本后验预测检查"
)
表18.5: 消费券模型的基本后验预测检查
statistic observed replicated_q025 replicated_median replicated_q975 upper_tail_probability
mean 59.855 56.32 59.846 63.436 0.495
sd 23.615 20.98 23.635 26.790 0.507
coupon_heterogeneity 6.038 -12.16 -1.531 9.085 0.073

后验预测尾部概率衡量复制统计量超过观测统计量的比例，并不是频率学派 \(p\) 值。接近 0 或 1 是模型难以再现相应特征的信号，不是机械的拒绝界线。若均值和总体标准差能够再现，而效应差位于复制分布尾部，就应回到交互项和群体结构，而不是仅增加 Gibbs 迭代次数。第19 章将把这种思想扩展到更系统的模型检查与比较。

18.5 贝叶斯逻辑回归：模型、近似与解释

18.5.1 概率模型与后验

许多经济统计问题的因变量是二元变量，例如是否违约、是否就业、是否参加某项培训、是否接受电子支付。对于二元响应变量，可以使用逻辑回归：

\[ Y_i\mid\mathbf x_i,\boldsymbol\beta \overset{\mathrm{ind}}{\sim}\operatorname{Bernoulli}(p_i), \qquad i=1,\ldots,n, \]

\[ \log\frac{p_i}{1-p_i} = \mathbf x_i^\top\boldsymbol\beta. \]

其中 \(Y_i\) 取 0 或 1，其观测实现记为 \(y_i\)，整组二元观测记为\(\mathbf y=(y_1,\ldots,y_n)^\top\)；\(p_i=P(Y_i=1\mid\mathbf x_i,\boldsymbol\beta)\)，\(\mathbf x_i\) 是解释变量向量。若给 \(\boldsymbol\beta\) 正态先验

\[ \boldsymbol\beta\sim N_p(\mathbf 0,\mathbf V_0), \]

后验密度正比于

\[ p(\boldsymbol\beta\mid \mathbf y) \propto \prod_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i} \times \exp\left( -\frac{1}{2}\boldsymbol\beta^\top\mathbf V_0^{-1}\boldsymbol\beta \right). \]

这个后验一般没有像线性回归那样的闭式正态形式。可以使用第 17 章介绍的随机游走Metropolis；Hamiltonian Monte Carlo 和变分近似也是常用扩展方法，但本书不展开其推导。本章采用MAP 加 Laplace 近似。

18.5.2 企业违约数据与参数尺度

下面构造一个模拟的小微企业贷款数据。变量包括资产负债率、现金流指数、是否年轻企业和近三个月逾期次数。目标是估计企业违约概率。数据仍然是模拟数据，不对应任何真实银行数据。

set.seed(1804)
n_firm <- 450

debt_ratio <- runif(n_firm, min = 0.1, max = 0.9)
cashflow <- rnorm(n_firm)
young <- rbinom(n_firm, size = 1, prob = 0.35)
late_pay <- rpois(n_firm,
                  lambda = exp(-1.1 + 1.5 * debt_ratio - 0.5 * cashflow))

eta_true <- -3.0 + 2.6 * debt_ratio - 0.9 * cashflow +
  0.7 * young + 0.55 * late_pay
default <- rbinom(n_firm, size = 1, prob = plogis(eta_true))

firm_dat <- data.frame(
  default = default,
  debt_ratio = debt_ratio,
  cashflow = cashflow,
  young = young,
  late_pay = late_pay
)

firm_dat$debt_10 <- (firm_dat$debt_ratio - mean(firm_dat$debt_ratio)) / 0.10
firm_dat$cashflow_z <- (firm_dat$cashflow - mean(firm_dat$cashflow)) /
  sd(firm_dat$cashflow)
firm_dat$late_pay_c <- firm_dat$late_pay - mean(firm_dat$late_pay)

mean(firm_dat$default)
#> [1] 0.3

18.5.3 MAP 与 Laplace 近似

为了数值稳定，先定义

\[ \log(1+\exp x) \]

的稳定计算函数。逻辑回归的对数似然可写为

\[ \ell(\boldsymbol\beta) = \sum_{i=1}^n \left[ y_i\eta_i-\log(1+\exp(\eta_i)) \right], \]

其中 \(\eta_i=\mathbf x_i^\top\boldsymbol\beta\)。再加上正态先验的对数密度，就得到对数后验。

记 \(\mathbf p=(p_1,\ldots,p_n)^\top\)，\(\mathbf D=\operatorname{diag}(s_0^2,\ldots,s_{p-1}^2)\) 为先验协方差，\(\mathbf W=\operatorname{diag}\{p_i(1-p_i)\}\)。零均值正态先验下，对数后验的梯度和负Hessian 分别是

\[ \nabla\log p(\boldsymbol\beta\mid\mathbf y) =\mathbf X^\top(\mathbf y-\mathbf p)-\mathbf D^{-1}\boldsymbol\beta, \]

\[ -\nabla^2\log p(\boldsymbol\beta\mid\mathbf y) =\mathbf X^\top\mathbf W\mathbf X+\mathbf D^{-1}. \]

这两个量既可加速优化，也可核验程序；负 Hessian 在 MAP 处还是 Laplace 协方差的精度矩阵。

本例把资产负债率的单位定义为 10 个百分点，把现金流变为一个样本标准差，并中心化逾期次数。相应地，截距先验标准差设为 2.5，各斜率先验标准差设为 1。这样的先验不是“无信息”：它允许一个单位变化使赔率显著改变，同时排除数量级极端的系数。优化之前，先把先验传播到观测到的企业画像上，查看它对违约概率意味着什么。

X_logit <- model.matrix(default ~ debt_10 + cashflow_z + young + late_pay_c,
                        data = firm_dat)
y_logit <- firm_dat$default
prior_sd_logit <- c(2.5, 1, 1, 1, 1)

set.seed(1803)
S_logit_prior <- 5000
beta_logit_prior <- sweep(
  matrix(rnorm(S_logit_prior * ncol(X_logit)), nrow = S_logit_prior),
  2, prior_sd_logit, FUN = "*"
)
prior_profile_id <- sample(
  seq_len(nrow(X_logit)), S_logit_prior, replace = TRUE
)
prior_default_prob <- plogis(rowSums(
  beta_logit_prior * X_logit[prior_profile_id, , drop = FALSE]
))

logit_prior_check <- data.frame(
  median = median(prior_default_prob),
  q025 = quantile(prior_default_prob, 0.025),
  q975 = quantile(prior_default_prob, 0.975),
  P_below_001 = mean(prior_default_prob < 0.01),
  P_above_099 = mean(prior_default_prob > 0.99)
)

knitr::kable(
  logit_prior_check, row.names = FALSE, digits = 3,
  caption = "企业画像上的先验预测违约概率"
)
表18.6: 企业画像上的先验预测违约概率
median q025 q975 P_below_001 P_above_099
0.508 0.001 0.999 0.107 0.099

先验预测若几乎全部挤在 0 和 1 附近，通常表示多个系数叠加后过于极端；若领域知识明确认为违约是低概率事件，则还应相应移动截距先验，而不是依赖数据事后纠正。这里使用样本协变量只是规定要在哪些企业画像上检查先验，没有使用响应变量 default。


log1pexp <- function(x) {
  pmax(x, 0) + log1p(exp(-abs(x)))
}

logpost_logit <- function(beta) {
  eta <- drop(X_logit %*% beta)
  loglik <- sum(y_logit * eta - log1pexp(eta))
  logprior <- sum(dnorm(beta, mean = 0, sd = prior_sd_logit, log = TRUE))
  loglik + logprior
}

gradient_logpost_logit <- function(beta) {
  eta <- drop(X_logit %*% beta)
  drop(crossprod(X_logit, y_logit - plogis(eta))) -
    beta / prior_sd_logit^2
}

negative_hessian_logpost_logit <- function(beta) {
  prob <- plogis(drop(X_logit %*% beta))
  weight <- prob * (1 - prob)
  crossprod(X_logit * sqrt(weight)) +
    diag(1 / prior_sd_logit^2)
}

logit_starts <- rbind(
  rep(0, ncol(X_logit)),
  rep(0.5, ncol(X_logit)),
  rep(-0.5, ncol(X_logit))
)

logit_fits <- lapply(seq_len(nrow(logit_starts)), function(j) {
  optim(
    logit_starts[j, ],
    fn = function(b) -logpost_logit(b),
    gr = function(b) -gradient_logpost_logit(b),
    method = "BFGS",
    control = list(reltol = 1e-11)
  )
})

logit_optim_check <- data.frame(
  start = seq_along(logit_fits),
  convergence = vapply(logit_fits, `[[`, integer(1), "convergence"),
  objective = vapply(logit_fits, `[[`, numeric(1), "value"),
  max_abs_gradient = vapply(
    logit_fits,
    function(f) max(abs(gradient_logpost_logit(f$par))),
    numeric(1)
  )
)

knitr::kable(
  logit_optim_check, row.names = FALSE, digits = 7,
  caption = "逻辑回归 MAP 的多初值优化检查"
)
表18.7: 逻辑回归 MAP 的多初值优化检查
start convergence objective max_abs_gradient
1 0 224.6 0.0000066
2 0 224.6 0.0000106
3 0 224.6 0.0001239

三个初值应到达相同目标值，且梯度接近零。convergence = 0 不能替代这些检查。正态先验还使后验在完全分离或近似分离时保持为适当分布，但先验尺度必须结合变量单位解释：debt_10 增加 1 表示资产负债率增加 10 个百分点，late_pay_c 增加 1 表示多一次逾期。

在 \(\widehat{\boldsymbol\beta}_{\mathrm{MAP}}\) 附近作二阶展开，得到

\[ \boldsymbol\beta\mid\mathbf y \ \dot\sim\ N_p\left(\widehat{\boldsymbol\beta}_{\mathrm{MAP}}, \left[-\nabla^2\log p(\widehat{\boldsymbol\beta}_{\mathrm{MAP}}\mid\mathbf y)\right]^{-1} \right). \]

fit_logit <- logit_fits[[1L]]
beta_map <- fit_logit$par
precision_laplace <- negative_hessian_logpost_logit(beta_map)
stopifnot(min(eigen(precision_laplace, symmetric = TRUE,
                    only.values = TRUE)$values) > 0)
cov_laplace <- chol2inv(chol(precision_laplace))

set.seed(1805)
beta_logit_draw <- rmvnorm_chol(8000, beta_map, cov_laplace)
colnames(beta_logit_draw) <- colnames(X_logit)

logit_summary <- data.frame(
  term = colnames(X_logit),
  mean = colMeans(beta_logit_draw),
  q025 = apply(beta_logit_draw, 2, quantile, 0.025),
  q975 = apply(beta_logit_draw, 2, quantile, 0.975)
)

knitr::kable(logit_summary, row.names = FALSE, digits = 3)
term mean q025 q975
(Intercept) -1.160 -1.457 -0.852
debt_10 0.250 0.141 0.358
cashflow_z -0.912 -1.187 -0.628
young 0.271 -0.210 0.749
late_pay_c 0.352 0.109 0.594

18.5.4 用重要性抽样检查局部近似

Laplace 近似只使用 MAP 附近的二阶形状。为检查它是否遗漏明显的偏斜或尾部，本例再用放大 1.5 倍标准差的多元正态作为建议分布。重要性抽样仍只需要未归一化后验；若权重 ESS 足够大、最大权重不突出，而且校正前后的摘要接近，才有理由认为 Laplace 在本例中工作良好。

log_dmvnorm <- function(x, mean, Sigma) {
  x <- as.matrix(x)
  centered <- sweep(x, 2, mean)
  R <- chol(Sigma)
  z <- t(backsolve(R, t(centered), transpose = TRUE))
  -0.5 * (
    ncol(x) * log(2 * pi) + 2 * sum(log(diag(R))) +
      rowSums(z^2)
  )
}

weighted_quantile_reg <- function(x, w, probs) {
  ord <- order(x)
  x <- x[ord]
  w <- w[ord] / sum(w)
  cw <- cumsum(w)
  vapply(probs, function(p) x[which(cw >= p)[1L]], numeric(1))
}

proposal_cov <- 1.5^2 * cov_laplace
set.seed(1807)
S_logit_is <- 12000
beta_logit_is <- rmvnorm_chol(S_logit_is, beta_map, proposal_cov)
colnames(beta_logit_is) <- colnames(X_logit)

log_target_is <- apply(beta_logit_is, 1, logpost_logit)
log_q_is <- log_dmvnorm(beta_logit_is, beta_map, proposal_cov)
log_weight_is <- log_target_is - log_q_is
weight_is <- exp(log_weight_is - max(log_weight_is))
weight_is <- weight_is / sum(weight_is)

is_ess <- 1 / sum(weight_is^2)
logit_is_diagnostic <- data.frame(
  draws = S_logit_is,
  weight_ess = is_ess,
  ess_fraction = is_ess / S_logit_is,
  max_weight = max(weight_is)
)

knitr::kable(
  logit_is_diagnostic, row.names = FALSE, digits = 4,
  caption = "逻辑回归 Laplace 建议分布的重要性权重诊断"
)
表18.8: 逻辑回归 Laplace 建议分布的重要性权重诊断
draws weight_ess ess_fraction max_weight
12000 4847 0.4039 6e-04 
logit_approximation_check <- data.frame(
  term = colnames(X_logit),
  Laplace_mean = colMeans(beta_logit_draw),
  importance_mean = colSums(beta_logit_is * weight_is)
)
logit_approximation_check$absolute_difference <- abs(
  logit_approximation_check$Laplace_mean -
    logit_approximation_check$importance_mean
)

knitr::kable(
  logit_approximation_check, row.names = FALSE, digits = 4,
  caption = "Laplace 摘要与重要性校正结果的比较"
)
表18.9: Laplace 摘要与重要性校正结果的比较
term Laplace_mean importance_mean absolute_difference
(Intercept) -1.1598 -1.1758 0.0160
debt_10 0.2500 0.2532 0.0033
cashflow_z -0.9115 -0.9305 0.0190
young 0.2710 0.2754 0.0044
late_pay_c 0.3520 0.3597 0.0077

重要性抽样在这里是低维核验工具，不是高维逻辑回归的通用答案。维数增加后，权重可能迅速退化；此时应转向第 17 章的抽样思想，或使用实现了更适合高维问题之抽样算法的成熟软件。

18.5.5 从系数到违约概率

逻辑回归系数表示对数赔率的变化。为了便于解释，常把系数指数化为赔率比。若某变量的系数为\(\beta_j\)，则该变量增加 1 个单位时，违约赔率乘以 \(\exp(\beta_j)\)，其他变量保持不变。

odds_draw_is <- exp(beta_logit_is[, -1, drop = FALSE])
odds_summary <- data.frame(
  term = colnames(X_logit)[-1],
  odds_ratio_mean = colSums(odds_draw_is * weight_is),
  odds_ratio_q025 = apply(
    odds_draw_is, 2, weighted_quantile_reg,
    w = weight_is, probs = 0.025
  ),
  odds_ratio_q975 = apply(
    odds_draw_is, 2, weighted_quantile_reg,
    w = weight_is, probs = 0.975
  )
)

knitr::kable(odds_summary, row.names = FALSE, digits = 3)
term odds_ratio_mean odds_ratio_q025 odds_ratio_q975
debt_10 1.290 1.156 1.439
cashflow_z 0.398 0.296 0.521
young 1.358 0.816 2.150
late_pay_c 1.444 1.129 1.839

表中的 debt_10 是资产负债率增加 10 个百分点的赔率比，cashflow_z 是现金流指数增加一个样本标准差的赔率比，late_pay_c 增加 1 则表示多一次逾期。赔率比不是风险比：当基准违约概率较高时，二者数值可能相差很大。

最后计算两个企业画像的违约概率后验分布：一家财务状况较稳健的企业，和一家负债率高、现金流偏弱且近期多次逾期的企业。

risk_dat <- data.frame(
  debt_ratio = c(0.35, 0.78),
  cashflow = c(0.6, -0.8),
  young = c(0, 1),
  late_pay = c(0, 3)
)

risk_dat$debt_10 <- (risk_dat$debt_ratio - mean(firm_dat$debt_ratio)) / 0.10
risk_dat$cashflow_z <- (risk_dat$cashflow - mean(firm_dat$cashflow)) /
  sd(firm_dat$cashflow)
risk_dat$late_pay_c <- risk_dat$late_pay - mean(firm_dat$late_pay)

X_risk <- model.matrix(~ debt_10 + cashflow_z + young + late_pay_c,
                       data = risk_dat)
risk_draw <- plogis(beta_logit_is %*% t(X_risk))

risk_summary <- data.frame(
  profile = c("财务稳健企业", "财务承压企业"),
  mean_default_prob = colSums(risk_draw * weight_is),
  q025 = apply(
    risk_draw, 2, weighted_quantile_reg,
    w = weight_is, probs = 0.025
  ),
  q975 = apply(
    risk_draw, 2, weighted_quantile_reg,
    w = weight_is, probs = 0.975
  ),
  prob_default_gt_20pct = colSums((risk_draw > 0.20) * weight_is)
)

knitr::kable(risk_summary, row.names = FALSE, digits = 3)
profile mean_default_prob q025 q975 prob_default_gt_20pct
财务稳健企业 0.085 0.053 0.126 0
财务承压企业 0.788 0.664 0.883 1

这个结果比单纯给出“违约/不违约”分类更有信息。信贷风控中，概率本身往往比硬分类更有用，因为银行可以结合风险偏好、贷款利率、抵押品和监管要求设定不同阈值。

这里的区间描述企业条件违约概率的不确定性。对下一家具有相同画像的企业，违约结果仍然只有 0 或 1，其后验预测违约概率等于表中的 mean_default_prob。此外，这些数据是观察性特征，系数和风险差异首先是条件关联；不能把“降低资产负债率一个单位”直接解释成由该系数保证的因果效应。

两个画像的结果也不是模型有效性的证明。投入实际风控前，还应检查不同预测风险组中的实际违约率、复制数据能否再现违约数量与关键分组差异，并在未参与拟合的数据上评价预测。后验区间只传播给定模型内部的不确定性，不会自动包含变量遗漏、样本选择变化或未来经济环境变化带来的误差。

18.6 从研究问题到可报告结果

贝叶斯回归不是在经典回归拟合之后再附加一个先验，而是一条从研究问题到概率结论的完整推理链。每一环都产生下一环所需的对象，也各有不能被后续计算补救的失败方式。

环节 需要明确的问题 本章中的检查
建模 响应分布、线性预测子和目标效应是什么 设计矩阵、交互项和条件解释
先验 系数的一单位变化代表什么 变量尺度、先验预测和敏感性分析
计算 后验是否可解析，近似或抽样是否可信 解析基准、梯度、Hessian、多链诊断或权重诊断
解释 哪个后验函数直接回答研究问题 政策概率、赔率比、画像风险和预测区间
模型检查 模型遗漏和分布假设会怎样影响结论 残差尺度、预测表现及下一章的后验预测检查

这一顺序不能任意颠倒。设计矩阵的列不满秩，不能靠增加 MCMC 迭代解决；先验在响应尺度上明显荒谬，也不能因为后验均值看似合理就忽略。反过来，模型和先验合理并不保证计算正确，因此应在有解析答案时建立基准，在非共轭模型中检查优化、局部曲率或抽样诊断。最后报告的不应只是系数表，还应包括研究问题对应的概率或预测、先验敏感性、数值误差以及不能由数据识别的因果限制。

18.7 本章小结

贝叶斯回归用后验分布而不是单一估计值描述系数不确定性。已知方差的高斯线性模型具有解析正态后验，其后验精度是数据精度与先验精度之和；误差方差未知时，可以用 Gibbs 抽样联合传播系数、方差和预测不确定性。变量中心化和尺度定义决定了系数先验的实际含义，先验预测则把这些判断还原到响应尺度。

交互项决定“效应”针对哪个群体，预测区间与条件均值可信区间回答不同问题。对于逻辑回归等非共轭模型，MAP 与 Hessian 给出 Laplace 近似，但仍需用多初值、曲率和重要性权重检查计算。高斯系数先验还连接到第 14 章的 Ridge 型 MAP；贝叶斯分析进一步保留了后验不确定性。

18.8 思考题

  1. 为什么同一个 \(N(0,10^2)\) 先验用于“收入（元）”和“收入（万元）”会表达完全不同的信息？
  2. 已知 \(\sigma\) 的解析后验与把估计得到的 \(\hat\sigma\) 当作已知量有什么区别？后者遗漏了什么？
  3. 为什么主效应模型中的消费券系数不能回答低收入组与高收入组的效应是否相同？
  4. 后验标准差、预测标准差与 MCMC 的 MCSE 分别描述哪一种不确定性？
  5. 正态系数先验为什么能缓解逻辑回归中的完全分离？这是否意味着任意宽的先验都同样有效？
  6. 赔率比为 2 是否意味着违约概率也变为原来的 2 倍？请用两个不同的基准概率说明。
  7. Laplace 近似与重要性校正结果接近，能否证明后验近似完全正确？还应检查哪些量？
  8. 为什么预测准确、后验概率很高或计算诊断良好，都不能单独建立因果解释？

18.9 上机实验（Lab）

  1. Lab 1：先验预测。 分别把消费券系数的先验标准差设为 5、20 和 100，比较先验预测范围、负消费概率以及后验政策概率，说明先验尺度怎样传播到观测尺度。

  2. Lab 2：异质效应。 把低收入阈值改为 \(-1\)、\(0\)，重新构造交互项。报告两个群体的样本量、效应区间和\(P(\beta_1+\beta_{1L}>15\mid\mathbf y)\)，讨论阈值选择与不确定性的关系。

  3. Lab 3：多链 Gibbs。 改变四条链的初值和烧入长度，比较轨迹、秩标准化 \(\widehat R\)、ESS 与 MCSE。不得只根据后验均值相近宣布收敛。

  4. Lab 4：预测区间。 为消费券案例构造三个新居民画像，分别报告条件均值可信区间与未来消费预测区间，并解释两者宽度不同的原因。

  5. Lab 5：逻辑回归尺度。 改用未经变换的 debt_ratio 拟合模型。相应调整先验，使它与 debt_10 模型表达相同信息，再比较两个参数化下的风险预测。

  6. Lab 6：近似失效。 把违约样本量减小到 40，比较 Laplace 与重要性校正后的系数和风险概率，同时检查权重 ESS 与最大权重。

  7. 拓展 Lab：决策阈值。 假设漏判一次违约的损失是误拒一家企业的 8 倍，推导相应概率阈值，并比较它与固定 20% 阈值下的决策。

18.10 延伸阅读

贝叶斯回归的建模、先验预测和后验解释可参见 Gelman et al. (2013)。关于后验计算、Laplace 近似与Monte Carlo 误差，可参见 Robert and Casella (2004) 和 Monahan (2011)。下一章将从参数解释转向后验预测检查与模型比较，进一步判断回归模型能否再现数据中的重要特征。

参考文献

Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2013. Bayesian Data Analysis. 3rd ed. Boca Raton, FL: CRC Press.
Monahan, John F. 2011. Numerical Methods of Statistics. 2nd ed. Cambridge: Cambridge University Press.
Robert, Christian P., and George Casella. 2004. Monte Carlo Statistical Methods. 2nd ed. New York: Springer.