Computando o valor p usando bootstrap com R

28

Eu uso o pacote "boot" para calcular um valor p aproximado de inicialização de dois lados, mas o resultado está muito longe do valor p de usar t.test. Não consigo descobrir o que fiz de errado no meu código R. Alguém pode me dar uma dica para isso

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 

O valor de p inicializado com dois lados (pvalor) = 0,4804, mas o valor de p com dois lados de t.test é 0,04342. Ambos os valores de p têm cerca de 11 vezes a diferença. Como isso pode acontecer?

Tu.2
fonte
como vem b3 $ t0 tem duas entradas?
Xian
1
é um nome de colônia!
Elvis
2
Você está calculando um valor incorretamente. A documentação diz que t 0 é a estatística observada, não a distribuição nula, como a notação sugere. Você precisa apresentar uma estimativa da amostra dist-n sob o valor nulo. Veja minha resposta para mais informações. Tente fazer um teste sem correção de viés. pt0 0mean(abs(b3$t0) < abs(b3$t-mean(b3$t)))
AdamO 5/05

Respostas:

31

Você está usando o bootstrap para gerar dados sob a distribuição empírica dos dados observados. Isso pode ser útil para fornecer um intervalo de confiança na diferença entre os dois meios:

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

Para obter um valor- , você precisa gerar permutações sob a hipótese nula. Isso pode ser feito, por exemplo:p

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

Nesta solução, o tamanho dos grupos não é fixo, você reatribui aleatoriamente um grupo para cada indivíduo, inicializando a partir do conjunto de grupos inicial. Parece-me legítimo, no entanto, uma solução mais clássica é fixar o número de indivíduos de cada grupo; portanto, você apenas permuta os grupos em vez de iniciar (isso geralmente é motivado pelo design do experimento, onde os tamanhos dos grupos são previamente fixados). ):

> 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
Elvis
fonte
5
Tecnicamente, este é o teste de permutação, não um valor p de autoinicialização.
AdamO 3/17/17
@AdamO Concordo que o que é apresentado nesta resposta é o teste de permutação (e sua variante levemente modificada); isso ocorre porque durante a reamostragem, os grupos são agrupados. Por outro lado, em um teste baseado em autoinicialização, os valores para cada grupo devem ser amostrados usando apenas os dados para o mesmo grupo. Aqui está uma resposta que explica como fazê-lo: stats.stackexchange.com/a/187630/28666 .
ameba diz Restabelecer Monica
@amoeba Acho que a resposta que você vincula também é um teste de permutação, relacionado ao bootstrap apenas na medida em que envolvem reamostragem. É perfeitamente bom relatar, mas para relatar são dois métodos que estão sendo usados. Um bootstrap não paramétrico é tecnicamente incapaz de gerar dados sob uma hipótese nula. Veja minha resposta para saber como os valores-p são gerados a partir de uma distribuição de inicialização.
AdamO 3/17/17
@AdamO Acho que é uma questão de terminologia, mas não vejo como o procedimento descrito na resposta vinculada pode ser chamado de teste de "permutação" porque nada é permutado lá: valores reamostrados para cada grupo são gerados usando os dados desse somente grupo.
Ameba diz Restabelecer Monica
1
Elvis, acho que o primeiro código na sua resposta também é o teste de permutação. Quando você reamostrar, você agrupa os grupos! É isso que define o teste de permutação.
Ameba diz Reinstate Monica
25

A resposta de Elvis depende de permutações, mas, na minha opinião, não deixa claro o que há de errado com a abordagem de autoinicialização original. Deixe-me discutir uma solução baseada apenas no bootstrap.

O problema crucial da sua simulação original é que o bootstrap sempre fornece a distribuição VERDADEIRA da estatística de teste. No entanto, ao calcular o valor-p, você deve comparar o valor obtido da estatística de teste com sua distribuição SOB H0, ou seja, não com a distribuição verdadeira!

[Vamos deixar claro. Por exemplo, sabe-se que a estatística T do teste t clássico tem a distribuição t "central" clássica sob H0 e uma distribuição não central em geral. No entanto, todos estão familiarizados com o fato de que o valor observado de T é comparado com a distribuição t "central" clássica, ou seja, não se tenta obter a verdadeira distribuição t [não-central] para fazer a comparação com T.]

Seu valor p 0,4804 é muito grande, porque o valor observado "t0" da estatística do teste Mean [1] -Mean [2] fica muito próximo ao centro da amostra inicializada "t". É natural e normalmente é sempre assim [ou seja, independentemente da validade de H0], porque a amostra inicializada "t" emula a distribuição REAL de Mean [1] -Mean [2]. Mas, como observado acima [e também por Elvis], o que você realmente precisa é da distribuição de Mean [1] -Mean [2] SOB H0. É obvio que

1) em H0, a distribuição de Média [1] -Média [2] será centrada em torno de 0,

2) sua forma não depende da validade de H0.

Esses dois pontos sugerem que a distribuição de Média [1] -Média [2] sob H0 pode ser emulada pela amostra inicializada "t" SHIFTED, de modo que seja centralizada em torno de 0. Em R:

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

e o valor p correspondente será:

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

que fornece um valor "muito bom" de 0,0232. :-)

Deixe-me observar que o ponto "2)" mencionado acima é chamado de "equivalência de tradução" da estatística de teste e NÃO precisa se sustentar em geral! Ou seja, para algumas estatísticas de teste, o deslocamento do "t" com inicialização inicial não fornece uma estimativa válida da distribuição da estatística de teste em HO! Dê uma olhada nesta discussão e especialmente na resposta de P. Dalgaard: http://tolstoy.newcastle.edu.au/R/e6/help/09/04/11096.html

Seu problema de teste produz uma distribuição perfeitamente simétrica da estatística de teste, mas lembre-se de que existem alguns problemas com a obtenção de valores p de DOIS LADOS, no caso de uma distribuição com distorção de inicialização da estatística de teste. Mais uma vez, leia o link acima.

[E, finalmente, eu usaria o teste de permutação "puro" em sua situação; ou seja, a segunda metade da resposta de Elvis. :-)]

jan s.
fonte
17

Existem várias maneiras de calcular os ICs de inicialização e os valores de p. A questão principal é que é impossível para o bootstrap gerar dados sob uma hipótese nula. O teste de permutação é uma alternativa viável baseada em reamostragem para isso. Para usar um bootstrap adequado, você deve fazer algumas suposições sobre a distribuição de amostragem da estatística de teste.

β0 0=β^-β^β0 0=β^-β^

inicialização normal

Uma abordagem é um bootstrap normal em que você obtém a média e o desvio padrão da distribuição do bootstrap, calcula a distribuição de amostragem sob o nulo deslocando a distribuição e usando os percentis normais da distribuição nula no ponto da estimativa na amostra original do bootstrap . Essa é uma abordagem razoável quando a distribuição do bootstrap é normal, a inspeção visual geralmente é suficiente aqui. Os resultados dessa abordagem geralmente são muito próximos da estimativa de erro robusta ou baseada em sanduíche, robusta contra a heterocedasticidade e / ou suposições de variância finita da amostra. A suposição de uma estatística de teste normal é uma condição mais forte das suposições no próximo teste de autoinicialização que discutirei.

porcentagem de autoinicialização

F0 02×min(F0 0(β^),1-F0 0(β^))

Bootstrap estudantilizado

p

Exemplo de programação

Como exemplo, usarei os citydados no pacote de inicialização. Os intervalos de confiança da inicialização são calculados com este código:

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

e produza esta saída:

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

O IC de 95% para o bootstrap normal é obtido calculando:

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

O valor p é assim obtido:

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

O que concorda que o IC normal de 95% não inclui o valor da razão nula de 1.

O IC do percentil é obtido (com algumas diferenças devido aos métodos de empate):

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

E o valor p para o bootstrap de percentil é:

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

Dá um p de 0,035, que também concorda com o intervalo de confiança em termos da exclusão de 1 do valor. Em geral, não podemos observar que, enquanto a largura do IC do percentil é quase tão ampla quanto o IC normal e que o IC do percentil está mais distante do nulo que o IC do percentil deve fornecer valores p inferiores. Isso ocorre porque o formato da distribuição amostral subjacente ao IC para o método do percentil não é normal.

AdamO
fonte
É uma resposta muito interessante @AdamO, mas você poderia fornecer alguns exemplos? No R, você pode usar a função boot.cie o argumento "type" para escolher um IC estudado (você também pode escolher um IC BCA). No entanto, como você pode calcular os valores de p? Você está usando a estimativa ou a estatística do teste? Eu tive uma pergunta semelhante que resposta seria muito apreciada.
precisa saber é o seguinte
1
+1 para uma explicação clara dos benefícios da inicialização estudantil.
Eric_kernfeld
@KevinOunet Dei dois exemplos de replicação de valores p de ICs no pacote de inicialização. Isso ajuda?
7118 AdamO
1
Obrigado @AdamO, isso realmente ajuda! Você poderia fornecer um último exemplo para a inicialização estudada?
21818 Kevin Zarca