P-Werte im Permutationstest gleich 0

Ich habe zwei Datensätze und möchte wissen, ob sie sich erheblich unterscheiden oder nicht (dies kommt von " Zwei Gruppen unterscheiden sich erheblich? Test zu verwenden ").

Ich habe mich für einen Permutationstest entschieden und in R Folgendes durchgeführt:

permutation.test <- function(coding, lncrna) {
    coding <- coding[,1] # dataset1
    lncrna <- lncrna[,1] # dataset2

    ### Under null hyphotesis, both datasets would be the same. So:
    d <- c(coding, lncrna)

    # Observed difference
    diff.observed = mean(coding) - mean(lncrna)
    number_of_permutations = 5000
    diff.random = NULL

    for (i in 1:number_of_permutations) {
        # Sample from the combined dataset
        a.random = sample (d, length(coding), TRUE)
        b.random = sample (d, length(lncrna), TRUE)
        # Null (permuated) difference
        diff.random[i] = mean(b.random) - mean(a.random)
    }

    # P-value is the fraction of how many times the permuted difference is equal or more extreme than the observed difference
    pvalue = sum(abs(diff.random) >= abs(diff.observed)) / number_of_permutations
    pvalue
}

Trotzdem sollten p-Werte laut diesem Artikel nicht 0 sein: http://www.statsci.org/smyth/pubs/permp.pdf

Was empfehlen Sie mir zu tun? So berechnen Sie den p-Wert:

pvalue = sum(abs(diff.random) >= abs(diff.observed)) / number_of_permutations

ein guter Weg? Oder ist es besser, Folgendes zu tun?

pvalue = sum(abs(diff.random) >= abs(diff.observed)) + 1 / number_of_permutations + 1

Diskussion

Ein Permutationstest generiert alle relevanten Permutationen eines Datensatzes, berechnet für jede dieser Permutationen eine festgelegte Teststatistik und bewertet die tatsächliche Teststatistik im Kontext der resultierenden Permutationsverteilung der Statistiken. Eine übliche Methode zur Bewertung besteht darin, den Anteil der Statistiken zu melden, der (in gewissem Sinne) "als oder extremer" als die tatsächliche Statistik ist. Dies wird oft als "p-Wert" bezeichnet.

Da es sich bei dem tatsächlichen Datensatz um eine dieser Permutationen handelt, gehört seine Statistik zwangsläufig zu denjenigen, die in der Permutationsverteilung zu finden sind. Daher kann der p-Wert niemals Null sein.

Sofern der Datensatz nicht sehr klein ist (normalerweise weniger als 20-30 Gesamtzahlen) oder die Teststatistik eine besonders schöne mathematische Form hat, ist es nicht praktikabel, alle Permutationen zu generieren. (Ein Beispiel, in dem alle Permutationen generiert werden, wird unter Permutationstest in R angezeigt .) Daher werden Computerimplementierungen von Permutationstests normalerweise aus der Permutationsverteilung entnommen . Sie erzeugen dazu einige unabhängige zufällige Permutationen und hoffen, dass die Ergebnisse eine repräsentative Stichprobe aller Permutationen sind.

Daher sind alle von einer solchen Stichprobe abgeleiteten Zahlen (wie z. B. ein "p-Wert") nur Schätzer für die Eigenschaften der Permutationsverteilung. Es ist durchaus möglich - und häufig bei großen Effekten -, dass der geschätzte p-Wert Null ist. Daran ist nichts auszusetzen, aber es wirft sofort die bisher vernachlässigte Frage auf, inwieweit der geschätzte p-Wert vom richtigen Wert abweichen kann. Da die Stichprobenverteilung eines Anteils (z. B. ein geschätzter p-Wert) binomisch ist, kann dieser Unsicherheit mit einem binomischen Konfidenzintervall begegnet werden .

Die Architektur

Eine gut aufgebaute Implementierung wird die Diskussion in jeder Hinsicht genau verfolgen. Es würde mit einer Routine beginnen, die Teststatistik zu berechnen, da diese die Mittelwerte von zwei Gruppen vergleicht:

diff.means <- function(control, treatment) mean(treatment) - mean(control)

Schreiben Sie eine weitere Routine, um eine zufällige Permutation des Datensatzes zu generieren und die Teststatistik anzuwenden. Die Schnittstelle zu dieser ermöglicht es dem Aufrufer, die Teststatistik als Argument anzugeben. Es werden die ersten mElemente eines Arrays (vermutlich eine Referenzgruppe) mit den verbleibenden Elementen (der "Behandlungs" -Gruppe) verglichen .

f <- function(..., sample, m, statistic) {
  s <- sample(sample)
  statistic(s[1:m], s[-(1:m)])
}

Die Permutation Test wird bestimmt, indem die Statistik für die eigentlichen Daten zuerst durchgeführt (hier angenommen , in zwei Arrays gespeichert werden controlund treatment) und dann Statistiken für viele unabhängigen Zufälle Permutationen davon zu finden:

z <- stat(control, treatment) # Test statistic for the observed data
sim<- sapply(1:1e4, f, sample=c(control,treatment), m=length(control), statistic=diff.means)

Berechnen Sie nun die Binomialschätzung des p-Wertes und ein Konfidenzintervall dafür. Eine Methode verwendet die binconfin das HMiscPaket integrierte Prozedur :

require(Hmisc)                                    # Exports `binconf`
k <- sum(abs(sim) >= abs(z))                      # Two-tailed test
zapsmall(binconf(k, length(sim), method='exact')) # 95% CI by default

Es ist keine schlechte Idee, das Ergebnis mit einem anderen Test zu vergleichen, auch wenn bekannt ist, dass dies nicht ganz zutreffend ist: Zumindest könnte man eine Größenordnung erkennen, wo das Ergebnis liegen sollte. In diesem Beispiel (zum Vergleichen der Mittelwerte) liefert ein Student-T-Test normalerweise trotzdem ein gutes Ergebnis:

t.test(treatment, control)

Diese Architektur wird in einer komplexeren Situation mit RArbeitscode unter Testen, ob Variablen der gleichen Verteilung folgen dargestellt .

Beispiel

$10$$0$$20$$1.5$

set.seed(17)
control <- rnorm(10)
treatment <- rnorm(20, 1.5)

Nachdem ich mit dem obigen Code einen Permutationstest durchgeführt hatte, zeichnete ich die Stichprobe der Permutationsverteilung zusammen mit einer vertikalen roten Linie auf, um die tatsächliche Statistik zu markieren:

h <- hist(c(z, sim), plot=FALSE)
hist(sim, breaks=h$breaks)
abline(v = stat(control, treatment), col="Red")

Die Berechnung der Binomialvertrauensgrenze ergab

 PointEst Lower        Upper
        0     0 0.0003688199

$0$$0.00037$3.16e-05$0.00037$$0.00037$$0.05$$0.01$$0.001$

Bemerkungen

$k$$N$ $k/N$$(k+1)/(N+1)$$N$

$10$$10^2=100$$0.000005$$1.6$$11.7$parts per million: etwas kleiner als der angegebene Student-T-Test. Obwohl die Daten mit normalen Zufallszahlengeneratoren generiert wurden, was die Verwendung des Student-T-Tests rechtfertigen würde, weichen die Ergebnisse des Permutationstests von den Ergebnissen des Student-T-Tests ab, da die Verteilungen innerhalb der einzelnen Beobachtungsgruppen nicht völlig normal sind.

$\frac{B}{M}$$\frac{B+1}{M+1}$

(B ist die Anzahl der zufälligen Permutationen, bei denen eine Statistik erhalten wird, die größer oder gleich der beobachteten ist, und M ist die Gesamtzahl der abgetasteten zufälligen Permutationen).

$\frac{B}{M}$