Usando o bootstrap em H0 para executar um teste para a diferença de duas médias: substituição nos grupos ou na amostra agrupada

18

Suponha que eu tenha dados com dois grupos independentes:

g1.lengths <- c (112.64, 97.10, 84.18, 106.96, 98.42, 101.66)

g2.lengths <- c (84.44, 82.10, 83.26, 81.02, 81.86, 86.80, 
                     85.84, 97.08, 79.64, 83.32, 91.04, 85.92,
                     73.52, 85.58, 97.70, 89.72, 88.92, 103.72,
                     105.02, 99.48, 89.50, 81.74)

group = rep (c ("g1", "g2"), c (length (g1.lengths), length (g2.lengths)))

lengths = data.frame( lengths = c(g1.lengths, g2.lengths), group)

É evidente que o tamanho da amostra por grupo é enviesado, onde g1 tem 6 observações e g2 tem 22 . A ANOVA tradicional sugere que os grupos têm médias diferentes quando o valor crítico é definido como 0,05 (o valor de p é 0,0044 ).

summary (aov (lengths~group, data = lengths))  

Dado que meu objetivo é comparar a diferença média, esses dados amostrados desequilibrados e pequenos podem gerar resultados inadequados com a abordagem tradicional. Portanto, quero executar o teste de permutação e a inicialização.

ENSAIO DE PERMUTAÇÃO

A hipótese nula (H0) afirma que as médias do grupo são as mesmas. Essa suposição no teste de permutação é justificada pelo agrupamento de grupos em uma amostra. Isso garante que as amostras para dois grupos foram retiradas da distribuição idêntica. Por amostragem repetida (ou mais precisamente - reorganização) dos dados agrupados, as observações são realocadas (embaralhadas) para as amostras de uma nova maneira, e a estatística de teste é calculada. Realizando isso n vezes, fornecerá distribuição de amostragem das estatísticas de teste sob a suposição de que H0 é TRUE. No final, sob H0, o valor de p é a probabilidade de que a estatística do teste seja igual ou superior ao valor observado.

s.size.g1 <- length (g1.lengths)
s.size.g2 <- length (g2.lengths)

pool <- lengths$lengths
obs.diff.p <- mean (g1.lengths) - mean (g2.lengths)
iterations <- 10000
sampl.dist.p <- NULL

set.seed (5)
for (i in 1 : iterations) {
        resample <- sample (c(1:length (pool)), length(pool))

        g1.perm = pool[resample][1 : s.size.g1]
        g2.perm = pool[resample][(s.size.g1+1) : length(pool)]
        sampl.dist.p[i] = mean (g1.perm) - mean (g2.perm) 
}
p.permute <- (sum (abs (sampl.dist.p) >= abs(obs.diff.p)) + 1)/ (iterations+1)

O valor p reportado do teste de permutação é 0,0053 . OK, se eu fiz corretamente, permutações e ANOVA paramétrica dão resultados quase idênticos.

BOOTSTRAP

Antes de tudo, sei que o bootstrap não pode ajudar quando o tamanho da amostra é muito pequeno. Este post mostrou que pode ser ainda pior e enganador . Além disso, o segundo destacou que o teste de permutação geralmente é melhor do que o bootstrap quando o teste de hipóteses é o objetivo principal. No entanto, este ótimo post aborda diferenças importantes entre os métodos intensivos em computador. No entanto, aqui quero levantar (acredito) uma questão diferente.

Deixe-me apresentar primeiro a abordagem de bootstrap mais comum (Bootstrap1: reamostragem na amostra em pool ):

s.size.g1 <- length (g1.lengths)
s.size.g2 <- length (g2.lengths)

pool <- lengths$lengths
obs.diff.b1 <- mean (g1.lengths) - mean (g2.lengths)
iterations <- 10000
sampl.dist.b1 <- NULL

set.seed (5)
for (i in 1 : iterations) {
        resample <- sample (c(1:length (pool)), length(pool), replace = TRUE) 
        # "replace = TRUE" is the only difference between bootstrap and permutations

        g1.perm = pool[resample][1 : s.size.g1]
        g2.perm = pool[resample][(s.size.g1+1) : length(pool)]
        sampl.dist.b1[i] = mean (g1.perm) - mean (g2.perm) 
}
p.boot1 <- (sum (abs (sampl.dist.b1) >= obs.diff.b1) + 1)/ (iterations+1)

O valor p do bootstrap realizado dessa maneira é 0,005 . Mesmo que isso pareça razoável e quase idêntico ao teste paramétrico de ANOVA e permutação, é apropriado justificar H0 nesse bootstrap com base no fato de termos acabado de reunir amostras das quais extraímos amostras subsequentes?

Abordagem diferente que encontrei em vários trabalhos científicos. Especificamente, vi que os pesquisadores modificam os dados para atender ao H0 antes do bootstrap. Pesquisando, encontrei um post muito interessante no CV, onde @ jan.s explicou resultados incomuns do bootstrap na pergunta do post, onde o objetivo era comparar dois meios. No entanto, nesse post, não é abordado como executar a inicialização quando os dados são modificados antes da inicialização. A abordagem na qual os dados são modificados antes da inicialização é assim:

  1. H0 afirma que as médias de dois grupos são as mesmas
  2. H0 é verdadeiro se subtrairmos as observações individuais da média da amostra agrupada

Nesse caso, a modificação dos dados deve afetar a média dos grupos e, portanto, sua diferença, mas não a variação dentro (e entre) dos grupos.

  1. Os dados modificados servirão de base para mais bootstrap, com ressalvas de que a amostragem é realizada dentro de cada grupo separadamente .
  2. A diferença entre a média inicializada de g1 e g2 é calculada e comparada com a diferença observada (não modificada) entre os grupos.
  3. A proporção de valores iguais ou mais extremos do que o observado, dividido pelo número de iterações, fornecerá o valor p.

Aqui está o código (Bootstrap2: reamostragem nos grupos após a modificação que H0 é VERDADEIRO ):

s.size.g1 <- length (g1.lengths)
s.size.g2 <- length (g2.lengths)

pool <- lengths$lengths
obs.diff.b2 <- mean (g1.lengths) - mean (g2.lengths)

# make H0 to be true (no difference between means of two groups)
H0 <- pool - mean (pool)

# g1 from H0 
g1.H0 <- H0[1:s.size.g1] 

# g2 from H0
g2.H0 <- H0[(s.size.g1+1):length(pool)]

iterations <- 10000
sampl.dist.b2 <- NULL

set.seed (5)
for (i in 1 : iterations) {
        # Sample with replacement in g1
        g1.boot = sample (g1.H0, replace = T)

        # Sample with replacement in g2
        g2.boot = sample (g2.H0, replace = T)

        # bootstrapped difference
        sampl.dist.b2[i] <- mean (g1.boot) - mean (g2.boot)  
}
p.boot2 <- (sum (abs (sampl.dist.b2) >= obs.diff.b2) + 1)/ (iterations+1)

Esse bootstrap executado dará um valor de p de 0,514, que é tremendamente diferente em comparação com testes anteriores. Acredito que isso tenha a ver com a explicação do @ jan.s , mas não consigo descobrir onde está a chave ...

Newbie_R
fonte
1
Problema interessante e bem apresentado. O bootstrap tem problemas quando o tamanho da amostra é muito pequeno apenas porque é mais provável que a amostra original não represente muito bem a população. O tamanho da amostra não precisa ser muito grande para o bootstrap funcionar. Os tamanhos de amostra de 6 e 22 podem não ser tão ruins. Em um artigo de Efron (1983), o bootstrap foi comparado a uma forma de CV para estimar taxas de erro para funções discriminantes lineares em problemas de classificação com 2 classes em que cada tamanho de amostra de treinamento é menor que 10. Eu cubro isso no meu livro Bootstrap Methods ( 2007).
Michael R. Chernick
2
Meu livro também possui um capítulo sobre quando o bootstrap falha e inclui uma discussão sobre a questão do tamanho pequeno da amostra.
Michael R. Chernick
A única linha que você precisa corrigir em seu bootstrap # 2 abordagem para fazer o trabalho é este: H0 <- pool - mean (pool). Ele precisa ser substituído por H0 <- c(g1.lengths - mean(g1.lengths), g2.lengths - mean(g2.lengths)). Então você obterá o valor-p de 0,0023. (É a mesma coisa que o Zenit explicou em sua resposta.) Isso é tudo o que existe, apenas um simples bug no código. CC para @MichaelChernick
ameba diz Reinstate Monica
esses métodos poderiam ser dominados? Quero dizer se eles poderiam detectar QUALQUER diferença como significativa, quando os grupos são muito grandes: pool> 43k.
Alex Alvarez Pérez

Respostas:

17

Aqui está minha opinião, com base no capítulo 16 da introdução de Efron e Tibshirani, uma introdução à inicialização (página 220-224). O resumo é que seu segundo algoritmo de inicialização foi implementado incorretamente, mas a ideia geral está correta.

Ao realizar testes de autoinicialização, é necessário garantir que o método de re-amostragem gere dados que correspondam à hipótese nula. Usarei os dados do sono em R para ilustrar este post. Observe que estou usando a estatística do teste estudado em vez da diferença de médias, recomendada pelo livro didático.

O teste t clássico, que usa um resultado analítico para obter informações sobre a distribuição amostral da estatística t, produz o seguinte resultado:

x <- sleep$extra[sleep$group==1] y <- sleep$extra[sleep$group==2]
t.test(x,y)
t = -1.8608, df = 17.776, p-value = 0.07939

n1n2

# pooled sample, assumes equal variance
pooled <- c(x,y)
for (i in 1:10000){
  sample.index <- sample(c(1:length(pooled)),replace=TRUE)
  sample.x <- pooled[sample.index][1:length(x)]
  sample.y <- pooled[sample.index][-c(1:length(y))]
  boot.t[i] <- t.test(sample.x,sample.y)$statistic
}
p.pooled <-  (1 + sum(abs(boot.t) > abs(t.test(x,y)$statistic))) / (10000+1) 
p.pooled
[1] 0.07929207

H0 0H0 0H0 0z¯

x~Eu=xEu-x¯+z¯
y~Eu=yEu-y¯+z¯

x~/y~z¯H0 0

# sample from H0 separately, no assumption about equal variance
xt <- x - mean(x) + mean(sleep$extra) #
yt <- y - mean(y) + mean(sleep$extra)

boot.t <- c(1:10000)
for (i in 1:10000){
  sample.x <- sample(xt,replace=TRUE)
  sample.y <- sample(yt,replace=TRUE)
  boot.t[i] <- t.test(sample.x,sample.y)$statistic
}
p.h0 <-  (1 + sum(abs(boot.t) > abs(t.test(x,y)$statistic))) / (10000+1)  # 
p.h0
[1] 0.08049195

Desta vez, acabamos com valores-p semelhantes para as três abordagens. Espero que isto ajude!

Zenit
fonte
1
você seria gentil e explicaria por que '1' foi adicionado ao seguinte: em (1 + sum(abs(boot.t) > abs(t.test(x,y)$statistic))) / (10000+1)vez de algo como isto: mean(abs(boot.t) > abs(t.test(x,y)$statistic))Obrigado pelo seu tempo.
TG_Montana
+1. Essa abordagem de inicialização de dados modificados para amostra de H0 tem um nome específico?
ameba diz Restabelecer Monica
3
H0 0p-vumaeuvocêe=número de vezes {t>tobs}BB
@amoeba: Não tenho certeza se esse procedimento tem um nome formal ou acordado, mas acho que pode ser descrito como um teste de inicialização para a igualdade de meios , em vez de distribuições . A página que mostra o procedimento completo está ausente no google books , mas sua motivação aparece na página 223. Outra descrição pode ser encontrada nessas notas, na página 13 ( galton.uchicago.edu/~eichler/stat24600/Handouts/bootstrap. pdf ).
Zenit
(+1) Excelente resposta. Você poderia explicar por que "esse algoritmo [reamostrando os dados eles mesmos sem centralizar] está realmente testando se a distribuição de xey é idêntica"? Entendo que essa estratégia de reamostragem garante que as distribuições sejam iguais, mas por que testar se são iguais?
passe metade de