7.4 ANOVA

7.4.1 问题

你想要使用 ANOVA 比较多组之间的差异。

7.4.2 方案

假设这是你的数据：

  1. data <- read.table(header = TRUE, text = "
  2.  subject sex   age before after
  3.      1   F   old    9.5   7.1
  4.      2   M   old   10.3  11.0
  5.      3   M   old    7.5   5.8
  6.      4   F   old   12.4   8.8
  7.      5   M   old   10.2   8.6
  8.      6   M   old   11.0   8.0
  9.      7   M young    9.1   3.0
  10.      8   F young    7.9   5.2
  11.      9   F   old    6.6   3.4
  12.     10   M young    7.7   4.0
  13.     11   M young    9.4   5.3
  14.     12   M   old   11.6  11.3
  15.     13   M young    9.9   4.6
  16.     14   F young    8.6   6.4
  17.     15   F young   14.3  13.5
  18.     16   F   old    9.2   4.7
  19.     17   M young    9.8   5.1
  20.     18   F   old    9.9   7.3
  21.     19   F young   13.0   9.5
  22.     20   M young   10.2   5.4
  23.     21   M young    9.0   3.7
  24.     22   F young    7.9   6.2
  25.     23   M   old   10.1  10.0
  26.     24   M young    9.0   1.7
  27.     25   M young    8.6   2.9
  28.     26   M young    9.4   3.2
  29.     27   M young    9.7   4.7
  30.     28   M young    9.3   4.9
  31.     29   F young   10.7   9.8
  32.     30   M   old    9.3   9.4
  33. ")
  35. # 确保 subject
  36. # 列是一个因子变量，这样不会当作连续变量对待
  37. data$subject <- factor(data$subject)

7.4.2.1 单因素 ANOVA 分析

  1. # 单因素： 独立变量: sex 依赖变量: before
  2. aov1 <- aov(before ~ sex, data = data)
  3. summary(aov1)
  4. #>             Df Sum Sq Mean Sq F value Pr(>F)
  5. #> sex          1    1.5    1.53    0.57   0.46
  6. #> Residuals   28   74.7    2.67
  8. # 显示均值
  9. model.tables(aov1, "means")
  10. #> Tables of means
  11. #> Grand mean
  12. #>       
  13. #> 9.703 
  14. #> 
  15. #>  sex 
  16. #>      F      M
  17. #>     10  9.532
  18. #> rep 11 19.000

7.4.2.2 双因素 ANOVA 分析

  1. # 2x2: 独立变量: sex 独立变量: age 依赖变量: after
  2. # 下面两种调用方式等价：
  3. aov2 <- aov(after ~ sex * age, data = data)
  4. aov2 <- aov(after ~ sex + age + sex:age, data = data)
  5. summary(aov2)
  6. #>             Df Sum Sq Mean Sq F value  Pr(>F)    
  7. #> sex          1   16.1    16.1    4.04  0.0550 .  
  8. #> age          1   39.0    39.0    9.79  0.0043 ** 
  9. #> sex:age      1   89.6    89.6   22.51 6.6e-05 ***
  10. #> Residuals   26  103.5     4.0                    
  11. #> ---
  12. #> Signif. codes:  
  13. #> 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  15. # 显示均值
  16. model.tables(aov2, "means")
  17. #> Tables of means
  18. #> Grand mean
  19. #>       
  20. #> 6.483 
  21. #> 
  22. #>  sex 
  23. #>          F      M
  24. #>      7.445  5.926
  25. #> rep 11.000 19.000
  26. #> 
  27. #>  age 
  28. #>        old  young
  29. #>      7.874  5.556
  30. #> rep 12.000 18.000
  31. #> 
  32. #>  sex:age 
  33. #>      age
  34. #> sex   old    young 
  35. #>   F    6.260  8.433
  36. #>   rep  5.000  6.000
  37. #>   M    9.157  4.042
  38. #>   rep  7.000 12.000

7.4.2.3 Tukey HSD post-hoc 检验

  1. TukeyHSD(aov2)
  2. #>   Tukey multiple comparisons of means
  3. #>     95% family-wise confidence level
  4. #> 
  5. #> Fit: aov(formula = after ~ sex + age + sex:age, data = data)
  6. #> 
  7. #> $sex
  8. #>       diff    lwr     upr p adj
  9. #> M-F -1.519 -3.073 0.03475 0.055
  10. #> 
  11. #> $age
  12. #>             diff    lwr     upr  p adj
  13. #> young-old -2.318 -3.846 -0.7893 0.0044
  14. #> 
  15. #> $`sex:age`
  16. #>                    diff    lwr     upr  p adj
  17. #> M:old-F:old      2.8971 -0.308  6.1022 0.0870
  18. #> F:young-F:old    2.1733 -1.141  5.4878 0.2966
  19. #> M:young-F:old   -2.2183 -5.132  0.6953 0.1833
  20. #> F:young-M:old   -0.7238 -3.769  2.3215 0.9139
  21. #> M:young-M:old   -5.1155 -7.719 -2.5122 0.0001
  22. #> M:young-F:young -4.3917 -7.129 -1.6548 0.0009

7.4.3 有受试内变量的 ANOVAs

对于有受试内变量的 ANOVA 分析，数据必须满足为长格式。上面提到的数据都是宽格式，所以我们需要先转换数据格式（参见长宽格式互相转换）。

同样地，有受试内变量的 ANOVA 分析需要一个识别列。当前数据里是 subject 列。识别变量必须是一个因子，如果是数值类型，函数会解析错误导致不能正常工作。

  1. library(tidyr)
  3. # 原始数据 subject sex age before after 1 F old 9.5 7.1
  4. # 2 M old 10.3 11.0 3 M old 7.5 5.8
  6. # 转换为长格式
  7. data_long <- gather(data, time, value, before:after)
  8. # Look at first few rows
  9. head(data_long)
  10. #>   subject sex age   time value
  11. #> 1       1   F old before   9.5
  12. #> 2       2   M old before  10.3
  13. #> 3       3   M old before   7.5
  14. #> 4       4   F old before  12.4
  15. #> 5       5   M old before  10.2
  16. #> 6       6   M old before  11.0
  18. # 确保subject列是一个因子
  19. data_long$subject <- factor(data_long$subject)

7.4.3.1 单因素被试内 ANOVA

首先，像上面展示的一样将数据从宽格式转换到长格式并确保 subject 列是因子变量。如果 subject 是数值向量，而不是因子，结果将会出错。

  1. # 独立变量 (被试内): time 依赖变量: value
  2. aov_time <- aov(value ~ time + Error(subject/time), data = data_long)
  3. summary(aov_time)
  4. #> 
  5. #> Error: subject
  6. #>           Df Sum Sq Mean Sq F value Pr(>F)
  7. #> Residuals 29    261    9.01               
  8. #> 
  9. #> Error: subject:time
  10. #>           Df Sum Sq Mean Sq F value  Pr(>F)    
  11. #> time       1  155.5   155.5    71.4 2.6e-09 ***
  12. #> Residuals 29   63.1     2.2                    
  13. #> ---
  14. #> Signif. codes:  
  15. #> 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  17. # 因为一些原因，下面的代码不工作
  18. model.tables(aov_time, "means")
  19. #> Tables of means
  20. #> Grand mean
  21. #>       
  22. #> 8.093 
  23. #> 
  24. #>  time 
  25. #> time
  26. #>  after before 
  27. #>  6.483  9.703

7.4.3.2 混合设计 ANOVA

首先，像上面展示的一样将数据从宽格式转换到长格式并确保 subject 列是因子变量。

  1. # 2x2 mixed: 独立变量（被试间） : age 独立变量（被试内）
  2. # : time 依赖变量: value
  3. aov_age_time <- aov(value ~ age * time + Error(subject/time), 
  4.   data = data_long)
  5. summary(aov_age_time)
  6. #> 
  7. #> Error: subject
  8. #>           Df Sum Sq Mean Sq F value Pr(>F)
  9. #> age        1   24.4   24.44    2.89    0.1
  10. #> Residuals 28  236.8    8.46               
  11. #> 
  12. #> Error: subject:time
  13. #>           Df Sum Sq Mean Sq F value  Pr(>F)    
  14. #> time       1  155.5   155.5    98.1 1.2e-10 ***
  15. #> age:time   1   18.8    18.8    11.8  0.0018 ** 
  16. #> Residuals 28   44.4     1.6                    
  17. #> ---
  18. #> Signif. codes:  
  19. #> 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  21. # 因为数据不平衡，下面代码不会工作
  22. model.tables(aov_age_time, "means")
  23. #> Tables of means
  24. #> Grand mean
  25. #>       
  26. #> 8.093 
  27. #> 
  28. #>  age 
  29. #>        old  young
  30. #>      8.875  7.572
  31. #> rep 24.000 36.000
  32. #> 
  33. #>  time 
  34. #>      after before
  35. #>      6.483  9.703
  36. #> rep 30.000 30.000
  37. #> 
  38. #>  age:time 
  39. #>        time
  40. #> age     after  before
  41. #>   old    7.950  9.800
  42. #>   rep   12.000 12.000
  43. #>   young  5.506  9.639
  44. #>   rep   18.000 18.000

7.4.3.3 更多被试内变量的 ANOVA

下面这些例子使用的不是上面的数据，但可以解释怎么进行相应的处理。首先，像上面展示的一样将数据从宽格式转换到长格式并确保 subject 列是因子变量。

  1. # # 两个被试内变量 aov.ww <- aov(y ~ w1*w2 +
  2. # Error(subject/(w1*w2)), data=data_long) #
  3. # 1个被试间变量，两个被试内变量 aov.bww <- aov(y ~
  4. # b1*w1*w2 + Error(subject/(w1*w2)) + b1,
  5. # data=data_long) # 两个被试间变量，一个被试内变量
  6. # aov.bww <- aov(y ~ b1*b2*w1 + Error(subject/(w1)) +
  7. # b1*b2, data=data_long)