7.3 频率检验

7.3.1 问题

你有分类数据然后想要检验是否这些数据值的频数分布是否与预期不符，或者是否组间的频数分布有（显著）差异。

7.3.2 方案

频数检验通常解决两类问题：

  • 频数分布与预期或者理论的分布（比如 50% 的 yes，50% 的 no ）符合吗？（拟合优度检验）
  • 两组或多组之间的频率分布有差异吗？（独立检验） 通常用于解决这样问题的统计检验方法，分为精确检验与近似检验两种。
期望分布比较组别
精确精确二项检验Fisher精确检验
近似卡方拟合优度独立卡方检验

注意：精确二项检验仅能用于有两个水平的单变量。Fisher 精确检验仅能用于二维列联表（比如，当存在一个独立变量和一个非独立变量时它可以使用；但不能用于两个独立变量和一个非独立变量的情况）。

想要检验配对或被试内效应，我们可以使用 McNemar 检验。使用该检验必须满足存在两个水平的独立变量和两个水平的非独立变量。

想要检验有重复测量的两个变量独立性，我们可以使用 Cochran-Mantel-Haenszel 检验。

假设你有下面的数据，其中每一行代表一个记录：

  1. data <- read.table(header = TRUE, text = "
  2.  condition result
  3.    control      0
  4.    control      0
  5.    control      0
  6.    control      0
  7.  treatment      1
  8.    control      0
  9.    control      0
  10.  treatment      0
  11.  treatment      1
  12.    control      1
  13.  treatment      1
  14.  treatment      1
  15.  treatment      1
  16.  treatment      1
  17.  treatment      0
  18.    control      0
  19.    control      1
  20.    control      0
  21.    control      1
  22.  treatment      0
  23.  treatment      1
  24.  treatment      0
  25.  treatment      0
  26.    control      0
  27.  treatment      1
  28.    control      0
  29.    control      0
  30.  treatment      1
  31.  treatment      0
  32.  treatment      1
  33. ")

相比于以记录的数据框存储，你的数据可能是计数的数据框，或者是一个列联表。

7.3.2.1 拟合优度检验 （期望频率）

7.3.2.1.1 卡方检验

想要检验假设：结果列 result（忽略条件 condition ）中的两个值在总体中几乎相等（50% - 50%）。

  1. # 　为result列创建列联表，包含 0 和 1 两个值 　注意，'0'
  2. # 和 '1' 是列名而不是实际的值
  3. ct <- table(data$result)
  4. ct
  5. #> 
  6. #>  0  1 
  7. #> 17 13
  9. # 也可以手动创建表格
  10. ct <- matrix(c(17, 13), ncol = 2)
  11. colnames(ct) <- c("0", "1")
  13. # 执行卡方检验
  14. chisq.test(ct)
  15. #> 
  16. #>  Chi-squared test for given probabilities
  17. #> 
  18. #> data:  ct
  19. #> X-squared = 0.53, df = 1, p-value = 0.5

想要检验有不同期望频率的样本（比如下面一个 0.75，一个 0.25 ）：

  1. # 概率表 —— 和必须为 1
  2. pt <- c(0.75, 0.25)
  3. chisq.test(ct, p = pt)
  4. #> 
  5. #>  Chi-squared test for given probabilities
  6. #> 
  7. #> data:  ct
  8. #> X-squared = 5.4, df = 1, p-value = 0.02

如果你想要从检验结果中提取信息，可以将结果保存进一个变量，然后用 str() 函数查看变量信息，接着把你想要的部分取出来。例如：

  1. chi_res <- chisq.test(ct, p = pt)
  2. # 查看所有组分
  3. str(chi_res)
  4. #> List of 9
  5. #>  $ statistic: Named num 5.38
  6. #>   ..- attr(*, "names")= chr "X-squared"
  7. #>  $ parameter: Named num 1
  8. #>   ..- attr(*, "names")= chr "df"
  9. #>  $ p.value  : num 0.0204
  10. #>  $ method   : chr "Chi-squared test for given probabilities"
  11. #>  $ data.name: chr "ct"
  12. #>  $ observed : num [1:2] 17 13
  13. #>  $ expected : num [1:2] 22.5 7.5
  14. #>  $ residuals: num [1:2] -1.16 2.01
  15. #>  $ stdres   : num [1:2] -2.32 2.32
  16. #>  - attr(*, "class")= chr "htest"
  18. # 获取卡方值
  19. chi_res$statistic
  20. #> X-squared 
  21. #>     5.378
  23. # 获取p值
  24. chi_res$p.value
  25. #> [1] 0.02039
7.3.2.1.2 精确二项检验

精确二项检验仅能用于存在两个值的单变量数据。

  1. ct <- table(data$result)
  2. ct
  3. #> 
  4. #>  0  1 
  5. #> 17 13
  7. binom.test(ct, p = 0.5)
  8. #> 
  9. #>  Exact binomial test
  10. #> 
  11. #> data:  ct
  12. #> number of successes = 17, number of trials = 30,
  13. #> p-value = 0.6
  14. #> alternative hypothesis: true probability of success is not equal to 0.5
  15. #> 95 percent confidence interval:
  16. #>  0.3743 0.7454
  17. #> sample estimates:
  18. #> probability of success 
  19. #>                 0.5667
  21. # 使用 75% 的期望概率——注意 1 在第二列，所以只需要令 p
  22. # = 0.25
  23. binom.test(ct, p = 0.25)
  24. #> 
  25. #>  Exact binomial test
  26. #> 
  27. #> data:  ct
  28. #> number of successes = 17, number of trials = 30,
  29. #> p-value = 2e-04
  30. #> alternative hypothesis: true probability of success is not equal to 0.25
  31. #> 95 percent confidence interval:
  32. #>  0.3743 0.7454
  33. #> sample estimates:
  34. #> probability of success 
  35. #>                 0.5667

如果你想要从检验结果中提取信息，可以将结果保存进一个变量，然后用 str() 函数查看变量信息，接着把你想要的部分取出来。例如：

  1. bin_res <- binom.test(ct, p = 0.25)
  2. # 字符串格式化后查看信息
  3. str(bin_res)
  4. #> List of 9
  5. #>  $ statistic  : Named num 17
  6. #>   ..- attr(*, "names")= chr "number of successes"
  7. #>  $ parameter  : Named num 30
  8. #>   ..- attr(*, "names")= chr "number of trials"
  9. #>  $ p.value    : Named num 0.000216
  10. #>   ..- attr(*, "names")= chr "0"
  11. #>  $ conf.int   : num [1:2] 0.374 0.745
  12. #>   ..- attr(*, "conf.level")= num 0.95
  13. #>  $ estimate   : Named num 0.567
  14. #>   ..- attr(*, "names")= chr "probability of success"
  15. #>  $ null.value : Named num 0.25
  16. #>   ..- attr(*, "names")= chr "probability of success"
  17. #>  $ alternative: chr "two.sided"
  18. #>  $ method     : chr "Exact binomial test"
  19. #>  $ data.name  : chr "ct"
  20. #>  - attr(*, "class")= chr "htest"
  22. # 获取 p 值
  23. bin_res$p.value
  24. #>         0 
  25. #> 0.0002157
  27. # 获取 95% 置信区间
  28. bin_res$conf.int
  29. #> [1] 0.3743 0.7454
  30. #> attr(,"conf.level")
  31. #> [1] 0.95

7.3.2.2 独立检验（比较组间）

7.3.2.2.1 卡方检验

想要检验控制和处理组结果的频数差异，使用二维列联表。

  1. ct <- table(data$condition, data$result)
  2. ct
  3. #>            
  4. #>              0  1
  5. #>   control   11  3
  6. #>   treatment  6 10
  8. chisq.test(ct)
  9. #> 
  10. #>  Pearson's Chi-squared test with Yates'
  11. #>  continuity correction
  12. #> 
  13. #> data:  ct
  14. #> X-squared = 3.6, df = 1, p-value = 0.06
  16. chisq.test(ct, correct = FALSE)
  17. #> 
  18. #>  Pearson's Chi-squared test
  19. #> 
  20. #> data:  ct
  21. #> X-squared = 5.1, df = 1, p-value = 0.02
7.3.2.2.2 Fisher 精确检验

对于小样本而言 Fisher 精确检验更为适合。小样本的 2x2 列表非常典型，样本更多、更复杂的列表计算强度非常大。当然，用R进行比较复杂的计算也是没有太大问题的。

  1. ct <- table(data$condition, data$result)
  2. ct
  3. #>            
  4. #>              0  1
  5. #>   control   11  3
  6. #>   treatment  6 10
  8. fisher.test(ct)
  9. #> 
  10. #>  Fisher's Exact Test for Count Data
  11. #> 
  12. #> data:  ct
  13. #> p-value = 0.03
  14. #> alternative hypothesis: true odds ratio is not equal to 1
  15. #> 95 percent confidence interval:
  16. #>   0.9669 45.5550
  17. #> sample estimates:
  18. #> odds ratio 
  19. #>      5.714
7.3.2.2.3 Cochran-Mantel-Haenszel 检验

Cochran-Mantel-Haenszel 检验 （或称为 Mantel-Haenszel 检验)）用于检验重复测量两离散变量的独立性。通常使用 2x2xK 列表表示，K是测量条件的次数。比如你想要指导是否一个处理（C vs. D）是否影响了恢复的概率（yes or no）。假设该处理一天监控测量三次——早上、中午和晚上，而你想要你的检验能够控制它。那么你可以使用 CMH 检验对 2x2x3 列联表进行操作，第三个变量是你想要控制的变量。

R 中的 CMH 检验可以处理比 2x2xK 维度更高的数据，例如你处理 3x3xK 列联表。

在接下来的例子里有三个变量：Location、Allele 和 Habitat。问题是——当控制 location 变量时，Allel（94 或非 94）和 Habitat（marine 或 estuarine）两个变量是否独立。

  1. fish <- read.table(header = TRUE, text = "
  2.   Location Allele   Habitat Count
  3.  tillamook     94    marine    56
  4.  tillamook     94 estuarine    69
  5.  tillamook non-94    marine    40
  6.  tillamook non-94 estuarine    77
  7.    yaquina     94    marine    61
  8.    yaquina     94 estuarine   257
  9.    yaquina non-94    marine    57
  10.    yaquina non-94 estuarine   301
  11.    alsea     94    marine    73
  12.    alsea     94 estuarine    65
  13.    alsea non-94    marine    71
  14.    alsea non-94 estuarine    79
  15.   umpqua     94    marine    71
  16.   umpqua     94 estuarine    48
  17.   umpqua non-94    marine    55
  18.   umpqua non-94 estuarine    48
  19. ")

注意上面的数据是计数的数据框，而不是像之前的例子是记录的数据框。这里我们使用 xtabs() 函数将它转换为列联表。

  1. # 制造一个三维的列联表，最后一个变量时要控制的 Location
  2. # 变量
  3. ct <- xtabs(Count ~ Allele + Habitat + Location, data = fish)
  4. ct
  5. #> , , Location = alsea
  6. #> 
  7. #>         Habitat
  8. #> Allele   estuarine marine
  9. #>   94            65     73
  10. #>   non-94        79     71
  11. #> 
  12. #> , , Location = tillamook
  13. #> 
  14. #>         Habitat
  15. #> Allele   estuarine marine
  16. #>   94            69     56
  17. #>   non-94        77     40
  18. #> 
  19. #> , , Location = umpqua
  20. #> 
  21. #>         Habitat
  22. #> Allele   estuarine marine
  23. #>   94            48     71
  24. #>   non-94        48     55
  25. #> 
  26. #> , , Location = yaquina
  27. #> 
  28. #>         Habitat
  29. #> Allele   estuarine marine
  30. #>   94           257     61
  31. #>   non-94       301     57
  33. # 以扁平化显示
  34. ftable(ct)
  35. #>                  Location alsea tillamook umpqua yaquina
  36. #> Allele Habitat                                          
  37. #> 94     estuarine             65        69     48     257
  38. #>        marine                73        56     71      61
  39. #> non-94 estuarine             79        77     48     301
  40. #>        marine                71        40     55      57
  42. # 按指定方式进行变量输出
  43. ftable(ct, row.vars = c("Location", "Allele"), col.vars = "Habitat")
  44. #>                  Habitat estuarine marine
  45. #> Location  Allele                         
  46. #> alsea     94                    65     73
  47. #>           non-94                79     71
  48. #> tillamook 94                    69     56
  49. #>           non-94                77     40
  50. #> umpqua    94                    48     71
  51. #>           non-94                48     55
  52. #> yaquina   94                   257     61
  53. #>           non-94               301     57
  55. mantelhaen.test(ct)
  56. #> 
  57. #>  Mantel-Haenszel chi-squared test with
  58. #>  continuity correction
  59. #> 
  60. #> data:  ct
  61. #> Mantel-Haenszel X-squared = 5, df = 1, p-value =
  62. #> 0.02
  63. #> alternative hypothesis: true common odds ratio is not equal to 1
  64. #> 95 percent confidence interval:
  65. #>  0.6006 0.9593
  66. #> sample estimates:
  67. #> common odds ratio 
  68. #>             0.759

根据检验结果，当控制 Location 变量时 Allele 与 Habitat 变量存在相关（p = 0.025）。

注意列联表的前两个维度处理是一致的，所以前后顺序变化都不会影响结果。而最后一个变量变化会导致结果的不同，下面是一个实例。

  1. # 下面两个看似不同的列联表，实际检验结果相同
  2. ct.1 <- xtabs(Count ~ Habitat + Allele + Location, data = fish)
  3. ct.2 <- xtabs(Count ~ Allele + Habitat + Location, data = fish)
  4. mantelhaen.test(ct.1)
  5. #> 
  6. #>  Mantel-Haenszel chi-squared test with
  7. #>  continuity correction
  8. #> 
  9. #> data:  ct.1
  10. #> Mantel-Haenszel X-squared = 5, df = 1, p-value =
  11. #> 0.02
  12. #> alternative hypothesis: true common odds ratio is not equal to 1
  13. #> 95 percent confidence interval:
  14. #>  0.6006 0.9593
  15. #> sample estimates:
  16. #> common odds ratio 
  17. #>             0.759
  19. mantelhaen.test(ct.2)
  20. #> 
  21. #>  Mantel-Haenszel chi-squared test with
  22. #>  continuity correction
  23. #> 
  24. #> data:  ct.2
  25. #> Mantel-Haenszel X-squared = 5, df = 1, p-value =
  26. #> 0.02
  27. #> alternative hypothesis: true common odds ratio is not equal to 1
  28. #> 95 percent confidence interval:
  29. #>  0.6006 0.9593
  30. #> sample estimates:
  31. #> common odds ratio 
  32. #>             0.759
  34. # 把 Allele 放到最后，结果不同了
  35. ct.3 <- xtabs(Count ~ Location + Habitat + Allele, data = fish)
  36. ct.4 <- xtabs(Count ~ Habitat + Location + Allele, data = fish)
  37. mantelhaen.test(ct.3)
  38. #> 
  39. #>  Cochran-Mantel-Haenszel test
  40. #> 
  41. #> data:  ct.3
  42. #> Cochran-Mantel-Haenszel M^2 = 168, df = 3,
  43. #> p-value <2e-16
  45. mantelhaen.test(ct.4)
  46. #> 
  47. #>  Cochran-Mantel-Haenszel test
  48. #> 
  49. #> data:  ct.4
  50. #> Cochran-Mantel-Haenszel M^2 = 168, df = 3,
  51. #> p-value <2e-16
  53. # 把 Habitat 放最后，结果也不同
  54. ct.5 <- xtabs(Count ~ Allele + Location + Habitat, data = fish)
  55. ct.6 <- xtabs(Count ~ Location + Allele + Habitat, data = fish)
  56. mantelhaen.test(ct.5)
  57. #> 
  58. #>  Cochran-Mantel-Haenszel test
  59. #> 
  60. #> data:  ct.5
  61. #> Cochran-Mantel-Haenszel M^2 = 2, df = 3, p-value
  62. #> = 0.6
  64. mantelhaen.test(ct.6)
  65. #> 
  66. #>  Cochran-Mantel-Haenszel test
  67. #> 
  68. #> data:  ct.6
  69. #> Cochran-Mantel-Haenszel M^2 = 2, df = 3, p-value
  70. #> = 0.6

7.3.2.3 McNemar 检验

McNemar 检验概念上是频数数据的一个被试内检验。例如，假设你想要检验是否一个处理增加了一个人对某个问题反应「yes」的概率，而且你只有每个人处理前和处理后的数据。标准的卡方检验将不合适，因为它假设了组别是独立的。取而代之，我们可以使用 McNemar 检验。该检验仅适用于当存在一个独立变量的两次测量时。用于 McNemar 的列联表与用于卡方检验的非常相似，但结构上是不同的。

假设你有下面的数据。每个对象有处理前和后的反应。

  1. data <- read.table(header = TRUE, text = "
  2.  subject time result
  3.      1  pre      0
  4.      1 post      1
  5.      2  pre      1
  6.      2 post      1
  7.      3  pre      0
  8.      3 post      1
  9.      4  pre      1
  10.      4 post      0
  11.      5  pre      1
  12.      5 post      1
  13.      6  pre      0
  14.      6 post      1
  15.      7  pre      0
  16.      7 post      1
  17.      8  pre      0
  18.      8 post      1
  19.      9  pre      0
  20.      9 post      1
  21.     10  pre      1
  22.     10 post      1
  23.     11  pre      0
  24.     11 post      0
  25.     12  pre      1
  26.     12 post      1
  27.     13  pre      0
  28.     13 post      1
  29.     14  pre      0
  30.     14 post      0
  31.     15  pre      0
  32.     15 post      1
  33. ")
  1. library(tidyr)
  2. data_wide <- spread(data, time, result)
  3. data_wide
  4. #>    subject post pre
  5. #> 1        1    1   0
  6. #> 2        2    1   1
  7. #> 3        3    1   0
  8. #> 4        4    0   1
  9. #> 5        5    1   1
  10. #> 6        6    1   0
  11. #> 7        7    1   0
  12. #> 8        8    1   0
  13. #> 9        9    1   0
  14. #> 10      10    1   1
  15. #> 11      11    0   0
  16. #> 12      12    1   1
  17. #> 13      13    1   0
  18. #> 14      14    0   0
  19. #> 15      15    1   0

接下来从数据框的 pre 和 post 列生成列联表：

  1. ct <- table(data_wide[, c("pre", "post")])
  2. ct
  3. #>    post
  4. #> pre 0 1
  5. #>   0 2 8
  6. #>   1 1 4
  8. # 下面是用于标准卡方检验的列联表，注意差别
  9. table(data[, c("time", "result")])
  10. #>       result
  11. #> time    0  1
  12. #>   post  3 12
  13. #>   pre  10  5

执行检验：

  1. mcnemar.test(ct)
  2. #> 
  3. #>  McNemar's Chi-squared test with continuity
  4. #>  correction
  5. #> 
  6. #> data:  ct
  7. #> McNemar's chi-squared = 4, df = 1, p-value =
  8. #> 0.05

对于小样本，它会使用连续校正。我们可以使用精确校正的 McNemar 检验替换这种校正方式，前者更加的精确，可通过 exact2x2 包获取。

  1. library(exact2x2)
  2. #> Loading required package: exactci
  3. #> Loading required package: ssanv
  5. mcnemar.exact(ct)
  6. #> 
  7. #>  Exact McNemar test (with central confidence
  8. #>  intervals)
  9. #> 
  10. #> data:  ct
  11. #> b = 8, c = 1, p-value = 0.04
  12. #> alternative hypothesis: true odds ratio is not equal to 1
  13. #> 95 percent confidence interval:
  14. #>    1.073 354.981
  15. #> sample estimates:
  16. #> odds ratio 
  17. #>          8