lvzhaopeng的数字花园

R语言那些事 统计系列（十二） rstatix One-Way MANOVA 多个响应变量方差分析+详细假设检验

17分钟阅读

缘由

统计是生信的基石之一

碎碎念

要爱一个具体的人，不要爱一个抽象的人； 要爱生活，不要爱生活的意义。 ——《卡拉马佐夫兄弟》

Ref

One-Way MANOVA in R: The Ultimate Practical Guide - Datanovia

MANOVA（Multivariate Analysis of Variance）是对多个响应变量同时检验组间是否存在差异的扩展版 ANOVA。与分别做多个单变量 ANOVA（并做多重比较校正）相比，MANOVA 考虑了响应变量间的相关性，通过联合信息提高检验效率并控制 Type I error 的累积。

举例：若 Sepal.Length 与 Petal.Length 都受 Species 影响，但两变量相关，MANOVA 会同时考虑两者的协方差结构判断 Species 是否对“这两个变量的联合向量”有显著影响。

多变量情况 (MANOVA)：检验一个自变量的不同组别在多个响应变量同时考虑（即构成一个均值向量） 时，是否有显著差异。

零假设 (H₀)：所有组的均值向量都相等。 ([µ₁₁, µ₁₂] = [µ₂₁, µ₂₂] = [µ₃₁, µ₃₂])

例子：想比较3种教学方法（自变量）对学生“数学成绩”和“逻辑思维能力得分”（两个响应变量）的联合影响。MANOVA检验的不是数学成绩或逻辑思维单独来看是否有差异，而是检验 “[数学成绩， 逻辑思维]”这个整体组合在不同教学方法间是否存在差异。

代码

library(tidyverse)
library(rstatix)
library(ggpubr)
iris2 <- iris %>%
  select(Sepal.Length, Petal.Length, Species) %>%
  add_column(id = 1:nrow(iris), .before = 1)
 
head(iris2)

假设检验

ggboxplot(
  iris2, x = "Species", y = c("Sepal.Length", "Petal.Length"), 
  merge = T, palette = "jco"
)
iris2 %>%
  group_by(Species) %>%
  get_summary_stats(Sepal.Length, Petal.Length, type = "mean_sd")

说明：按 Species 分组，计算 Sepal.Length 与 Petal.Length 的均值和标准差（mean_sd）。用于查看每组的中心和离散程度，帮助判断是否需要方差齐性检验或方差不齐时选择替代方法。

#样本数据量
iris2 %>%
  group_by(Species) %>%
  summarise(N = n())

说明：显示每组样本量。MANOVA 对样本量平衡（或接近平衡）更稳健，但不要求严格平衡；若组大小非常不平衡，某些检验（如 Box’s M）会更敏感。

#单变量异常值
iris2 %>%
  group_by(Species) %>%
  identify_outliers(Sepal.Length)
 
iris2 %>%
  group_by(Species) %>%
  identify_outliers(Petal.Length)

说明：逐变量（univariate）识别异常值（基于 IQR 或基于 rstatix 的定义）。单变量异常值是后续判断多变量异常值的基础 — 有时一个变量的极端值会影响协方差结构，从而影响 MANOVA。

#多变量异常值
iris2 %>%
  group_by(Species) %>%
  mahalanobis_distance(-id) %>%
  filter(is.outlier == TRUE) %>%
  as.data.frame()

说明：使用马氏距离（Mahalanobis distance）按 Species 识别多变量异常值。

意义：多变量异常值会对协方差矩阵估计和 MANOVA 的判别产生较大影响，发现后应考虑是否剔除或单独分析。

#单变量正态性检测
iris2 %>%
  group_by(Species) %>%
  shapiro_test(Sepal.Length, Petal.Length) %>%
  arrange(variable)

说明：对每个 Species 的每个变量做 Shapiro-Wilk 单变量正态检验。输出会显示每个变量、每个组的 W 值和 p 值。
注意：Shapiro-Wilk 对样本量非常敏感（当 N 很大时，即使微小偏离也会显著）；因此常结合 Q-Q 图进行视觉判断。

#对于每个组，所有点都大致沿参考线落下。 所以我们可以假设数据是正态性的。
#如果您的样本量大于 50，则首选正态 QQ 图，因为在较大的样本量下，Shapiro-Wilk 检验即使对与正态性的微小偏差也变得非常敏感。
ggqqplot(iris2, "Sepal.Length", facet.by = "Species",
         ylab = "Sepal Length", ggtheme = theme_bw())

说明：按组绘制 Q-Q 图，用于视觉判断正态性。Q-Q 图比单纯的显著性检验更能反映数据是否“足够近似正态”以满足分析假设。

#多变量正态性检测
iris2 %>%
  select(Sepal.Length, Petal.Length) %>%
  mshapiro_test()

说明：执行多变量（multivariate）Shapiro-Wilk 检验（mshapiro_test()）。注意事项：rstatix::mshapiro_test() 实际上是基于 mvShapiroTest 的实现（或 stats::shapiro.test 扩展） —— 有常见的限制：通常要求样本量不大（例如 R 中某些实现需要 n <= 50）；如果你的总样本量 > 50，则该检验可能无法使用或结果不可靠。对于较大样本，常用方法是：

  • 分组分别检查单变量正态性 + Q-Q 图（如你已做）
  • 检查多变量的散点矩阵/椭圆，或使用偏态/峰度指标
  • 或使用稳健方法（如 Pillai’s trace 对偏离正态较稳健）
#识别多重共线性
iris2 %>% cor_test(Sepal.Length, Petal.Length)

说明：计算两变量之间的相关系数（Pearson）及 p 值。若高度相关（接近 ±1），说明存在多重共线性问题（在多变量比较中会影响判别与解释）。

#线性假设
library(GGally)
results <- iris2 %>%
  select(Sepal.Length, Petal.Length, Species) %>%
  group_by(Species) %>%
  doo(~ggpairs(.) + theme_bw(), result = "plots")
results
results$plots

说明：ggpairs() 绘制变量间成对散点图矩阵（含密度/直方图与相关系数），按 Species 分组后可以检查 线性关系 假设（MANOVA 假设变量之间在每组内部近似线性关系）。doo() 是 rstatix 的便捷函数用于对每个分组执行操作并把结果保存在 results$plots。

#协方差同质性 ，如果您采用平衡设计（即大小相似的组），则无需过多担心违反方差-协方差矩阵的同质性，您可以继续分析。
box_m(iris2[, c("Sepal.Length", "Petal.Length")], iris2$Species)

说明：Box’s M 检验用于检验各组的协方差矩阵是否相等（homogeneity of covariance matrices）。MANOVA 的一个关键假设就是各组的协方差矩阵相近；若 Box’s M 显著（p < .05）说明拒绝同质性假设，需谨慎。
注意：Box’s M 对偏离多元正态非常敏感并且对样本量不平衡敏感；当样本量较大或不平衡时，即使轻微差异也可导致显著。实务上：

  • 若 Box’s M 不显著：可以继续使用经典 MANOVA（例如 Pillai’s trace）。
  • 若 Box’s M 显著但样本量接近且组大小相似：Pillai’s trace 通常更稳健，可继续；若组大小严重不平衡/偏离严重，考虑稳健替代或对协方差进行修正。
#方差的同质性
iris2 %>% 
  gather(key = "variable", value = "value", Sepal.Length, Petal.Length) %>%
  group_by(variable) %>%
  levene_test(value ~ Species)

说明：Levene 检验用于逐变量检验组间方差是否齐（homogeneity of variance）。若 Levene 显著，说明该变量在不同组间方差不齐，进而可能影响后续单变量 ANOVA 的类型 I/II/III 解读，且在 MANOVA 中也反映为协方差不同（Box’s M 也会发现）。Levene 更适合替代 Bartlett（后者对正态性敏感）。

model <- lm(cbind(Sepal.Length, Petal.Length) ~ Species, iris2)
Manova(model, test.statistic = "Pillai")

说明：

  • 使用 lm(cbind(... ) ~ Species, data) 来建立多变量线性模型（每个响应变量的联合模型）。这是在 R 中进行 MANOVA 的常见做法（stats::lm() 支持 cbind() 多响应）。
  • Manova()（注意大小写，通常来自 car 包；若你加载了 car，Manova() 可以计算多种检验统计量）对 lm 对象执行 MANOVA，test.statistic = "Pillai" 指定使用 Pillai’s trace。
    Pillai’s trace 通常被认为在违反某些假设（如协方差略有不齐或偏离正态）时更为稳健，特别是在组大小不平衡时比 Wilks 更稳健。
#事后测试
grouped.data <- iris2 %>%
  gather(key = "variable", value = "value", Sepal.Length, Petal.Length) %>%
  group_by(variable)

说明：将数据转换成长格式（long），为后续对每个变量分别进行单变量检验或事后比较做准备。

#单变量单因素方差分析 是否不同
 
# Do welch one way anova test  #当违反方差同质性假设
grouped.data %>% welch_anova_test(value ~ Species)
 
# or do Kruskal-Wallis test #单因素方差分析检验的非参数替代方案
grouped.data %>% kruskal_test(value ~ Species)
 
# or use aov() #当满足方差的正态性和同质性假设时可以使用
grouped.data %>% anova_test(value ~ Species)

说明：对每个单变量分别进行单因素检验，选择依据：

  • anova_test()（常规 ANOVA）：当满足正态性与方差齐时使用。
  • welch_anova_test()：当方差不齐时使用（Welch 更稳健）。
  • kruskal_test()：非参数替代（不需要正态假设），但检验的是中位数差异而非均值差异。
    这三者是互斥选项，按你的检验结果（Shapiro/Levene/BoxM）选择适当方法。
#Compute multiple pairwise comparisons 回答哪些组不同
 
#当方差不齐的时候
pwc <- iris2 %>%
  gather(key = "variables", value = "value", Sepal.Length, Petal.Length) %>%
  group_by(variables) %>%
  games_howell_test(value ~ Species) %>%
  select(-estimate, -conf.low, -conf.high) # Remove details
 
pwc

说明：games_howell_test() 是用于两两比较的事后检验，当存在方差不齐且/或样本量不等时 推荐使用（比 Tukey 更稳健）。你在这里对两个变量分别做了 Games-Howell 的成对比较，并去掉了一些列以简化显示。

#多个成对比较
 
pwc <- pwc %>% add_xy_position(x = "Species")
 
test.label <- create_test_label(
  description = "MANOVA", statistic.text = quote(italic("F")),
  statistic = 71.83, p= "<0.0001", parameter = "4,294",
  type = "expression", detailed = TRUE
)
 
ggboxplot(
  iris2, x = "Species", y = c("Sepal.Length", "Petal.Length"), 
  merge = TRUE, palette = "jco"
) + 
  stat_pvalue_manual(
    pwc, hide.ns = TRUE, tip.length = 0, 
    step.increase = 0.1, step.group.by = "variables",
    color = "variables"
  ) +
  labs(
    subtitle = test.label,
    caption = get_pwc_label(pwc, type = "expression")
  )

为什么要检验这些假设、违反后如何处理

1) 多变量正态性（Multivariate normality）

  • 定义：在每个组内，响应变量的联合分布应接近多变量正态分布（joint normal）。这比单变量正态更严格。

  • 为什么重要：MANOVA 的 p 值推断（Wilks, Pillai, Hotelling-Lawley 等统计量）基于正态性假设来推导其分布；严重偏离可能会影响 Type I error。

  • 检验方法：

    • mshapiro_test()（多元 Shapiro-Wilk）：直接检验（注意样本量限制，许多实现要求 n <= 50 或对大样本不适用）。
    • 视觉检查：多元散点图、主成分投影下的 Q-Q 图、或每个变量的单变量 Q-Q 图。
    • 检查偏度/kurtosis（skewness/kurtosis）也是常用手段。

  • 若违背：

    • 对样本量较大（>50），即便轻微偏离，检验显著也不必惊慌——Pillai’s trace 对轻度偏离较稳健。
    • 使用非参数或稳健方法（例如基于排列检验的 MANOVA、或使用稳健协方差估计）。
    • 对变量做变换（log、Box-Cox）以减轻偏态。

2) 协方差矩阵同质性（Homogeneity of covariance matrices across groups）

  • 定义：每个组内响应变量的协方差矩阵应相同。

  • 为什么重要：若协方差不同，MANOVA 的统计量（尤其 Wilks’ Lambda）假设被破坏，可能导致误判。

  • 检验方法：Box’s M 检验（box_m()）。

  • 若违背：

    • Pillai’s trace 较 Wilks 更稳健，通常在 Manova(..., test.statistic="Pillai") 推荐使用。
    • 若组样本量相近且不是极端不平衡，Pillai 仍可使用。
    • 若协方差差异很大且组大小不平衡，建议使用稳健/非参数方法或对单变量做 Welch / Games-Howell 这类对方差不齐稳健的事后检验。
    • 也可使用广义线性模型或混合效应模型并采用群内协方差结构建模。

3) 线性关系（Linearity among dependent variables within each group）

  • 定义：对于 MANOVA 来说，多个响应变量在每个组内应当近似呈线性关系（pairwise 线性），否则结果解释会复杂。
  • 为什么重要：MANOVA 基于线性判别与线性模型，若变量间关系强烈非线性，线性模型可能无法捕获差异。
  • 检验方法：散点图、ggpairs() 可视化。

4) 无多重共线性（Absence of multicollinearity / singularity）

  • 定义：响应变量不能高度相关（接近完全线性关系），否则协方差矩阵接近奇异，逆矩阵计算不稳定。
  • 为什么重要：若高度共线性，MANOVA 的统计推断变得不可靠，某些统计量可能无法计算（矩阵不可逆）。
  • 处理：删除/合并高度相关变量，或使用主成分分析（PCA）提取少数不相关主成分再进行 MANOVA。

5) 独立性（Independence of observations）

  • 定义：观测值之间要独立（没有重复测量、配对或群内相关性）。
  • 为什么重要：若不独立（例如重复测量）, 则应使用重复测量 MANOVA 或混合效应模型。
  • 检验/处理：设计阶段保证随机采样；若数据存在层次结构，用线性混合模型（LMM）或广义 estimating equations（GEE）。

6) 样本量和自由度

  • MANOVA 需要每组样本量足够大以估计协方差矩阵；通常，样本量应大于响应变量的数量（p），并且每组内应至少有 p+1 个观测值。更稳健的方法要求每组样本远大于变量数。

总结

  • 若响应变量互相关联（通常是这样），MANOVA 优于多个独立 ANOVA，因为它考虑相关性并减少多重检验问题。

  • 若 MANOVA 显著：说明在联合响应空间里不同组之间有差异。接下来通常进行事后分析来找出哪些变量或哪些组对差异贡献最大：

    1. 检查各个单变量 ANOVA（或 Welch/Kruskal）看哪些变量显著；
    2. 进行成对比较（Tukey / Games-Howell / Bonferroni 等），并结合效应量（Cohen’s d 或多变量效果量如 partial eta
    3. 如果变量高度相关，考虑用判别分析（LDA）看哪些线性组合最能区分组。

关系图谱

反向链接