Какую реализацию теста перестановки в R использовать вместо t-тестов (парных и непарных)?

У меня есть данные из эксперимента, которые я проанализировал с помощью t-тестов. Зависимая переменная масштабируется по интервалу, а данные либо непарные (т. Е. 2 ​​группы), либо парные (т. Е. Внутри-субъекты). Например (в рамках предметов):

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)

Тем не менее, данные не являются нормальными, поэтому один рецензент попросил нас использовать что-то кроме t-критерия. Однако, как можно легко видеть, данные не только не распределены нормально, но и неравномерно распределены между условиями:

Поэтому обычные непараметрические критерии, U-критерий Манна-Уитни (непарный) и критерий Уилкоксона (парный), нельзя использовать, поскольку они требуют равного распределения между условиями. Поэтому я решил, что какой-то тест с передискретизацией или перестановкой будет лучше.

Теперь я ищу R-реализацию эквивалента t-критерия на основе перестановок или любой другой совет о том, что делать с данными.

Я знаю, что есть некоторые R-пакеты, которые могут сделать это для меня (например, Coin, Perm, SparkRestTest и т. Д.), Но я не знаю, какой из них выбрать. Так что, если кто-то, имеющий некоторый опыт использования этих тестов, может дать мне толчок, это будет Ubercool.

ОБНОВЛЕНИЕ: было бы идеально, если бы вы могли предоставить пример того, как сообщить о результатах этого теста.

Это не должно иметь большого значения, так как тестовая статистика всегда будет иметь различие в средствах (или что-то эквивалентное). Небольшие различия могут быть связаны с применением методов Монте-Карло. Попытка трех пакетов с вашими данными с односторонним тестом для двух независимых переменных:

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

Чтобы проверить точное значение p с помощью ручного расчета всех перестановок, я ограничу данные первыми 9 значениями.

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

coinи exactRankTestsоба от одного автора, но, coinкажется, являются более общими и обширными - также с точки зрения документации. exactRankTestsне активно развивается больше. Поэтому я бы выбрал coin(также из-за информативных функций, таких как support()), если вы не любите иметь дело с объектами S4.

РЕДАКТИРОВАТЬ: для двух зависимых переменных, синтаксис

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

Несколько комментариев, я полагаю, в порядке.

1) Я бы посоветовал вам попробовать несколько визуальных отображений ваших данных, потому что они могут фиксировать вещи, которые потеряны (например, графики), и я также настоятельно рекомендую рисовать на осях рядом друг с другом. В этом случае я не верю, что гистограммы отлично справляются с передачей основных характеристик ваших данных. Например, посмотрите на боксы бок о бок:

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

Или даже рядом стрип-чарты:

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) Вы не очень подробно объяснили, откуда берутся ваши данные и как они были измерены, но эта информация очень важна, когда приходит время выбирать статистическую процедуру. Ваши два образца выше независимы? Есть ли основания полагать, что предельные распределения двух образцов должны быть одинаковыми (например, за исключением различий в местоположении)? Каковы были соображения до исследования, которые заставили вас искать доказательства различия между двумя группами?

$z$$p$$p$

$p$

5) По моему мнению, эти данные являются идеальным (?) Примером того, что правильно выбранная картина стоит 1000 проверок гипотез. Нам не нужна статистика, чтобы определить разницу между карандашом и сараем. На мой взгляд, подходящим утверждением для этих данных будет «Эти данные демонстрируют заметные различия в отношении местоположения, масштаба и формы». Вы можете использовать (надежную) описательную статистику для каждого из них, чтобы количественно оценить различия и объяснить, что эти различия означают в контексте вашего первоначального исследования.

$p$$p$

Этот ответ намного длиннее, чем я предполагал. Прости за это.

Мои комментарии не о реализации теста перестановки, а о более общих проблемах, поднятых этими данными и обсуждении их, в частности поста Дж. Керна.

Эти два распределения на самом деле выглядят очень похоже на меня, КРОМЕ для группы 0 в Y1, которые сильно отличаются от других наблюдений в этой выборке (следующее наименьшее - около 50 по шкале 0-100), а также всех в X1. Сначала я бы выяснил, было ли что-то другое в этих наблюдениях.

Во-вторых, если предположить, что эти 0 действительно принадлежат анализу, то сказать, что тест на перестановку недействителен, потому что распределения, кажется, отличаются, напрашивается вопрос. Если бы значение NULL было истинным (распределения были идентичными), могли бы вы (с разумной вероятностью) получить распределения, выглядящие такими же разными, как эти два? Отвечая, в этом весь смысл теста, не так ли? Возможно, в этом случае некоторые сочтут ответ очевидным без запуска теста, но с этими небольшими, своеобразными дистрибутивами, я не думаю, что смогу.

Когда этот вопрос снова возник, я могу добавить еще один ответ, вдохновленный недавним сообщением в блоге через R-Bloggers от Robert Kabacoff, автора Quick-R и R in Action, использующего lmPermпакет.

Однако этот метод дает резко контрастирующие (и очень нестабильные) результаты по сравнению с результатами, полученными coinпакетом в ответе @caracakl (p-значение анализа внутри объекта 0.008). Анализ также основан на ответе @ 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)))

производит:

> 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 

Если вы запустите это несколько раз, значения p будут колебаться между ~ .05 и ~ .1.

Хотя это и есть ответ на вопрос, позвольте мне в конце поставить вопрос (я могу перенести его на новый вопрос, если это необходимо):
любые идеи о том, почему этот анализ настолько нестабилен и действительно дает такие расходящиеся значения р анализ монет? Я сделал что-то не так?

Являются ли эти оценки пропорциями? Если это так, то вам определенно не следует использовать гауссовский параметрический тест, и, хотя вы можете использовать непараметрический подход, такой как тест на перестановку или загрузку средств, я бы посоветовал вам получить больше статистической мощности используя подходящий негауссовский параметрический подход. В частности, в любое время, когда вы можете вычислить меру пропорции в интересующей единице (например, участник эксперимента), вы можете и, вероятно, должны использовать модель смешанных эффектов, которая определяет наблюдения с биномиально распределенной ошибкой. Смотри Диксон 2004 .

Просто добавить еще один подход, ezPermв ezупаковке:

> # 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     *

Это , как представляется, соответствует к oneway_testв coinупаковке:

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

Однако обратите внимание, что это не тот пример, предоставленный @caracal . В его примере он включает alternative="greater", следовательно, разницу в значениях р ~0.008против ~0.016.

aovpПакет предлагается в одном из ответов производят подозрительно снижающий р-значение, и работаю подозрительно быстро , даже когда я пытаюсь высокими значениями для Iter, Caи maxIterаргументов:

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

Тем не менее, аргументы, похоже, немного уменьшают вариации значений p от ~.03и ~.1(я получил больший диапазон, чем сообщаемый @Henrik), до 0.03и 0.07.