R Language
Выполнение теста перестановки

Достаточно общая функция

Мы будем использовать встроенный набор данных о росте зуба . Нас интересует, существует ли статистически значимая разница в росте зубов, когда морским свинкам дают витамин С против апельсинового сока.

Вот полный пример:

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 функции будет интерпретировать этот вектор как 1 для TRUE и 0 для FALSE , и дать нам процент 1 в смеси, то есть количество раз, когда наша перетасованная векторная средняя разница превысила или сравнила то, что мы наблюдали.

Наконец, умножим на 2, потому что распределение нашей тестовой статистики очень симметрично, и мы действительно хотим знать, какие результаты «более экстремальные», чем наш наблюдаемый результат.

0.06093939 только вывести значение p, которое оказалось равным 0.06093939 . Интерпретация этого значения субъективна, но я бы сказал, что это похоже на то, что витамин С способствует росту зубов намного больше, чем апельсиновый сок.