Wartości P równe 0 w teście permutacji

15

Mam dwa zestawy danych i chciałbym wiedzieć, czy są one znacząco różne, czy nie (pochodzi od „ Dwie grupy są znacząco różne? Test do użycia ”).

Zdecydowałem się użyć testu permutacji, wykonując następujące czynności w języku R:

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
}

Niemniej jednak wartości p nie powinny wynosić 0 zgodnie z tym artykułem: http://www.statsci.org/smyth/pubs/permp.pdf

Co polecasz mi zrobić? Czy w ten sposób obliczasz wartość p:

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

dobra droga? A może lepiej wykonać następujące czynności?

pvalue = sum(abs(diff.random) >= abs(diff.observed)) + 1 / number_of_permutations + 1
użytkownik 2886545
źródło
(1) Ostatni wiersz pytania jest błędny, ponieważ nie zawiera nawiasów niezbędnych do wykonania zamierzonego obliczenia. (Gwarantowane jest uzyskanie wyników większych niż , co nie jest możliwe dla żadnej wartości p.) (2) W rzeczywistości nie przeprowadzasz testu permutacji: dwie próbki i rzadko zawierają losową partycję danych, ale zwykle nakładają się w zasadzie. Zamiast tego oblicz jako uzupełnienie unii i . 1a.randomb.randomb.randoma.randomcodinglncrna
whuber
Ponieważ wartość p jest zbiorem wartości co najmniej tak ekstremalnych jak obserwowane, jeśli ocenia się rozkład permutacji, obserwowana statystyka znajduje się w zliczonych „permutacjach”. Podczas przeprowadzania randomizacji często zalicza się obserwowaną statystykę do rozważanych statystyk permutacji (z podobnych powodów).
Glen_b

Odpowiedzi:

15

Dyskusja

Test permutacji generuje wszystkie odpowiednie permutacje zestawu danych, oblicza wyznaczoną statystykę testową dla każdej takiej permutacji i ocenia rzeczywistą statystykę testową w kontekście wynikowego rozkładu permutacji statystyk. Częstym sposobem oceny jest zgłaszanie odsetka statystyk, które są (w pewnym sensie) „bardziej lub bardziej ekstremalne” niż statystyki rzeczywiste. Jest to często nazywane „wartością p”.

Ponieważ rzeczywisty zestaw danych jest jedną z tych permutacji, jego statystyki z pewnością będą należeć do tych znalezionych w rozkładzie permutacji. Dlatego wartość p nigdy nie może wynosić zero.

O ile zestaw danych nie jest bardzo mały (zwykle mniej niż około 20-30 liczb ogółem) lub statystyka testowa ma szczególnie ładną formę matematyczną, nie jest możliwe wygenerowanie wszystkich permutacji. (Przykład generowania wszystkich permutacji pojawia się w teście permutacji w R. ) Dlatego komputerowe implementacje testów permutacji zazwyczaj próbkują z rozkładu permutacji. Robią to, generując niektóre niezależne losowe permutacje i licząc, że wyniki będą reprezentatywną próbką wszystkich permutacji.

Dlatego dowolne liczby (takie jak „wartość p”) uzyskane z takiej próbki są jedynie estymatorami właściwości rozkładu permutacji. Jest całkiem możliwe - i często zdarza się, gdy efekty są duże - że oszacowana wartość p wynosi zero. Nie ma w tym nic złego, ale natychmiast rodzi to zaniedbywane dotąd pytanie, o ile szacunkowa wartość p może różnić się od poprawnej? Ponieważ rozkład próbkowania proporcji (taki jak szacowana wartość p) jest dwumianowy, niepewność tę można rozwiązać za pomocą przedziału ufności dwumianowej .


Architektura

Dobrze skonstruowane wdrożenie będzie ściśle śledzić dyskusję pod każdym względem. Zacząłby się od rutyny obliczania statystyki testowej, ponieważ ta służy do porównania średnich dwóch grup:

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

Napisz inną procedurę, aby wygenerować losową permutację zestawu danych i zastosować statystyki testowe. Interfejs do tego pozwala wywołującemu dostarczyć statystyki testowe jako argument. Porówna pierwsze melementy tablicy (przypuszczalnie, że jest grupą odniesienia) z pozostałymi elementami (grupą „leczenia”).

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

Test permutacji jest przeprowadzany najpierw przez znalezienie statystyki dla rzeczywistych danych (zakładanych tutaj, aby były przechowywane w dwóch tablicach controli treatment), a następnie znalezienie statystyki dla wielu niezależnych losowych permutacji:

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)

Teraz obliczyć dwumianową estymatę wartości p i przedział ufności dla niej. Jedna metoda wykorzystuje wbudowaną binconfprocedurę w HMiscpakiecie:

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

Nie jest złym pomysłem porównywanie wyniku z innym testem, nawet jeśli wiadomo, że nie jest to do końca odpowiednie: przynajmniej możesz uzyskać poczucie wielkości rzędu, w którym wynik powinien się znajdować. W tym przykładzie (porównania średnich) test t-Studenta zazwyczaj daje dobry wynik:

t.test(treatment, control)

Ta architektura jest zilustrowana w bardziej złożonej sytuacji, z działającym Rkodem, w Test, czy zmienne mają tę samą dystrybucję .


Przykład

100201.5

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

Po użyciu poprzedniego kodu do przeprowadzenia testu permutacji narysowałem próbkę rozkładu permutacji wraz z pionową czerwoną linią, aby zaznaczyć rzeczywistą statystykę:

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

Postać

Wynikało z obliczenia dwumianowego limitu ufności

 PointEst Lower        Upper
        0     0 0.0003688199

00.000373.16e-050.000370.000370.050.010.001


Komentarze

kN k/N(k+1)/(N+1)N

10102=1000.0000051.611.7części na milion: nieco mniej niż w teście t-Studenta. Chociaż dane zostały wygenerowane za pomocą normalnych generatorów liczb losowych, co uzasadniałoby zastosowanie testu t-Studenta, wyniki testu permutacji różnią się od wyników testu t-Studenta, ponieważ rozkłady w każdej grupie obserwacji nie są całkowicie normalne.

Whuber
źródło
Cytowany powyżej artykuł Smyth & Phipson wyraźnie pokazuje, dlaczego k / N jest złym wyborem dla estymatora wartości p. W skrócie, dla odpowiednich poziomów istotności, takich jak alfa = 0,05, P ((k / N) <alfa | H0) może być zaskakująco większy niż alfa. Oznacza to, że losowy test permutacji z wykorzystaniem k / N jako estymatora wartości p i 0,05 jako progu odrzucenia odrzuci hipotezę zerową więcej niż 5% razy! Zerowa wartość p jest skrajnym przypadkiem tego problemu - przy kryterium alfa = 0 spodziewamy się, że nigdy nie odrzucimy wartości zerowej, ale b / m może równać się zeru pod zerą, co prowadzi do fałszywego odrzucenia.
Trisoloriansunscreen
1
@Tal „Zły wybór” do określonego celu. To, co nas wyróżnia jako statystów od innych, to nasze rozumienie roli zmienności w analizie danych i podejmowaniu decyzji, a także nasza zdolność do odpowiedniego oszacowania tej zmienności. Takie podejście jest przedstawione (i domyślnie zalecane) w mojej odpowiedzi tutaj. Gdy jest przeprowadzany, nie ma takiego problemu, jaki opisujesz, ponieważ użytkownik procedury permutacji jest zmuszony zrozumieć swoje ograniczenia i mocne strony i będzie miał swobodę działania zgodnie ze swoimi celami.
whuber
13

BMB+1M+1

(B to liczba losowych permutacji, w których uzyskana jest statystyka większa lub równa niż zaobserwowana, a M to całkowita liczba losowych permutacji próbkowanych).

BM

Trisoloriansunscreen
źródło
1
+1 To dobre podsumowanie głównej kwestii artykułu. Szczególnie doceniam waszą uwagę na rozróżnienie między oszacowaną wartością p a rzeczywistą wartością p permutacji.
whuber