P-значения, равные 0 в тесте перестановки

15

У меня есть два набора данных, и я хотел бы знать, значительно ли они различаются или нет (это происходит из « Две группы существенно различаются? Тест для использования »).

Я решил использовать тест перестановки, выполнив в 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
}

Тем не менее, p-значения не должны быть 0 в соответствии с этим документом: http://www.statsci.org/smyth/pubs/permp.pdf

Что вы мне порекомендуете сделать? Это способ расчета р-значения:

pvalue = sum(abs(diff.random) >= abs(diff.observed)) / number_of_permutations

Хороший путь? Или лучше делать следующее?

pvalue = sum(abs(diff.random) >= abs(diff.observed)) + 1 / number_of_permutations + 1
user2886545
источник
(1) Последняя строка в вопросе ошибочна, поскольку она не содержит скобок, необходимых для выполнения предполагаемого вычисления. (Это гарантировано результатов производят больше , чем , что невозможно для любого значения р.) (2) Вы на самом деле не проводите тест перестановки: два образца и редко будет включать в себя случайное разбиение данных , но , как правило , будет перекрывать по существу. Вместо этого вычислите как дополнение внутри объединения и . 1a.randomb.randomb.randoma.randomcodinglncrna
whuber
Поскольку p-значение - это набор значений, по крайней мере, таких же экстремальных, как и наблюдаемые, если оценивать распределение перестановок, наблюдаемая статистика находится в подсчитанных «перестановках». При выполнении рандомизации обычно учитывают наблюдаемую статистику среди рассматриваемой статистики перестановок (по аналогичным причинам).
Glen_b

Ответы:

15

обсуждение

Тест перестановки генерирует все соответствующие перестановки набора данных, вычисляет назначенную статистику теста для каждой такой перестановки и оценивает фактическую статистику теста в контексте результирующего распределения перестановок статистики. Распространенным способом оценки является отчет о доле статистики, которая (в некотором смысле) «как или более экстремальная», чем фактическая статистика. Это часто называют «р-значением».

Поскольку фактический набор данных является одной из этих перестановок, его статистика обязательно будет среди тех, которые находятся в распределении перестановок. Следовательно, значение p никогда не может быть нулевым.

Если набор данных не очень мал (обычно менее 20-30 общих чисел) или тестовая статистика имеет особенно хорошую математическую форму, практически невозможно создать все перестановки. (Пример, где все перестановки генерируются, появляется в Тесте перестановки в R ). Поэтому компьютерные реализации тестов перестановки обычно выбирают из распределения перестановок. Они делают это путем генерации некоторых независимых случайных перестановок и надеются, что результаты являются репрезентативной выборкой всех перестановок.

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


Архитектура

Хорошо продуманная реализация будет внимательно следить за обсуждением во всех отношениях. Он будет начинаться с подпрограммы для вычисления статистики теста, так как эта будет сравнивать средние значения двух групп:

diff.means <- function(control, treatment) mean(treatment) - mean(control)

Напишите другую подпрограмму для генерации случайной перестановки набора данных и примените тестовую статистику. Интерфейс этого позволяет вызывающей стороне предоставлять статистику теста в качестве аргумента. Он будет сравнивать первые mэлементы массива (предположительно эталонную группу) с остальными элементами (группа «обработка»).

f <- function(..., sample, m, statistic) {
  s <- sample(sample)
  statistic(s[1:m], s[-(1:m)])
}

Тест перестановки выполняется сначала путем нахождения статистики для фактических данных (предполагается, что они здесь сохраняются в двух массивах controlи treatment), а затем нахождения статистики для многих независимых случайных перестановок:

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)

Теперь вычислите биномиальную оценку значения p и доверительный интервал для него. Один метод использует встроенную binconfпроцедуру в HMiscпакете:

require(Hmisc)                                    # Exports `binconf`
k <- sum(abs(sim) >= abs(z))                      # Two-tailed test
zapsmall(binconf(k, length(sim), method='exact')) # 95% CI by default

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

t.test(treatment, control)

Эта архитектура проиллюстрирована в более сложной ситуации с рабочим Rкодом в разделе «Тестируют ли переменные одно и то же распределение» .


пример

100201,5

set.seed(17)
control <- rnorm(10)
treatment <- rnorm(20, 1.5)

После использования предыдущего кода для запуска теста перестановки я нанес на график образец распределения перестановок вместе с вертикальной красной линией, чтобы отметить фактическую статистику:

h <- hist(c(z, sim), plot=FALSE)
hist(sim, breaks=h$breaks)
abline(v = stat(control, treatment), col="Red")

фигура

Биноминальный расчет доверительного интервала привел к

 PointEst Lower        Upper
        0     0 0.0003688199

00,000373.16e-050,000370,000370,050,010,001


Комментарии

КN К/N(К+1)/(N+1)N

10102знак равно1000.0000051,611,7частей на миллион: немного меньше, чем сообщил t-критерий Стьюдента. Хотя данные были получены с помощью обычных генераторов случайных чисел, что оправдывало бы использование t-критерия Стьюдента, результаты теста на перестановку отличаются от результатов t-критерия Стьюдента, поскольку распределения в каждой группе наблюдений не совсем нормальные.

Whuber
источник
Приведенная выше работа Smyth & Phipson ясно показывает, почему k / N является плохим выбором для оценки p-значения. В двух словах, для соответствующих уровней значимости, таких как альфа = 0,05, P ((k / N) <альфа | H0) может быть удивительно больше, чем альфа. Это означает, что тест случайной перестановки, использующий k / N в качестве оценки p-значения и 0,05 в качестве порога отклонения, отклонит нулевую гипотезу более чем в 5% случаев! Нулевое p-значение является крайним случаем этой проблемы - с критерием альфа = 0 мы ожидаем, что никогда не отклоним ноль, но b / m может равняться нулю при нуле, что приводит к ложному отклонению.
Trisoloriansunscreen
1
@ Тал "Плохой выбор" для определенной цели. Что отличает нас как статистиков от других, так это наше понимание роли изменчивости в анализе данных и принятии решений, а также наша способность количественно определять эту изменчивость. Этот подход проиллюстрирован (и неявно защищен) в моем ответе здесь. Когда это выполняется, такой проблемы, как вы описываете, нет, потому что пользователь процедуры перестановки должен понимать ее ограничения и сильные стороны и будет иметь свободу действий в соответствии со своими целями.
whuber
13

ВMВ+1M+1

(B - количество случайных перестановок, в которых получена статистика, большая или равная наблюдаемой, а M - общее количество выборок случайных перестановок).

ВM

Trisoloriansunscreen
источник
1
+1 Это хорошее резюме основного пункта статьи. Я особенно ценю ваше внимание на различие между оцененным p-значением и истинным p-значением перестановки.
whuber