# 第 27 章 Tidy Statistics

library(tidyverse)

## 27.1 从一个案例开始

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

wages %>%
knitr::kable()
earn height sex race ed age
79571 73.89 male white 16 49
96397 66.23 female white 16 62
48711 63.77 female white 16 33
80478 63.22 female other 16 95
82089 63.08 female white 17 43
15313 64.53 female white 15 30

## 27.2 单因素方差分析

t.test(earn ~ sex, data = wages)
##
## 	Welch Two Sample t-test
##
## data:  earn by sex
## t = -12, df = 768, p-value <2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -25324 -18171
## sample estimates:
## mean in group female   mean in group male
##                24246                45993
lm(earn ~ sex, data = wages) %>%
summary()
##
## Call:
## lm(formula = earn ~ sex, data = wages)
##
## Residuals:
##    Min     1Q Median     3Q    Max
## -46092 -20516  -4639  11722 271956
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    24246       1004    24.1   <2e-16 ***
## sexmale        21748       1636    13.3   <2e-16 ***
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 29400 on 1377 degrees of freedom
## Multiple R-squared:  0.114,	Adjusted R-squared:  0.113
## F-statistic:  177 on 1 and 1377 DF,  p-value: <2e-16
aov(earn ~ sex, data = wages) %>%
summary()
##               Df   Sum Sq  Mean Sq F value Pr(>F)
## sex            1 1.53e+11 1.53e+11     177 <2e-16 ***
## Residuals   1377 1.19e+12 8.66e+08
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 27.3 双因素方差分析

• len : 牙齿长度
• supp : 两种喂食方法 (橙汁和维生素C)
• dose : 抗坏血酸剂量 (0.5, 1, and 2 mg)
library("ggpubr")
my_data <- ToothGrowth %>%
mutate_at(vars(supp, dose), ~ as_factor(.))

my_data %>% head()
##    len supp dose
## 1  4.2   VC  0.5
## 2 11.5   VC  0.5
## 3  7.3   VC  0.5
## 4  5.8   VC  0.5
## 5  6.4   VC  0.5
## 6 10.0   VC  0.5
my_data %>%
ggplot(aes(x = supp, y = len, fill = supp)) +
geom_boxplot(position = position_dodge()) +
facet_wrap(vars(dose)) +
labs(title = "VC剂量和摄入方式对豚鼠牙齿的影响")

aov(len ~ supp + dose, data = my_data) %>%
broom::tidy()
## # A tibble: 3 x 6
##   term         df sumsq meansq statistic   p.value
##   <chr>     <dbl> <dbl>  <dbl>     <dbl>     <dbl>
## 1 supp          1  205.  205.       14.0  4.29e- 4
## 2 dose          2 2426. 1213.       82.8  1.87e-17
## 3 Residuals    56  820.   14.7      NA   NA

aov(len ~ supp + dose, data = my_data) %>%
TukeyHSD(which = "dose") %>%
broom::tidy()
## # A tibble: 3 x 7
##   term  contrast null.value estimate conf.low conf.high
##   <chr> <chr>         <dbl>    <dbl>    <dbl>     <dbl>
## 1 dose  1-0.5             0     9.13     6.22     12.0
## 2 dose  2-0.5             0    15.5     12.6      18.4
## 3 dose  2-1               0     6.37     3.45      9.28
## # ... with 1 more variable: adj.p.value <dbl>
aov(len ~ supp + dose, data = my_data) %>%
TukeyHSD(which = "supp") %>%
broom::tidy()
## # A tibble: 1 x 7
##   term  contrast null.value estimate conf.low conf.high
##   <chr> <chr>         <dbl>    <dbl>    <dbl>     <dbl>
## 1 supp  VC-OJ             0     -3.7    -5.68     -1.72
## # ... with 1 more variable: adj.p.value <dbl>

aov(len ~ supp * dose, data = my_data) %>%
broom::tidy()
## # A tibble: 4 x 6
##   term         df sumsq meansq statistic   p.value
##   <chr>     <dbl> <dbl>  <dbl>     <dbl>     <dbl>
## 1 supp          1  205.  205.      15.6   2.31e- 4
## 2 dose          2 2426. 1213.      92.0   4.05e-18
## 3 supp:dose     2  108.   54.2      4.11  2.19e- 2
## 4 Residuals    54  712.   13.2     NA    NA