R语言多重插补数据集下glmnet岭回归模型系数合并方案咨询

pool()函数仅支持自带标准误提取方法的模型对象（如lm、glm类），glmnet类的正则化模型默认不输出系数标准误，因此无法直接适配。

步骤1：提取每个插补模型的系数与标准误 #

可以用parameters包的model_parameters()函数直接提取glmnet模型的系数和对应标准误，代码如下：

# 加载需要的依赖包
library(parameters)
library(purrr)
library(dplyr)
library(jtools)

# 遍历所有插补模型，提取系数和标准误
coef_list <- map(mods3, function(mod) {
  mod_params <- model_parameters(mod, ci_method = "normal")
  as.data.frame(mod_params) %>% 
    select(term = Parameter, estimate = Coefficient, std.error = SE)
})

步骤2：按Rubin法则合并多重插补结果 #

多重插补的结果合并遵循Rubin规则，直接手动计算即可：

m <- 10 # 你的插补次数

# 合并所有插补的系数结果
all_coef <- bind_rows(coef_list, .id = "imp_id")

# 按变量分组计算合并后的系数、标准误
pooled_result <- all_coef %>%
  group_by(term) %>%
  summarise(
    # 合并系数为所有插补估计值的均值
    estimate = mean(estimate),
    # 插补内方差均值
    ubar = mean(std.error^2),
    # 插补间方差
    b = var(estimate),
    # 合并后的总方差
    total_var = ubar + (1 + 1/m)*b,
    # 合并后的标准误
    std.error = sqrt(total_var)
  ) %>%
  # 可选补充统计量，适配绘图需求
  mutate(
    z_score = estimate / std.error,
    p_value = 2 * pnorm(abs(z_score), lower.tail = F),
    conf.low = estimate - qnorm(0.975) * std.error,
    conf.high = estimate + qnorm(0.975) * std.error
  )

步骤3：调用plot_summs绘图 #

直接将合并后的结果传入plot_summs，指定对应列即可：

plot_summs(
  pooled_result,
  coefs = "estimate",
  se = "std.error",
  coef.names = pooled_result$term,
  ci_level = 0.95,
  # 可自定义其他绘图参数，比如点颜色、线条样式等
  point.color = "#2c3e50"
)

• 若对parameters包提取的glmnet标准误可靠性有要求，可在每个插补数据集拟合模型时，通过bootstrap自助法计算系数标准误，替换上述代码中的std.error列即可。
• 小样本场景下，可补充计算Rubin法则校正的自由度，替换上述的正态近似计算p值和置信区间，结果会更严谨。

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