第 27 章 Tidy Statistics

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

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

想归想,但本章还是打算介绍基本的方差分析内容。

library(tidyverse)

27.1 从一个案例开始

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

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

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

wages %>% 
  head() %>% 
  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 双因素方差分析

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

  • 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剂量和摄入方式对豚鼠牙齿的影响")

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

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

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

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

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

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