mice多重插补后面板固定/随机效应模型拟合与检验问题咨询

报错和broom.mixed包无关，完全是with()函数传参位置错误导致的。mice提供的with.mids()方法仅会执行表达式位置的建模语句，你将plm()所需的data、index、effect、model参数放在了with()的外层参数位，这些参数不会传递给plm()，最终函数拿到的是未拆解的mids类插补对象，才会触发No tidy method for objects of class mids报错。
另外你当前的插补流程没有适配面板数据结构，存在结果偏差风险，下面给出完整修正方案。

普通PMM插补未考虑面板的个体聚类特征，会导致插补值不符合个体内序列规律，需要把个体标识作为聚类变量纳入预测矩阵，不要将其完全排除在插补流程外：

library(mice)
library(plm)

# 修正预测矩阵：保留id作为聚类分组变量，不插补id、Name、Year
predictormatrix <- quickpred(
  data,
  include = c("a", "b", "c", "x", "y", "id"), # 将id加入预测变量集合
  exclude = c("Name", "Year"), # 不插补Name/Year，若需要控制时间趋势可保留Year
  mincor = 0.1
)
# 强制id变量不被插补
predictormatrix[, c("id")] <- 0

# 重新运行插补，可适当提高m和maxit提升结果稳定性：m建议至少设为核心变量缺失率的100倍（x缺失35%时设m=35更严谨，最低不要小于10）
imp <- mice(
  data,
  predictorMatrix = predictormatrix,
  m = 10,
  maxit = 10,
  meth = "pmm",
  seed = 123 # 设置随机种子保证结果可复现
)

所有plm()的参数必须全部写在with()的表达式内部，不要放在外层。with.mids()会自动依次调用每个插补完成的完整数据集作为建模数据源，只需要在plm()里指定data = .data（代表当前调用的单个插补数据集）即可：

# 拟合个体固定效应模型（within模型）
fitimp.fe <- with(
  imp,
  plm(
    y ~ x + a + b + c,
    data = .data, # .data是mice自动传入的当前插补完整数据集
    index = c("id", "Year"), # 用id做个体标识比Name更稳妥，可避免重名问题
    effect = "individual",
    model = "within"
  )
)
# 池化合并固定效应结果
fe_res <- summary(pool(fitimp.fe))

# 拟合随机效应模型
fitimp.re <- with(
  imp,
  plm(
    y ~ x + a + b + c,
    data = .data,
    index = c("id", "Year"),
    effect = "individual",
    model = "random"
  )
)
# 池化合并随机效应结果
re_res <- summary(pool(fitimp.re))

运行上述代码不会再出现类报错，2020年后更新的mice和plm版本原生支持plm对象的池化合并，不需要额外加载broom.mixed包。

不能直接将池化后的固定/随机效应结果套入原生phtest()函数，因为phtest()需要单个模型的协方差矩阵。正确做法是对每一个插补数据集单独做Hausman检验，再按照Rubin规则合并检验统计量：

# 遍历每个插补数据集，单独完成Hausman检验
hausman_list <- lapply(1:imp$m, function(i) {
  # 提取第i个插补数据集
  current_data <- complete(imp, i)
  # 分别拟合当前数据集的FE、RE模型
  fe_current <- plm(y ~ x + a + b + c, data = current_data, index = c("id","Year"), model = "within")
  re_current <- plm(y ~ x + a + b + c, data = current_data, index = c("id","Year"), model = "random")
  # 返回当前数据集的Hausman检验结果
  phtest(fe_current, re_current)
})

# 提取m个检验的卡方统计量和自由度，按Rubin规则合并
chi_stats <- sapply(hausman_list, function(x) x$statistic)
df <- hausman_list[[1]]$parameter
# 计算合并后的卡方统计量和对应p值
pooled_chi <- mean(chi_stats)
p_value <- pchisq(pooled_chi, df = df, lower.tail = FALSE)

# 输出检验结果
cat("合并后Hausman检验卡方统计量：", round(pooled_chi,3), "，自由度：", df, "，p值：", round(p_value,4), "\n")

若检验p值小于0.05则选择固定效应模型，反之选择随机效应模型。

• 不要在插补前手动对面板数据做组内去均值处理，会人为扭曲变量协方差结构，造成插补偏差，直接使用原始数据完成插补即可。
• 如果需要控制时间固定效应，可将plm()的effect参数改为"twoways"，插补环节也可将Year转化为虚拟变量纳入预测矩阵，提升插补准确性。
• 你之前设置的m=5、maxit=5对于核心自变量缺失率35%的数据集来说，插补次数和迭代次数都不足，结果稳定性较差，建议至少设置m=20、maxit=15，插补完成后确认变量分布、迭代链收敛性达标后再开展后续建模。

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