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促进牙齿生长的作用比橙汁大得多。