r - 调查加权和多重估算数据的汇总平均边际效应

Question

除了使用 mouse() 估算的缺失数据外，我还正在处理调查数据及其相关权重。我最终运行的模型包含变量之间复杂的相互作用，我希望获得平均边际效应。

这个任务在 STATA 中似乎微不足道，但我宁愿留在 R 中，因为这是我最了解的。检索每个单独的估算数据集的 AME 并对估计值进行平均似乎很容易。但是，我需要使用 pool() (来自小鼠)来确保我得到正确的标准错误。

这是一个可重现的示例:

library(tidyverse)
library(survey)
library(mice)
library(margins)

df <- tibble(y = c(0, 5, 0, 4, 0, 1, 2, 3, 1, 12), region = c(1, 1, 1, 1, 1, 3, 3, 3, 3, 3), 
             weight = c(7213, 2142, 1331, 4342, 9843, 1231, 1235, 2131, 7548, 2348), 
             x1 = c(1.14, 2.42, -0.34, 0.12, -0.9, -1.2, 0.67, 1.24, 0.25, -0.3),
             x2 = c(12, NA, 10, NA, NA, 12, 11, 8, 9, 9))

在简单(非多个)svyglm 上使用 margins() 可以顺利工作。使用 which() 对每个插补运行 svyglm 并汇集结果也效果很好。

m <- with(surv_obj, svyglm(y ~ x1 * x2))
pool(m)

但是，将 margins() 包装到 which() 中会返回错误“.svycheck(design) 中的错误:缺少参数“design”，没有默认值”

with(surv_obj, margins(svyglm(y ~ x1 * x2), design = surv_obj))

如果我在 svyglm 调用中指定设计，则会收到“UseMethod("svyglm", design) 中的错误:没有适用于“svyglm”的方法应用于类“svyimputationList”的对象”

with(surv_obj, margins(svyglm(y ~ x1 * x2, design = surv_obj), design = surv_obj))

如果我放弃调查图层，并简单地尝试在每个估算集上运行边距，然后进行池化，我会收到警告:“get.dfcom(object, dfcom) 中的警告:假定样本大小无限大。”。

m1 <- with(imputed_df, margins(lm(y ~ x1 * x2)))
pool(m1)

考虑到 pool() 可能在计算中使用样本大小，这让我很担心。

有谁知道有什么方法可以(a)使用which()，margins()和pool()来检索汇集的平均边际效应或(b)知道我应该将margins()的哪些元素传递给pool( ) (或 pool.scalar()) 来达到预期的结果？

Answer 1

根据 Vincent 的评论进行更新

想要根据 Vincent 的评论和相关包 margineffects() 更新这篇文章，最终解决了我的问题。希望这对遇到类似问题的其他人有所帮助。

我实现了文森特评论中链接的小插图中的代码，添加了一些允许调查加权和建模的步骤。值得注意的是 svydesign() 将删除聚类/加权变量上丢失的任何观察结果，因此 margineffects() 无法将值预测回原始“dat”数据，并会抛出错误。汇集我的实际数据仍然会出现“假设无限样本大小”，这(如前所述)应该没问题，但我仍在研究修复。

library(tidyverse)
library(survey)
library(mice)
library(marginaleffects)

fit_reg <- function(dat) {
  
    svy <- svydesign(ids = ~ 1, cluster = ~ region, weight = ~weight, data = dat)
    mod <- svyglm(y ~ x1 + x2*factor(x3), design = svy)
    out <- marginaleffects(mod, newdata = dat)
    
    class(out) <- c("custom", class(out))
    return(out)
}

tidy.custom <- function(x, ...) {
    out <- marginaleffects:::tidy.marginaleffects(x, ...)
    out$term <- paste(out$term, out$contrast)
    return(out)
}

df <- tibble(y = c(0, 5, 0, 4, 0, 1, 2, 3, 1, 12), region = c(1, 1, 1, 1, 1, 3, 3, 3, 3, 3), 
             weight = c(7213, 2142, 1331, 4342, 9843, 1231, 1235, 2131, 7548, 2348), 
             x1 = c(1.14, 2.42, -0.34, 0.12, -0.9, -1.2, 0.67, 1.24, 0.25, -0.3),
             x2 = c(12, NA, 10, NA, NA, 12, 11, 8, 9, 9),
             x3 = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2))

imputed_df <- mice(df, m = 2, seed = 123)

dat_mice <- complete(imputed_df, "all")
mod_imputation <- lapply(dat_mice, fit_reg)
mod_imputation <- pool(mod_imputation)

summary(mod_imputation)