tidystats_aov.Rmd

# Tidy Statistics {#tidystats-aov}

```{r, include=FALSE}
knitr::opts_chunk$set(
   echo         = TRUE, 
   warning      = FALSE, 
   message      = FALSE,
   fig.showtext = TRUE
)
```

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


对非统计专业的初学者（比如我）感觉 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个对象关于收入、身高、教育水平等信息的数据集，数据在课件首页下载。

首先，我们下载后导入数据
```{r tidystats-1, message = FALSE, warning = FALSE}
library(tidyverse)
wages <- read_csv("./demo_data/wages.csv")

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

我们的问题：男性是否就比女性挣的多？
<!-- ## 单样本的t检验 -->
<!-- <http://www.biye5u.com/article/R/2019/6398.html> -->

## 单因素方差分析

```{r tidystats-2}
t.test(earn ~ sex, data = wages)
```


```{r tidystats-3}
lm(earn ~ sex, data = wages) %>% 
  summary()
```


```{r tidystats-4}
aov(earn ~ sex, data = wages) %>% 
  summary()

```


## 双因素方差分析

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

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


```{r tidystats-5, message=FALSE, warning=FALSE}
library(ggpubr)

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

my_data %>% head()
```



```{r tidystats-6}
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)

```{r tidystats-7}
aov(len ~ supp + dose, data = my_data) %>%
  broom::tidy()
```


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

```{r tidystats-8}
aov(len ~ supp + dose, data = my_data) %>%
  TukeyHSD(which = "dose") %>%
  broom::tidy()
```




```{r tidystats-9}
aov(len ~ supp + dose, data = my_data) %>%
  TukeyHSD(which = "supp") %>%
  broom::tidy()
```

思考：交互效应是否显著？
```{r tidystats-10}
aov(len ~ supp * dose, data = my_data) %>%
  broom::tidy()
```

## 在tidyverse中的应用

我们也可以配合强大的tidyverse函数，完成不同分组下的方差分析，比如
```{r tidystats-11}
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)
```

更多使用可参考第 \@ref(tidyverse-colwise) 章。

## 使用rstatix包

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

```{r}
library(rstatix)

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


```{r}
# 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)

```

```{r}
# One-way ANOVA test

my_data %>% 
  anova_test(len ~ dose)
```


```{r}
# Two-way ANOVA test

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



```{r}
# 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))
```





```{r tidystats-12, echo = F}
# remove the objects
# rm(list=ls())
rm(my_data, wages)
```

```{r tidystats-13, echo = F, message = F, warning = F, results = "hide"}
pacman::p_unload(pacman::p_loaded(), character.only = TRUE)
```