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

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:

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.

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

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

Lub nawet paski wykresów obok siebie:

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

$x1$$y1$$x1$$y1$$x1$$y1$$y1$

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?

$z$$p$$p$

$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.

$p$$p$

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

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ł.

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?

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 .

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.