plot_slopes迭代函数中无法计算预测值的问题排查与优化

循环内外过滤的差异原因 #

核心问题出在过滤的粒度和针对性不同：

• 循环外全局过滤：通常是一次性移除所有因变量（或任意变量）含NA的行，这会导致每个因变量的有效样本被过度压缩——比如因变量y1的NA在row1，y2的NA在row2，全局过滤会同时去掉row1和row2，而实际上处理y1时row2是有效样本，处理y2时row1是有效样本。当结合subset_var_list拆分分组后，某些分组下的有效样本量可能被不合理地减少，甚至出现模型拟合后，marginaleffects无法在该子集内生成预测值的情况（比如样本量不足、变量水平缺失）。
• 循环内针对性过滤：针对当前迭代的因变量，仅移除该因变量含NA的行，保留其他变量的有效样本。这样每个因变量在对应的分组子集里，都是基于自身的完整有效样本拟合模型，和你手动单独运行每个回归的数据集完全一致，因此marginaleffects可以正常计算预测值。

另外，marginaleffects在计算预测值时会依赖模型拟合时的原始数据集结构，如果循环外过滤后的数据集，与手动拟合单个模型时的数据集（仅过滤当前因变量NA）结构不一致（比如变量的观测数、因子水平的出现频率不同），也会触发无法计算预测值的错误。

循环内过滤前提下的提速方案 #

1. 提前预拆分分组数据集 #

先根据subset_var_list把原始数据拆分为分组数据集列表，避免每次迭代都重复拆分：

library(dplyr)
# 假设subset_var_list是分组变量名的向量
grouped_data <- df %>%
  group_by(across(all_of(subset_var_list))) %>%
  group_split()

后续迭代直接基于这个列表操作，减少重复的分组拆分开销。

2. 并行化迭代 #

用furrr（purrr的并行版）代替mapply，利用多核CPU加速：

library(furrr)
plan(multisession, workers = 4) # 根据CPU核心数设置进程数

# 定义单步处理函数
process_one <- function(data, outcome_var) {
  data_filtered <- data %>% filter(!is.na(.data[[outcome_var]]))
  model <- lm(paste(outcome_var, "~ x * z", sep = " "), data = data_filtered)
  plot_slopes(model, variables = "x", condition = "z")
}

# 并行迭代所有分组和因变量组合
plots <- future_map2(grouped_data, outcome_vars, process_one)

3. 预生成因变量过滤索引 #

提前为每个因变量生成无NA的行索引，循环内直接用索引过滤，避免重复调用is.na和filter：

# 预生成每个因变量的有效行索引
outcome_na_indices <- lapply(outcome_vars, function(var) {
  which(!is.na(df[[var]]))
})
names(outcome_na_indices) <- outcome_vars

# 循环内直接用索引过滤
process_one <- function(data, outcome_var) {
  idx <- outcome_na_indices[[outcome_var]]
  data_filtered <- data[idx, ]
  # 后续模型拟合和绘图操作
}

4. 优化plot_slopes计算 #

• 指定fixed参数固定非交互变量的取值（比如取均值或中位数），减少预测值的计算量：

  plot_slopes(model, variables = "x", condition = "z", fixed = list(other_var = median))
  

• 如果不需要置信区间，设置conf_level = NULL关闭计算：

  plot_slopes(model, variables = "x", condition = "z", conf_level = NULL)
  

• 先计算斜率结果再批量绘图：用slopes()先获取数值结果，再统一用ggplot2绘图，减少plot_slopes的绘图开销。

5. 缓存重复模型 #

如果多个因变量共享同一个分组的模型结构（比如仅因变量不同，自变量相同），可以缓存拟合好的模型框架，仅更换因变量重新拟合，减少重复的模型初始化开销：

# 预定义模型公式模板
formula_template <- function(outcome_var) {
  as.formula(paste(outcome_var, "~ x * z", sep = " "))
}

# 缓存分组数据的模型框架
grouped_models <- lapply(grouped_data, function(data) {
  lm(~ x * z, data = data) # 先拟合空框架，后续替换因变量
})

# 循环内替换因变量重新拟合
process_one <- function(model_frame, outcome_var, data) {
  data_filtered <- data %>% filter(!is.na(.data[[outcome_var]]))
  model_frame$y <- data_filtered[[outcome_var]]
  model <- update(model_frame, formula_template(outcome_var))
  plot_slopes(model, variables = "x", condition = "z")
}

内容的提问来源于stack exchange，提问作者flâneur