使用R合并混合效应模型的方差-协方差矩阵

针对多个插补数据集的方差-协方差矩阵合并，不能仅对10个矩阵逐元素取均值，需要严格遵循Rubin法则，结合组内变异和组间变异两部分计算：

• 步骤1：计算组内方差-协方差矩阵（$\bar{V}$）
  对10个数据集各自输出的方差-协方差矩阵，逐元素取算术平均值，这部分反映单个插补数据集中的抽样变异。

• 步骤2：计算组间方差-协方差矩阵（$B$）
  用每个数据集的系数估计值$\hat{\beta}_m$，减去你已经合并得到的总系数估计值$\bar{\beta}$，得到偏差向量；对每个偏差向量计算外积（向量转置后相乘得到矩阵），再将这10个外积矩阵取平均（通常用$m-1$做分母得到无偏估计，$m=10$是数据集数量），这部分反映不同插补数据集之间的系数变异。

• 步骤3：合并得到总方差-协方差矩阵（$T$）
  用组内矩阵加上校正后的组间矩阵，公式为：
  $$T = \bar{V} + \left(1 + \frac{1}{m}\right)B$$
  其中$\left(1 + \frac{1}{m}\right)$是校正因子，用于调整组间方差的权重，弥补插补数量有限带来的偏差。

针对你问题的具体解答 #

1. 仅对矩阵取均值只完成了组内部分，缺少组间变异的校正，结果不准确；
2. 不需要提取随机效应估计值，直接用每个数据集的系数估计值和合并后的总系数即可计算组间部分；
3. 协方差的处理逻辑和方差完全一致：组内协方差是各矩阵对应协方差元素的均值，组间协方差是偏差向量外积的对应元素均值，最后同样代入上述合并公式得到总协方差。

R代码示例 #

假设你有一个包含10个拟合模型的列表fit_list，每个模型支持coef()提取系数、vcov()提取方差-协方差矩阵：

m <- length(fit_list)

# 已完成的合并系数估计
beta_combined <- colMeans(sapply(fit_list, coef))

# 计算组内方差-协方差矩阵
vcov_within <- Reduce("+", lapply(fit_list, vcov)) / m

# 计算组间方差-协方差矩阵
beta_mat <- sapply(fit_list, coef)
vcov_between <- (beta_mat - beta_combined) %*% t(beta_mat - beta_combined) / (m - 1)

# 合并总方差-协方差矩阵
vcov_total <- vcov_within + (1 + 1/m) * vcov_between

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