R Language
Przeprowadzanie testu permutacji

Dość ogólna funkcja

Użyjemy wbudowanego zestawu danych dotyczących wzrostu zębów . Interesuje nas, czy istnieje statystycznie istotna różnica we wzroście zębów, gdy świnkom morskim podaje się witaminę C w porównaniu z sokiem pomarańczowym.

Oto pełny przykład:

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

Po przeczytaniu w CSV definiujemy funkcję

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
}

Ta funkcja pobiera dwa wektory i tasuje ich zawartość, a następnie wykonuje testStat na testStat wektorach. Wynik teststat jest dodawany do trials , co jest wartością zwracaną.

Robi to N = 10^5 razy. Zauważ, że wartość N mogłaby równie dobrze być parametrem funkcji.

Pozostaje nam nowy zestaw danych, trials , zestaw środków, które mogą powstać, jeśli naprawdę nie ma związku między tymi dwiema zmiennymi.

Teraz, aby zdefiniować naszą statystykę testową:

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

Wykonaj test:

result = permutationTest(vec1, vec2, subtractMeans)

Oblicz naszą rzeczywistą zaobserwowaną średnią różnicę:

observedMeanDifference = subtractMeans(vec1, vec2)

Zobaczmy, jak wygląda nasza obserwacja na histogramie naszej statystyki testowej.

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

Nie wygląda na to, aby nasz zaobserwowany wynik pojawił się przypadkowo ...

Chcemy obliczyć wartość p, prawdopodobieństwo oryginalnego zaobserwowanego wyniku, jeśli nie ma związku między dwiema zmiennymi.

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

Rozbijmy to trochę:

result >= (observedMeanDifference)

Stworzy wektor boolowski, taki jak:

FALSE TRUE FALSE FALSE TRUE FALSE ...

Przy pomocy TRUE każdym razem wartość result jest większa lub równa observedMean wartości observedMean .

mean funkcji zinterpretuje ten wektor jako 1 dla TRUE i 0 dla FALSE , i da nam procent 1 w mieszance, tj. Liczbę przypadków, gdy nasza tasowana średnia różnica wektora przekroczyła lub wyrównała to, co zaobserwowaliśmy.

Na koniec mnożymy przez 2, ponieważ rozkład naszej statystyki testowej jest wysoce symetryczny i naprawdę chcemy wiedzieć, które wyniki są „bardziej ekstremalne” niż nasz zaobserwowany wynik.

Pozostało tylko wyprowadzenie wartości p, która okazuje się wynosić 0.06093939 . Interpretacja tej wartości jest subiektywna, ale powiedziałbym, że wygląda na to, że witamina C promuje wzrost zębów znacznie bardziej niż sok pomarańczowy.