Welche Permutationstest-Implementierung in R soll anstelle von t-Tests (gepaart und ungepaart) verwendet werden?

Ich habe Daten aus einem Experiment, das ich mit t-Tests analysiert habe. Die abhängige Variable ist intervallskaliert und die Daten sind entweder ungepaart (dh 2 Gruppen) oder gepaart (dh innerhalb von Subjekten). ZB (innerhalb von Fächern):

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)

Da die Daten jedoch nicht normal sind, hat uns ein Rezensent gebeten, etwas anderes als den T-Test zu verwenden. Wie man jedoch leicht sehen kann, sind die Daten nicht nur nicht normal verteilt, sondern die Verteilungen sind zwischen den Bedingungen nicht gleich:

Daher können die üblichen nichtparametrischen Tests, der Mann-Whitney-U-Test (ungepaart) und der Wilcoxon-Test (gepaart), nicht verwendet werden, da sie gleiche Verteilungen zwischen den Bedingungen erfordern. Daher entschied ich mich für einen Resampling- oder Permutationstest.

Jetzt suche ich nach einer R-Implementierung eines permutationsbasierten Äquivalents des t-Tests oder einem anderen Rat, was mit den Daten zu tun ist.

Ich weiß, dass es einige R-Pakete gibt, die dies für mich tun können (z. B. Coin, Perm, exactRankTest usw.), aber ich weiß nicht, welches ich auswählen soll. Also, wenn jemand mit etwas Erfahrung mit diesen Tests mir einen Kick-Start geben könnte, wäre das überaus cool.

UPDATE: Es wäre ideal, wenn Sie ein Beispiel dafür liefern könnten, wie die Ergebnisse dieses Tests gemeldet werden.

Es sollte nicht so wichtig sein, da die Teststatistik immer die Mittelwertdifferenz (oder etwas Äquivalentes) ist. Kleine Unterschiede können sich aus der Implementierung von Monte-Carlo-Methoden ergeben. Probieren Sie die drei Pakete mit Ihren Daten mit einem einseitigen Test für zwei unabhängige Variablen aus:

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

Um den genauen p-Wert mit einer manuellen Berechnung aller Permutationen zu überprüfen, beschränke ich die Daten auf die ersten 9 Werte.

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

coinund exactRankTestssind beide vom selben Autor, coinscheinen aber allgemeiner und umfangreicher zu sein - auch in Bezug auf die Dokumentation. exactRankTestsist nicht mehr aktiv entwickelt. Ich würde daher wählen coin(auch wegen informativer Funktionen wie support()), es sei denn, Sie möchten sich nicht mit S4-Objekten befassen.

EDIT: Für zwei abhängige Variablen lautet die Syntax

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

Ein paar Bemerkungen sind meiner Meinung nach angebracht.

1) Ich würde Sie dazu ermutigen, mehrere visuelle Darstellungen Ihrer Daten zu versuchen, da diese Dinge erfassen können, die durch (Diagramme wie) Histogramme verloren gegangen sind, und ich empfehle Ihnen nachdrücklich, auf nebeneinander liegenden Achsen zu zeichnen. In diesem Fall glaube ich nicht, dass die Histogramme die herausragenden Merkmale Ihrer Daten sehr gut vermitteln. Schauen Sie sich beispielsweise nebeneinander liegende Boxplots an:

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

Oder auch Side-by-Side-Stripcharts:

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

Schauen Sie sich die Zentren, Ausbreitungen und Formen an! Etwa drei Viertel der Daten liegen weit über dem Median der Daten. Die Ausbreitung von ist winzig, während die Ausbreitung von groß ist. Sowohl als auch sind stark nach links geneigt, jedoch auf unterschiedliche Weise. Zum Beispiel hat fünf (!) Wiederholte Werte von Null.y 1 x 1 y 1 x 1 y 1 y $x1$$y1$$x1$$y1$$x1$$y1$$y1$

2) Sie haben weder ausführlich erklärt, woher Ihre Daten stammen, noch wie sie gemessen wurden, aber diese Informationen sind sehr wichtig, wenn es darum geht, ein statistisches Verfahren auszuwählen. Sind Ihre beiden obigen Proben unabhängig? Gibt es Gründe zu der Annahme, dass die Grenzverteilungen der beiden Stichproben gleich sein sollten (außer zum Beispiel bei unterschiedlichen Standorten)? Welche Überlegungen haben Sie vor der Studie veranlasst, nach Belegen für einen Unterschied zwischen den beiden Gruppen zu suchen?

3) Der t-Test ist für diese Daten nicht geeignet, da die Randverteilungen deutlich nicht normal sind und in beiden Proben Extremwerte aufweisen. Wenn Sie möchten, können Sie auf die CLT ansprechen (aufgrund Ihrer mäßig großen Probe) ein verwenden -Test (die zu einem z-Test für große Proben ähnlich sein würden), aber die Schiefe (in beiden Variablen) gegeben von Ihre Angaben würde ich als solche Berufung nicht sehr überzeugend beurteilen. Natürlich können Sie damit auch einen Wert berechnen , aber was bringt das für Sie? Wenn die Annahmen nicht erfüllt sind, ist ein Wert nur eine Statistik. es sagt nichts darüber aus, was Sie (vermutlich) wissen wollen: ob es Beweise dafür gibt, dass die beiden Samples aus unterschiedlichen Distributionen stammen.p $z$$p$$p$

4) Ein Permutationstest wäre für diese Daten ebenfalls ungeeignet. Die einzige und häufig übersehene Annahme für Permutationstests ist, dass die beiden Proben unter der Nullhypothese austauschbar sind . Das würde bedeuten, dass sie identische Randverteilungen haben (unter der Null). Sie stecken jedoch in Schwierigkeiten, da die Grafiken darauf hindeuten, dass sich die Verteilungen sowohl in der Position als auch im Maßstab (und auch in der Form) unterscheiden. Sie können also nicht (gültig) auf Standortunterschiede testen, weil sich die Maßstäbe unterscheiden, und Sie können nicht (gültig) auf Maßstabunterschiede testen, weil sich die Standorte unterscheiden. Hoppla. Sie können den Test trotzdem durchführen und einen Wert erhalten, aber was nun? Was hast du wirklich erreicht?$p$

5) Meiner Meinung nach sind diese Daten ein perfektes (?) Beispiel dafür, dass ein gut ausgewähltes Bild 1000 Hypothesentests wert ist. Wir brauchen keine Statistiken, um den Unterschied zwischen einem Bleistift und einer Scheune zu erkennen. Für diese Daten würde aus meiner Sicht die entsprechende Aussage lauten: "Diese Daten weisen deutliche Unterschiede in Bezug auf Position, Maßstab und Form auf." Anschließend können Sie eine (aussagekräftige) deskriptive Statistik für jede dieser Statistiken erstellen, um die Unterschiede zu quantifizieren und zu erläutern, was die Unterschiede im Kontext Ihrer ursprünglichen Studie bedeuten.

6) Ihr Rezensent wird wahrscheinlich (und leider) auf einer Art Wert als Voraussetzung für die Veröffentlichung bestehen. Seufzer! Wenn ich es wäre, würde ich angesichts der Unterschiede in Bezug auf alles wahrscheinlich einen nichtparametrischen Kolmogorov-Smirnov-Test verwenden, um einen Wert auszuspucken , der zeigt, dass die Verteilungen unterschiedlich sind, und dann mit der oben beschriebenen Statistik fortzufahren. Sie müssten den beiden Samples etwas Rauschen hinzufügen, um die Krawatten loszuwerden. (Und dies setzt natürlich voraus, dass Ihre Stichproben unabhängig sind, was Sie nicht ausdrücklich angegeben haben.)$p$$p$

Diese Antwort ist viel länger, als ich ursprünglich beabsichtigt hatte. Das tut mir leid.

In meinen Kommentaren geht es nicht um die Implementierung des Permutationstests, sondern um die allgemeineren Fragen, die durch diese Daten aufgeworfen werden, und die Diskussion darüber, insbesondere um den Beitrag von G. Jay Kerns.

Die beiden Verteilungen sehen mir mit AUSNAHME der Gruppe der Nullen in Y1 sehr ähnlich, was sich von den anderen Beobachtungen in dieser Stichprobe (die nächstkleinere ist etwa 50 auf der Skala von 0 bis 100) sowie von allen in X1 stark unterscheidet. Ich würde zuerst untersuchen, ob an diesen Beobachtungen etwas anderes war.

Zweitens stellt sich die Frage, wenn angenommen wird, dass diese Nullen in die Analyse gehören und der Permutationstest nicht gültig ist, weil die Verteilungen unterschiedlich zu sein scheinen. Wenn die Null wahr wäre (Verteilungen sind identisch), könnten Sie (mit angemessener Wahrscheinlichkeit) Verteilungen erhalten, die so unterschiedlich aussehen wie diese beiden? Antworten, das ist der springende Punkt des Tests, nicht wahr? Vielleicht werden in diesem Fall einige die Antwort als offensichtlich betrachten, ohne den Test auszuführen, aber mit diesen kleinen, eigenartigen Verteilungen glaube ich nicht, dass ich das tun würde.

Da diese Frage erneut auftauchte, kann ich eine weitere Antwort hinzufügen, die von einem kürzlich veröffentlichten Blogbeitrag über R-Blogger von Robert Kabacoff, dem Autor von Quick-R und R in Action , der das lmPermPaket verwendet, inspiriert wurde .

Diese Methode erzeugt jedoch sehr kontrastreiche (und sehr instabile) Ergebnisse zu denen, die von der coinPackung in der Antwort von @caracakl erzeugt wurden (der p-Wert der Analyse innerhalb der Subjekte ist 0.008). Die Analyse bezieht die Datenaufbereitung auch aus der Antwort von @ 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)))

produziert:

> 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 

Wenn Sie dies mehrmals ausführen, bewegen sich die p-Werte zwischen ~ .05 und ~ .1.

Obwohl es sich um eine Antwort auf die Frage handelt, lassen Sie mich am Ende eine Frage stellen (ich kann diese auf Wunsch auf eine neue Frage verschieben):
Irgendwelche Ideen, warum diese Analyse so instabil ist und zu so unterschiedlichen p-Werten führt die Münzanalyse? Habe ich etwas falsch gemacht?

Sind diese Punkte Proportionen? In diesem Fall sollten Sie auf keinen Fall einen Gauß-Parametertest verwenden, und obwohl Sie einen nicht-parametrischen Ansatz wie einen Permutationstest oder einen Bootstrap der Mittel anwenden könnten, würde ich vorschlagen, dass Sie mehr statistische Leistung erhalten Anwendung eines geeigneten nicht-gaußschen parametrischen Ansatzes. Insbesondere können und sollten Sie jedes Mal, wenn Sie ein Proportionsmaß innerhalb einer interessierenden Einheit (z. B. Teilnehmer an einem Experiment) berechnen können, ein Modell mit gemischten Effekten verwenden, das Beobachtungen mit binomial verteiltem Fehler angibt. Siehe Dixon 2004 .

Hinzufügen eines anderen Ansatzes ezPermdes ezPakets:

> # 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     *

Dies scheint im Einklang mit oneway_testdem coinPaket zu stehen:

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

Beachten Sie jedoch, dass dies nicht dasselbe Beispiel ist, das von @caracal bereitgestellt wird . In seinem Beispiel enthält er alternative="greater"daher den Unterschied in der p-Werten ~0.008vs ~0.016.

Das aovpPaket vorgeschlagen in eine der Antworten erzeugt verdächtig niedrige p-Werte und läuft verdächtig schnell , auch wenn ich hohe Werte für die versuchen Iter, Caund maxIterArgumente:

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))

Das heißt, die Argumente scheinen die Variationen der p-Werte von ~.03und ~.1(ich habe einen größeren Bereich als den von @Henrik berichteten) zu 0.03und etwas zu verringern 0.07.