如何为quantreg::rq中分类变量获取全局p值

针对你用quantreg::rq拟合的中位数回归模型，要获取4水平分类变量Education level的全局p值，以下是几种可行方案：

方法1：使用quantreg内置的秩检验 #

quantreg包的anova.rq函数原生支持rq对象的全局效应检验，默认的秩检验适配分位数回归的特性，通过对比全模型与仅含截距的零模型即可得到全局p值：

library(quantreg)

# 拟合全模型
qr <- rq(Total ~ `Education level`, data = df, tau = 0.5)
# 拟合仅含截距的零模型
qr_null <- rq(Total ~ 1, data = df, tau = 0.5)

# 执行秩检验，输出全局p值
anova(qr_null, qr, test = "rank")

结果中的Pr(>|tau|)即为Education level的全局p值。

方法2：Wald检验（借助car::linearHypothesis） #

虽然car::Anova不兼容rq对象，但car::linearHypothesis可直接对rq模型的系数做联合检验，构造假设为分类变量所有非基准水平的系数同时为0：

library(car)

# 拟合模型
qr <- rq(Total ~ `Education level`, data = df, tau = 0.5)

# 构造联合假设：分类变量3个非基准水平系数全为0
hypotheses <- paste0("`Education level`", levels(df$`Education level`)[-1], " = 0", collapse = " & ")

# 执行Wald检验
linearHypothesis(qr, hypothesis.matrix = hypotheses)

如果需要用自助法（bootstrap）的标准误计算Wald检验，可先通过summary.rq获取boot协方差矩阵再传入：

# 获取带boot标准误的模型协方差矩阵
qr_boot_sum <- summary(qr, se = "boot", bsmethod = "xy")

# 用boot协方差矩阵做Wald检验
linearHypothesis(qr, hypothesis.matrix = hypotheses, vcov = qr_boot_sum$cov)

方法3：自助法（Bootstrap）自定义全局检验 #

若需要更灵活的检验逻辑，可手动实现自助抽样：每次重抽样后拟合模型并计算检验统计量，最后通过统计量的分布推导p值：

library(boot)

# 定义自助函数：返回Wald检验统计量
boot_wald <- function(data, indices) {
  d <- data[indices, ]
  qr_boot <- rq(Total ~ `Education level`, data = d, tau = 0.5)
  # 构造排除截距的系数矩阵
  coef_sub <- coef(qr_boot)[-1]
  vcov_sub <- vcov(qr_boot)[-1, -1]
  # 计算Wald统计量
  wald_stat <- t(coef_sub) %*% solve(vcov_sub) %*% coef_sub
  return(wald_stat)
}

# 执行1000次自助抽样
set.seed(111)
boot_result <- boot(data = df, statistic = boot_wald, R = 1000)

# 计算原模型Wald统计量大于自助样本统计量的比例，即为p值
original_wald <- t(coef(qr)[-1]) %*% solve(vcov(qr)[-1, -1]) %*% coef(qr)[-1]
p_value <- mean(boot_result$t >= original_wald)
cat("Bootstrap全局p值：", p_value, "\n")

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