R Language
Een permutatietest uitvoeren

Een vrij algemene functie

We zullen de ingebouwde gegevensset voor tandgroei gebruiken . We zijn geïnteresseerd in de vraag of er een statistisch significant verschil is in tandgroei wanneer de cavia's vitamine C versus sinaasappelsap krijgen.

Hier is het volledige voorbeeld:

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

Nadat we de CSV hebben gelezen, definiëren we de functie

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
}

Deze functie neemt twee vectoren en schudt hun inhoud samen en voert vervolgens de functie testStat op de testStat vectoren. Het resultaat van teststat wordt toegevoegd aan trials , wat de retourwaarde is.

Het doet deze N = 10^5 keer. Merk op dat de waarde N heel goed een parameter voor de functie had kunnen zijn.

Dit laat ons met een nieuwe set gegevens, trials , de set middelen die zou kunnen resulteren als er echt geen verband is tussen de twee variabelen.

Om nu onze teststatistiek te definiëren:

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

Voer de test uit:

result = permutationTest(vec1, vec2, subtractMeans)

Bereken ons feitelijk waargenomen gemiddelde verschil:

observedMeanDifference = subtractMeans(vec1, vec2)

Laten we eens kijken hoe onze observatie eruit ziet op een histogram van onze teststatistiek.

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

Het ziet er niet uit als onze waargenomen resultaat is zeer waarschijnlijk ontstaan door toeval ...

We willen de p-waarde berekenen, de waarschijnlijkheid van het oorspronkelijke waargenomen resultaat als er geen verband is tussen de twee variabelen.

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

Laten we dat een beetje opsplitsen:

result >= (observedMeanDifference)

Zal een booleaanse vector maken, zoals:

FALSE TRUE FALSE FALSE TRUE FALSE ...

Met TRUE elke keer de waarde van het result groter dan of gelijk aan de observedMean .

De functie mean zal deze vector interpreteren als 1 voor TRUE en 0 voor FALSE , en geef ons het percentage van 1 's in de mix, dat wil zeggen het aantal keren dat onze geschudde vector gemiddelde verschil overtroffen of geëvenaard wat we waargenomen.

Ten slotte vermenigvuldigen we met 2 omdat de verdeling van onze teststatistiek zeer symmetrisch is en we echt willen weten welke resultaten "extremer" zijn dan ons waargenomen resultaat.

Het enige dat overblijft is om de p-waarde uit te voeren, die 0.06093939 blijkt te zijn. De interpretatie van deze waarde is subjectief, maar ik zou zeggen dat het lijkt alsof vitamine C de tandgroei aanzienlijk meer bevordert dan sinaasappelsap.