Quelle implémentation de test de permutation dans R utiliser au lieu de tests t (appariés et non appariés)?

56

J'ai des données provenant d'une expérience que j'ai analysée à l'aide de tests t. La variable dépendante est mise à l'échelle par intervalles et les données sont soit non appariées (c'est-à-dire 2 groupes), soit appariées (c'est-à-dire intra-sujets). Par exemple (au sein des sujets):

x1 <- c(99, 99.5, 65, 100, 99, 99.5, 99, 99.5, 99.5, 57, 100, 99.5, 
        99.5, 99, 99, 99.5, 89.5, 99.5, 100, 99.5)
y1 <- c(99, 99.5, 99.5, 0, 50, 100, 99.5, 99.5, 0, 99.5, 99.5, 90, 
        80, 0, 99, 0, 74.5, 0, 100, 49.5)

Cependant, les données n'étant pas normales, un critique nous a demandé d'utiliser autre chose que le test t. Cependant, comme on peut facilement le constater, les données ne sont pas non plus normalement distribuées, mais les distributions ne sont pas égales entre les conditions: texte alternatif

Par conséquent, les tests non paramétriques habituels, le test de Mann-Whitney-U (non apparié) et le test de Wilcoxon (appariés), ne peuvent pas être utilisés car ils nécessitent des distributions égales entre les conditions. J'ai donc décidé qu'un test de rééchantillonnage ou de permutation serait préférable.

Maintenant, je recherche une implémentation R d'un équivalent du test t basé sur la permutation, ou tout autre conseil sur l'utilisation des données.

Je sais qu'il existe des packages R qui peuvent le faire pour moi (par exemple, une pièce de monnaie, une permanente, un test exact de test de pile, etc.), mais je ne sais pas lequel choisir. Donc, si quelqu'un avec une certaine expérience dans l'utilisation de ces tests pouvait me donner un élan, ce serait ubercool.

MISE À JOUR: L’ idéal serait que vous donniez un exemple de la manière de rapporter les résultats de ce test.

Henrik
la source

Réponses:

43

Cela ne devrait pas avoir beaucoup d'importance, car la statistique de test sera toujours la différence de moyen (ou quelque chose d'équivalent). La mise en œuvre des méthodes de Monte-Carlo peut présenter de petites différences. Essayer les trois packages avec vos données avec un test unilatéral pour deux variables indépendantes:

DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
library(coin)                    # for oneway_test(), pvalue()
pvalue(oneway_test(DV ~ IV, alternative="greater", 
                   distribution=approximate(B=9999)))
[1] 0.00330033

library(perm)                    # for permTS()
permTS(DV ~ IV, alternative="greater", method="exact.mc", 
       control=permControl(nmc=10^4-1))$p.value
[1] 0.003

library(exactRankTests)          # for perm.test()
perm.test(DV ~ IV, paired=FALSE, alternative="greater", exact=TRUE)$p.value
[1] 0.003171822

Pour vérifier la valeur p exacte avec un calcul manuel de toutes les permutations, je limiterai les données aux 9 premières valeurs.

x1 <- x1[1:9]
y1 <- y1[1:9]
DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
pvalue(oneway_test(DV ~ IV, alternative="greater", distribution="exact"))
[1] 0.0945907

permTS(DV ~ IV, alternative="greater", exact=TRUE)$p.value
[1] 0.0945907

# perm.test() gives different result due to rounding of input values
perm.test(DV ~ IV, paired=FALSE, alternative="greater", exact=TRUE)$p.value
[1] 0.1029412

# manual exact permutation test
idx  <- seq(along=DV)                 # indices to permute
idxA <- combn(idx, length(x1))        # all possibilities for different groups

# function to calculate difference in group means given index vector for group A
getDiffM <- function(x) { mean(DV[x]) - mean(DV[!(idx %in% x)]) }
resDM    <- apply(idxA, 2, getDiffM)  # difference in means for all permutations
diffM    <- mean(x1) - mean(y1)       # empirical differencen in group means

# p-value: proportion of group means at least as extreme as observed one
(pVal <- sum(resDM >= diffM) / length(resDM))
[1] 0.0945907

coinet exactRankTestssont tous deux du même auteur, mais coinsemble être plus général et plus complet - également en termes de documentation. exactRankTestsn'est plus développé activement. Je choisirais donc coin(notamment à cause de fonctions informatives support()), à moins que vous n'aimiez pas traiter avec des objets S4.

EDIT: pour deux variables dépendantes, la syntaxe est la suivante:

id <- factor(rep(1:length(x1), 2))    # factor for participant
pvalue(oneway_test(DV ~ IV | id, alternative="greater",
                   distribution=approximate(B=9999)))
[1] 0.00810081
caracal
la source
Merci pour votre bonne réponse! Deux autres questions: Votre deuxième exemple signifie-t-il que cette pièce offre en fait toutes les permutations possibles et constitue un test exact? Y at-il un avantage à ne pas fournir un test exact dans mon cas?
Henrik
10
(+1) Il n’est pas surprenant que le test t (non apparié) donne essentiellement la même valeur p, 0,000349. Malgré les propos de l'examinateur, le test t est applicable à ces données. La raison en est que les distributions d'échantillonnage des moyennes sont à peu près normales, même si les distributions des données ne le sont pas. De plus, comme vous pouvez le constater, le test t est en réalité plus conservateur que le test de permutation. (Cela signifie qu'un résultat significatif avec le test t indique que le test de permutation sera probablement aussi significatif.)
whuber
2
@Henrik Pour certaines situations (test choisi et complexité numérique), coinon peut en effet calculer exactement la distribution de permutation (sans passer par toutes les permutations, il existe des algorithmes plus élégants que cela). Compte tenu de ce choix, la distribution exacte semble préférable, mais la différence par rapport à une approximation de Monte-Carlo avec un nombre élevé de répétitions devrait être faible.
Caracal
1
@Caracal Merci pour la clarification. Une question demeure: les données que j'ai présentées sont appariées. Par conséquent, j'ai besoin de l'équivalent du test t apparié. La oneway_testfonction est -elle précise? Et si oui, lequel est le correct pour les données non appariées?
Henrik
2
@Henrik L' coinauteur m'a écrit que je oneway_test()ne peux pas calculer la distribution exacte pour le cas dépendant, vous devez utiliser l'approximation MC (cela ne wilcoxsign_test()convient que pour le test exact). Je ne le savais pas et préférerais une erreur dans ce cas, mais MC devrait être suffisamment précis avec un nombre élevé de répétitions.
caracal
29

Je pense que quelques commentaires sont recevables.

1) Je vous encourage à essayer plusieurs affichages visuels de vos données, car ils peuvent capturer des éléments perdus par des histogrammes (graphiques, par exemple), et je vous recommande également fortement de tracer sur des axes côte à côte. Dans ce cas, je ne pense pas que les histogrammes communiquent très bien les principales caractéristiques de vos données. Par exemple, examinons les boîtes à moustaches côte à côte:

boxplot(x1, y1, names = c("x1", "y1"))

texte alternatif

Ou même des graphiques à bandes côte à côte:

stripchart(c(x1,y1) ~ rep(1:2, each = 20), method = "jitter", group.names = c("x1","y1"), xlab = "")

texte alternatif

Regardez les centres, les spreads et les formes de ceux-ci! Environ les trois quarts des données situent bien au-dessus de la médiane des données . La propagation de est minime, alors que celle de est énorme. Les deux et sont fortement biaisé à gauche, mais de différentes manières. Par exemple, a cinq valeurs (!) Répétées de zéro.y 1 x 1 y 1 x 1 y 1 y 1x1y1x1y1x1y1y1

2) Vous n'avez pas expliqué en détail l'origine de vos données, ni comment elles ont été mesurées, mais ces informations sont très importantes pour le choix d'une procédure statistique. Vos deux échantillons ci-dessus sont-ils indépendants? Y a-t-il des raisons de croire que les distributions marginales des deux échantillons devraient être identiques (à l'exception d'une différence de localisation, par exemple)? Quelles étaient les considérations préalables à l'étude qui vous ont amené à rechercher des preuves d'une différence entre les deux groupes?

3) Le test t n'est pas approprié pour ces données car les distributions marginales sont nettement non normales, avec des valeurs extrêmes dans les deux échantillons. Si vous le souhaitez, vous pouvez faire appel au CLT (en raison de votre échantillon de taille moyenne) pour utiliser un test (qui serait similaire à un test z pour les grands échantillons), mais étant donné l'asymétrie (dans les deux variables) de vos données, je ne jugerais pas un tel appel très convaincant. Bien sûr, vous pouvez quand même l'utiliser pour calculer une valeur de , mais qu'est-ce que cela fait pour vous? Si les hypothèses ne sont pas satisfaites, alors une n'est qu'une statistique; il ne dit pas ce que vous voulez (probablement) vouloir savoir: s'il existe une preuve que les deux échantillons proviennent de distributions différentes.p pzpp

4) Un test de permutation serait également inapproprié pour ces données. L'hypothèse unique et souvent négligée pour les tests de permutation est que les deux échantillons sont échangeables sous l'hypothèse nulle. Cela voudrait dire qu'ils ont des distributions marginales identiques (sous le zéro). Mais vous êtes en difficulté, car les graphiques suggèrent que les distributions diffèrent tant par leur emplacement que par leur échelle (et leur forme également). Ainsi, vous ne pouvez pas (valablement) tester une différence d'emplacement car les échelles sont différentes et vous ne pouvez pas (valablement) rechercher une différence d'échelle car les emplacements sont différents. Oops. Encore une fois, vous pouvez quand même faire le test et obtenir une valeur de , mais alors quoi? Qu'as-tu vraiment accompli?p

5) À mon avis, ces données illustrent parfaitement (?) Le fait qu'une image bien choisie vaut 1 000 tests d'hypothèses. Nous n'avons pas besoin de statistiques pour faire la différence entre un crayon et une grange. À mon avis, l'énoncé approprié pour ces données serait "Ces données présentent des différences marquées en ce qui concerne l'emplacement, l'échelle et la forme". Vous pouvez ensuite utiliser des statistiques descriptives (robustes) pour chacune d’elles afin de quantifier les différences et expliquer leur signification dans le contexte de votre étude initiale.

6) Votre critique va probablement (et malheureusement) insister sur une sorte de valeur comme condition préalable à la publication. Soupir! Si c’était moi, étant donné les différences par rapport à tout, j’utiliserais probablement un test non paramétrique de Kolmogorov-Smirnov pour cracher une valeur qui démontre que les distributions sont différentes, puis procéder à la statistique descriptive comme ci-dessus. Vous devrez ajouter du bruit aux deux échantillons pour vous débarrasser des liens. (Et bien sûr, tout cela suppose que vos échantillons sont indépendants, ce que vous n'avez pas explicitement déclaré.)ppp

Cette réponse est beaucoup plus longue que ce que j'avais initialement prévu. Désolé pour ça.


la source
Je suis curieux de savoir si l’approche suivante quasi visualisée est la suivante: estimez par bootstrap les moments des deux groupes (moyennes, variances et moments les plus hauts, si vous le souhaitez), puis tracez ces estimations et leurs intervalles de confiance en regardant pour le degré de chevauchement entre les groupes à chaque instant. Cela vous permet de parler des différences potentielles entre les différentes caractéristiques de distribution. Si les données sont appariées, calculez les scores de différence et amorcez les moments de cette distribution unique. Pensées?
Mike Lawrence
2
(+1) Bonne analyse. Vous avez parfaitement raison de dire que les résultats sont évidents et qu'il n'est pas nécessaire d'appuyer sur le point avec une valeur p. Vous pouvez être un peu extrême dans votre énoncé de (3), car le test t ne nécessite pas de données distribuées normalement. Si vous êtes concerné, il existe des ajustements pour l'asymétrie (par exemple, la variante de Chen): vous pouvez voir si la valeur p du test ajusté change la réponse. Sinon, vous allez probablement bien. Dans cette situation particulière, avec ces données (fortement asymétriques), le test t fonctionne correctement.
whuber
(+1) belle prise! Et de très bons commentaires.
chl
Nous semblons accepter l'idée que la distribution sous-jacente est "similaire" à l'instanciation aléatoire. On ne peut donc pas poser la question: sont-ils tous deux de la version bêta (0,25, 0,25) et ensuite, le test serait de savoir s'ils ont le même paramètre de (non) centralité. Et cela ne justifierait-il pas l’utilisation d’un test de permutation ou de Wilcoxon?
DWin
4
Il y a beaucoup de bonnes informations ici, mais elles sont très fortement formulées et mélangées avec certaines choses qui n'ont pas beaucoup de sens pour moi. Par exemple, au point 3, ma compréhension de la relation entre le z-test et le t-test est qu’elles sont essentiellement les mêmes, mais z est utilisé lorsque le SD est connu a priori & t est utilisé lorsque le SD est estimé à partir des données. Je ne vois pas comment, si le CLT garantit la normalité des distorsions d'échantillonnage, cela autorise le test z tout en laissant le test t non valide. Je ne pense pas non plus que les SD annulent t si le CLT vous couvre, tant qu’un ajustement (par exemple, Welch – Satterthwaite) est utilisé.
Gay - Rétablir Monica
5

Mes commentaires ne concernent pas la mise en œuvre du test de permutation, mais les problèmes plus généraux soulevés par ces données et leur discussion, en particulier le message de G. Jay Kerns.

Les deux distributions me ressemblent beaucoup, SAUF le groupe des 0 de Y1, qui sont très différents des autres observations de cet échantillon (la plus petite suivante est environ 50 sur l’échelle 0-100) ainsi que toutes celles de X1. J'examinerais d'abord s'il y avait quelque chose de différent à propos de ces observations.

Deuxièmement, en supposant que ces 0 appartiennent à l'analyse, affirmant que le test de permutation n'est pas valide, car les distributions semblent différer soulèvent la question. Si la valeur null était vraie (les distributions sont identiques), pourriez-vous (avec une probabilité raisonnable) avoir des distributions aussi différentes que ces deux-là? Répondre, c'est le but du test, n'est-ce pas? Peut-être que dans ce cas, certains jugeront la réponse évidente sans exécuter le test, mais avec ces distributions minuscules et particulières, je ne pense pas que je le ferais.

Andrew Taylor
la source
Il semble que cela devrait être un ou plusieurs commentaires, pas une réponse. Si vous cliquez sur le petit "ajouter un commentaire" gris, vous pouvez placer vos pensées dans la conversation au-dessous de la question ou d'une réponse particulière, à laquelle elles appartiennent. Vous soulevez des points de fond ici, mais il n’est pas clair que ce n’est pas l’endroit approprié pour eux.
gung - Réintégrer Monica
1
@gung Il faut un peu de réputation pour pouvoir poster un commentaire ;-).
whuber
4
C'est un bon point concernant l'applicabilité du test de permutation. Quant à l'évidence de la différence entre les groupes, c'est peut-être une question d'expérience :-). Pour l'intuition, comme il est évident que la principale différence réside dans les petites valeurs, nous pourrions nous renseigner sur les chances que les sept valeurs les plus petites d'un ensemble de 40 valeurs tombent dans un sous-ensemble aléatoire de 20. En gros, chacune a environ 1/2 chance d'être dans le sous-ensemble ou son complément, donc tous les sept seront dans le même groupe avec des chances autour de . Cette arithmétique mentale fournit une orientation initiale rapide. 2(1/2)7.01
whuber
4

Comme cette question est réapparue , je pourrais ajouter une autre réponse inspirée par un récent article de blog de RAB Bloggers de Robert Kabacoff, auteur de Quick-R et R in Action utilisant le lmPermpackage.

Cependant, cette méthode produit des résultats très contrastés (et très instables) par rapport à celui produit par le coinpackage dans la réponse de @caracakl (la valeur p de l'analyse intra-sujets est 0.008). L'analyse prend la préparation des données de la réponse de @ caracal:

x1 <- c(99, 99.5, 65, 100, 99, 99.5, 99, 99.5, 99.5, 57, 100, 99.5, 
        99.5, 99, 99, 99.5, 89.5, 99.5, 100, 99.5)
y1 <- c(99, 99.5, 99.5, 0, 50, 100, 99.5, 99.5, 0, 99.5, 99.5, 90, 
        80, 0, 99, 0, 74.5, 0, 100, 49.5)

DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
id <- factor(rep(1:length(x1), 2)) 

library(lmPerm)

summary(aovp( DV ~ IV + Error(id)))

produit:

> summary(aovp( DV ~ IV + Error(id)))
[1] "Settings:  unique SS "

Error: id
Component 1 :
          Df R Sum Sq R Mean Sq
Residuals 19    15946       839


Error: Within
Component 1 :
          Df R Sum Sq R Mean Sq Iter Pr(Prob)  
IV         1     7924      7924 1004    0.091 .
Residuals 19    21124      1112                
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 

Si vous l'exécutez plusieurs fois, les valeurs de valeurs oscillent entre ~ 0,05 et ~ 0,1.

Bien qu'il soit une réponse à la question me laisse permets de poser une question à la fin (je peux me déplacer ce à une nouvelle question si on le souhaite):
Toute idée de pourquoi cette analyse est si instable et ne produit une sorte divergeant des valeurs de p à l'analyse des pièces? Est-ce que j'ai fait quelque chose de mal?

Henrik
la source
2
Il peut être préférable de poser cette question séparément, s’il s’agit vraiment d’une question à laquelle vous souhaitez répondre, plutôt que d’une autre solution possible que vous souhaitez énumérer avec le reste. Je remarque que vous spécifiez une strate d'erreur, mais pas @caracal; ce serait ma première hypothèse sur la différence b / t de cette sortie & sa. En outre, lors de la simulation, les valeurs sautent généralement. pour la reproductibilité, vous spécifiez la graine, par exemple set.seed(1); pour plus de précision dans l'estimation MC, vous augmentez le nombre d'itérations; Je ne sais pas si l'une ou l'autre de ces réponses est la «bonne» réponse à votre question, mais elles sont probablement pertinentes.
gung - Réintégrer Monica
2
Encore une fois, je suggère de comparer les résultats de la MC avec les calculs manuels en utilisant le test de permutation complète (re-randomisation). Voir le code de votre exemple pour une comparaison de oneway_anova()(toujours proche du résultat correct) et aovp()(généralement loin du résultat correct). Je ne sais pas pourquoi aovp()donne des résultats extrêmement variables, mais au moins pour ce cas-ci, ils sont invraisemblables. @gung le dernier appel oneway_test(DV ~ IV | id, ...)dans ma réponse originale a spécifié les strates d'erreur pour le cas dépendant (syntaxe différente de celle utilisée par aov()).
Caracal
@ caracal, vous avez raison. Je n'ai pas regardé le dernier bloc de code après l'édition. Je regardais le bloc de code supérieur - bâclé de ma part.
Gay - Rétablir Monica
Je n'ai pas vraiment besoin de réponse. C’est une autre possibilité qui mérite d’être mentionnée ici. Malheureusement, il est loin des autres résultats que je mérite également de noter.
Henrik
1
@Henrik lance aovp avec maxExact = 1000. Si cela prend trop de temps, définissez iter = 1000000 et Ca = 0,001. Le calcul se termine lorsque l'erreur type estimée de p est inférieure à Ca * p. (Les valeurs
basses
1

Ces scores sont-ils des proportions? Si tel est le cas, vous ne devriez certainement pas utiliser un test paramétrique gaussien, et même si vous pouvez utiliser une approche non paramétrique comme un test de permutation ou un bootstrap des moyens, je suggérerais que vous obteniez plus de puissance statistique employant une approche paramétrique non gaussienne appropriée. Plus précisément, chaque fois que vous pouvez calculer une mesure de proportion au sein d’une unité d’intérêt (par exemple un participant à une expérience), vous pouvez et devriez probablement utiliser un modèle à effets mixtes qui spécifie des observations avec une erreur distribuée de façon binomiale. Voir Dixon 2004 .

Mike Lawrence
la source
Les scores ne sont pas des proportions mais des estimations des participants sur une échelle de 0 à 100 (les données présentées sont des moyennes d'estimations pour plusieurs éléments de cette échelle).
Henrik
Les méthodes non paramétriques semblent alors constituer la méthode traditionnelle. Cela dit, je me suis demandé si de telles données d'échelle pourraient utilement être déduites d'un processus binomial et donc analysées comme telles. C'est-à-dire que vous dites que chaque score est la moyenne de plusieurs éléments et que chaque élément correspond à une échelle de 10 points. Dans ce cas, je représenterais une réponse de "8" par exemple qui ont la valeur 1 et deux dont la valeur 0, tous étiquetés avec la même étiquette dans une variable "item". Avec ces données étendues / binomiales, vous pouvez ensuite calculer le modèle binomial à effets mixtes.
Mike Lawrence
Suite à mon précédent commentaire, je devrais noter que, dans les données développées / binomialisées, vous pouvez modéliser la variable "item" comme un effet fixe ou aléatoire. Je pense que je préférerais modéliser cet effet comme un effet fixe, car il est probable que vous voudriez peut-être non seulement expliquer, mais aussi évaluer les différences entre éléments et toute interaction possible entre l'élément et les autres variables prédictives.
Mike Lawrence
0

Juste en ajoutant une autre approche, ezPermdu ezpaquet:

> # preparing the data
> DV <- c(x1, y1)
> IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
> id <- factor(rep(1:length(x1), 2))
> df <- data.frame(id=id,DV=DV,IV=IV)
>
> library(ez)
> ezPerm( data = df, dv = DV, wid = id, within = IV, perms = 1000)
|=========================|100%              Completed after 17 s 
  Effect     p p<.05
1     IV 0.016     *

Cela semble être cohérent avec oneway_testle coinpaquet:

> library(coin)
> pvalue(oneway_test(DV ~ IV | id,  distribution=approximate(B=999999)))
[1] 0.01608002
99 percent confidence interval:
 0.01575782 0.01640682

Cependant, notez qu'il ne s'agit pas du même exemple fourni par @caracal . Dans son exemple, il inclut alternative="greater"donc la différence entre les valeurs p ~0.008et ~0.016.

Le aovppaquet suggéré dans l'une des réponses produit des p-values ​​étrangement inférieures et est incroyablement rapide, même lorsque j'essaie des valeurs élevées pour les arguments Iter, Caet maxIter:

library(lmPerm)
summary(aovp(DV ~ IV + Error(id/IV), data=df,  maxIter = 1000000000))
summary(aovp(DV ~ IV + Error(id/IV), data=df,  Iter = 1000000000))
summary(aovp(DV ~ IV + Error(id/IV), data=df,  Ca = 0.00000000001))

Cela dit, les arguments semblent réduire légèrement les variations des valeurs de p de ( ~.03et ~.1j'ai (une plage plus grande que celle rapportée par @Henrik), de 0.03et 0.07.

toto_tico
la source