使用R中quantreg::rq时p值不稳定的问题及解决咨询

问题背景 #

使用quantreg::rq对时间序列回测的特征做分位数回归，基于p值≤5%筛选特征。多次运行程序时，特征的beta系数始终一致，但p值却不稳定——例如某特征第一次p值为0.049被选中，第二次变为0.0536被剔除。已确认输入数据无变化，设置set.seed()后问题仍存在，怀疑与显著性水平计算时的随机抽样有关。

可复现代码：

require(quantreg)
data(engel)
set.seed(12345)

QR_taus <- c(.10, .2, 0.33, .50, 0.66, .80, .90)

for (i in 1:5){
  mod <- rq(foodexp ~ income,
            tau = QR_taus,
            data = engel)
  summ <- summary(mod, se = "boot")
  
  Residuals_QR <- summ[[1]][["coefficients"]]
  assign(paste("Residuals_QR_",i,sep=""),Residuals_QR)
}

Residuals_QR_1
Residuals_QR_2
Residuals_QR_3
Residuals_QR_4
Residuals_QR_5

对比结果可见，系数一致，但标准误、t值和p值每次运行都存在差异。

问题原因 #

核心原因是**summary.rq中se="boot"参数采用了自助法（bootstrap）抽样，而随机种子的设置未覆盖自助抽样的随机过程**：

• rq()本身是确定性算法，因此回归系数始终稳定；
• 自助法需要反复从原始数据中随机抽取样本以估计标准误，这个过程的随机性如果未被固定，每次运行都会产生不同的标准误，进而导致p值波动；
• 全局设置的set.seed(12345)仅在循环前生效一次，而每次调用summary(mod, se="boot")时，自助抽样会重新启动随机序列，因此即使设置了全局种子，结果仍不稳定。

解决方法 #

方法1：将set.seed()放入循环内部，固定每次自助抽样的随机种子 #

在每次调用summary.rq前设置相同的随机种子，确保每次自助抽样的随机过程完全一致：

require(quantreg)
data(engel)

QR_taus <- c(.10, .2, 0.33, .50, 0.66, .80, .90)

for (i in 1:5){
  set.seed(12345)  # 循环内每次都固定种子
  mod <- rq(foodexp ~ income,
            tau = QR_taus,
            data = engel)
  summ <- summary(mod, se = "boot")
  
  Residuals_QR <- summ[[1]][["coefficients"]]
  assign(paste("Residuals_QR_",i,sep=""),Residuals_QR)
}

# 此时所有Residuals_QR_*的结果完全一致

方法2：使用非自助法的标准误计算方式 #

如果不需要自助法的标准误，可选择quantreg支持的确定性标准误计算方法，避免随机波动：

• se="iid"：基于独立同分布假设的渐近标准误；
• se="nid"：基于正态独立同分布假设的渐近标准误；
• se="ker"：核估计的渐近标准误。

示例代码：

require(quantreg)
data(engel)
set.seed(12345)

QR_taus <- c(.10, .2, 0.33, .50, 0.66, .80, .90)

for (i in 1:5){
  mod <- rq(foodexp ~ income,
            tau = QR_taus,
            data = engel)
  summ <- summary(mod, se = "iid")  # 使用确定性的iid标准误
  
  Residuals_QR <- summ[[1]][["coefficients"]]
  assign(paste("Residuals_QR_",i,sep=""),Residuals_QR)
}

# 所有运行结果完全一致

方法3：优化自助法参数（需结合种子固定） #

若必须使用自助法，可通过以下方式提升结果稳定性：

1. 增加自助抽样次数（B参数）：次数越多，标准误和p值的波动越小，但计算时间会增加；
2. 选择更稳定的自助类型（boot.type参数），如"perc"（百分位自助法）；
3. 必须配合循环内的种子固定，才能保证结果完全可复现。

示例代码：

require(quantreg)
data(engel)

QR_taus <- c(.10, .2, 0.33, .50, 0.66, .80, .90)

for (i in 1:5){
  set.seed(12345)
  mod <- rq(foodexp ~ income,
            tau = QR_taus,
            data = engel)
  summ <- summary(mod, se = "boot", 
                  boot.type = "perc",  # 选择百分位自助法
                  B = 1000)  # 增加抽样次数，默认100
  
  Residuals_QR <- summ[[1]][["coefficients"]]
  assign(paste("Residuals_QR_",i,sep=""),Residuals_QR)
}

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