R Language
Utför ett permutationstest

En ganska allmän funktion

Vi kommer att använda den inbyggda datatillväxten . Vi är intresserade av om det finns en statistiskt signifikant skillnad i tandtillväxt när marsvinen ges C-vitamin kontra apelsinjuice.

Här är det fullständiga exemplet:

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

När vi läst i CSV definierar vi funktionen

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
}

Den här funktionen tar två vektorer, och blandar deras innehåll tillsammans, då utför funktionen testStat på blandade vektorer. Resultatet av teststat läggs till i trials , vilket är returvärdet.

Det gör detta N = 10^5 gånger. Observera att värdet N mycket väl kunde ha varit en parameter för funktionen.

Detta ger oss en ny uppsättning data, trials , den uppsättning medel som kan uppstå om det verkligen inte finns något samband mellan de två variablerna.

Nu för att definiera vår teststatistik:

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

Utför testet:

result = permutationTest(vec1, vec2, subtractMeans)

Beräkna vår faktiska observerade genomsnittliga skillnad:

observedMeanDifference = subtractMeans(vec1, vec2)

Låt oss se hur vår observation ser ut på ett histogram i vår teststatistik.

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

Det ser inte ut som vårt observerade resultat mycket troligtvis kommer att inträffa av slumpmässig chans ...

Vi vill beräkna p-värdet, sannolikheten för det ursprungliga observerade resultatet om det inte är något samband mellan de två variablerna.

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

Låt oss dela upp det lite:

result >= (observedMeanDifference)

Skapar en boolesk vektor, som:

FALSE TRUE FALSE FALSE TRUE FALSE ...

Med TRUE varje gång värdet på result är större än eller lika med den observedMean .

Funktionen mean kommer att tolka denna vektor som 1 för TRUE och 0 för FALSE , och ge oss andelen 1 är i mixen, det vill säga hur många gånger vår SHUFFLE vektor genomsnittliga skillnaden överträffade eller motsvarade vad vi observerat.

Slutligen multiplicerar vi med 2 eftersom fördelningen av vår teststatistik är mycket symmetrisk, och vi vill verkligen veta vilka resultat som är "mer extrema" än vårt observerade resultat.

Allt som återstår är att mata ut p-värdet, vilket visar sig vara 0.06093939 . Tolkning av detta värde är subjektivt, men jag skulle säga att det ser ut som C-vitamin främjar tandtillväxt ganska mycket mer än apelsinjuice gör.