方差分析

ANOVA 多组比较

统计建模篇

一个事实是，用统计的，往往不是学统计的。

对非统计专业的初学者（比如我）感觉 t-test， ANOVAs， Chi-Square test等太不友好了，每次用的时候，我都要去翻书看我用对了没有，还要担心p-value是否徘徊在0.05附近。或许，从t-test等统计检验方法开始学统计是个错误的开始。我有时候在想，我们是不是应该更关心模型的理解，或者模型背后的理论呢。（如果和我的想法一样，就跳过本章吧）

想归想，但同学们对这方面的需求很大，所以还是打算介绍基本的方差分析内容。

方法的区分

比较几组数据之间是否有显著性差异，最常用的方法有以下几种

X变量类型	X组别数量	Y变量类型	分析方法	R语法
定类	2组或者多组	定量	方差	`aov()`
定类	仅仅2组	定量	t检验	`t.test()`
定类	2组或者多组	定类	卡方	`chisq.test()`

根据X变量的个数，方差分析又分为单因素方差分析和多因素方差分析，当X的个数（不是组别数量）为1个时，我们称之为单因素方差；X的个数为2个时，则为双因素方差。

从一个案例开始

从这是一份1994年收集1379个对象关于收入、身高、教育水平等信息的数据集，数据在课件首页下载。

首先，我们下载后导入数据

library(tidyverse)
wages <- read_csv("./demo_data/wages.csv")

wages %>% 
  head() %>% 
  knitr::kable()

我们的问题：男性是否就比女性挣的多？

单因素方差分析

t.test(earn ~ sex, data = wages)

lm(earn ~ sex, data = wages) %>% 
  summary()

aov(earn ~ sex, data = wages) %>% 
  summary()

双因素方差分析

我们采用ggpubr宏包下的ToothGrowth来说明，这个数据集包含60个样本，记录着每10只豚鼠在不同的喂食方法和不同的药物剂量下，牙齿的生长情况.

len : 牙齿长度
supp : 两种喂食方法 (橙汁和维生素C)
dose : 抗坏血酸剂量 (0.5, 1, and 2 mg)

library(ggpubr)

my_data <- ToothGrowth %>%
  mutate(
    across(c(supp, dose), ~ as_factor(.x))
  )

my_data %>% head()

my_data %>%
  ggplot(aes(x = supp, y = len, fill = supp)) +
  geom_boxplot(position = position_dodge()) +
  facet_wrap(vars(dose)) +
  labs(title = "Effects of VC dose and intake mode on teeth of guinea pigs")

问题：豚鼠牙齿的长度是否与药物的食用方法和剂量有关？

线性回归时，我们是通过独立变量来预测响应变量，但现在我们关注的重点会从预测转向不同组别差异之间的分析，这即为方差分析（ANOVA）。

这里是两个解释变量，所以问题需要双因素方差分析 (ANOVA)

aov(len ~ supp + dose, data = my_data) %>%
  broom::tidy()

检验表明不同类型之间存在显著差异，但是并没有告诉我们具体谁与谁之间的不同。需要多重比较帮助我们解决这个问题。使用TurkeyHSD函数

aov(len ~ supp + dose, data = my_data) %>%
  TukeyHSD(which = "dose") %>%
  broom::tidy()

aov(len ~ supp + dose, data = my_data) %>%
  TukeyHSD(which = "supp") %>%
  broom::tidy()

思考：交互效应是否显著？

aov(len ~ supp * dose, data = my_data) %>%
  broom::tidy()

在tidyverse中的应用

我们也可以配合强大的tidyverse函数，完成不同分组下的方差分析，比如

mtcars %>%
  group_by(cyl) %>%
  summarise(
    broom::tidy(aov(mpg ~ gear, data = cur_data())),
    .groups = "keep"
  ) %>% 
  select(term, statistic, p.value) %>% 
  filter(term != "Residuals") %>% 
  arrange(p.value)

更多使用可参考相关章节。

使用rstatix包

rstatix包吸收了tidyverse的设计哲学，让熟悉dplyr语法的用户能更方便的完成t-test, Wilcoxon test, ANOVA, Kruskal-Wallis等基础统计检验，同时也增强了代码的可读性，在实际应用中还是挺受用户欢迎的，比如这本书。下面，就ToothGrowth数据，列举rstatix包的一些用法。

library(rstatix)

my_data %>% 
  group_by(dose) %>%
  t_test(len ~ 1, mu = 0)

# T-test
stat.test <- my_data %>% 
  t_test(len ~ supp, paired = FALSE) 

# Create a box plot
p <- ggboxplot(
  my_data, x = "supp", y = "len", 
  color = "supp", palette = "jco", ylim = c(0,40)
  )

# Add the p-value manually
p + 
  stat_pvalue_manual(stat.test, label = "p", y.position = 35) +
  stat_pvalue_manual(stat.test, label = "T-test, p = {p}", 
                     y.position = 36)

# One-way ANOVA test

my_data %>% 
  anova_test(len ~ dose)

# Two-way ANOVA test

my_data %>% 
  anova_test(len ~ supp*dose)

# Two-way repeated measures ANOVA

my_data %>% 
  mutate(id = rep(1:10, 6) ) %>%  # Add individuals id
  anova_test(dv = len, wid = id, within = c(supp, dose))

# remove the objects
# rm(list=ls())
rm(my_data, wages)

pacman::p_unload(pacman::p_loaded(), character.only = TRUE)