一个相当普遍的功能

我们将使用内置的牙齿生长数据集 。我们感兴趣的是，当给予豚鼠维生素 C 与橙汁时，牙齿生长是否存在统计学上的显着差异。

这是完整的例子：

teethVC = ToothGrowth[ToothGrowth$supp == 'VC',]
teethOJ = ToothGrowth[ToothGrowth$supp == 'OJ',]

permutationTest = function(vectorA, vectorB, testStat){
  N = 10^5
  fullSet = c(vectorA, vectorB)
  lengthA = length(vectorA)
  lengthB = length(vectorB)
  trials <- replicate(N, 
                      {index <- sample(lengthB + lengthA, size = lengthA, replace = FALSE)
                      testStat((fullSet[index]), fullSet[-index])  } )
  trials
}
vec1 =teethVC$len;
vec2 =teethOJ$len;
subtractMeans = function(a, b){ return (mean(a) - mean(b))}
result = permutationTest(vec1, vec2, subtractMeans)
observedMeanDifference = subtractMeans(vec1, vec2)
result = c(result, observedMeanDifference)
hist(result)
abline(v=observedMeanDifference, col = "blue")
pValue = 2*mean(result <= (observedMeanDifference))
pValue

在我们读取 CSV 之后，我们定义了该函数

permutationTest = function(vectorA, vectorB, testStat){
  N = 10^5
  fullSet = c(vectorA, vectorB)
  lengthA = length(vectorA)
  lengthB = length(vectorB)
  trials <- replicate(N, 
                      {index <- sample(lengthB + lengthA, size = lengthA, replace = FALSE)
                      testStat((fullSet[index]), fullSet[-index])  } )
  trials
}

该函数接受两个向量，并将它们的内容混合在一起，然后在混洗向量上执行函数 testStat。teststat 的结果被添加到 trials，这是返回值。

这样做了 4 次。请注意，值 N 很可能是函数的参数。

这给我们留下了一组新的数据，trials，如果两个变量之间确实没有关系，可能会产生一系列均值。

现在定义我们的测试统计：

subtractMeans = function(a, b){ return (mean(a) - mean(b))}

执行测试：

result = permutationTest(vec1, vec2, subtractMeans)

计算我们实际观察到的平均差异：

observedMeanDifference = subtractMeans(vec1, vec2)

让我们看看我们的测试统计数据的直方图上的观察结果。

hist(result)
abline(v=observedMeanDifference, col = "blue")

它看起来不像我们观察到的结果很可能是随机发生的……

我们想要计算 p 值，即原始观察结果的可能性，如果它们在两个变量之间没有关系的话。

pValue = 2*mean(result >= (observedMeanDifference))

让我们稍微打破一下：

result >= (observedMeanDifference)

将创建一个布尔向量，如：

FALSE TRUE FALSE FALSE TRUE FALSE ...

每当 result 的值大于或等于 observedMean 时，使用 TRUE。

mean 函数将此向量解释为 TRUE 的 1 和 FALSE 的 0，并给出混合中 1 的百分比，即我们的混洗向量平均差异超过或等于我们观察到的次数。

最后，我们乘以 2，因为我们的测试统计量的分布是高度对称的，我们真的想知道哪些结果比我们观察到的结果更极端。

剩下的就是输出 p 值，结果是 0.06093939。对这个值的解释是主观的，但我会说维生素 C 看起来比橙汁更能促进牙齿生长。