Вычисление p-значения с помощью начальной загрузки с R

28

Я использую пакет «boot» для вычисления приблизительного 2-стороннего загрузочного p-значения, но результат слишком далек от p-значения при использовании t.test. Я не могу понять, что я сделал неправильно в моем коде R. Может кто-нибудь, пожалуйста, дайте мне подсказку для этого

time = c(14,18,11,13,18,17,21,9,16,17,14,15,
         12,12,14,13,6,18,14,16,10,7,15,10)
group=c(rep(1:2, each=12))
sleep = data.frame(time, group)

require(boot)
diff = function(d1,i){
    d = d1[i,]
    Mean= tapply(X=d$time, INDEX=d$group, mean)
    Diff = Mean[1]-Mean[2]
    Diff
}

set.seed(1234)
b3 = boot(data = sleep, statistic = diff, R = 5000, strata=sleep$group)

pvalue = mean(abs(b3$t) > abs(b3$t0))
pvalue 

Двустороннее загрузочное p-значение (pvalue) = 0,4804, но двустороннее p-значение t.test составляет 0,04342. Оба значения р примерно в 11 раз больше. Как это может случиться?

Tu.2
источник
почему у b3 $ t0 есть две записи?
Сиань
1
это имя!
Элвис
2
Вы неправильно вычисляете значение. В документации говорится, что t 0 - наблюдаемая статистика, а не нулевое распределение, как можно предположить в обозначениях. Вам необходимо составить оценку выборки dist-n под нулевым значением. Смотрите мой ответ для получения дополнительной информации. Попробуйте для исправления несмещенного теста. пT0mean(abs(b3$t0) < abs(b3$t-mean(b3$t)))
AdamO

Ответы:

31

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

> quantile(b3$t,c(0.025,0.975))
     2.5%     97.5% 
0.4166667 5.5833333 

Чтобы получить значение, вам нужно сгенерировать перестановки по нулевой гипотезе. Это можно сделать, например, так:п

diff2 = function(d1,i){
    d = d1; 
    d$group <- d$group[i];  # randomly re-assign groups
    Mean= tapply(X=d$time, INDEX=d$group, mean)
    Diff = Mean[1]-Mean[2]
    Diff
}

> set.seed(1234)
> b4 = boot(data = sleep, statistic = diff2, R = 5000)
> mean(abs(b4$t) > abs(b4$t0))
[1] 0.046

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

> R <- 10000; d <- sleep
> b5 <- numeric(R); for(i in 1:R) { 
+    d$group <- sample(d$group, length(d$group)); 
+    b5[i] <- mean(d$time[d$group==1])-mean(d$time[d$group==2]); 
+ }
> mean(abs(b5) > 3)
[1] 0.0372
Элвис
источник
5
Технически это тест перестановки, а не загрузочное p-значение.
AdamO
@AdamO Я согласен, что в этом ответе представлен тест перестановки (и его слегка измененный вариант); это потому, что во время повторной выборки группы объединяются. Напротив, в тесте на основе начальной загрузки значения для каждой группы должны выбираться с использованием только данных для этой же группы. Вот один ответ, объясняющий, как это сделать: stats.stackexchange.com/a/187630/28666 .
говорит амеба: восстанови
@amoeba Я думаю, что ответ, на который вы ссылаетесь, также является тестом перестановки, связанным с начальной загрузкой только постольку, поскольку он включает повторную выборку. Совершенно нормально сообщать, но для составления отчета используются два метода. Непараметрический загрузчик технически неспособен генерировать данные при нулевой гипотезе. Смотрите мой ответ для некоторых из того, как p-значения генерируются из дистрибутива начальной загрузки.
AdamO
@ AdamO Я полагаю, что это вопрос терминологии, но я не понимаю, как процедуру, описанную в связанном ответе, можно назвать тестом «перестановки», потому что там ничего не переставляется: пересчитанные значения для каждой группы генерируются с использованием данных из этого только группа.
говорит амеба, восстанови Монику
1
Элвис, я думаю, что первый фрагмент кода в твоем ответе тоже тест на перестановку. При повторной выборке вы объединяете группы вместе! Это то, что определяет тест перестановки.
говорит амеба, восстанови Монику
25

Ответ Элвиса основан на перестановках, но, по моему мнению, он не дает четкого представления о том, что не так с оригинальным подходом начальной загрузки. Позвольте мне обсудить решение, основанное исключительно на начальной загрузке.

Основная проблема вашей оригинальной симуляции заключается в том, что bootstrap всегда предоставляет вам ИСТИННОЕ распределение статистики теста. Однако при вычислении p-значения необходимо сравнивать полученное значение тестовой статистики с ее распределением ПОД H0, т.е. не с истинным распределением!

[Давайте сделаем это ясно. Например, известно, что тестовая статистика T классического t-теста имеет классическое «центральное» t-распределение при H0 и нецентральное распределение в целом. Однако всем знаком тот факт, что наблюдаемое значение T сравнивается с классическим «центральным» t-распределением, то есть никто не пытается получить истинное [нецентральное] t-распределение для сравнения с T.]

Ваше значение p 0,4804 настолько велико, потому что наблюдаемое значение «t0» тестовой статистики Среднее [1] -Среднее [2] лежит очень близко к центру начальной загрузки образца «t». Это естественно и, как правило, так всегда [т.е. независимо от действительности H0], потому что загруженный образец «t» эмулирует АКТУАЛЬНОЕ распределение Mean [1] -Mean [2]. Но, как отмечено выше [и также Элвисом], вам действительно нужно распределение Mean [1] -Mean [2] ПОД H0. Очевидно, что

1) при H0 распределение Mean [1] -Mean [2] будет сосредоточено вокруг 0,

2) его форма не зависит от справедливости H0.

Эти две точки означают, что распределение Mean [1] -Mean [2] в H0 может эмулироваться загрузочным образцом «t» SHIFTED, так что он центрируется около 0. В R:

b3.under.H0 <- b3$t - mean(b3$t)

и соответствующее значение p будет:

mean(abs(b3.under.H0) > abs(b3$t0))

что дает вам «очень хорошее» значение 0,0232. :-)

Позвольте мне заметить, что упомянутая выше точка «2) называется« переводной эквивалентностью »тестовой статистики, и она НЕ должна выполняться в целом! Т.е. для некоторых тестовых статистик смещение начальной загрузки "t" не дает вам достоверной оценки распределения тестовой статистики в соответствии с HO! Посмотрите на эту дискуссию и особенно на ответ П. Далгаарда: http://tolstoy.newcastle.edu.au/R/e6/help/09/04/11096.html

Ваша задача тестирования дает абсолютно симметричное распределение статистики теста, но имейте в виду, что существуют некоторые проблемы с получением ДВУСТОРОННИХ p-значений в случае перекошенного начального распределения статистики теста. Опять же, прочитайте приведенную выше ссылку.

[И, наконец, я бы использовал «чистый» тест перестановки в вашей ситуации; т.е. вторая половина Элвиса ответит. :-)]

Ян с.
источник
17

Существует множество способов расчета CI начальной загрузки и p-значений. Основная проблема заключается в том, что загрузчик не может генерировать данные при нулевой гипотезе. Тест на перестановку является жизнеспособной альтернативой этому. Чтобы использовать правильный загрузчик, вы должны сделать некоторые предположения о распределении выборки статистики теста.

β0*знак равноβ^-β^*β0*знак равноβ^*-β^

нормальный бутстрап

Одним из подходов является нормальная начальная загрузка, где вы берете среднее и стандартное отклонение распределения начальной загрузки, вычисляете распределение выборки под нулевым значением, сдвигая распределение и используя нормальные процентили из нулевого распределения в точке оценки в исходной выборке начальной загрузки. , Это разумный подход, когда распределение начальной загрузки нормальное, обычно здесь достаточно визуального осмотра. Результаты, использующие этот подход, обычно очень близки к надежной или основанной на сэндвичах оценке ошибок, которая устойчива к предположениям о гетероскедастичности и / или конечной выборочной дисперсии. Допущение нормальной тестовой статистики является более сильным условием допущений в следующем тесте начальной загрузки, который я буду обсуждать.

процентиль бутстрап

F0*) вычитая распределение начальной загрузки из оценок, а затем вычисляя, какой процент DSNH является «более экстремальным», чем ваша оценка, используя 2×мин(F0*(β^),1-F0*(β^))

Студенческий бутстрап

п

Пример программирования

В качестве примера я буду использовать cityданные в пакете начальной загрузки. Доверительные интервалы начальной загрузки рассчитываются с помощью этого кода:

ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)
city.boot <- boot(city, ratio, R = 999, stype = "w", sim = "ordinary")
boot.ci(city.boot, conf = c(0.90, 0.95),
        type = c("norm", "basic", "perc", "bca"))

и произвести этот вывод:

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates

CALL : 
boot.ci(boot.out = city.boot, conf = c(0.9, 0.95), type = c("norm", 
    "basic", "perc", "bca"))

Intervals : 
Level      Normal              Basic         
90%   ( 1.111,  1.837 )   ( 1.030,  1.750 )   
95%   ( 1.042,  1.906 )   ( 0.895,  1.790 )  

Level     Percentile            BCa          
90%   ( 1.291,  2.011 )   ( 1.292,  2.023 )   
95%   ( 1.251,  2.146 )   ( 1.255,  2.155 )  
Calculations and Intervals on Original Scale

95% -й доверительный интервал для нормальной начальной загрузки получается путем вычисления:

with(city.boot, 2*t0 - mean(t) + qnorm(c(0.025, 0.975)) %o% sqrt(var(t)[1,1]))

Таким образом, p-значение получается:

> with(city.boot, pnorm(abs((2*t0 - mean(t) - 1) / sqrt(var(t)[1,1])), lower.tail=F)*2)
[1] 0.0315

Который соглашается, что 95% -й нормальный CI не включает нулевое значение отношения 1.

Процент CI получается (с некоторыми отличиями в связи с методами для связей):

quantile(city.boot$t, c(0.025, 0.975))

И значение р для процентиля начальной загрузки:

cvs <- quantile(city.boot$t0 - city.boot$t + 1, c(0.025, 0.975))
mean(city.boot$t > cvs[1] & city.boot$t < cvs[2])

Дает значение 0,035, которое также согласуется с доверительным интервалом с точки зрения исключения 1 из значения. В общем, мы не можем заметить, что, хотя ширина процентиля CI почти равна ширине нормального CI, и что процентиль CI дальше от нуля, что процентиль CI должен обеспечивать более низкие значения p. Это связано с тем, что форма распределения выборки, лежащего в основе CI для метода процентилей, является ненормальной.

Adamo
источник
Это очень интересный ответ @AdamO, но не могли бы вы привести несколько примеров? На R вы можете использовать функцию boot.ciи использовать аргумент «тип» для выбора студенческого CI (вы также можете выбрать BCA CI). Тем не менее, как вы можете рассчитать р-значения? Вы используете оценку или статистику теста? У меня был похожий вопрос, ответ на который был бы очень признателен.
Кевин Зарка,
1
+1 за четкое объяснение преимуществ студенческого бутстрапа.
eric_kernfeld
@KevinOunet Я привел два примера репликации p-значений из CI в загрузочном пакете. Это помогает?
AdamO
1
Спасибо @AdamO, это действительно помогает! Не могли бы вы привести последний пример для начальной загрузки?
Кевин Зарка,