使用marginaleffects聚类Bootstrap标准误遇数据掩码过时错误的解决方法

我针对不同因变量拟合了多个线性模型，并用marginaleffects包计算边际效应，基础代码如下：

library(palmerpenguins)
library(marginaleffects)
library(sandwich)
library(tidyr)
library(dplyr)

long_pengs = penguins |>
  pivot_longer(cols = c(body_mass_g, flipper_length_mm),
               names_to = 'outcome',
              values_to = 'vals') |>
  drop_na(sex) |>
  summarise(mods = list(lm(vals ~ sex * bill_length_mm, data = pick(everything()))), .by = outcome)

comps = long_pengs |>
  rowwise(outcome) |>
  reframe(avg_comparisons(mods,
                          variables = 'sex',
                          subset(sex == 'female')))

但尝试用聚类Bootstrap计算标准误时，出现了以下错误：

# 这段代码可以运行
long_pengs |> 
  rowwise(outcome) |>
  reframe(avg_comparisons(mods,
    variables = 'sex',
    subset(sex == 'female')) |>
      inferences(method = 'rsample')) 

# 这段代码报错
long_pengs |>
rowwise(outcome) |>
reframe(avg_comparisons(mods,
                        variables = 'sex',
                        subset(sex == 'female'),
                        vcov = vcovBS(mods, cluster = ~species)))
#> Error in `reframe()`:
#> ℹ In argument: `avg_comparisons(...)`.
#> ℹ In row 1.
#> Caused by error:
#> ! Obsolete data mask.
#> ✖ Too late to resolve `species` after the end of `dplyr::summarise()`.
#> ℹ Did you save an object that uses `species` lazily in a column in the
#>   `dplyr::summarise()` expression ?

# 这段代码也报错
long_pengs |> 
  rowwise(outcome) |>
  reframe(avg_comparisons(mods,
    variables = 'sex',
    subset(sex == 'female')) |>
      inferences(method = 'rsample', strata = species))
#> Error in `reframe()`:
#> ℹ In argument: `inferences(...)`.
#> ℹ In row 1.
#> Caused by error:
#> ! object 'species' not found

期望得到的输出是类似手动拟合每个模型再合并的结果：

## 期望输出
m1 = lm(body_mass_g ~ sex * bill_length_mm, data = penguins)

c1 = avg_comparisons(m1, variables = 'sex',
                     subset(sex == 'female'),
                    vcov = vcovBS(m1, cluster = ~species))

m2 = lm(flipper_length_mm ~ sex * bill_length_mm, data = penguins)

c2 = avg_comparisons(m2, variables = 'sex',
                     subset(sex == 'female'),
                     vcov = vcovBS(m2, cluster = ~species))

rbind(c1, c2)
#> 
#>  Estimate Std. Error     z Pr(>|z|)   S   2.5 % 97.5 %
#>   420.487     293.10 1.435    0.151 2.7 -153.97  994.9
#>     0.392       5.16 0.076    0.939 0.1   -9.73   10.5
#> 
#> Term: sex
#> Type:  response 
#> Comparison: mean(male) - mean(female)
#> Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted

解决方案 #

方法1：放弃非标准评估，直接循环处理

既然不局限于非标准评估，直接遍历因变量列表，逐个拟合模型并计算边际效应是最直观的方式，避免tidyverse数据掩码的问题：

library(palmerpenguins)
library(marginaleffects)
library(sandwich)
library(dplyr)

# 定义要分析的因变量
outcomes = c("body_mass_g", "flipper_length_mm")

# 循环处理每个因变量
results = lapply(outcomes, function(outcome) {
  # 拟合模型
  mod = lm(paste(outcome, "~ sex * bill_length_mm"), data = penguins |> drop_na(sex))
  # 计算边际效应并指定聚类标准误
  avg_comparisons(mod, 
                  variables = 'sex',
                  subset = sex == 'female',
                  vcov = vcovBS(mod, cluster = ~species)) |>
    # 添加因变量标识列
    mutate(outcome = outcome)
})

# 合并结果
bind_rows(results)

方法2：调整tidyverse流程，保留聚类变量的上下文

如果想继续用tidyverse的分组处理，需要确保species变量在数据掩码中可用，或者在拟合模型时提前指定聚类信息：

long_pengs = penguins |>
  pivot_longer(cols = c(body_mass_g, flipper_length_mm),
               names_to = 'outcome',
               values_to = 'vals') |>
  drop_na(sex) |>
  # 额外保留species列，确保后续能访问到
  summarise(mods = list(lm(vals ~ sex * bill_length_mm, data = pick(everything()))),
            species_data = list(species), .by = outcome)

# 计算边际效应时，从保留的species_data中提取聚类变量
long_pengs |>
  rowwise(outcome) |>
  reframe(
    avg_comparisons(mods,
                    variables = 'sex',
                    subset = sex == 'female',
                    vcov = vcovBS(mods, cluster = species_data[[1]])) |>
      mutate(outcome = outcome)
  )

两种方法都能得到期望的输出结果，解决数据掩码和变量找不到的问题。

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