R语言分位数回归偏差编码下分类变量缺失水平的Beta系数与95%CI获取

背景 #

用quantreg包做中位数回归（tau=0.5）时，对含A、B、C三组的分类变量group采用偏差编码（contr.sum），模型仅输出A、B组的系数，C组系数需通过β_C = -(β_A + β_B)手动计算，下面给出两种获取C组95%置信区间的实用方法。

示例基础代码 #

先构造数据并拟合模型：

set.seed(123)
dat <- data.frame(
  y = rnorm(100),
  group = factor(sample(c("A", "B", "C"), 100, replace = TRUE))
)

# 设置偏差编码
contrasts(dat$group) <- contr.sum(3)
colnames(contrasts(dat$group)) <- c("A", "B")

# 拟合中位数回归
library(quantreg)
qr_model <- rq(y ~ group, data = dat, tau = 0.5)
summary(qr_model)

模型输出示例：

Call: rq(formula = y ~ group, tau = 0.5, data = dat)

tau: [1] 0.5

Coefficients:
            coefficients lower bd upper bd
(Intercept)  0.04791     -0.1772  0.2558  
groupA       0.20614     -0.1234  0.5357  
groupB      -0.18058     -0.5099  0.1487  

方法一：基于渐近正态性的计算 #

利用模型的方差-协方差矩阵，通过线性组合的方差公式快速计算置信区间：

# 提取系数和方差-协方差矩阵
qr_coef <- coef(qr_model)
qr_vcov <- vcov(qr_model)

# 定义线性组合权重：β_C = -β_A -β_B
weights <- c(0, -1, -1)

# 计算C组系数
beta_C <- t(weights) %*% qr_coef

# 计算标准误
se_C <- sqrt(t(weights) %*% qr_vcov %*% weights)

# 95%置信区间
ci_C_lower <- beta_C - 1.96 * se_C
ci_C_upper <- beta_C + 1.96 * se_C

# 输出结果
cat("C组系数：", round(beta_C, 4), "\n")
cat("95%置信区间（渐近正态）：", round(ci_C_lower, 4), "至", round(ci_C_upper, 4), "\n")

方法二：自助法（更稳健） #

通过重复抽样估计C组系数的分布，适合小样本或对结果稳健性要求高的场景：

# 定义自助抽样的统计量函数
boot_qr_C <- function(data, indices) {
  d <- data[indices, ]
  contrasts(d$group) <- contr.sum(3)
  colnames(contrasts(d$group)) <- c("A", "B")
  model <- rq(y ~ group, data = d, tau = 0.5)
  coefs <- coef(model)
  return(-(coefs["groupA"] + coefs["groupB"]))
}

# 运行自助抽样（1000次）
library(boot)
set.seed(456)
boot_result <- boot(data = dat, statistic = boot_qr_C, R = 1000)

# 提取百分位数法的95%置信区间
ci_C_boot <- boot.ci(boot_result, type = "perc")$percent[4:5]

# 输出结果
cat("C组系数（自助均值）：", round(mean(boot_result$t), 4), "\n")
cat("95%置信区间（自助百分位数）：", round(ci_C_boot[1], 4), "至", round(ci_C_boot[2], 4), "\n")

关键说明 #

• 偏差编码的核心约束是三组系数之和为0，因此C组系数必然是A、B组系数的相反数之和。
• 渐近正态法计算高效，但依赖大样本假设；自助法无需分布假设，结果更稳健。

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