R Language
Durchführen eines Permutationstests

Eine ziemlich allgemeine Funktion

Wir werden den integrierten Zahnwachstumsdatensatz verwenden . Wir sind daran interessiert, ob es einen statistisch signifikanten Unterschied im Zahnwachstum gibt, wenn den Meerschweinchen Vitamin C gegen Orangensaft gegeben wird.

Hier ist das vollständige Beispiel:

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

Nachdem wir die CSV eingelesen haben, definieren wir die Funktion

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
}

Diese Funktion nimmt zwei Vektoren, mischt ihren Inhalt zusammen und führt dann die Funktion testStat für die gemischten Vektoren aus. Das Ergebnis von teststat wird zu den trials hinzugefügt. Dies ist der Rückgabewert.

Das macht N = 10^5 mal. Beachten Sie, dass der Wert N durchaus ein Parameter für die Funktion gewesen sein könnte.

Dadurch bleiben neue Daten, trials , die Mittel, die sich ergeben könnten, wenn wirklich keine Beziehung zwischen den beiden Variablen besteht.

Nun definieren Sie unsere Teststatistik:

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

Führen Sie den Test durch:

result = permutationTest(vec1, vec2, subtractMeans)

Berechnen Sie unsere tatsächlich beobachtete mittlere Differenz:

observedMeanDifference = subtractMeans(vec1, vec2)

Mal sehen, wie unsere Beobachtung in einem Histogramm unserer Teststatistik aussieht.

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

Es sieht nicht so aus, als ob unser beobachtetes Ergebnis sehr wahrscheinlich zufällig auftritt ...

Wir wollen den p-Wert berechnen, die Wahrscheinlichkeit des ursprünglichen beobachteten Ergebnisses, wenn zwischen den beiden Variablen keine Beziehung besteht.

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

Lassen Sie uns das ein bisschen brechen:

result >= (observedMeanDifference)

Erzeugt einen booleschen Vektor wie:

FALSE TRUE FALSE FALSE TRUE FALSE ...

Bei TRUE ist der result jedes Mal größer oder gleich dem observedMean .

Die Funktion mean wird dieses Vektors als interpretieren 1 für TRUE und 0 für FALSE , und geben Sie uns den Prozentsatz von 1 ‚s in der Mischung, dh die Anzahl der Male , unsere neu gemischt Vektor mittlere Differenz übertroffen oder erreicht , was wir beobachtet.

Schließlich multiplizieren wir mit 2, da die Verteilung unserer Teststatistik hochsymmetrisch ist und wir wirklich wissen wollen, welche Ergebnisse "extremer" sind als unser beobachtetes Ergebnis.

Alles, was übrig bleibt, ist die Ausgabe des p-Werts, der 0.06093939 . Die Interpretation dieses Wertes ist subjektiv, aber ich würde sagen, dass Vitamin C das Zahnwachstum wesentlich mehr fördert als Orangensaft.