双因素方差分析之类检验同样要求数据满足正态性、方差齐性，但是如果数据不满足这些条件，目前来讲没有很好的解决方案（对于方差分析来讲），因此，就像大佬说的一样，只能在检验之前默默祈祷，方法满足齐性条件或者能够通过数据转换解决，实在不行就换种方法吧。

来源：http://www.biostathandbook.com/homoscedasticity.html

这里先提一下上面讲的Ⅲ型平方和、Ⅱ型平方和与Ⅰ型平方和，详细的计算过程可参考 Daniel Wollschläger. Sum of Squares Type I, II, III: the underlying hypotheses, model comparisons, and their calculation in R.

• Ⅰ型平方和在进行双因素方差分析时，因子型变量的顺序会影响计算结果；
• Ⅲ型平方和则计算的是简单主效应，模型中每个因子在进行检验时，会控制另一个因子的水平及交互作用相同；
• Ⅱ型平方和与Ⅲ型平方和相似，但是不会控制交互作用

具体选择建议：在两因素交互作用显著时，建议使用Ⅲ型平方和，而在交互作用不显著时，建议使用Ⅱ型平方和，如果是平衡的实验设计（因子水平相同），那恭喜你不用考虑这个问题，三种平方和计算结果相同，可进一步查看 湿地生态环境研究 公众号的相关推文。

下面进行代码展示

方法1（仍然使用rstatix包）：

data("jobsatisfaction", package = "datarium")
#描述性统计
jobsatisfaction %>%
  group_by(gender, education_level) %>%
  get_summary_stats(score, type = "mean_sd")
#可视化
(bxp <- ggboxplot(
  jobsatisfaction, x = "gender", y = "score",
  color = "education_level", palette = "jco"
))
#前提条件检验
# Build the linear model
model  <- lm(score ~ gender*education_level,
             data = jobsatisfaction)
# Create a QQ plot of residuals
ggqqplot(residuals(model))
# 正态性检验
shapiro_test(residuals(model))
#或者
jobsatisfaction %>%
  group_by(gender, education_level) %>%
  shapiro_test(score)
#方差齐性检验
jobsatisfaction %>% levene_test(score ~ gender*education_level)

数据残差正态性、方差齐性都满足

#双因素方差分析
options(contrasts = c("contr.sum", "contr.poly"))#使用III型平方和时需先运行这个，对其他类型平方和没影响
(res.aov <- jobsatisfaction %>% anova_test(score ~ gender * education_level,type=3))#type=3表示使用Ⅲ型平方和

交互作用显著，只能先进行简单主效应分析

# Group the data by gender and fit  anova
model <- lm(score ~ gender * education_level, data = jobsatisfaction)
jobsatisfaction %>%
  group_by(gender) %>%
  anova_test(score ~ education_level, error = model)
#The argument error is used to specify the ANOVA model from which the pooled error sum of squares and degrees of freedom are to be calculated.

简单主效应仍然显著，进行多重比较

#简单主效应显著后，进行多重比较
(pwc <- jobsatisfaction %>% 
  group_by(gender) %>%
  emmeans_test(score ~ education_level, p.adjust.method = "bonferroni") )

#可视化
pwc <- pwc %>% add_xy_position(x = "gender")
bxp +
  stat_pvalue_manual(pwc) +
  labs(
    subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc)
  )

方法2（使用car包）：

#方法2：
model1  <- lm(score ~ gender*education_level,
             data = jobsatisfaction)
#残差正态性检验
shapiro.test(resid(model1))
#方差齐性检验
library(car)
leveneTest(score ~ gender*education_level,
           data = jobsatisfaction)
#显著性检验
options(contrasts = c("contr.sum", "contr.poly"))#使用III型平方和时需先运行这个，对其他类型平方和没影响
Anova(model1,type=3)
summary(model1)
#交互作用显著，先做简单主效应分析
joint_tests(model1,by="gender")
joint_tests(model1,by="education_level")
#多重比较
library(emmeans)
#控制education_level，比较gender各水平
emmeans(model1,pairwise~gender|education_level,adjust="bonferroni")
#控制gender，比较education_level各水平
emmeans(model1,pairwise~education_level|gender,adjust="bonferroni")

---

• 如交互作用不显著（大同小异，不作赘述）

# -------------------------------------------------------------------------
#若交互作用不显著时，换成II型平方和，仍用上述数据仅做示例
(res.aov <- jobsatisfaction %>% 
   anova_test(score ~ gender * education_level,type=2))
#或者
Anova(model1,type=2)
#对应的多重比较的方法：
jobsatisfaction %>% 
  emmeans_test(
    score ~ education_level, p.adjust.method = "bonferroni",
    model = model1
  )
#或者
jobsatisfaction %>%
  pairwise_t_test(
    score ~ education_level, 
    p.adjust.method = "bonferroni"
  )
#或者
emmeans(model1,pairwise~education_level,adjust="bonferroni")
#或者
TukeyHSD(aov(model1))$education_level

参考资料：

[1].https://www-dwoll-de.translate.goog/r/ssTypes.php?_x_tr_sch=http&_x_tr_sl=en&_x_tr_tl=zh-CN&_x_tr_hl=zh-CN&_x_tr_pto=op,sc

[2].https://mp.weixin.qq.com/s/edRFYK8tRwS9Gn_PdUccLA

[3].http://www.biostathandbook.com/homoscedasticity.html

[4].https://www.datanovia.com/en/lessons/anova-in-r/