在R中为MICE插补+CBPS加权数据计算聚类标准误遇问题求助

问题概述 #

通过mice完成多重插补，再用MatchThem::weightthem做CBPS加权后，尝试三种方法计算以schoolID为聚类的标准误均报错：

• estimatr包：Error in eval_tidy(mfargs[[da]], data = data) : object 'schoolID' not found（schoolID在无聚类模型中可正常调用）
• miceadds包：Error in as.data.frame.default(data) : cannot coerce class "wimids" to a data.frame
• sandwich+lmtest包：Error in UseMethod("estfun") : no applicable method for 'estfun' applied to an object of class "c('mimira', 'mira')"

核心原因分析 #

1. schoolID找不到：mice默认仅对含缺失值的变量进行插补并保留，若schoolID无缺失，会被排除在插补后的数据集之外。
2. wimids类无法转成数据框：weightthem返回的wimids是加权多重插补专用对象，不能直接作为普通数据框传入miceadds的函数。
3. mira对象不支持estfun：sandwich包默认没有为多重插补模型的mira/mimira类实现对应的方差提取方法。

可行解决步骤与代码示例 #

步骤1：确保插补数据集保留schoolID #

修改mice调用，强制保留schoolID（即使无缺失）：

# 定义插补方法：对schoolID不插补但保留，其他变量用pmm
meth <- make.method(d)
meth["schoolID"] <- ""  # 空字符串表示跳过该变量的插补，直接保留

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

# 重新生成加权数据
weighted_data <- weightthem(trtmnt ~ x1 + x2 + x3,
                            data = tempdata,
                            method = "cbps",
                            estimand = "ATT")

步骤2：提取加权插补数据集并拟合聚类模型 #

先将wimids对象转为单个插补数据集的列表，再对每个数据集拟合带权重和聚类的模型：

library(estimatr)
library(purrr)

# 提取所有m个加权插补数据集
weighted_datasets <- complete(weighted_data, "all")

# 定义拟合函数：替换为你的结局变量和协变量
fit_cluster_model <- function(data) {
  lm_robust(
    outcome ~ trtmnt + x1 + x2 + x3,  # 替换为实际结局变量与协变量
    data = data,
    weights = weights,  # weightthem生成的权重变量名为weights
    clusters = schoolID,
    se_type = "stata"  # 采用Stata风格的聚类标准误
  )
}

# 批量拟合模型
models <- map(weighted_datasets, fit_cluster_model)

步骤3：用Rubin规则手动合并聚类稳健结果 #

由于mice的pool()函数默认不支持聚类稳健方差的合并，需手动计算：

# 提取每个模型的系数、方差、自由度
coefs <- map_dfr(models, ~as.data.frame(t(coef(.x))))
vars <- map(models, ~vcov(.x))
df_vec <- map_dbl(models, ~.x$df.residual)

m <- length(models)  # 插补次数

# 计算合并系数
coef_pooled <- colMeans(coefs)

# 插补内方差：各模型方差的均值
var_within <- colMeans(do.call(rbind, vars))

# 插补间方差：系数的方差，乘以Rubin调整因子(m+1)/m
var_between <- apply(coefs, 2, var) * (m + 1)/m

# 合并方差（含插补内、插补间及调整项）
var_pooled <- var_within + var_between + var_between/m

# 计算标准误、t值与p值
se_pooled <- sqrt(var_pooled)
t_stats <- coef_pooled / se_pooled
df_pooled <- (m - 1) * (1 + var_within/(var_between + var_between/m))^2
p_vals <- 2 * pt(abs(t_stats), df = df_pooled, lower.tail = FALSE)

# 输出结果
results <- data.frame(
  Variable = names(coef_pooled),
  Coefficient = coef_pooled,
  Clustered_SE = se_pooled,
  t_value = t_stats,
  df = df_pooled,
  p_value = p_vals
)
print(results, row.names = FALSE)

替代方案：使用mitml包简化流程 #

mitml专门针对多重插补数据的模型拟合与结果合并，支持加权和聚类：

library(mitml)
library(lme4)

# 将wimids对象转为mitml长格式数据
mitml_data <- mitmlComplete(weighted_data, "long")

# 拟合带权重和聚类的混合效应模型
fit <- with(mitml_data, lmer(
  outcome ~ trtmnt + x1 + x2 + x3 + (1|schoolID),
  weights = weights
))

# 合并结果并计算聚类稳健标准误
pooled_results <- testEstimates(fit, var.comp = "cluster", cluster = "schoolID")
summary(pooled_results)

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