R 中估算和加权数据的聚类标准误差

Question

我正在尝试使用估算（MICE）和加权（CBPS）的数据来获得聚类SE（在我的数据中的学校级别）。我尝试了几种不同的方法，但都引发了不同的错误。

这是我必须开始的，效果很好：

library(tidyverse)
library(mice)
library(MatchThem)
library(CBPS)

tempdata <- mice(d, m = 10, maxit = 50, meth = "pmm", seed = 99)

weighted_data <- weightthem(trtmnt ~ x1 + x2 + x3,
                            data = tempdata,
                            method = "cbps",
                            estimand = "ATT")

使用此（https://www.r-bloggers.com/2021/05/clustered-standard-errors-with-r/）作为指导，我尝试了所有 3 个，这都导致了各种类型的错误消息。

我的数据位于受限制的服务器中，因此不幸的是我无法将其带入此处以准确地重现内容，尽管如果它有用，我可以尝试重新创建一些示例数据。

因此，首先尝试使用

estimatr

，我收到此错误：

m1 <- estimatr::lm_robust(outcome ~ trtmnt + x1 + x2 + x3,
                          clusters = schoolID,
                          data = weighted_data)

Error in eval_tidy(mfargs[[da]], data = data) :
  object 'schoolID' not found

我不知道 schoolID 变量会在哪里丢失/不被识别。它不是加权过程的一部分，但它仍然应该在数据框中......如果我将它用作标准模型中的协变量而不进行聚类，它就在那里。

我也尝试使用

miceadds

并收到此错误：

m2 <- miceadds::lm.cluster(outcome ~ trtmnt + x1 + x2 + x3,
                           cluster = "schoolID",
                           data = weighted_data)

Error in as.data.frame.default(data) : 
  cannot coerce class `"wimids"` to a data.frame

最后，用

sandwich

和

lmtest

：

library(sandwich)
library(lmtest)

m3 <- weighted_models <- with(weighted_data,
                              exp=lm(outcome ~ trtmnt + x1 + x2 + x3))

msandwich <- coeftest(m3, vcov = vcovCL, cluster = ~schoolID)

Error in UseMethod("estfun") :
   no applicable method for `estfun` applied to an object of class "c(`mimira`, `mira`)"

对上述任何方法有什么想法，或者下一步该去哪里？

Answer 1

更新（2025/01/13）：

您现在可以在

lm_weightit()

内使用

with()

来拟合线性结果模型，该模型考虑了权重的估计并允许聚类标准误差。执行此操作的代码如下：

weighted_models <- with(weighted_data,
                        lm_weightit(outcome ~ trtmnt + x1 + x2 + x3,
                                    cluster = ~schoolID))

---

你们真的很亲密。您需要使用

with(weighted_data, .)

来拟合加权数据集中的模型，并且需要使用

estimatr::lm_robust()

来获取聚类标准误差。所以请尝试以下操作：

weighted_models <- with(weighted_data,
                        estimatr::lm_robust(outcome ~ trtmnt + x1 + x2 + x3,
                                            cluster = schoolID))

您的第一种和第二种方法是不正确的，因为您将

weighted_data

提供给单个模型，就好像它是一个数据框一样，但事实并非如此；这是一个复杂的

wimids

对象。您需要使用

with()

基础设施将模型拟合到估算的加权数据。

您的第三种方法很接近，但是

coeftest()

需要在单个模型上使用，而不是

mimira

对象，其中包含适合估算数据集的所有模型。虽然您可以在

coeftest()

内部将

with()

与

mira

对象一起使用，但不能对

mimira

中的

MatchThem

对象执行此操作。这就是

estimatr::lm_robust()

的用武之地，因为它能够在每个估算数据集中应用聚类。

我还建议您查看这篇博客文章，了解使用多重插补数据加权后估计治疗效果。您的情况与帖子中提供的代码的唯一区别是，无论您使用哪个函数，您都可以将

vcov = "HC3"

更改为

vcov = ~schoolID

。