如何从MIcombine获取组内组间方差及λ比值？

问题背景 #

在对多重插补数据应用raking加权后，需要合并估计值并报告以下参数以描述插补变异性的影响：

• U-hat：估计值的平均组内方差
• B：组间方差（缺失数据导致的额外方差）
• λ：缺失数据导致的总方差占比（B与总方差的比值）

由于使用mitools包的MIcombine替代了mice的pool()，已能提取组内标准误，但需要手动计算B和λ。

解决方案步骤 #

1. 提取各插补数据集的估计值与方差 #

从raking后的调查设计对象中，提取每个插补数据集的BMI均值估计值及其方差：

# 提取每个插补集的raking后估计结果
rake_results <- with(small.svy.rake, svymean(~bmi))

# 提取所有插补集的估计值向量
estimates <- sapply(rake_results, function(x) coef(x))
# 提取所有插补集的组内方差向量
variances <- sapply(rake_results, function(x) vcov(x)[1,1])

2. 计算U-hat（平均组内方差） #

U-hat是各插补数据集方差的平均值，直接对variances取均值即可：

U_hat <- mean(variances)

3. 计算B（组间方差） #

组间方差反映插补带来的估计值变异，计算公式为：
[ B = \frac{m}{m-1} \times \text{var}(\text{estimates}) ]
其中m是插补次数（本例中m=5）：

m <- length(estimates)
B <- var(estimates) * m / (m - 1)

4. 计算总方差T与λ #

总方差是组内方差、组间方差的合并值，公式为：
[ T = \hat{U} + B + \frac{B}{m} ]
λ则是组间方差占总方差的比例：

T <- U_hat + B + B/m
lambda <- B / T

5. 完整验证代码 #

整合以上步骤，直接运行即可得到所有参数：

# 提取结果
rake_results <- with(small.svy.rake, svymean(~bmi))
estimates <- sapply(rake_results, function(x) coef(x))
variances <- sapply(rake_results, function(x) vcov(x)[1,1])

# 计算各参数
m <- length(estimates)
U_hat <- mean(variances)
B <- var(estimates) * m / (m - 1)
T <- U_hat + B + B/m
lambda <- B / T

# 输出结果
cat("U-hat（平均组内方差）：", U_hat, "\n")
cat("B（组间方差）：", B, "\n")
cat("λ（缺失数据总方差占比）：", lambda, "\n")

结果说明 #

• 计算得到的T的平方根，与summary(MIcombine(rake_results))中的标准误一致，可验证计算正确性。
• λ值越大，说明缺失数据导致的额外方差占比越高，插补变异性对估计结果的影响越大。

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