Jak tworzymy przedział ufności dla parametru testu permutacji?

Testy permutacyjne to testy istotności oparte na próbkach permutacyjnych losowo pobranych z oryginalnych danych. Próbki permutacji są rysowane bez zamiany, w przeciwieństwie do próbek bootstrap, które są rysowane z zamianą. Oto przykład, który zrobiłem w R prostego testu permutacji. (Twoje komentarze są mile widziane)

Testy permutacyjne mają ogromne zalety. Nie wymagają określonych kształtów populacji, takich jak normalność. Odnoszą się one do różnych statystyk, nie tylko do statystyk, które mają prosty rozkład pod hipotezą zerową. Mogą dawać bardzo dokładne wartości p, niezależnie od kształtu i wielkości populacji (jeśli zastosuje się wystarczającą liczbę permutacji).

Przeczytałem również, że często przydatne jest podanie przedziału ufności wraz z testem, który jest tworzony przy użyciu ponownego próbkowania bootstrap zamiast ponownego próbkowania permutacyjnego.

Czy możesz wyjaśnić (lub po prostu podać kod R), jak konstruowany jest przedział ufności (tj. Dla różnicy między średnimi dwóch próbek w powyższym przykładzie)?

EDYTOWAĆ

Po pewnym googlowaniu znalazłem tę ciekawą lekturę .

Można użyć ponownego próbkowania permutacji. To naprawdę zależy od wielu czynników. Jeśli twoje permutacje są względnie niskie, twoje oszacowanie przedziału ufności nie jest tak dobre w przypadku permutacji. Twoje permutacje znajdują się w nieco szarym obszarze i prawdopodobnie są w porządku.

Jedyna różnica w porównaniu z poprzednim kodem polega na tym, że próbki są generowane losowo zamiast z permutacjami. I wygenerowałbyś ich więcej, powiedzmy na przykład 1000. Uzyskaj wyniki różnic dla 1000 powtórzeń eksperymentu. Weź granice dla środkowego 950 (95%). To twój przedział ufności. Wypada bezpośrednio z paska startowego.

Zrobiłeś już większość tego w swoim przykładzie. różnica leczenia wynosi 462 pozycji. Dlatego potrzebujesz dolnej wartości granicznej 2,5% i górnej 2,5% (około 11 pozycji na każdym końcu).

Używanie kodu sprzed ...

y <- sort(dif.treat)
ci.lo <- y[11]
ci.hi <- y[462-11]

Powiedziałbym, że 462 jest trochę za niski, ale znajdziesz bootstrap do 10.000, który ma trochę inne wyniki (prawdopodobnie bliższe średniej).

Pomyślałem, że dodam również prosty kod wymagający biblioteki rozruchowej (na podstawie twojego poprzedniego kodu).

diff <- function(x,i) mean(x[i[6:11]]) - mean(x[i[1:5]])
b <- boot(total, diff, R = 1000)
boot.ci(b)

Ponieważ test permutacji jest testem dokładnym , dającym dokładną wartość p. Uruchomienie testu permutacji nie ma sensu.

Poza tym określenie przedziału ufności wokół statystyki testowej również nie ma sensu, ponieważ jest obliczane na podstawie próbki, a nie oszacowania. Przedziały ufności ustalane są wokół szacunków takich jak średnie i podobne, ale nie wokół statystyk testowych.

Testów permutacyjnych nie należy używać w zestawach danych, które są tak duże, że nie można już obliczyć wszystkich możliwych permutacji. W takim przypadku użyj procedury ładowania początkowego, aby ustalić wartość graniczną dla używanej statystyki testowej. Ale znowu nie ma to nic wspólnego z 95% przedziałem ufności.

Przykład: używam tutaj klasycznej statystyki T, ale używam prostego podejścia do ładowania początkowego w celu obliczenia rozkładu empirycznego mojej statystyki. Na tej podstawie obliczam empiryczną wartość p:

x <- c(11.4,25.3,29.9,16.5,21.1)
y <- c(23.7,26.6,28.5,14.2,17.9,24.3)

t.sample <- t.test(x,y)$statistic
t.dist <- apply(
      replicate(1000,sample(c(x,y),11,replace=F)),2,
      function(i){t.test(i[1:5],i[6:11])$statistic})

# two sided testing
center <- mean(t.dist)
t.sample <-abs(t.sample-center)
t.dist <- abs(t.dist - center)
p.value <- sum( t.sample < t.dist ) / length(t.dist)
p.value

Weź pod uwagę, że to dwustronne testowanie działa tylko w przypadku rozkładów symetrycznych. Rozkłady niesymetryczne są zwykle testowane tylko jednostronnie.

EDYTOWAĆ :

OK, źle zrozumiałem pytanie. Jeśli chcesz obliczyć przedział ufności dla oszacowania różnicy, możesz użyć kodu wymienionego tutaj do ładowania początkowego w każdej próbce. Pamiętaj, że jest to tendencyjne oszacowanie: generalnie daje to zbyt niski CI. Zobacz także podany tam przykład jako powód, dla którego musisz zastosować inne podejście do przedziału ufności i wartości p.

Z kodu Jorisa Meysa w odpowiedziach, ale z modyfikacją umożliwiającą zastosowanie go w więcej niż jednej sytuacji:

Próbowałem edytować drugi, ale nie miałem czasu na dokończenie iz jakiegoś powodu nie mogę komentować (może dlatego, że jest to stare pytanie).

x <- c(11.4,25.3,29.9,16.5,21.1)
y <- c(23.7,26.6,28.5,14.2,17.9,24.3)

t.sample <- t.test(x,y)$statistic

t.dist <- apply(
          replicate(1000,sample(c(x,y),length(c(x,y)),replace=F)), 2,
          function(i){t.test(i[1:length(x)],i[length(x)+1:length(c(x,y))])$statistic})

# two sided testing
center <- mean(t.dist)
t.sample <-abs(t.sample-center)
t.dist <- abs(t.dist - center)
p.value <- sum( t.sample < t.dist ) / length(t.dist)
p.value