我试图用插补(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`)"
对上述任何一种方法有什么想法,或者下一步该怎么做?
1条答案
按热度按时间z4iuyo4d1#
您已经非常接近了,您需要使用
with(weighted_data, .)
来拟合加权数据集中的模型,并且需要使用estimatr::lm_robust()
来获得聚类标准误差,因此请尝试以下操作:您的第一种和第二种方法是不正确的,因为您将
weighted_data
提供给单个模型,就好像它是一个 Dataframe 一样,但它不是;这是一个复杂的wimids
对象,您需要使用with()
基础设施来使模型适合估算的加权数据。您的第三种方法很接近,但是
coeftest()
需要用于单个模型,而不是mimira
对象,后者包含适合插补数据集的所有模型。对于MatchThem
中的mimira
对象,您不能这样做。这就是estimatr::lm_robust()
的用武之地,因为它能够在每个估算的数据集中应用聚类。我还建议你看看这个blog post,在用多重插补数据加权后估计治疗效果。你的情况与文章中给出的代码的唯一区别是,你会在使用任何函数时将
vcov = "HC3"
改为vcov = ~schoolID
。