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

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

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

$10$$0$$20$$1.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")

Wynikało z obliczenia dwumianowego limitu ufności

 PointEst Lower        Upper
        0     0 0.0003688199

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

Komentarze

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

$10$$10^2=100$$0.000005$$1.6$$11.7$częś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.

$\frac{B}{M}$$\frac{B+1}{M+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).

$\frac{B}{M}$