R Language
Effectuer un test de permutation

Une fonction assez générale

Nous utiliserons le jeu de données de croissance des dents intégré. Nous nous demandons s'il existe une différence statistiquement significative dans la croissance des dents lorsque les cobayes reçoivent de la vitamine C par rapport au jus d'orange.

Voici l'exemple complet:

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

Après avoir lu dans le fichier CSV, nous définissons la fonction

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
}

Cette fonction prend deux vecteurs et mélange leurs contenus, puis effectue la fonction testStat sur les vecteurs mélangés. Le résultat de teststat est ajouté à trials , qui correspond à la valeur teststat .

Cela fait N = 10^5 fois. Notez que la valeur N pourrait très bien avoir été un paramètre pour la fonction.

Cela nous laisse un nouvel ensemble de données, d’ trials , l’ensemble des moyens qui pourraient résulter de l’absence de relation entre les deux variables.

Maintenant, pour définir notre statistique de test:

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

Effectuez le test:

result = permutationTest(vec1, vec2, subtractMeans)

Calculez notre différence moyenne observée réelle:

observedMeanDifference = subtractMeans(vec1, vec2)

Voyons à quoi ressemble notre observation sur un histogramme de notre statistique de test.

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

Il ne ressemble pas à notre résultat observé est très susceptible de se produire par hasard ...

Nous voulons calculer la valeur de p, la vraisemblance du résultat original observé s’il n’ya pas de relation entre les deux variables.

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

Brisons un peu ça:

result >= (observedMeanDifference)

Va créer un vecteur booléen, comme:

FALSE TRUE FALSE FALSE TRUE FALSE ...

Avec TRUE chaque fois que la valeur du result est supérieure ou égale à la valeur observedMean .

La mean fonctions interprétera ce vecteur comme étant 1 pour TRUE et 0 pour FALSE , et nous donnera le pourcentage de 1 dans le mélange, c’est-à-dire le nombre de fois que notre différence moyenne

Enfin, nous multiplions par 2 car la distribution de notre statistique de test est très symétrique et nous voulons vraiment savoir quels résultats sont "plus extrêmes" que notre résultat observé.

Il ne reste plus qu'à afficher la valeur p, qui s'avère être 0.06093939 . L'interprétation de cette valeur est subjective, mais je dirais que la vitamine C favorise beaucoup plus la croissance des dents que le jus d'orange.