Sprawdź, czy zmienne mają ten sam rozkład

Jeśli chcesz sprawdzić, czy dwie zmienne mają ten sam rozkład, czy dobrym pomysłem byłoby posortowanie obu zmiennych, a następnie sprawdzenie ich korelacji? Jeśli jest wysoki (co najmniej 0,9?), Wówczas zmienne najprawdopodobniej pochodzą z tego samego rozkładu.

Z rozkładem tutaj rozumiem „normalny”, „chi-kwadrat”, „gamma” itp.

Sprawdźmy, czy to dobry test, czy nie. Jest o wiele więcej niż tylko twierdzenie, że jest źle lub pokazanie w jednym przypadku, że nie działa dobrze. Większość testów działa w niektórych okolicznościach źle, dlatego często mamy do czynienia z określeniem okoliczności, w których każdy proponowany test może być dobrym wyborem.

Opis testu

Jak każdy test hipotezy, ta składa się z (a) hipotezy zerowej i alternatywnej oraz (b) statystyki testowej (współczynnika korelacji) mającego na celu rozróżnienie hipotez.

Hipotezą zerową jest to, że dwie zmienne pochodzą z tego samego rozkładu. Mówiąc ściślej, nazwijmy zmienne i Y i załóżmy, że zaobserwowaliśmy n x wystąpień X , zwanych x i = ( x 1 , x 2 , … , x n x ) , oraz n y wystąpień Y , zwanych y i . Hipoteza zerowa jest taka, że ​​wszystkie wystąpienia X i Y są niezależne i identycznie rozmieszczone (iid).$X$$Y$$n_x$$X$$x_i = (x_1, x_2, \ldots, x_{n_x})$$n_y$$Y$$y_i$$X$$Y$

Przyjmijmy jako alternatywną hipotezę, że (a) wszystkie wystąpienia są iid zgodnie z pewnym podstawowym rozkładem F X i (b) wszystkie wystąpienia Y są iid zgodnie z pewnym podstawowym rozkładem F Y, ale (c) F X różni się od F Y . (W ten sposób nie będziemy szukać korelacji wśród osób x i korelacje wśród y í , korelacje między x i i y j , czy różnice rozkładu między x „y lub $X$$F_X$$Y$$F_Y$$F_X$$F_Y$$x_i$$y_i$$x_i$$y_j$$x$$y$jest osobno: zakłada się, że nie jest to prawdopodobne).

Proponowany statystyka testu zakłada się, że (połączenie to wspólna wartość n ) i oblicza się współczynnik korelacji, ( X [ I ] , Y [ i ] ) (w którym, jak zwykle, [ I ] Wyznacza i p najmniejsza danych). Nazwij to t ( x , y ) .$n_x = n_y$$n$$(x_{[i]}, y_{[i]})$$[i]$$i^\text{th}$$t(x,y)$

Testy permutacyjne

W tej sytuacji - bez względu na proponowaną statystykę - zawsze możemy przeprowadzić test permutacyjny. Zgodnie z hipotezą zerową prawdopodobieństwo danych ( ( x 1 , x 2 , … , x n ) , ( y 1 , y 2 , … , y n ) ) jest takie samo, jak prawdopodobieństwo permutacji 2 n wartości danych. Innymi słowy, przypisanie połowy danych do X, a drugiej połowy do $t$$\left((x_1, x_2, \ldots, x_n), (y_1, y_2, \ldots, y_n)\right)$$2n$$X$$Y$jest czystym przypadkowym zbiegiem okoliczności. Jest to prosty, bezpośrednią konsekwencją IID założeń i hipoteza zerowa, że .$F_X=F_Y$

W związku z tym, rozmieszczenie próbkowania , uzależnione od obserwacji x I i y I jest rozmieszczenie wszystkich wartości t osiągnięte dla wszystkich ( 2 N ) ! permutacje danych. Interesuje nas to, ponieważ dla dowolnego zamierzonego rozmiaru testu α , takiego jak α = 0,05 (co odpowiada 95 % ufności), zbudujemy dwustronny obszar krytyczny z rozkładu próbkowania t : składa się z najbardziej ekstremalnego$t(x,y)$$x_i$$y_i$$t$$(2n)!$$\alpha$$\alpha = .05$$95$$t$ % możliwych wartości t (po stronie wysokiej, ponieważ wysoka korelacja jest zgodna z podobnymi rozkładami, a niska korelacja nie). W ten sposób ustalamy, jak duży musi być współczynnik korelacji, aby zdecydować, że dane pochodzą z różnych rozkładów.$100\alpha$$t$

Symulowanie zerowego rozkładu próbkowania

Ponieważ (lub, jeśli chcesz, ( 2 $(2n)!$, który zlicza liczbę sposobów podziału2ndanych na dwie części o rozmiarzen), staje się duży nawet dla małychn, nie jest możliwe dokładne obliczenie rozkładu próbkowania, więc próbkujemy go za pomocą symulacji. (Na przykład, gdyn=16, ( 2$\binom{2n}{n}/2$$2n$$n$$n$$n=16$i(2n)! ≈2,63×1035). Często wystarcza około tysiąca próbek (i na pewno będą to poszukiwania, które zamierzamy podjąć).$\binom{2n}{n}/2 = 300\ 540\ 195$$(2n)! \approx 2.63\times 10^{35}$

Musimy się dowiedzieć dwóch rzeczy: po pierwsze, jak wygląda rozkład próbkowania pod hipotezą zerową. Po drugie, jak dobrze ten test rozróżnia różne dystrybucje?

Istnieje komplikacja: rozkład próbkowania zależy od charakteru danych. Wszystko, co możemy zrobić, to spojrzeć na realistyczne dane, stworzone w celu naśladowania tego, czym jesteśmy zainteresowani badaniem, i mamy nadzieję, że to, czego uczymy się z symulacji, będzie miało zastosowanie do naszej sytuacji.

Realizacja

Aby to zilustrować, przeprowadziłem tę pracę w R. Naturalnie dzieli się na trzy części.

1. Funkcja do obliczenia statystyki testowej . Ponieważ chcę być trochę bardziej ogólny, moja wersja obsługuje zestawy danych o różnych rozmiarach ( n x ≠ n y ), interpolując liniowo wartości w (posortowanym) większym zbiorze danych, aby utworzyć dopasowania z (posortowanym) mniejszym zestawem danych. Ponieważ jest to już zrobione przez funkcję , po prostu biorę jej wyniki:$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))
   }

   Mały zwrot - niepotrzebny, ale pomocny w wizualizacji - ponownie wyraża współczynnik korelacji w sposób, który sprawi, że rozkład statystyki zerowej będzie w przybliżeniu symetryczny. Tak się transformdzieje.

2. Symulacja rozkładu próbkowania. Dla danych wejściowych funkcja akceptuje liczbę iteracji n.iterwraz z dwoma zestawami danych w tablicach xi y. Wyprowadza tablicę n.iterwartości statystyki testowej. Jego wewnętrzne działanie powinno być przezroczyste, nawet dla osób niebędących Rużytkownikami:

   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. Chociaż to wszystko, czego potrzebujemy do przeprowadzenia testu, aby go przestudiować, będziemy chcieli powtórzyć test wiele razy. Tak więc przeprowadzamy test raz i owijamy ten kod w trzecią warstwę funkcjonalną, ftutaj ogólnie nazywaną tutaj, którą możemy wielokrotnie wywoływać. Aby uczynić go wystarczająco ogólnym dla szerokiego badania, dla danych wejściowych przyjmuje rozmiary zestawów danych do symulacji ( n.xi n.y), liczbę iteracji dla każdego testu permutacji ( n.iter), odniesienie do funkcji testobliczania statystyki testu (zobaczysz chwilowo dlaczego nie chcieliśmy tego kodować na stałe) i dwie funkcje do generowania iid wartości losowych, jedna dla ( ) i jedna dla Y ( ). Opcja$X$dist.x$Y$dist.yplot.it jest przydatny, aby zobaczyć, co się dzieje.

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

   Dane wyjściowe to symulowana „wartość p”: odsetek symulacji dających statystykę, która wygląda na bardziej ekstremalną niż obliczona dla danych.

Części (2) i (3) są niezwykle ogólne: możesz przeprowadzić badanie takie jak ten dla innego testu, po prostu zastępując test.statisticje innymi obliczeniami. Robimy to poniżej.

Pierwsze wyniki

Domyślnie nasz kod porównuje dane pochodzące z dwóch jednolitych dystrybucji. Pozwalam to zrobić (dla , które są dość małymi zestawami danych i dlatego przedstawiają średnio trudny przypadek testowy), a następnie powtarzam to dla porównania równomiernie-normalnie i równomiernie wykładniczo. (Rozkłady jednolite nie są łatwe do odróżnienia od rozkładów normalnych, chyba że masz nieco więcej niż 16 wartości, ale rozkłady wykładnicze - mające wysoką skośność i długi prawy ogon - zwykle łatwo odróżnić od rozkładów jednolitych.)$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$niezależne wartości od każdego. To dość niska moc. Ale może jest to nieuniknione, więc kontynuujmy.

Wykresy po prawej stronie podobnie testują rozkład równomierny względem wykładniczego. Ten wynik jest dziwny. W tym teście najczęściej stwierdza się, że jednolite dane i dane wykładnicze wyglądają tak samo. Wydaje się „myśleć”, że zmienne jednolite i wykładnicze są bardziej podobne niż dwie zmienne jednolite! Co tu się dzieje?

Problem polega na tym, że dane z rozkładu wykładniczego będą miały zwykle kilka bardzo wysokich wartości. Kiedy utworzysz wykres rozproszenia tych względem równomiernie rozłożonych wartości, w prawym górnym rogu wszystkich pozostałych będzie kilka punktów. Odpowiada to bardzo wysokiemu współczynnikowi korelacji. Tak więc, ilekroć którykolwiek z rozkładów generuje kilka skrajnych wartości, współczynnik korelacji jest okropnym wyborem do pomiaru, jak różne są rozkłady. Prowadzi to do jeszcze jednego, jeszcze gorszego problemu: wraz ze wzrostem rozmiarów zbioru danych zwiększa się szansa na uzyskanie kilku ekstremalnych obserwacji. Dlatego możemy oczekiwać, że ten test będzie działał coraz gorzej wraz ze wzrostem ilości danych. Jak bardzo okropne ...

Lepszy test

$y=x$

Oto Rimplementacja:

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

Zgadza się: jest wbudowany w oprogramowanie, więc musimy tylko to nazwać. Ale poczekaj! Jeśli uważnie przeczytać instrukcję, dowiesz się, że (a) dostawy testowe p-wartość, ale (b) wartość p jest (rażąco) niepoprawne gdy oba xi ysą zbiory danych. Jest przeznaczony do użycia, gdy uważasz, że wiesz dokładnie, z jakiej dystrybucji xpochodzą dane i chcesz sprawdzić, czy to prawda. Zatem test nie uwzględnia w odpowiedni sposób niepewności co do rozkładu, z którego pochodzą dane y.

Nie ma problemu! Struktura testu permutacji jest nadal równie ważna. Dokonując poprzedniej zmiany test.statistic, wszystko, co musimy zrobić, to ponownie uruchomić poprzednie badanie bez zmian. Oto wyniki.

$p=0.20$

$70$$0.05$$11$

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

Wnioski

Tak więc problemy z testem korelacji nie wynikają z pewnych nieodłącznych trudności w tym ustawieniu. Test korelacji nie tylko działa bardzo źle, ale jest zły w porównaniu do powszechnie znanego i dostępnego testu. (Sądzę, że jest niedopuszczalny, co oznacza, że ​​zawsze będzie działał średnio gorzej niż wersja permutacyjna testu KS, co oznacza, że ​​nie ma powodu, aby go używać).

Nie, korelacja nie jest dobrym sprawdzianem tego.

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

Nie znam dobrego testu, który porówna, czy np. Dwie dystrybucje są normalne, ale być może mają inną średnią i sd. Pośrednio, możesz przetestować normalność każdej z nich osobno, a jeśli obie wydają się normalne, zgadnij, że obie były.

Jeśli istnieje wystarczająco duża liczba zmiennych, może to wykazywać większą korelację z wartościami uporządkowanymi według wielkości. Nie wydaje się to jednak szczególnie przydatną metodą, zwłaszcza dlatego, że zapewnia niewiele środków do oszacowania pewności, że mogą skorzystać z tego samego modelu.

Problem, którego możesz doświadczyć, polega na tym, że masz modele o podobnej średniej i skośności, ale z różnicą w kurtozie, ponieważ umiarkowana liczba pomiarów może zmieścić się wystarczająco dobrze, aby wyglądać dość dobrze skorelowana.

Bardziej uzasadnione wydaje się modelowanie obu zmiennych dla różnych rozkładów, aby zobaczyć, która jest najbardziej prawdopodobna dla każdej z nich, i porównać wyniki.

Normalizacja obu wartości, sortowanie i kreślenie każdej z nich może mieć pewne zalety - pozwoli to zobaczyć, jak pasują do siebie - i można również wykreślić możliwy model dla obu, który byłby związany z tym, co zasugerowałeś, ale nie oczekiwanie konkretnej odpowiedzi, tylko wizualny pomysł na bliskość dystrybucji.