LM模型残差bootstrap回归系数置信区间覆盖概率计算咨询

两种计算方案的本质差异 #

两个方案结果差异极大的核心原因是二者的统计定义完全不同，不存在等价关系：

• 第一种方案（区间完全落入原始HC3置信区间的比例）：这个指标衡量的是bootstrap生成的区间和单次原始样本估计值的一致性，和覆盖概率没有任何关联。原始样本计算得到的HC3区间本身就是带抽样误差的单次估计，将其作为刚性标尺去卡所有bootstrap生成的区间，得到的7.6%、40.5%的低数值，本质是在反映bootstrap区间和单次估计的重叠范围，不是对覆盖能力的度量。由于两个区间宽度相近，大部分时候bootstrap区间会略微超出原始区间的边界，这个低结果是完全正常的。
• 第二种方案（区间包含已知真实参数的比例）：这个思路完全匹配覆盖概率的标准定义，是正确的计算方向。

正确的覆盖概率计算流程 #

覆盖概率的核心定义是：针对给定的置信区间构造方法，通过重复抽样生成大量独立的置信区间，计算其中包含真实总体参数值的区间占总重复次数的比例，用来验证方法是否达到标称的置信水平。
结合你最初采用的QR分解高效实现思路，可以完全抛弃每个bootstrap样本都调用lm()的低效操作，复用预计算的QR矩阵大幅提升计算速度：

1. 预计算公共矩阵：对固定的设计矩阵X做QR分解，在所有bootstrap步骤中复用，不需要重复计算
2. 残差预处理：对原始模型的残差做中心化处理（残差bootstrap的标准操作，避免人为引入分布偏差）
3. 循环重抽：每次从中心化残差中有放回抽样，生成bootstrap样本的y值，用预计算的QR矩阵快速求解回归系数
4. 置信区间构造：基于所有bootstrap得到的系数分布，用分位数法/偏差校正法构造对应参数的置信区间
5. 覆盖统计：对每一组构造的置信区间，判断其是否包含已知的真实参数值，计算总覆盖比例

常见认知误区 #

固定单个原始样本做残差重抽得到的覆盖比例，不能直接等价为方法的真实覆盖概率。真实的覆盖概率衡量的是方法在所有可能的样本上的长期表现，需要通过以下流程计算：

1. 从已知参数的总体中，重复生成大量独立的n=100的原始样本
2. 对每个原始样本，执行完整的残差bootstrap流程，生成对应参数的置信区间
3. 统计所有区间中包含真实参数的比例，才是方法的实际覆盖水平
   你当前固定单个样本的计算结果，仅能反映该特定样本下的bootstrap分布特征，会因为原始样本的抽样波动出现和真实覆盖水平的偏差。

当前代码的结果偏差来源 #

你基于lm()计算得到的截距覆盖概率0.842、斜率覆盖概率0.945，和标称95%的置信水平存在偏差，核心原因是：

• 你在每个bootstrap生成的样本上单独计算了同方差假设下的OLS t置信区间，没有利用bootstrap本身的分布信息。标准的残差bootstrap区间应该基于bootstrap得到的系数经验分布构造，而非每个重抽样本独立做t检验。
• 残差重抽过程没有做自由度调整，小样本下区间宽度会偏窄，导致覆盖概率低于标称值。

高效实现参考代码 #

set.seed(42)
n <- 100
true_beta <- c(Intercept=50, Slope=25)
x <- rnorm(n)
e <- rnorm(n)
y <- true_beta[1] + true_beta[2]*x + e
X <- cbind(Intercept=1, x)

# 预计算QR分解，所有bootstrap步骤复用
qrX <- qr(X, LAPACK=TRUE)
# 原始样本系数与中心化残差
beta_hat <- qr.coef(qrX, y)
resid <- y - X %*% beta_hat
resid_centered <- resid - mean(resid)

# 单次bootstrap计算
boot_once <- function() {
  y_boot <- X %*% beta_hat + sample(resid_centered, n, replace=TRUE)
  qr.coef(qrX, y_boot)
}
# 生成1000次bootstrap系数
R <- 1000
boot_coef <- t(replicate(R, boot_once()))

# 构造95%分位数置信区间
ci_intercept <- quantile(boot_coef[,1], c(0.025, 0.975))
ci_slope <- quantile(boot_coef[,2], c(0.025, 0.975))

# 计算单样本下的快速覆盖验证值
cover_intercept <- sum(apply(boot_coef, 1, function(b) {
  se <- sd(boot_coef[,1])
  lower <- b[1] - 1.96*se
  upper <- b[1] + 1.96*se
  lower <= true_beta[1] & true_beta[1] <= upper
}))/R
cover_slope <- sum(apply(boot_coef, 1, function(b) {
  se <- sd(boot_coef[,2])
  lower <- b[2] - 1.96*se
  upper <- b[2] + 1.96*se
  lower <= true_beta[2] & true_beta[2] <= upper
}))/R

内容的提问来源于stack exchange，提问作者cliu