如何用R中pool.scalar计算多重插补模型系数的p值？

一、获取p值的方法 #

由于pool()函数不支持分位数回归（rq）模型，使用pool.scalar()实现Rubin法则后，可通过Barnard-Rubin调整自由度的t分布计算p值，具体步骤如下：

1. 从每个插补数据集的rq模型中提取目标系数的估计值和对应方差（分位数回归的方差可通过summary(rq_model, se="nid")或se="boot"等方法获取）
2. 用pool.scalar()合并得到核心统计量：qbar（合并系数）、U（插补内方差均值）、B（插补间方差）、T（总方差，公式为U + (1 + 1/m)*B，m为插补次数）
3. 计算调整后的自由度：

   m <- length(imputed_datasets) # 插补次数
   df <- (m - 1) * (1 + U / ((1 + 1/m) * B))^2
   

4. 计算t统计量并得到双侧p值：

   t_val <- qbar / sqrt(T)
   p_val <- 2 * pt(abs(t_val), df = df, lower.tail = FALSE)
   

示例代码片段：

library(quantreg)
library(mice)

# 假设imp是已完成多重插补的mice对象
m <- imp$m
coef_list <- list()
var_list <- list()

# 遍历插补数据集，拟合分位数回归并提取系数与方差
for (i in 1:m) {
  dat <- complete(imp, i)
  model <- rq(y ~ x1 + x2, data = dat, tau = 0.5)
  mod_sum <- summary(model, se = "nid") # 用nid方法估计方差
  coef_list[[i]] <- mod_sum$coefficients[, "Value"]
  var_list[[i]] <- diag(mod_sum$cov) # 提取系数的方差
}

# 以第一个自变量系数为例，执行Rubin合并
target_coef_idx <- 2 # x1对应的系数位置
q_list <- sapply(coef_list, function(x) x[target_coef_idx])
u_list <- sapply(var_list, function(x) x[target_coef_idx])
pool_result <- pool.scalar(q = q_list, u = u_list)

# 计算调整自由度与p值
U <- pool_result$u
B <- pool_result$b
T_total <- pool_result$t
df <- (m - 1) * (1 + U / ((1 + 1/m)*B))^2
t_val <- pool_result$qbar / sqrt(T_total)
p_val <- 2 * pt(abs(t_val), df = df, lower.tail = FALSE)

二、置信区间的其他计算方式 #

除了你已使用的qbar ± qt(0.975, df)*sqrt(T)（基于调整t分布），还有两种常用替代方式：

1. 正态近似法 #

当插补次数较多（m≥20）时，可直接使用标准正态分布分位数计算，计算效率更高：

ci_norm <- pool_result$qbar + qnorm(c(0.025, 0.975)) * sqrt(T_total)

2. Bootstrap合并法 #

对每个插补数据集做Bootstrap抽样拟合分位数回归，再合并所有Bootstrap结果，无需依赖Rubin法则的方差假设，更适配分位数回归的非正态特性：

# 定义单轮Bootstrap的分位数回归系数提取函数
boot_rq_coef <- function(dat, tau=0.5, target_idx) {
  boot_dat <- dat[sample(nrow(dat), replace=TRUE), ]
  model <- rq(y ~ x1 + x2, data=boot_dat, tau=tau)
  coef(model)[target_idx]
}

# 每个插补数据集执行1000次Bootstrap
B_boot <- 1000
boot_coefs <- list()
for (i in 1:m) {
  dat <- complete(imp, i)
  boot_coefs[[i]] <- replicate(B_boot, boot_rq_coef(dat, target_idx=target_coef_idx))
}

# 合并所有Bootstrap系数，计算95%置信区间
all_boot_coefs <- unlist(boot_coefs)
ci_boot <- quantile(all_boot_coefs, c(0.025, 0.975))

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