R一个相当普通的功能
示例
我们将使用内置的牙齿生长数据集。我们对豚鼠服用维生素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，即返回值。
这N=10^5一次。注意，该值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 ...
随着TRUE的值每次result大于或等于observedMean。
函数mean将按1forTRUE和0for解释此向量FALSE，并提供1混合中的的百分比，即我们混洗后的向量均值差超过或等于我们观察到的次数。
最后，我们乘以2是因为测试统计量的分布是高度对称的，我们真的想知道哪些结果比我们观察到的结果“更极端”。
剩下的就是输出p值，结果为0.06093939。该值的解释是主观的，但我想说维生素C促进牙齿生长的作用比橙汁大得多。