12.2　用coin包做置换检验

对于独立性问题，coin包提供了一个进行置换检验的一般性框架。通过该包，你可以回答如下问题。

  • 响应值与组的分配独立吗？

  • 两个数值变量独立吗？

  • 两个类别型变量独立吗？

使用包中提供的函数（见表12-2），我们可以很便捷地做置换检验，它们与第7章的大部分传统统计检验是等价的。 表12-2　相对于传统检验，提供可选置换检验的coin函数

检　　验 coin函数
两样本和K样本置换检验 oneway_test(y ~ A)
含一个分层（区组）因子的两样本和K样本置换检验 oneway_test(y ~ A | C)
Wilcoxon-Mann-Whitney秩和检验 wilcox_test(y ~ A)
Kruskal-Wallis检验 kruskal_test(y ~ A)
Person卡方检验 chisq_test(A ~ B)
Cochran-Mantel-Haenszel检验 cmh_test(A ~ B | C)
线性关联检验 lbl_test(D ~ E)
Spearman检验 spearman_test(y ~ x)
Friedman检验 friedman_test(y ~ A | C)
Wilcoxon符号秩检验 wilcoxsign_test(y1 ~ y2)

* 在coin函数中，y和x是数值变量，A和B是分类因子，C是类别型区组变量，D和E是有序因子，y1和y2是相匹配的数值变量。

表12-2列出来的每个函数都是如下形式：

  1. function_name( formula, data, distribution= )

其中：

  • formula描述的是要检验变量间的关系。示例可参见表12-2；

  • data是一个数据框；

  • distribution指定经验分布在零假设条件下的形式，可能值有exact，asymptotic和approximate。

若distribution = "exact"，那么在零假设条件下，分布的计算是精确的（即依据所有可能的排列组合）。当然，也可以根据它的渐进分布（distribution = "asymptotic"）或蒙特卡洛重抽样（distribution = "approxiamate(B = #)"）来做近似计算，其中#指所需重复的次数。distribution = "exact"当前仅可用于两样本问题。

注意　在coin包中，类别型变量和序数变量必须分别转化为因子和有序因子。另外，数据要以数据框形式存储。

在本节余下部分，我们将把表12-2中的一些置换检验应用到在先前章节出现的问题中，这样，你可以对传统的参数方法和非参数方法进行比较。本节最后，我们将通过一些高级拓展应用对coin包进行讨论。

12.2.1　独立两样本和K样本检验

首先，根据表12-2的虚拟数据，我们对独立样本t检验和单因素精确检验进行比较。结果见代码清单12-1。

代码清单12-1　虚拟数据中的t检验与单因素置换检验

  1. > library(coin)
  2. > score <- c(40, 57, 45, 55, 58, 57, 64, 55, 62, 65)
  3. > treatment <- factor(c(rep("A",5), rep("B",5)))
  4. > mydata <- data.frame(treatment, score)
  5. > t.test(score~treatment, data=mydata, var.equal=TRUE)
  7.         Two Sample t-test
  9. data:  score by treatment
  10. t = -2.3, df = 8, p-value = 0.04705
  11. alternative hypothesis: true difference in means is not equal to 0
  12. 95 percent confidence interval:
  13.  -19.04  -0.16
  14. sample estimates:
  15. mean in group A mean in group B
  16.              51              61
  18. > oneway_test(score~treatment, data=mydata, distribution="exact")
  20.         Exact 2-Sample Permutation Test
  22. data:  score by treatment (A, B)
  23. Z = -1.9, p-value = 0.07143
  24. alternative hypothesis: true mu is not equal to 0

传统t检验表明存在显著性差异（p < 0.05），而精确检验却表明差异并不显著（p > 0.072）。由于只有10个观测，我更倾向于相信置换检验的结果，在做出最后结论之前，还要多收集些数据。

现在来看Wilcoxon-Mann-Whitney U检验。第7章中，我们用wilcox.test()函数检验了美国南部监禁概率与非南部间的差异。现使用Wilcoxon秩和检验，可得：

  1. > library(MASS)
  2. > UScrime <- transform(UScrime, So = factor(So))
  3. > wilcox_test(Prob ~ So, data=UScrime, distribution="exact")
  5.         Exact Wilcoxon Mann-Whitney Rank Sum Test
  7. data:  Prob by So (0, 1)
  8. Z = -3.7, p-value = 8.488e-05
  9. alternative hypothesis: true mu is not equal to 0

结果表明监禁在南部可能更多。注意在上面的代码中，数值变量So被转化为因子，因为coin包规定所有的类别型变量都必须以因子形式编码。另外，聪明的读者可能会发现此处结果与第7章wilcox.test()计算结果一样，这是因为wilcox.test()默认计算的也是精确分布。

最后，探究下K样本检验问题。在第9章，对于50个患者的样本，我们使用了单因素方差分析来评价五种药物疗法对降低胆固醇的效果。下面代码对其做了近似的K样本置换检验：

  1. > library(multcomp)
  2. > set.seed(1234)
  3. > oneway_test(response~trt, data=cholesterol,
  4.   distribution=approximate(B=9999))
  6.         Approximative K-Sample Permutation Test
  8. data:  response by
  9.          trt (1time, 2times, 4times, drugD, drugE)
  10. maxT = 4.7623, p-value < 2.2e-16

此处，参考分布得自于数据9999次的置换。设定随机数种子可让结果重现。结果表明各组间病人的响应值显著不同。

12.2.2　列联表中的独立性

通过chisq_test()或cmh_test()函数，我们可用置换检验判断两类别型变量的独立性。当数据可根据第三个类别型变量进行分层时，需要使用后一个函数。若变量都是有序型，可使用lbl_test()函数来检验是否存在线性趋势。

第7章中，我们用卡方检验评价了关节炎的治疗（treatment）与效果（improvement）间的关系。治疗有两个水平（安慰剂、治疗），效果有三个水平（无、部分、显著），变量Improved以有序因子形式编码。

若想实施卡方检验的置换版本，可用如下代码：

  1. > library(coin)
  2. > library(vcd)
  3. > Arthritis <- transform(Arthritis,
  4.   Improved=as.factor(as.numeric(Improved)))
  5. > set.seed(1234)
  6. > chisq_test(Treatment~Improved, data=Arthritis,
  7.              distribution=approximate(B=9999))
  9.                 Approximative Pearson's Chi-Squared Test
  11. data:  Treatment by Improved (1, 2, 3)
  12. chi-squared = 13.055, p-value = 0.0018

此处经过9999次的置换，可获得一个近似的卡方检验。你可能会有疑问，为什么需要把变量Improved从一个有序因子变成一个分类因子？（好问题！）这是因为，如果你用有序因子，coin()将会生成一个线性与线性趋势检验，而不是卡方检验。虽然趋势检验在本例中是一个不错的选择，但是此处使用卡方检验可以同第7章所得的结果进行比较。

12.2.3　数值变量间的独立性

spearman_test()函数提供了两数值变量的独立性置换检验。第7章中，我们检验了美国文盲率与谋杀率间的相关性。如下代码可进行相关性的置换检验：

  1. > states <- as.data.frame(state.x77)
  2. > set.seed(1234)
  3. > spearman_test(Illiteracy~Murder, data=states,
  4.                 distribution=approximate(B=9999))
  6.         Approximative Spearman Correlation Test
  8. data:  Illiteracy by Murder
  9. Z = 4.7065, p-value < 2.2e-16
  10. alternative hypothesis: true mu is not equal to 0

基于9999次重复的近似置换检验可知：独立性假设并不被满足。注意，state.x77是一个矩阵，在coin包中，必须将其转化为一个数据框。

12.2.4　两样本和K样本相关性检验

当处于不同组的观测已经被分配得当，或者使用了重复测量时，样本相关检验便可派上用场。对于两配对组的置换检验，可使用wilcoxsign_test()函数；多于两组时，使用friedman_test()函数。

第7章中，我们比较了城市男性中14～24年龄段（U1）与35～39年龄段（U2）间的失业率差异。由于两个变量对于美国50个州都有记录，你便有了一个两依赖组设计（state是匹配变量），可使用精确Wilcoxon符号秩检验来判断两个年龄段间的失业率是否相等：

  1. > library(coin)
  2. > library(MASS)
  3. > wilcoxsign_test(U1~U2, data=UScrime, distribution="exact")
  5.         Exact Wilcoxon-Signed-Rank Test
  7. data:  y by x (neg, pos)
  8.          stratified by block
  9. Z = 5.9691, p-value = 1.421e-14
  10. alternative hypothesis: true mu is not equal to 0

结果表明两者的失业率是不同的。

12.2.5　深入探究

coin包提供了一个置换检验的一般性框架，可以分析一组变量相对于其他任意变量，是否与第二组变量（可根据一个区组变量分层）相互独立。特别地，independence_test()函数可以让用户从置换角度来思考大部分传统检验，进而在面对无法用传统方法解决的问题时，使用户可以自己构建新的统计检验。当然，这种灵活性也是有门槛的：要正确使用该函数必须具备丰富的统计知识。更多函数细节请参阅包附带的文档（运行vignette("coin")即可）。

下一节，你将学习lmPerm包，它提供了线性模型的置换方法，包括回归和方差分析。