Por que os testes de hipótese em conjuntos de dados reamostrados rejeitam o nulo com muita frequência?

10

tl; dr: Começando com um conjunto de dados gerado sob o valor nulo, fiz uma nova amostragem de casos com substituição e conduzi um teste de hipótese em cada conjunto de dados reamostrado. Esses testes de hipótese rejeitam o nulo mais de 5% do tempo.

Na simulação abaixo, muito simples, eu conjuntos de dados com e um modelo OLS simples para cada um. Em seguida, para cada conjunto de dados, eu gero 1000 novos conjuntos de dados re-amostrando linhas do conjunto de dados original com substituição (um algoritmo especificamente descrito no texto clássico de Davison & Hinkley como apropriado para a regressão linear). Para cada um deles, eu me encaixo no mesmo modelo OLS. Por fim, cerca de 16% dos testes de hipótese nas amostras de bootstrap rejeitam o nulo , enquanto que devemos obter 5% (como fazemos nos conjuntos de dados originais).XN(0 0,1 1)⨿YN(0 0,1 1)

Eu suspeitava que isso tivesse algo a ver com observações repetidas, causando associações infladas; portanto, para comparação, tentei duas outras abordagens no código abaixo (comentado). No método 2, eu corrijo e substituo por resíduos reamostrados do modelo OLS no conjunto de dados original. No método 3, eu desenho uma subamostra aleatória sem substituição. Ambas as alternativas funcionam, ou seja, seus testes de hipótese rejeitam o nulo 5% das vezes.XY

Minha pergunta: Estou certo de que observações repetidas são as culpadas? Se sim, considerando que essa é uma abordagem padrão para o bootstrap, onde exatamente estamos violando a teoria padrão do bootstrap?

Atualização # 1: Mais simulações

Eu tentei um cenário ainda mais simples, um modelo de regressão somente interceptar para . O mesmo problema ocorre.Y

# note: simulation takes 5-10 min on my laptop; can reduce boot.reps
#  and n.sims.run if wanted
# set the number of cores: can change this to match your machine
library(doParallel)
registerDoParallel(cores=8)
boot.reps = 1000
n.sims.run = 1000

for ( j in 1:n.sims.run ) {

  # make initial dataset from which to bootstrap
  # generate under null
  d = data.frame( X1 = rnorm( n = 1000 ), Y1 = rnorm( n = 1000 ) )

  # fit OLS to original data
  mod.orig = lm( Y1 ~ X1, data = d )
  bhat = coef( mod.orig )[["X1"]]
  se = coef(summary(mod.orig))["X1",2]
  rej = coef(summary(mod.orig))["X1",4] < 0.05

  # run all bootstrap iterates
  parallel.time = system.time( {
    r = foreach( icount( boot.reps ), .combine=rbind ) %dopar% {

      # Algorithm 6.2: Resample entire cases - FAILS
      # residuals of this model are repeated, so not normal?
      ids = sample( 1:nrow(d), replace=TRUE )
      b = d[ ids, ]

      # # Method 2: Resample just the residuals themselves - WORKS
      # b = data.frame( X1 = d$X1, Y1 = sample(mod.orig$residuals, replace = TRUE) )

      # # Method 3: Subsampling without replacement - WORKS
      # ids = sample( 1:nrow(d), size = 500, replace=FALSE )
      # b = d[ ids, ]

      # save stats from bootstrap sample
      mod = lm( Y1 ~ X1, data = b ) 
      data.frame( bhat = coef( mod )[["X1"]],
                  se = coef(summary(mod))["X1",2],
                  rej = coef(summary(mod))["X1",4] < 0.05 )

    }
  } )[3]


  ###### Results for This Simulation Rep #####
  r = data.frame(r)
  names(r) = c( "bhat.bt", "se.bt", "rej.bt" )

  # return results of each bootstrap iterate
  new.rows = data.frame( bt.iterate = 1:boot.reps,
                         bhat.bt = r$bhat.bt,
                         se.bt = r$se.bt,
                         rej.bt = r$rej.bt )
  # along with results from original sample
  new.rows$bhat = bhat
  new.rows$se = se
  new.rows$rej = rej

  # add row to output file
  if ( j == 1 ) res = new.rows
  else res = rbind( res, new.rows )
  # res should have boot.reps rows per "j" in the for-loop

  # simulation rep counter
  d$sim.rep = j

}  # end loop over j simulation reps



##### Analyze results #####

# dataset with only one row per simulation
s = res[ res$bt.iterate == 1, ]

# prob of rejecting within each resample
# should be 0.05
mean(res$rej.bt); mean(s$rej)

Atualização # 2: A resposta

Várias possibilidades foram propostas nos comentários e respostas, e eu fiz mais simulações para testá-las empiricamente. Acontece que o JWalker está certo de que o problema é que precisamos centralizar as estatísticas de autoinicialização pela estimativa dos dados originais para obter a distribuição de amostragem correta em . No entanto, também acho que o comentário do whuber sobre a violação dos pressupostos dos testes paramétricos também está correto, embora neste caso realmente recebamos falsos positivos nominais quando solucionamos o problema do JWalker.H0 0

meia passagem
fonte
11
No bootstrap padrão, você consideraria apenas a distribuição do coeficiente de X1, não seus valores de p associados. Portanto, não é um problema do bootstrap. No entanto, sua observação é interessante e não intuitiva.
Michael M
11
@ MichaelM, isso é verdade. Mas como o CDF conjunto dos dados nas reamostragens deve convergir em n e o número de autoinicializações itera para o CDF verdadeiro que gerou os dados originais, eu não esperaria que os valores-p também diferissem.
meia-passagem
Direita. Estou certo de que o efeito advém de observações não independentes (como você disse), gerando erros padrão muito otimistas. Na sua simulação, parece ser a única suposição violada do modelo linear normal. Talvez possamos derivar o fator de deflação da variância correspondente.
Michael M
2
Uma coisa que está clara no método 1 é a violação da suposição de erro iid: ao reamostrar com substituição, os resíduos de qualquer valor são perfeitamente correlacionados em vez de independentes! Portanto, você não está inicializando corretamente, isso é tudo. Como demonstração, após a computação, substitua-os por, mas proceda exatamente como antes. Isso lida corretamente com as duplicatas (embora produza uma amostra menor). Você obterá uma distribuição uniforme dos valores-p. xidsids <- unique(ids)
whuber
2
@whuber. Eu vejo. E isso explicaria por que a reamostragem de resíduos com substituição funciona, apesar das observações repetidas: os resíduos desse modelo são mais uma vez independentes de X. Se você gostaria de transformar isso em uma resposta, ficaria feliz em aceitar.
passe metade

Respostas:

5

Quando você redimensiona o valor nulo, o valor esperado do coeficiente de regressão é zero. Quando você reamostrar alguns conjuntos de dados observados, o valor esperado é o coeficiente observado para esses dados. Não é um erro do tipo I se P <= 0,05 quando você reamostrar os dados observados. De fato, é um erro do tipo II se você tiver P> 0,05.

Você pode obter alguma intuição calculando a correlação entre o abs (b) e a média (P). Aqui está o código mais simples para replicar o que você fez, além de calcular a correlação entre o erro "tipo I" no conjunto de simulações

boot.reps = 1000
n.sims.run = 10
n <- 1000
b <- matrix(NA, nrow=boot.reps, ncol=n.sims.run)
p <- matrix(NA, nrow=boot.reps, ncol=n.sims.run)
for(sim_j in 1:n.sims.run){
  x <- rnorm(n)
  y <- rnorm(n)
  inc <- 1:n
  for(boot_i in 1:boot.reps){
    fit <- lm(y[inc] ~ x[inc])
    b[boot_i, sim_j] <- abs(coefficients(summary(fit))['x[inc]', 'Estimate'])
    p[boot_i, sim_j] <- coefficients(summary(fit))['x[inc]', 'Pr(>|t|)']
    inc <- sample(1:n, replace=TRUE)
  }
}
# note this is not really a type I error but whatever
type1 <- apply(p, 2, function(x) sum(x <= 0.05))/boot.reps
# correlation between b and "type I"
cor(b[1, ], type1)

atualizar a resposta por grand_chat não é a razão pela qual a frequência de P <= 0,05 é> 0,05. A resposta é muito simples e o que eu disse acima - o valor esperado da média de cada nova amostra é a média original observada. Esta é toda a base do bootstrap, que foi desenvolvido para gerar erros / limites de confiança padrão em uma média observada e não como um teste de hipótese. Como a expectativa não é zero, é claro que o "erro tipo I" será maior que alfa. E é por isso que haverá uma correlação entre a magnitude do coeficiente (a que distância de zero) e a magnitude do desvio do "erro tipo I" de alfa.

JWalker
fonte
H0 0:β=β^H0 0:β=0 0
H0 0:β=βˆ testa a equivalência e requer uma abordagem diferente do desenho do estudo. é usado quando o importante é garantir que suas diferenças observadas não sejam aleatórias, equivalência quando você deseja garantir que sua previsão esteja correta. Infelizmente, costuma ser visto como tamanho único, mas depende dos riscos em sua situação. É típico o uso de na pesquisa em estágio inicial para filtrar os flukes quando você não sabe o suficiente para definir uma hipótese alternativa e, quando mais se sabe, pode fazer sentido mudar para testar a correção do seu conhecimento. H0 0:β=0 0H0 0:β=0 0
ReneBt
2

Se você fizer uma amostra com a substituição da sua amostra normal original, a amostra de inicialização resultante não será normal . A distribuição conjunta da amostra de bootstrap segue uma distribuição de mistura gnarly que provavelmente contém registros duplicados, enquanto que valores duplicados têm probabilidade zero quando você coleta amostras de uma distribuição normal.

Como um exemplo simples, se sua amostra original for duas observações de uma distribuição normal univariada, uma amostra de autoinicialização com substituição consistirá na metade do tempo na amostra original e na metade do tempo consistirá em um dos valores originais duplicados. É claro que a variação da amostra de bootstrap será, em média, menor que a do original - na verdade, será metade do original.

A principal conseqüência é que a inferência que você está fazendo com base na teoria normal retorna os valores incorretos quando aplicados ao exemplo de autoinicialização. Em particular, a teoria normal produz regras de decisão anticonservadoras, porque sua amostra de bootstrap produzirá estatísticas cujos denominadores são menores do que o esperado na teoria normal, devido à presença de duplicatas. Como resultado, o teste de hipótese da teoria normal acaba rejeitando a hipótese nula mais do que o esperado.pt

grand_chat
fonte
Mas, se esse for o caso, não teríamos exatamente o mesmo problema ao reamostrar os resíduos com a substituição? No entanto, de fato, essa abordagem rejeita com probabilidade nominal.
passe metade
Além disso, um teste t com n = 1000 não deve ter problemas com dados não normais.
passe metade
0

Concordo totalmente com a resposta da @ JWalker.

Há outro aspecto desse problema. Isso está no seu processo de reamostragem. Você espera que o coeficiente de regressão seja centrado em torno de zero porque você assume Xe Yé independente. No entanto, na sua reamostragem, você faz

ids = sample( 1:nrow(d), replace=TRUE )
  b = d[ ids, ]

o que cria correlação porque você está amostrando Xe Yjuntos. Por exemplo, digamos que a primeira linha do conjunto de dados dseja (x1, y1): No conjunto de dados reamostrado P(Y = y1|X = x1) = 1, enquanto se Xe Yfor independente, P(Y|X = x1)segue uma distribuição normal.

Então, outra maneira de corrigir isso é usar

b = data.frame( X1 = rnorm( n = 1000 ), Y1 = rnorm( n = 1000 ) )

o mesmo código que você usa para gerar d, a fim de tornar X e Y independentes um do outro.

O mesmo motivo explica por que ele funciona com reamostragem residual (porque Xé independente do novo Y).

Tianxia Zhou
fonte
Por um tempo, também pensei que as observações reamostradas pudessem ser independentes, mas, pensando muito mais sobre isso, não acho que seja esse o caso: stats.stackexchange.com/questions/339237/…
half -pass
O problema que estou descrevendo acima é diferente do seu post. O que você se referiu é a independência de x's. A que me referi é a correlação entre Xs e Ys.
Tianxia Zhou
-1

O maior problema aqui é que os resultados do modelo são espúrios e, portanto, altamente instáveis, porque o modelo é apenas um ruído adequado. Num sentido muito literal. Y1 não é uma variável dependente devido à forma como os dados da amostra foram gerados.


Editar, em resposta aos comentários. Deixe-me tentar outra vez explicar o meu pensamento.

Com um OLS, a intenção geral é descobrir e quantificar os relacionamentos subjacentes nos dados. Com dados do mundo real, geralmente não sabemos exatamente.

Mas esta é uma situação de teste artificial. Conhecemos o mecanismo EXATO de geração de dados, está bem ali no código postado pelo OP É

X1 = rnorm (n = 1000), Y1 = rnorm (n = 1000)

Se expressarmos isso na forma familiar de uma regressão OLS, ou seja,

Y1 = interceptação + Beta1 * X1 + Erro
que se torna
Y1 = média (X1) + 0 (X1) + Erro

Então, na minha opinião, este é um modelo expresso em FORM linear, mas NÃO é realmente uma relação / modelo linear, porque não há inclinação. Beta1 = 0.000000.

Quando gerarmos os 1000 pontos de dados aleatórios, o gráfico de dispersão será semelhante ao spray circular de espingarda clássica. Pode haver alguma correlação entre X1 e Y1 na amostra específica de 1000 pontos aleatórios gerados, mas, nesse caso, é uma ocorrência aleatória. Se o OLS encontrar uma correlação, ou seja, rejeitar a hipótese nula de que não há inclinação, uma vez que sabemos definitivamente que realmente não há relação entre essas duas variáveis, o OLS encontrou literalmente um padrão no Componente de Erro. Eu caracterizei isso como "ajustando o barulho" e "espúrio".

Além disso, uma das suposições / requisitos padrão de um OLS é que (+/-) "o modelo de regressão linear é" linear em parâmetros ". Dados os dados, minha opinião é de que não satisfazemos essa suposição.Portanto, as estatísticas de teste subjacentes à significância são imprecisas.Minha convicção é de que a violação da suposição de linearidade é a causa direta dos resultados não intuitivos do bootstrap.

Quando eu li esse problema pela primeira vez, não percebi que o OP pretendia testar sob a hipótese [nula].

Mas os mesmos resultados não intuitivos aconteceriam se o conjunto de dados tivesse sido gerado como

X1 = rnorm (n = 1000), Y1 = X1 * .4 + rnorm (n = 1000)?

Doug Dame
fonte
4
Y1 1Y1 1
(-1) pelas mesmas razões que @whuber deu.
half-pass
11
Resposta à pergunta final da sua edição: sim, definitivamente. Tente você mesmo com uma simulação. (Mas ter cuidado com a interpretação, porque você tem que considerar que o nulo é e qual é o estado real das coisas é.)
whuber