ANOVA

问题

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

方案

假设这是你的数据：

  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. data$subject <- factor(data$subject)

单因素ANOVA分析

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

双因素ANOVA分析

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

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.519139 -3.073025 0.03474709 0.0549625
  10. #> 
  11. #> $age
  12. #>               diff       lwr        upr     p adj
  13. #> young-old -2.31785 -3.846349 -0.7893498 0.0044215
  14. #> 
  15. #> $`sex:age`
  16. #>                       diff        lwr        upr     p adj
  17. #> M:old-F:old      2.8971429 -0.3079526  6.1022384 0.0869856
  18. #> F:young-F:old    2.1733333 -1.1411824  5.4878491 0.2966111
  19. #> M:young-F:old   -2.2183333 -5.1319553  0.6952887 0.1832890
  20. #> F:young-M:old   -0.7238095 -3.7691188  2.3214997 0.9138789
  21. #> M:young-M:old   -5.1154762 -7.7187601 -2.5121923 0.0000676
  22. #> M:young-F:young -4.3916667 -7.1285380 -1.6547953 0.0008841

有受试内变量的ANOVAs

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

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

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

One-way within ANOVA

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

  1. # 独立变量 (within): time
  2. # 依赖变量:          value
  3. aov_time <- aov(value ~ time + Error(subject/time), data=data_long)
  4. summary(aov_time)
  5. #> 
  6. #> Error: subject
  7. #>           Df Sum Sq Mean Sq F value Pr(>F)
  8. #> Residuals 29  261.2   9.009               
  9. #> 
  10. #> Error: subject:time
  11. #>           Df Sum Sq Mean Sq F value   Pr(>F)    
  12. #> time       1 155.53  155.53   71.43 2.59e-09 ***
  13. #> Residuals 29  63.14    2.18                     
  14. #> ---
  15. #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  18. # This won't work here for some reason (?)
  19. model.tables(aov_time, "means")
  20. #> Tables of means
  21. #> Grand mean
  22. #>          
  23. #> 8.093333 
  24. #> 
  25. #>  time 
  26. #> time
  27. #>  after before 
  28. #>  6.483  9.703

混合设计 ANOVA

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

  1. # 2x2 mixed:
  2. # 独立变量（被试间） : age
  3. # 独立变量（被试内） :  time
  4. # 依赖变量:         value
  5. aov_age_time <- aov(value ~ age*time + Error(subject/time), data=data_long)
  6. summary(aov_age_time)
  7. #> 
  8. #> Error: subject
  9. #>           Df Sum Sq Mean Sq F value Pr(>F)
  10. #> age        1  24.44  24.440    2.89    0.1
  11. #> Residuals 28 236.81   8.457               
  12. #> 
  13. #> Error: subject:time
  14. #>           Df Sum Sq Mean Sq F value   Pr(>F)    
  15. #> time       1 155.53  155.53   98.14 1.18e-10 ***
  16. #> age:time   1  18.77   18.77   11.84  0.00184 ** 
  17. #> Residuals 28  44.37    1.58                     
  18. #> ---
  19. #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  21. # This won't work here because the data is unbalanced
  22. model.tables(aov_age_time, "means")
  23. #> Tables of means
  24. #> Grand mean
  25. #>          
  26. #> 8.093333 
  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

更多被试内变量的ANOVA

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

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

更多使用R进行ANOVA分析，参见：