R语言中mice多重插补后的MANOVA分析及结果合并求助

假设你已经用mice()生成了包含30个插补数据集的对象imp，以下是针对这些数据集运行MANOVA并合并结果的完整步骤：

1. 加载所需包 #

library(mice)

2. 在所有插补数据集上运行MANOVA #

使用mice的with()函数批量处理每个插补数据集，这里假设你的多因变量为y1、y2、y3，自变量为group，可根据实际情况修改公式：

# 定义MANOVA公式：多因变量用cbind()合并
manova_formula <- cbind(y1, y2, y3) ~ group

# 在每个插补数据集上执行MANOVA
manova_fits <- with(imp, manova(manova_formula))

3. 提取每个插补结果的关键统计量 #

我们需要提取每个MANOVA结果中Wilks' Lambda的检验信息（F值、自由度、p值），定义一个提取函数：

extract_manova_results <- function(fit) {
  # 提取Wilks' Lambda检验表
  wilks_table <- summary(fit, test = "Wilks")$stats
  # 提取自变量（这里是group）的统计量
  stats <- wilks_table["group", c("Wilks approx F", "num Df", "den Df", "Pr(>F)")]
  names(stats) <- c("F_value", "num_df", "den_df", "p_value")
  return(stats)
}

# 批量提取所有插补结果的统计量
all_stats <- lapply(manova_fits$analyses, extract_manova_results)
stats_df <- do.call(rbind, all_stats)

4. 用Rubin's规则合并结果 #

多重插补的结果合并需要遵循Rubin's规则，针对MANOVA的Wilks' Lambda统计量，我们通过其对数转换（近似正态分布）来合并：

# 提取每个插补结果的Wilks' Lambda对数
extract_log_wilks <- function(fit) {
  wilks_value <- summary(fit, test = "Wilks")$stats["group", "Wilks"]
  return(log(wilks_value))
}

log_wilks_list <- lapply(manova_fits$analyses, extract_log_wilks)
log_wilks <- unlist(log_wilks_list)

# 应用Rubin's规则合并
m <- length(log_wilks) # 插补数据集数量（30）
Q_bar <- mean(log_wilks) # 对数Wilks的均值
U_bar <- var(log_wilks) / m # 内插补方差
B <- var(log_wilks) # 间插补方差
T <- U_bar + (1 + 1/m) * B # 合并方差

# 转换回原始Wilks' Lambda，并计算合并的F值和p值
merged_wilks <- exp(Q_bar)
s <- 3 # 多因变量的数量，根据你的实际情况修改
df1 <- stats_df$num_df[1] # 自变量自由度
df2 <- stats_df$den_df[1] # 误差自由度

# Wilks' Lambda转F统计量的公式
merged_F <- ((1 - merged_wilks^(1/s)) / merged_wilks^(1/s)) * ((df2 - df1*s + s)/df1)
merged_p <- pf(merged_F, df1*s, df2 - df1*s + s, lower.tail = FALSE)

5. 输出合并后的结果 #

cat("=== 合并后的MANOVA结果 ===\n")
cat("Wilks' Lambda: ", round(merged_wilks, 4), "\n", sep = "")
cat("F统计量: ", round(merged_F, 4), "\n", sep = "")
cat("自由度: ", df1*s, " / ", df2 - df1*s + s, "\n", sep = "")
cat("p值: ", round(merged_p, 4), "\n", sep = "")

注意事项 #

• 如果你的模型包含多个自变量（如协变量），只需修改manova_formula，并在提取统计量时替换"group"为对应变量名即可。
• 运行MANOVA前需确保数据满足多元正态性、方差-协方差矩阵齐性等前提假设，插补后的数据集也需要验证这些假设。
• 若你偏好更简单的合并方式，可使用Fisher方法合并所有插补结果的p值：

fisher_stat <- -2 * sum(log(stats_df$p_value))
merged_p_fisher <- pchisq(fisher_stat, df = 2*m, lower.tail = FALSE)

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