Testen Sie, ob Variablen der gleichen Verteilung folgen

Wenn Sie testen möchten, ob zwei Variablen der gleichen Verteilung folgen, ist es ein guter Test, einfach beide Variablen zu sortieren und dann ihre Korrelation zu überprüfen? Wenn es hoch ist (mindestens 0,9?), Stammen die Variablen höchstwahrscheinlich aus derselben Verteilung.

Mit Verteilung meine ich hier "normal", "Chi-Quadrat", "Gamma" usw.

Lassen Sie uns herausfinden, ob dies ein guter Test ist oder nicht. Es steckt viel mehr dahinter, als nur zu behaupten, es sei schlecht oder zu zeigen, dass es nicht gut funktioniert. Die meisten Tests funktionieren unter bestimmten Umständen schlecht. Daher müssen wir häufig feststellen, unter welchen Umständen ein vorgeschlagener Test möglicherweise eine gute Wahl darstellt.

Beschreibung des Tests

Wie jeder Hypothesentest besteht dieser aus (a) einer Null- und einer Alternativhypothese und (b) einer Teststatistik (dem Korrelationskoeffizienten), mit der zwischen den Hypothesen unterschieden werden soll.

Die Nullhypothese lautet, dass die beiden Variablen aus derselben Verteilung stammen. Um genau zu sein, nennen wir die Variablen und Y und nehmen an, wir haben n x Instanzen von X mit der Bezeichnung x i = ( x 1 , x 2 , … , x n x ) und n y Instanzen von Y mit der Bezeichnung y beobachtet ich . Die Nullhypothese ist, dass alle Instanzen von $X$$Y$$n_x$$X$$x_i = (x_1, x_2, \ldots, x_{n_x})$$n_y$$Y$$y_i$$X$ und unabhängig und identisch verteilt sind (iid).$Y$

Nehmen wir als alternative Hypothese an, dass (a) alle Instanzen von gemäß einer zugrunde liegenden Verteilung F X sind und (b) alle Instanzen von Y gemäß einer zugrunde liegenden Verteilung F Y sind, aber (c) F X von F verschieden ist Y . (Daher werden wir nicht nach Korrelationen zwischen x i , Korrelationen zwischen y i , Korrelationen zwischen x i und y j oder Verteilungsunterschieden zwischen x 's oder y suchen$X$$F_X$$Y$$F_Y$$F_X$$F_Y$$x_i$$y_i$$x_i$$y_j$$x$$y$ist separat: das wird als nicht plausibel angenommen.)

Die vorgeschlagene Teststatistik wird davon ausgegangen , dass (nennen wir diesen gemeinsamen Wert n ) und berechnet den Korrelationskoeffizienten des ( x [ i ] , Y [ i ] ) (wobei, wie üblich, [ i ] bezeichnet den i - ten kleinsten der Daten). Nenne dies t ( x , y ) .$n_x = n_y$$n$$(x_{[i]}, y_{[i]})$$[i]$$i^\text{th}$$t(x,y)$

Permutationstests

In dieser Situation - egal welche Statistik vorgeschlagen wird - wir können immer ein Verhalten Permutationstest. Unter der Nullhypothese ist die Wahrscheinlichkeit der Daten ( ( x 1 , x 2 , … , x n ) , ( y 1 , y 2 , … , y n ) ) dieselbe wie die Wahrscheinlichkeit einer Permutation der 2 n Datenwerte. Mit anderen Worten die Zuordnung der Hälfte der Daten zu $t$$\left((x_1, x_2, \ldots, x_n), (y_1, y_2, \ldots, y_n)\right)$$2n$$X$ und der anderen Hälfte zu $Y$ist ein reiner Zufall. Dies ist eine einfache, direkte Folge der iid-Annahmen und der Nullhypothese, dass $F_X=F_Y$ .

Daher ist die Stichprobenverteilung von , bedingt auf den Beobachtungen x i und y i ist die Verteilung aller Werte von t für alle erreicht ( 2 n ) ! Permutationen der Daten. Dies interessiert uns, da wir für jede vorgesehene Testgröße α , wie α = 0,05 (entsprechend 95 % Konfidenz), einen zweiseitigen kritischen Bereich aus der Stichprobenverteilung von konstruieren$t(x,y)$$x_i$$y_i$$t$$(2n)!$$\alpha$$\alpha = .05$$95$: Er besteht aus dem extremsten$t$ % der möglichen Werte von t (auf der hohen Seite, weil eine hohe Korrelation mit ähnlichen Verteilungen konsistent ist und eine niedrige Korrelation nicht). So bestimmen wir, wie groß der Korrelationskoeffizient sein muss, um zu entscheiden, ob die Daten aus unterschiedlichen Verteilungen stammen.$100\alpha$$t$

Simulation der Null-Sampling-Verteilung

Weil (oder, wenn Sie möchten, ( 2 $(2n)!$, die die Anzahl der Arten der Spaltung zählt die2nDaten in zwei Stücke der Größen) wird großselbst für kleinen,es nicht praktikabel istdie Stichprobenverteilung genau zu berechnen, so dass wir sie abtasten eine Simulation. (Zum Beispiel, wennn=16, ( 2$\binom{2n}{n}/2$$2n$$n$$n$$n=16$und(2n)! ≤2,63×10$\binom{2n}{n}/2 = 300\ 540\ 195$$(2n)! \approx 2.63\times 10^{35}$ ) Ungefähr tausend Proben reichen oft aus (und werden es sicherlich für die Erkundungen tun, die wir in Kürze durchführen werden).

Wir müssen zwei Dinge herausfinden: Erstens, wie sieht die Stichprobenverteilung unter der Nullhypothese aus. Zweitens, wie gut unterscheidet dieser Test zwischen verschiedenen Verteilungen?

Es gibt eine Komplikation: Die Stichprobenverteilung hängt von der Art der Daten ab. Alles, was wir tun können, ist, realistische Daten zu betrachten, die erstellt wurden, um das zu emulieren, woran wir interessiert sind, und zu hoffen, dass das, was wir aus den Simulationen lernen, auf unsere eigene Situation zutrifft.

Implementierung

Zur Veranschaulichung habe ich diese Arbeit in ausgeführt R. Es fällt natürlich in drei Teile.

1. Eine Funktion zur Berechnung der Teststatistik . Da ich etwas allgemeiner sein möchte, verarbeitet meine Version Datensätze unterschiedlicher Größe ( n x ≠ n y ), indem sie linear zwischen den Werten im (sortierten) größeren Datensatz interpoliert, um Übereinstimmungen mit dem (sortierten) kleineren Datensatz zu erstellen. Da dies bereits von der Funktion erledigt wird , nehme ich nur die Ergebnisse:$t(x,y)$$n_x \ne n_y$Rqqplot

   test.statistic <- function(x, y) {
     transform <- function(z) -log(1-z^2)/2
     fit <- qqplot(x,y, plot.it=FALSE)
     transform(cor(fit$x, fit$y))
   }

   Eine kleine Drehung - unnötig, aber hilfreich für die Visualisierung - drückt den Korrelationskoeffizienten so um, dass die Verteilung der Nullstatistik annähernd symmetrisch wird. Das ist, was transformmacht.

2. Die Simulation der Stichprobenverteilung. Für die Eingabe akzeptiert diese Funktion die Anzahl der Iterationen n.iterzusammen mit den beiden Datensätzen in Arrays xund y. Es gibt ein Array von n.iterWerten der Teststatistik aus. Das Innenleben sollte auch für Nichtbenutzer transparent sein R:

   permutation.test <- function(n.iter, x, y) {
     z <- c(x,y)
     n.x <- length(x)
     n.y <- length(y)
     n <- length(z)
     k <- min(n.x, n.y)
     divide <- function() {
       i <- sample.int(n, size=k)
       test.statistic(z[i], z[-i])
     }
     replicate(n.iter, divide())
   }

3. Obwohl dies alles ist, was wir zur Durchführung des Tests benötigen, werden wir den Test zum Studieren viele Male wiederholen wollen. Also führen wir den Test einmal durch und binden diesen Code in eine dritte Funktionsschicht fein, die wir hier allgemein nennen und die wir wiederholt aufrufen können. Damit es für eine umfassende Studie ausreichend allgemein ist, akzeptiert es für die Eingabe die Größe der zu simulierenden Datensätze ( n.xund n.y), die Anzahl der Iterationen für jeden Permutationstest ( n.iter) und einen Verweis auf die Funktion testzur Berechnung der Teststatistik (siehe Abbildung) momentan, warum wir dies vielleicht nicht hart codieren wollen, und zwei Funktionen, um iid-Zufallswerte zu generieren, eine für ( ) und eine für Y ( ). Eine Option$X$dist.x$Y$dist.yplot.it ist nützlich, um zu sehen, was los ist.

   f <- function(n.x, n.y, n.iter, test=test.statistic, dist.x=runif, dist.y=runif, 
       plot.it=FALSE) {
     x <- dist.x(n.x)
     y <- dist.y(n.y)
     if(plot.it) qqplot(x,y)
   
     t0 <- test(x,y)
     sim <- permutation.test(n.iter, x, y)
     p <- mean(sim > t0) + mean(sim==t0)/2
     if(plot.it) {
       hist(sim, xlim=c(min(t0, min(sim)), max(t0, max(sim))), 
            main="Permutation distribution")
       abline(v=t0, col="Red", lwd=2)
     }
     return(p)
   }

   Der Ausgang ist ein simulierter „p-Wert“: Der Anteil der Simulationen eine Statistik wodurch man sieht , dass mehr extreme als die tatsächlich für die Daten berechnet.

Die Teile (2) und (3) sind äußerst allgemein gehalten: Sie können eine Studie wie diese für einen anderen Test durchführen, indem Sie sie einfach durch test.statisticeine andere Berechnung ersetzen . Das machen wir weiter unten.

Erste Ergebnisse

Standardmäßig vergleicht unser Code Daten aus zwei gleichmäßigen Verteilungen. Ich lasse es das tun (für , die ziemlich kleine Datensätze sind und daher einen mäßig schwierigen Testfall darstellen) und wiederhole es dann für einen Vergleich von Normal zu Normal und von Exponential zu Exponential. (Gleichmäßige Verteilungen sind nicht leicht von normalen Verteilungen zu unterscheiden, es sei denn, Sie haben etwas mehr als 16 Werte, aber exponentielle Verteilungen - mit einer hohen Schiefe und einem langen rechten Schwanz - sind normalerweise leicht von gleichmäßigen Verteilungen zu unterscheiden.)$n.x = n.y = 16$$16$

set.seed(17)             # Makes the results reproducible
n.per.rep <- 1000        # Number of iterations to compute each p-value
n.reps <- 1000           # Number of times to call `f`
n.x <- 16; n.y <- 16     # Dataset sizes

par(mfcol=c(2,3))        # Lay results out in three columns
null <- replicate(n.reps, f(n.x, n.y, n.per.rep))
hist(null, breaks=20)
plot(null)

normal <- replicate(n.reps, f(n.x, n.y, n.per.rep, dist.y=rnorm))
hist(normal, breaks=20)
plot(normal)

exponential <- replicate(n.reps, f(n.x, n.y, n.per.rep, dist.y=function(n) rgamma(n, 1)))
hist(exponential, breaks=20)
plot(exponential)

$X$$Y$$X$$Y$

$16$$x_i$$16$$y_i$f$0.05$$11$$16$unabhängige Werte von jedem. Das ist ziemlich wenig Strom. Aber vielleicht ist es unvermeidlich, also lasst uns fortfahren.

Die rechten Diagramme testen auf ähnliche Weise eine gleichmäßige Verteilung gegen eine exponentielle. Dieses Ergebnis ist bizarr. Dieser Test lässt meistens den Schluss zu, dass einheitliche Daten und exponentielle Daten gleich aussehen. Es scheint , zu „denken“ , dass einheitliche und exponentielles variates ist mehr ähnlich als zwei einheitliche Variablen! Was ist denn hier los?

Das Problem ist, dass Daten aus einer Exponentialverteilung dazu neigen, einige extrem hohe Werte zu haben. Wenn Sie ein Streudiagramm dieser Werte gegen gleichmäßig verteilte Werte erstellen, befinden sich im Übrigen einige Punkte weit rechts oben. Das entspricht einem sehr hohen Korrelationskoeffizienten. Somit , wenn entweder der Verteilungen einige Extremwerte erzeugt, wird der Korrelationskoeffizient eine schreckliche Wahl für die Messung , wie unterschiedlich die Verteilungen sind. Dies führt zu einem weiteren, noch schlimmeren Problem: Je größer der Datensatz wird, desto größer ist die Wahrscheinlichkeit, dass Sie einige extreme Beobachtungen erhalten. Daher können wir erwarten , dass dieser Test mit zunehmender Datenmenge immer schlechter wird. Wie schrecklich ...

Ein besserer Test

$y=x$

Hier ist eine RImplementierung:

test.statistic <- function(x, y) {
  ks.test(x,y)$statistic
}

Das ist richtig: Es ist in die Software integriert, also müssen wir es nur aufrufen. Aber warte! Wenn Sie das Handbuch sorgfältig lesen, werden Sie feststellen, dass (a) der Test einen p-Wert liefert, (b) dieser p-Wert jedoch (grob) falsch ist, wenn beide xund yDatensätze sind. Es ist zur Verwendung gedacht, wenn Sie glauben, genau zu wissen , von welcher Verteilung die Daten xstammen und ob dies zutrifft. Daher berücksichtigt der Test die Unsicherheit über die Verteilung, aus der die Daten ystammen, nicht richtig .

Kein Problem! Das Permutationstest-Framework ist immer noch genauso gültig. Nachdem wir die vorherige Änderung an vorgenommen test.statistichaben, müssen wir die vorherige Studie unverändert wiederholen. Hier sind die Ergebnisse.

$p=0.20$

$70$$0.05$$11$

$30$$\alpha=5$$50$$\alpha=10$$0.10$

Schlussfolgerungen

Somit sind die Probleme mit dem Korrelationstest nicht auf eine inhärente Schwierigkeit bei dieser Einstellung zurückzuführen. Der Korrelationstest ist nicht nur sehr schlecht, sondern auch im Vergleich zu einem allgemein bekannten und verfügbaren Test schlecht. (Ich würde vermuten, dass es unzulässig ist, was bedeutet, dass es im Durchschnitt immer schlechter abschneidet als die Permutationsversion des KS-Tests, was bedeutet, dass es keinen Grund gibt, es jemals zu verwenden.)

Nein, die Korrelation ist kein guter Test dafür.

x <- 1:100 #Uniform
y <- sort(rnorm(100)) #Normal
cor(x,y) #.98

Ich kenne keinen guten Test, der vergleicht, ob z. B. zwei Verteilungen normal sind, aber möglicherweise mit unterschiedlichem Mittelwert und indirekter SD-Verteilung. Sie könnten die Normalität von jeder einzeln testen, und wenn beide normal erscheinen, raten Sie, dass beide normal sind.

Wenn es eine ausreichend große Anzahl von Variablen gibt, kann dies eine stärkere Korrelation mit Größenordnungswerten anzeigen. Es scheint jedoch keine besonders nützliche Methode zu sein, nicht zuletzt, weil sie wenig Mittel bietet, um das Vertrauen zu schätzen, dass sie dasselbe Modell verwenden könnten.

Ein Problem, das Sie wahrscheinlich haben, ist, wenn Sie Modelle mit einem ähnlichen Mittelwert und einer ähnlichen Schiefe, aber einem Unterschied in der Kurtosis haben, da eine moderate Anzahl von Messungen ausreichend gut passt, um ziemlich gut korreliert auszusehen.

Es erscheint sinnvoller, beide Variablen anhand unterschiedlicher Verteilungen zu modellieren, um die jeweils wahrscheinlichste zu ermitteln und die Ergebnisse zu vergleichen.

Das Normalisieren beider Werte, das Sortieren und Zeichnen beider Werte kann von Vorteil sein. Auf diese Weise können Sie sehen, wie die Übereinstimmungen verglichen werden. Außerdem können Sie ein mögliches Modell für beide Werte zeichnen, das sich auf die von Ihnen vorgeschlagenen Werte bezieht Erwarten Sie eine konkrete Antwort, nur eine visuelle Vorstellung von der Nähe der Distributionen.