Jakiej implementacji testu permutacji w R użyć zamiast testów t (sparowanych i niesparowanych)?

56

Mam dane z eksperymentu, który przeanalizowałem za pomocą testów t. Zmienna zależna jest skalowana w odstępach czasu, a dane są niesparowane (tj. 2 grupy) lub sparowane (tj. W obrębie osobników). Np. (W ramach przedmiotów):

x1 <- c(99, 99.5, 65, 100, 99, 99.5, 99, 99.5, 99.5, 57, 100, 99.5, 
        99.5, 99, 99, 99.5, 89.5, 99.5, 100, 99.5)
y1 <- c(99, 99.5, 99.5, 0, 50, 100, 99.5, 99.5, 0, 99.5, 99.5, 90, 
        80, 0, 99, 0, 74.5, 0, 100, 49.5)

Jednak dane nie są normalne, więc jeden z recenzentów poprosił nas o użycie czegoś innego niż test t. Jednak, jak łatwo zauważyć, dane nie tylko nie są normalnie dystrybuowane, ale rozkłady nie są równe między warunkami: alternatywny tekst

Dlatego zwykłe testy nieparametryczne, test U Manna-Whitneya (niesparowany) i test Wilcoxona (sparowany) nie mogą być stosowane, ponieważ wymagają równych rozkładów między warunkami. Dlatego zdecydowałem, że najlepszy będzie test ponownego próbkowania lub permutacji.

Teraz szukam implementacji R ekwiwalentu testu t opartego na permutacji lub jakiejkolwiek innej porady, co zrobić z danymi.

Wiem, że istnieją pewne pakiety R, które mogą to dla mnie zrobić (np. Moneta, perm, exactRankTest itp.), Ale nie wiem, który wybrać. Tak więc, jeśli ktoś z pewnym doświadczeniem w korzystaniu z tych testów mógłby dać mi szansę na start, to byłby to ubercool.

AKTUALIZACJA: Byłoby idealnie, gdyby można podać przykład zgłaszania wyników z tego testu.

Henrik
źródło

Odpowiedzi:

43

Nie powinno to mieć większego znaczenia, ponieważ statystyki testowe zawsze będą różnicą w środkach (lub czymś równoważnym). Niewielkie różnice mogą wynikać z wdrożenia metod Monte-Carlo. Wypróbowanie trzech pakietów z danymi za pomocą jednostronnego testu dla dwóch niezależnych zmiennych:

DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
library(coin)                    # for oneway_test(), pvalue()
pvalue(oneway_test(DV ~ IV, alternative="greater", 
                   distribution=approximate(B=9999)))
[1] 0.00330033

library(perm)                    # for permTS()
permTS(DV ~ IV, alternative="greater", method="exact.mc", 
       control=permControl(nmc=10^4-1))$p.value
[1] 0.003

library(exactRankTests)          # for perm.test()
perm.test(DV ~ IV, paired=FALSE, alternative="greater", exact=TRUE)$p.value
[1] 0.003171822

Aby sprawdzić dokładną wartość p za pomocą ręcznego obliczenia wszystkich permutacji, ograniczę dane do pierwszych 9 wartości.

x1 <- x1[1:9]
y1 <- y1[1:9]
DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
pvalue(oneway_test(DV ~ IV, alternative="greater", distribution="exact"))
[1] 0.0945907

permTS(DV ~ IV, alternative="greater", exact=TRUE)$p.value
[1] 0.0945907

# perm.test() gives different result due to rounding of input values
perm.test(DV ~ IV, paired=FALSE, alternative="greater", exact=TRUE)$p.value
[1] 0.1029412

# manual exact permutation test
idx  <- seq(along=DV)                 # indices to permute
idxA <- combn(idx, length(x1))        # all possibilities for different groups

# function to calculate difference in group means given index vector for group A
getDiffM <- function(x) { mean(DV[x]) - mean(DV[!(idx %in% x)]) }
resDM    <- apply(idxA, 2, getDiffM)  # difference in means for all permutations
diffM    <- mean(x1) - mean(y1)       # empirical differencen in group means

# p-value: proportion of group means at least as extreme as observed one
(pVal <- sum(resDM >= diffM) / length(resDM))
[1] 0.0945907

coini exactRankTestsoba pochodzą od tego samego autora, ale coinwydaje się być bardziej ogólne i obszerne - także pod względem dokumentacji. exactRankTestsnie jest już aktywnie rozwijany. Dlatego wybrałbym coin(także z powodu funkcji informacyjnych, takich jak support()), chyba że nie lubisz zajmować się obiektami S4.

EDYCJA: dla dwóch zmiennych zależnych składnia to

id <- factor(rep(1:length(x1), 2))    # factor for participant
pvalue(oneway_test(DV ~ IV | id, alternative="greater",
                   distribution=approximate(B=9999)))
[1] 0.00810081
karakal
źródło
Dzięki za wspaniałą odpowiedź! 2 kolejne pytania: czy twój drugi przykład oznacza, że ​​ta moneta faktycznie zapewnia wszystkie możliwe kombinacje i jest dokładnym testem? Czy jest jakaś korzyść z nieprzeprowadzenia dokładnego testu w moim przypadku?
Henrik
10
(+1) Nie jest zaskoczeniem, że (niesparowany) test t daje zasadniczo tę samą wartość p, 0,000349. Pomimo tego, co powiedział recenzent, test t ma zastosowanie do tych danych. Powodem jest to, że rozkłady próbkowania średnich są w przybliżeniu normalne, nawet jeśli rozkłady danych nie są. Co więcej, jak widać z wyników, test t jest w rzeczywistości bardziej konserwatywny niż test permutacji. (Oznacza to, że znaczący wynik testu t wskazuje, że test permutacji również może być znaczący.)
whuber
2
@Henrik W niektórych sytuacjach (wybrany test i złożoność numeryczna) coinmoże rzeczywiście obliczyć dokładny rozkład permutacji (bez faktycznego przechodzenia przez wszystkie permutacje, istnieją bardziej eleganckie algorytmy niż to). Biorąc pod uwagę wybór, dokładny rozkład wydaje się preferowany, ale różnica w stosunku do przybliżenia Monte-Carlo przy dużej liczbie powtórzeń powinna być niewielka.
caracal
1
@ Caracal Dzięki za wyjaśnienie. Pozostaje jedno pytanie: dane, które przedstawiłem, są sparowane. Dlatego potrzebuję odpowiednika sparowanego testu t. Czy oneway_testdokładna funkcja? A jeśli tak, to który z nich jest odpowiedni dla danych niesparowanych?
Henrik,
2
@Henrik coinAutor napisał do mnie, że oneway_test()nie mogę obliczyć dokładnego rozkładu przypadku zależnego, musisz przejść do przybliżenia MC ( wilcoxsign_test()nadaje się tylko do dokładnego testu). Nie wiedziałem o tym i wolałbym w takim przypadku błąd, ale MC powinien być wystarczająco dokładny przy dużej liczbie powtórzeń.
caracal
29

Uważam, że kilka komentarzy jest w porządku.

1) Zachęcam do wypróbowania wielu wizualnych prezentacji danych, ponieważ mogą one uchwycić rzeczy, które zostały utracone przez histogramy (takie jak wykresy), a także zdecydowanie zalecam tworzenie wykresów na osiach obok siebie. W tym przypadku nie wydaje mi się, aby histogramy dobrze radziły sobie z przekazywaniem istotnych cech danych. Na przykład spójrz na równoległe wykresy pudełkowe:

boxplot(x1, y1, names = c("x1", "y1"))

alternatywny tekst

Lub nawet paski wykresów obok siebie:

stripchart(c(x1,y1) ~ rep(1:2, each = 20), method = "jitter", group.names = c("x1","y1"), xlab = "")

alternatywny tekst

x1y1x1y1x1y1y1

2) Nie wyjaśniłeś zbyt szczegółowo, skąd pochodzą twoje dane ani w jaki sposób zostały zmierzone, ale ta informacja jest bardzo ważna, jeśli chodzi o wybór procedury statystycznej. Czy twoje dwie próbki są powyżej niezależne? Czy istnieją powody, by sądzić, że rozkład krańcowy dwóch próbek powinien być taki sam (z wyjątkiem na przykład różnicy w lokalizacji)? Jakie były rozważania przed badaniem, które doprowadziły cię do poszukiwania dowodów na różnicę między tymi dwiema grupami?

zpp

p

5) Moim zdaniem dane te są doskonałym (?) Przykładem, że dobrze wybrany obraz jest wart 1000 testów hipotez. Nie potrzebujemy statystyk, aby odróżnić ołówek od stodoły. Moim zdaniem stosowne oświadczenie dla tych danych brzmiałoby: „Dane te wykazują wyraźne różnice w odniesieniu do lokalizacji, skali i kształtu”. Możesz śledzić (solidne) statystyki opisowe dla każdej z nich, aby obliczyć różnice i wyjaśnić, co oznaczają różnice w kontekście oryginalnego badania.

pp

Ta odpowiedź jest znacznie dłuższa niż pierwotnie zamierzałem. Przepraszam za to.


źródło
Jestem ciekawy, czy weźmiesz pod uwagę następujące odpowiednie quasi-wizualizowane podejście: szacunki bootstrap dla momentów dwóch grup (średnie, wariancje i wyższe momenty, jeśli chcesz), a następnie wykreśl te szacunki i ich przedziały ufności, patrząc dla stopnia nakładania się grup w każdej chwili. Pozwala to mówić o potencjalnych różnicach w różnych charakterystykach dystrybucji. Jeśli dane są sparowane, oblicz wyniki różnic i rozpocznij od początku momenty tej pojedynczej dystrybucji. Myśli?
Mike Lawrence
2
(+1) Dobra analiza. Masz całkowitą rację, że wyniki są oczywiste i nie trzeba naciskać punktu z wartością p. Możesz być trochę ekstremalny w swoim stwierdzeniu (3), ponieważ test t nie wymaga normalnie rozproszonych danych. Jeśli jesteś zaniepokojony, istnieją korekty skośności (np. Wariant Chena): możesz zobaczyć, czy wartość p dla skorygowanego testu zmienia odpowiedź. Jeśli nie, prawdopodobnie wszystko w porządku. W tej konkretnej sytuacji, przy tych (mocno wypaczonych) danych, test t działa dobrze.
whuber
(+1) Niezły chwyt! I bardzo dobre komentarze.
chl
Wydaje się, że akceptujemy pogląd, że podstawowy rozkład jest „podobny” do losowej instancji. Więc nie można postawić pytania: czy oba pochodzą z wersji beta (0,25; 0,25), a następnie sprawdzianem byłby, czy mają one ten sam parametr (nie) centralności. I czy nie uzasadniałoby to ani testu permutacji, ani Wilcoxona?
DW
4
5

Moje komentarze nie dotyczą wdrożenia testu permutacji, ale bardziej ogólnych kwestii poruszonych przez te dane i ich dyskusji, w szczególności posta G. Jaya Kernsa.

Te dwa rozkłady faktycznie wyglądają dość podobnie do mnie Z WYJĄTKIEM dla grupy zer w Y1, które znacznie różnią się od innych obserwacji w tej próbce (następna najmniejsza to około 50 w skali 0-100), a także wszystkie w X1. Najpierw zbadam, czy w tych obserwacjach było coś innego.

Po drugie, zakładając, że te 0 należą do analizy, stwierdzenie, że test permutacji nie jest prawidłowy, ponieważ rozkłady wydają się różne, rodzi pytanie. Jeśli zerowa wartość jest prawdą (rozkłady są identyczne), czy (z rozsądnym prawdopodobieństwem) można uzyskać rozkłady wyglądające tak odmiennie jak te dwa? Odpowiedź to sedno testu, prawda? Może w tym przypadku niektórzy uznają odpowiedź za oczywistą bez uruchamiania testu, ale przy tych niewielkich, osobliwych dystrybucjach nie sądzę, żebym to zrobił.

Andrew Taylor
źródło
Wygląda na to, że powinien to być jeden lub więcej komentarzy, a nie odpowiedź. Jeśli klikniesz mały szary „dodaj komentarz”, możesz umieścić swoje myśli w rozmowie pod pytaniem lub określoną odpowiedzią, tam gdzie one należą. Robisz tutaj istotne uwagi, ale nie jest jasne, czy nie jest to dla nich odpowiednie miejsce.
Gung - Przywróć Monikę
1
@ gung Potrzeba trochę reputacji, aby móc dodać komentarz ;-).
whuber
4
2(1/2)7.01
4

Ponieważ pytanie to pojawiło się ponownie, mogę dodać kolejną odpowiedź zainspirowaną najnowszym postem na blogu autorstwa R-Bloggerów autorstwa Roberta Kabacoffa, autora Quick-R i R w akcji korzystającego z lmPermpakietu.

Jednak metody te dają wyraźnie kontrastujące (i bardzo niestabilne) wyniki z wynikami otrzymanymi przez coinpakiet w odpowiedzi na @ caracakl (wartość p analizy wewnątrz badanych wynosi 0.008). Analiza bierze również przygotowanie danych z odpowiedzi @ caracal:

x1 <- c(99, 99.5, 65, 100, 99, 99.5, 99, 99.5, 99.5, 57, 100, 99.5, 
        99.5, 99, 99, 99.5, 89.5, 99.5, 100, 99.5)
y1 <- c(99, 99.5, 99.5, 0, 50, 100, 99.5, 99.5, 0, 99.5, 99.5, 90, 
        80, 0, 99, 0, 74.5, 0, 100, 49.5)

DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
id <- factor(rep(1:length(x1), 2)) 

library(lmPerm)

summary(aovp( DV ~ IV + Error(id)))

produkuje:

> summary(aovp( DV ~ IV + Error(id)))
[1] "Settings:  unique SS "

Error: id
Component 1 :
          Df R Sum Sq R Mean Sq
Residuals 19    15946       839


Error: Within
Component 1 :
          Df R Sum Sq R Mean Sq Iter Pr(Prob)  
IV         1     7924      7924 1004    0.091 .
Residuals 19    21124      1112                
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 

Jeśli uruchomisz to wiele razy, wartości p będą się zmieniać między ~ 0,05 a ~ .1.

Chociaż jest to odpowiedź na pytanie, pozwól mi zadać pytanie na końcu (w razie potrzeby mogę przenieść je do nowego pytania):
Wszelkie pomysły na to, dlaczego ta analiza jest tak niestabilna i daje tak różne wartości p, aby analiza monet? Czy zrobiłem coś złego?

Henrik
źródło
2
Lepiej zadać to jako osobne pytanie, jeśli naprawdę jest to pytanie, na które chcesz uzyskać odpowiedź, niż inne możliwe rozwiązanie, które chcesz wymienić z resztą. Zauważam, że określasz warstwy błędów, ale @caracal nie; to byłby mój pierwszy przypuszczenie o różnicy między tym wyjściem a jego. Ponadto podczas symulacji wartości zwykle przeskakują; dla odtwarzalności podajesz ziarno, np set.seed(1).; dla większej precyzji w oszacowaniu MC zwiększa się liczbę iteracji; Nie jestem pewien, czy którekolwiek z nich są „właściwą” odpowiedzią na twoje pytanie, ale prawdopodobnie są trafne.
gung - Przywróć Monikę
2
Ponownie sugeruję porównanie wyników MC z obliczeniami ręcznymi przy użyciu testu pełnej permutacji (ponownej randomizacji). Zobacz kod dla swojego przykładu, aby porównać oneway_anova()(zawsze blisko poprawnego wyniku) i aovp()(zwykle daleko od poprawnego wyniku). Nie wiem, dlaczego aovp()daje tak różne wyniki, ale przynajmniej w tym przypadku są nieprawdopodobne. @ ostatnie ostatnie wezwanie do oneway_test(DV ~ IV | id, ...)w mojej oryginalnej odpowiedzi określiło warstwy błędów dla zależnej wielkości liter (inna składnia niż używana przez aov()).
caracal
@caracal, masz rację. Po edycji nie patrzyłem na ostatni blok kodu. Patrzyłem na górny blok kodu - z mojej strony niechlujstwa.
gung - Przywróć Monikę
Naprawdę nie potrzebuję odpowiedzi. Jest to kolejna możliwość, o której warto tutaj wspomnieć. Niestety jest daleki od innych wyników, na które również warto zwrócić uwagę.
Henrik
1
@Henrik uruchom aovp z maxExact = 1000. Jeśli trwa to zbyt długo, ustaw iter = 1000000 i Ca = 0,001. Obliczenia kończą się, gdy szacowany błąd standardowy p jest mniejszy niż Ca * p. (Niższe wartości dają bardziej stabilne wyniki.)
xmjx,
1

Czy te proporcje wyników? Jeśli tak, to na pewno nie powinieneś używać gaussowskiego testu parametrycznego i chociaż możesz zastosować podejście nieparametryczne, takie jak test permutacji lub bootstrap środków, sugerowałbym, aby uzyskać większą moc statystyczną poprzez zastosowanie odpowiedniego niegaussowskiego podejścia parametrycznego. W szczególności, za każdym razem, gdy możesz obliczyć miarę proporcji w jednostce zainteresowania (np. Uczestnik eksperymentu), możesz i prawdopodobnie powinieneś użyć modelu efektów mieszanych, który określa obserwacje z błędem o rozkładzie dwumianowym. Zobacz Dixon 2004 .

Mike Lawrence
źródło
Oceny nie są proporcjami, ale szacunkami uczestników w skali od 0 do 100 (przedstawione dane są średnimi szacunkami dla kilku pozycji z tą skalą).
Henrik
Wtedy nieparametryczna wydaje się tradycyjna droga. To powiedziawszy, zastanawiałem się, czy takie dane w skali mogłyby być pożyteczne, by wywnioskować z dwumianowego procesu i tym samym przeanalizowane jako takie. Oznacza to, że mówisz, że każdy wynik jest średnią z kilku pozycji, i powiedzmy, że każdy element jest skalą 10-punktową, w którym to przypadku reprezentowałbym odpowiedź, powiedzmy, „8” jako serię prób, 8 z które mają wartość 1, a dwie z nich mają wartość 0, wszystkie oznaczone tą samą etykietą w zmiennej „item”. Dzięki tym rozszerzonym / dwumianowym danym możesz następnie obliczyć dwumianowy model efektów mieszanych.
Mike Lawrence
W związku z moim poprzednim komentarzem należy zauważyć, że w danych rozwiniętych / dwumianowych można modelować zmienną „item” jako efekt stały lub losowy. Myślę, że skłaniam się ku modelowaniu tego jako ustalonego efektu, ponieważ prawdopodobnie możesz być zainteresowany nie tylko rozliczaniem, ale także oceną różnic pozycji i wszelkich możliwych interakcji między pozycją a innymi zmiennymi predykcyjnymi.
Mike Lawrence
0

Wystarczy dodanie innego podejścia, ezPermz ezpakietu:

> # preparing the data
> DV <- c(x1, y1)
> IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
> id <- factor(rep(1:length(x1), 2))
> df <- data.frame(id=id,DV=DV,IV=IV)
>
> library(ez)
> ezPerm( data = df, dv = DV, wid = id, within = IV, perms = 1000)
|=========================|100%              Completed after 17 s 
  Effect     p p<.05
1     IV 0.016     *

Wydaje się to spójne oneway_testz coinpakietem:

> library(coin)
> pvalue(oneway_test(DV ~ IV | id,  distribution=approximate(B=999999)))
[1] 0.01608002
99 percent confidence interval:
 0.01575782 0.01640682

Zauważ jednak, że nie jest to ten sam przykład podany przez @caracal . W swoim przykładzie obejmuje alternative="greater"zatem różnicę wartości p ~0.008względem ~0.016.

aovpPakiet zasugerował w jednym z odpowiedziami produkować podejrzanie niższe wartości p i uruchamia podejrzanie szybko nawet przy próbie wysokie wartości dla Iter, Cai maxIterargumenty:

library(lmPerm)
summary(aovp(DV ~ IV + Error(id/IV), data=df,  maxIter = 1000000000))
summary(aovp(DV ~ IV + Error(id/IV), data=df,  Iter = 1000000000))
summary(aovp(DV ~ IV + Error(id/IV), data=df,  Ca = 0.00000000001))

To powiedziawszy, argumenty wydają się nieznacznie zmniejszać warianty wartości p od ~.03i ~.1(dostałem większy zakres niż raportowany przez @Henrik), do 0.03i 0.07.

toto_tico
źródło