R中带Bootstrap的逻辑回归报错：替换项长度不匹配

问题背景 #

需要基于含19个变量（1个二分类因变量）的数据集，在R中构建多元逻辑回归模型，通过Bootstrap结合向后逐步消元筛选变量，保留至少在70%Bootstrap模型中出现的变量，但运行代码时出现错误：

Error in t.star[r, ] <- res[[r]] : 
  number of items to replace is not a multiple of replacement length

错误原因 #

1. 返回结果长度不一致：boot包要求每次bootstrap迭代返回的结果（系数）长度必须相同，但向后逐步消元在不同bootstrap样本中选中的变量数量不同，导致返回的系数向量长度不一致，无法填充到统一的结果矩阵中。
2. 全局变量计数不可靠：直接在bootstrap函数中修改全局变量variable_counts，会因boot函数的环境隔离问题导致计数不准确，甚至在并行运行时出现冲突。
3. 多重插补数据处理疏漏：使用complete(Rubin_imp, "long")得到的长格式数据包含原始数据集（imp=0）和40个插补数据集，未区分插补组就直接bootstrap，可能引入数据混杂。

修正后的代码实现 #

# Step 1: 多重插补（MICE）
library(mice)
library(boot)
library(MASS)  # 替代car包的stepAIC，避免依赖冲突

# 生成40个插补数据集
set.seed(12345)
Rubin_imp <- mice(Final_model_above_5, m = 40, printFlag = FALSE)

# 定义所有候选变量（除因变量外）
all_vars <- setdiff(names(Final_model_above_5), "eaten_within_3_months")
full_formula <- as.formula(paste("eaten_within_3_months ~", paste(all_vars, collapse = "+")))

# 存储所有bootstrap迭代中选中的变量
selected_vars_list <- list()

# 定义Bootstrap函数
bootstrap_regression <- function(data, indices) {
  # 提取bootstrap样本
  boot_sample <- data[indices, ]
  
  # 拟合全模型后做向后逐步消元
  full_model <- glm(full_formula, data = boot_sample, family = "binomial")
  step_model <- stepAIC(full_model, direction = "backward", scope = list(lower = ~1, upper = full_formula), trace = FALSE)
  
  # 获取当前模型选中的变量（排除截距项）
  selected_vars <- setdiff(names(coef(step_model)), "(Intercept)")
  
  # 返回所有变量的系数，未选中的填NA
  coef_vec <- rep(NA, length(all_vars) + 1)  # +1是截距项
  names(coef_vec) <- c("(Intercept)", all_vars)
  coef_vec[names(coef(step_model))] <- coef(step_model)
  
  # 记录选中的变量
  selected_vars_list[[length(selected_vars_list) + 1]] <<- selected_vars
  
  return(coef_vec)
}

# 对每个插补数据集分别进行Bootstrap（符合多重插补分析规范）
# 初始化变量计数表
total_var_counts <- table(factor(all_vars, levels = all_vars))

# 循环处理每个插补数据集
for (i in 1:Rubin_imp$m) {
  imputed_data <- complete(Rubin_imp, i)
  # 清空当前插补集的选中变量列表
  selected_vars_list <- list()
  # 执行Bootstrap
  boot_result <- boot(data = imputed_data, statistic = bootstrap_regression, R = 100)
  # 统计当前插补集的变量频率，累加到总计数
  current_counts <- table(factor(unlist(selected_vars_list), levels = all_vars))
  total_var_counts <- total_var_counts + current_counts
}

# 计算总体出现频率
total_iterations <- Rubin_imp$m * 100
var_frequency <- total_var_counts / total_iterations

# 筛选出现频率≥70%的变量
threshold <- 0.7
final_selected <- names(var_frequency[var_frequency >= threshold])

# 输出结果
cat("变量出现频率:\n")
print(var_frequency)
cat("\n最终筛选出的变量:\n")
print(final_selected)

关键优化点 #

• 统一返回结果长度：预先定义包含所有候选变量（含截距）的系数向量，未被选中的变量系数填NA，确保每次迭代返回的结果长度一致。
• 可靠的变量计数：通过列表记录每次迭代选中的变量，最后用table统计频率，避免全局变量的环境问题；同时对每个插补数据集单独处理，符合多重插补的分析逻辑。
• 提升稳定性：增大Bootstrap迭代次数（R）到100以上，让频率统计更可靠；关闭stepAIC的trace输出，减少冗余信息。
• 规范多重插补处理：循环处理每个插补数据集，避免不同插补组数据混杂，保证结果的严谨性。

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