Comment créer un intervalle de confiance pour le paramètre d'un test de permutation?

9

Les tests de permutation sont des tests de signification basés sur des rééchantillons de permutation tirés au hasard à partir des données originales. Les rééchantillons de permutation sont dessinés sans remplacement, contrairement aux échantillons bootstrap, qui sont dessinés avec remplacement. Voici un exemple que j'ai fait dans R d'un test de permutation simple. (Vos commentaires sont les bienvenus)

Les tests de permutation présentent de grands avantages. Ils ne nécessitent pas de formes de population spécifiques telles que la normalité. Ils s'appliquent à une variété de statistiques, pas seulement aux statistiques qui ont une distribution simple sous l'hypothèse nulle. Ils peuvent donner des valeurs de p très précises, quelles que soient la forme et la taille de la population (si suffisamment de permutations sont utilisées).

J'ai également lu qu'il est souvent utile de donner un intervalle de confiance avec un test, qui est créé en utilisant le rééchantillonnage bootstrap plutôt que le rééchantillonnage par permutation.

Pourriez-vous expliquer (ou simplement donner le code R) comment un intervalle de confiance est construit (c'est-à-dire pour la différence entre les moyennes des deux échantillons dans l'exemple ci-dessus)?

ÉDITER

Après quelques recherches sur Google, j'ai trouvé cette lecture intéressante .

George Dontas
la source

Réponses:

7

Vous pouvez utiliser le rééchantillonnage par permutation. Cela dépend vraiment d'un certain nombre de facteurs. Si vos permutations sont un nombre relativement faible, votre estimation de votre intervalle de confiance n'est pas si grande avec les permutations. Vos permutations sont dans une sorte de zone grise et sont probablement très bien.

La seule différence par rapport à votre code précédent est que vous généreriez vos échantillons au hasard plutôt qu'avec des permutations. Et vous en généreriez plus, disons 1000 par exemple. Obtenez les scores de différence pour vos 1000 répétitions de votre expérience. Prenez les seuils pour le milieu 950 (95%). Voilà votre intervalle de confiance. Il tombe directement du bootstrap.

Vous avez déjà fait la plupart de cela dans votre exemple. dif.treat a une longueur de 462 articles. Par conséquent, vous avez besoin des seuils inférieurs de 2,5% et supérieurs de 2,5% (environ 11 éléments à chaque extrémité).

Utiliser votre code d'avant ...

y <- sort(dif.treat)
ci.lo <- y[11]
ci.hi <- y[462-11]

D'un côté, je dirais que 462 est un peu faible, mais vous trouverez un bootstrap à 10000 qui sort avec des scores peu différents (probablement plus proches de la moyenne).

Je pensais aussi ajouter du code simple nécessitant la bibliothèque de démarrage (basé sur votre code précédent).

diff <- function(x,i) mean(x[i[6:11]]) - mean(x[i[1:5]])
b <- boot(total, diff, R = 1000)
boot.ci(b)
John
la source
Je vous remercie. Est-il correct de générer les échantillons en utilisant sampleet replace=TRUE? Y a-t-il une raison d'utiliser un package comme boot?
George Dontas
En règle générale, il est effectué avec remplacement, vous devez donc définir ce paramètre sur TRUE. Quant à savoir pourquoi ... le paquet est optimisé, il devrait donc s'exécuter plus rapidement ... jamais chronométré. Cela peut être un problème si vous définissez R à grande. Et, comme vous pouvez le voir, le code est agréable et concis. Il possède également de nombreuses fonctionnalités que vous n'auriez pas facilement accès aux vôtres.
John
boot.ci renvoie l'intervalle de confiance. Y a-t-il une fonction (boot) qui donne la valeur p. (comme le rapport du nombre de différences au moins aussi élevé que celui observé, sur le nombre total d'échantillons générés)
George Dontas
ok, j'ai trouvé un moyen de l'obtenir:sum(b$t>=b$t0)/b$R
George Dontas
@ gd047: tenez compte du fait qu'il s'agit d'une valeur p unilatérale que vous calculez.
Joris Meys
4

Comme un test de permutation est un test exact , vous donnant une valeur de p exacte. L'amorçage d'un test de permutation n'a pas de sens.

À côté de cela, déterminer un intervalle de confiance autour d'une statistique de test n'a pas non plus de sens, car il est calculé en fonction de votre échantillon et non d'une estimation. Vous déterminez des intervalles de confiance autour d'estimations comme les moyennes et les goûts, mais pas autour des statistiques de test.

Les tests de permutation ne doivent pas être utilisés sur des ensembles de données si gros que vous ne pouvez plus calculer toutes les permutations possibles. Si tel est le cas, utilisez une procédure d'amorçage pour déterminer la coupure de la statistique de test que vous utilisez. Mais encore une fois, cela n'a pas grand-chose à voir avec un intervalle de confiance à 95%.

Un exemple: j'utilise ici la statistique T classique, mais j'utilise une approche simple du bootstrap pour le calcul de la distribution empirique de ma statistique. Sur cette base, je calcule une valeur de p empirique:

x <- c(11.4,25.3,29.9,16.5,21.1)
y <- c(23.7,26.6,28.5,14.2,17.9,24.3)

t.sample <- t.test(x,y)$statistic
t.dist <- apply(
      replicate(1000,sample(c(x,y),11,replace=F)),2,
      function(i){t.test(i[1:5],i[6:11])$statistic})

# two sided testing
center <- mean(t.dist)
t.sample <-abs(t.sample-center)
t.dist <- abs(t.dist - center)
p.value <- sum( t.sample < t.dist ) / length(t.dist)
p.value

Tenez compte du fait que ce test bilatéral ne fonctionne que pour les distributions symétriques. Les distributions non symétriques ne sont généralement testées que sur une seule face.

ÉDITER :

OK, j'ai mal compris la question. Si vous souhaitez calculer un intervalle de confiance sur l'estimation de la différence, vous pouvez utiliser le code mentionné ici pour l'amorçage au sein de chaque échantillon. Attention, il s'agit d'une estimation biaisée: cela donne généralement un IC trop petit. Voir également l'exemple donné comme une raison pour laquelle vous devez utiliser une approche différente pour l'intervalle de confiance et la valeur de p.

Joris Meys
la source
1
Pouvez-vous expliquer pourquoi les tests de permutation ne devraient pas être utilisés sur des ensembles de données pour lesquels vous ne pouvez pas calculer toutes les permutations possibles?
Andy W
@Andy W: Définissez d'abord "test de permutation". pour moi, les tests de permutation sont des tests exacts, utilisant toutes les permutations possibles. C'est impossible sur de plus grands ensembles de données. Les "tests de permutation approximatifs" sont en fait la méthode simple de Monte Carlo, et devraient être traités de cette façon. À côté de cela, le théorème central limite garantit dans la plupart des cas que les hypothèses concernant la distribution des statistiques de test sont respectées lors de l'utilisation de grands ensembles de données. Dans les tests complexes, l'utilisation de tests de permutation sur de grands ensembles de données rend les temps de calcul insupportablement longs sans ajouter de valeur significative. my2cents
Joris Meys
Je n'ai rien dit comme amorcer un test de permutation. Je suis entré dans cette question après avoir lu le dernier paragraphe de la [SECTION 14.5 | Résumé], dans le pdf lié.
George Dontas
@ gd047 Alors j'ai mal lu votre question. Mais vous devez vraiment garder les intervalles de confiance et les valeurs p strictement séparés. L'intervalle de confiance est estimé sur la base du bootstrap au sein de chaque échantillon (bien qu'il soit biaisé par définition), le test de permutation est effectué par permutations sur l'ensemble de données complet. Ce sont deux choses complètement différentes.
Joris Meys
@Kevin: Le code était sacrément juste. Relisez le code: le x[6:11]fait référence à l'argument xde la fonction anonyme dans l'appliquer. Peut-être déroutant, mais votre montage a donné de très mauvais résultats. Veuillez commenter ce que vous pensez qu'il devrait être avant de modifier le code. Me sauve un retour en arrière. Pour éviter toute confusion, j'ai changé cela xeni
Joris Meys
0

Du code Joris Meys dans les réponses mais avec des modifications pour lui permettre d'être appliqué dans plus d'une seule situation:

J'ai essayé de modifier l'autre mais je n'ai pas eu le temps de terminer et pour une raison quelconque, je ne peux pas commenter (peut-être parce que c'est une vieille question).

x <- c(11.4,25.3,29.9,16.5,21.1)
y <- c(23.7,26.6,28.5,14.2,17.9,24.3)

t.sample <- t.test(x,y)$statistic

t.dist <- apply(
          replicate(1000,sample(c(x,y),length(c(x,y)),replace=F)), 2,
          function(i){t.test(i[1:length(x)],i[length(x)+1:length(c(x,y))])$statistic})

# two sided testing
center <- mean(t.dist)
t.sample <-abs(t.sample-center)
t.dist <- abs(t.dist - center)
p.value <- sum( t.sample < t.dist ) / length(t.dist)
p.value
Kevin
la source