如何在插补数据集的回归模型中提取sigma值？

问题说明 #

常规线性模型（lm）可以直接用sigma()函数提取残差标准差：

library(mice)
nhanes$group <- as.factor(c(rep("A", 12), rep("B", 13)))
m_1 <- lm(bmi ~ group + age, nhanes)
sigma(m_1)

但通过mice包插补数据集后，用with()构建的mira类模型对象无法直接调用sigma()，会返回空值并报错；尝试with(m_2, sigma())或sigma(pool(m_2))也无法得到正确结果。

解决方法 #

1. 获取单个插补模型的sigma值 #

mira对象本质是包含多个插补模型的列表，可通过遍历每个模型提取对应sigma：

# 遍历所有插补模型，提取每个模型的sigma
sapply(m_2$analyses, sigma)

执行后会返回一个向量，每个元素对应一个插补数据集的残差标准差。

2. 计算合并后的加权平均sigma #

若需要得到整体的合并估计值，可基于每个模型的sigma和自由度进行加权平均（采用方差加权后开平方的方式，因为方差具有可加性）：

# 提取每个模型的sigma和自由度
sigma_list <- sapply(m_2$analyses, sigma)
df_list <- sapply(m_2$analyses, function(x) nobs(x) - length(coef(x)))

# 计算加权平均的sigma
pooled_sigma <- sqrt(weighted.mean(sigma_list^2, df_list))
pooled_sigma

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