一個相當普遍的功能

我們將使用內建的牙齒生長資料集 。我們感興趣的是，當給予豚鼠維生素 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 看起來比橙汁更能促進牙齒生長。