Avaliando a potência de um teste de normalidade (em R)

9

Quero avaliar a precisão dos testes de normalidade em diferentes tamanhos de amostra em R (percebo que os testes de normalidade podem ser enganosos ). Por exemplo, para examinar o teste de Shapiro-Wilk, estou realizando a seguinte simulação (bem como plotando os resultados) e esperaria que, à medida que o tamanho da amostra aumentasse, a probabilidade de rejeitar as reduções nulas:

n <- 1000
pvalue_mat <- matrix(NA, ncol = 1, nrow = n)

for(i in 10:n){
    x1 <- rnorm(i, mean = 0, sd = 1)
    pvalue_mat[i,] <- shapiro.test(x1)$p.value
}   

plot(pvalue_mat)

Meu pensamento seria que, à medida que o tamanho da amostra cresce, deve haver uma menor taxa de rejeição, no entanto, parece bastante uniforme. Acho que estou entendendo mal isso - todos e quaisquer pensamentos são bem-vindos.

user94759
fonte
2
Você pode dar uma olhada em: stats.stackexchange.com/questions/2492/…
nico

Respostas:

7

Você está simulando sob a hipótese nula (distribuição normal); portanto, a taxa de rejeição tenderá ao nível de significância conforme o esperado. Para avaliar o poder, você precisa simular sob qualquer distribuição não normal. Existem infinitas possibilidades / cenários (por exemplo, distribuições gama com assimetria crescente, distribuição t com df decrescente etc.) para escolher, dependendo do escopo do seu estudo.

Michael M
fonte
Obrigado pela resposta. Quando simulo sobre distribuições não normais, observo um padrão convexo em relação à origem - ou seja, à medida que o tamanho da amostra aumenta para qualquer distribuição não normal, aumenta a probabilidade de rejeitar o nulo de normalidade. No entanto, não entendo por que não é o inverso ao desenhar a partir de uma distribuição normal: por que a probabilidade de rejeitar o nulo diminui à medida que o tamanho da amostra aumenta? Graças
user94759
3
Porque a probabilidade de cometer esse erro tipo 1 é, por definição, igual ao nível de significância, que é constante. Em outras palavras, os valores de p são distribuídos uniformemente sob o valor nulo. Aliás, é aconselhável fazer muitas simulações por configuração, incluindo a escolha de n, não apenas uma como no seu código.
Michael M
7

A compreensão da análise de poder dos testes estatísticos de hipóteses pode ser aprimorada, realizando-se alguns e observando atentamente os resultados.


αα

O segundo critério exige que estipulemos de que maneira (s) e por quanto o nulo falha em ser verdadeiro. Nos casos de livros didáticos, isso é fácil, porque as alternativas são limitadas em escopo e claramente especificadas. Em testes de distribuição como o Shapiro-Wilk, as alternativas são muito mais vagas: são "não-normais". Ao escolher entre os testes de distribuição, é provável que o analista precise realizar seu próprio estudo de potência pontual para avaliar quão bem os testes funcionam em relação a hipóteses alternativas mais específicas que são motivo de preocupação no problema em questão.

ν1ν

αR

  • rdist, o nome de uma função para produzir uma amostra aleatória de alguma distribuição

  • n, o tamanho das amostras a serem solicitadas rdist

  • n.iter, o número dessas amostras para obter

  • ...rdistν

Os demais parâmetros controlam a exibição dos resultados; eles são incluídos principalmente como uma conveniência para gerar os números nesta resposta.

sim <- function(rdist, n, n.iter, prefix="",
                breaks=seq(0, 1, length.out=20), alpha=0.05,
                plot=TRUE, ...) {

  # The simulated P-values.
  # NB: The optional arguments "..." are passed to `rdist` to specify
  #     its parameters (if any).
  x <- apply(matrix(rdist(n*n.iter, ...), ncol=n.iter), 2, 
             function(y) shapiro.test(y)$p.value)

  # The histogram of P-values, if requested.
  if (plot) {
    power <- mean(x <= alpha)
    round.n <- 1+ceiling(log(1 + n.iter * power * (1-power), base=10) / 2)
    hist(x[x <= max(breaks)], xlab=paste("P value (n=", n, ")", sep=""), 
         breaks=breaks, 
         main=paste(prefix, "(power=", format(power, digits=round.n), ")", sep=""))
    # Specially color the "significant" part of the histogram
    hist(x[x <= alpha], breaks=breaks, col="#e0404080", add=TRUE)
  }

  # Return the array of P-values for any further processing.
  return(x)
}

5100.n520.

n.iter <- 10^5                 # Number of samples to generate
n.spec <- c(5, 10, 20)         # Sample sizes to study
par(mfrow=c(1,length(n.spec))) # Organize subsequent plots into a tableau
system.time(
  invisible(sapply(n.spec, function(n) sim(rnorm, n, n.iter, prefix="DF = Inf ")))
)

Depois de especificar os parâmetros, esse código também é apenas uma linha. Ele produz a seguinte saída:

Histogramas para o nulo

01α=0.05,.04810.0499

10.2

νν=100ν=11001000), o que não leva tempo. O código agora requer um loop duplo (e em situações mais complexas, geralmente precisamos de loops triplos ou quádruplos para acomodar todos os aspectos que precisamos variar): um para estudar como a potência varia com o tamanho da amostra e outro para estudar como ele varia com os graus de liberdade. Mais uma vez, porém, tudo é feito em apenas uma linha de código (a terceira e a final):

df.spec <- c(64, 16, 4, 2, 1)
par(mfrow=c(length(n.spec), length(df.spec)))
for (n in n.spec) 
  for (df in df.spec)
    tmp <- sim(rt, n, n.iter, prefix=paste("DF =", df, ""), df=df)

Histogramas para as alternativas

Um pequeno estudo desse quadro fornece uma boa intuição sobre o poder. Gostaria de chamar a atenção para seus aspectos mais importantes e úteis:

  • ν=64ν=1

  • n=5n=20

  • 1

  • 201110086.57=13%205%95%

  • αα=0.10

    α0.05020%α=0.01α=0.005breakssim


É divertido saber que muito pode ser coletado do que, na verdade, equivale a três linhas de código: uma para simular amostras de iid de uma distribuição especificada, uma para aplicá-la a uma matriz de distribuições nulas e a terceira para aplicá-la a uma matriz de distribuições alternativas. Estas são as três etapas que entram em qualquer análise de poder: o resto é apenas resumir e interpretar os resultados.

whuber
fonte
1

(Mais do que um comentário, talvez não seja uma resposta completa)

[I] esperaria que, à medida que o tamanho da amostra aumenta, a probabilidade de rejeitar o valor nulo diminua.

Deixando de lado as considerações dos testes tendenciosos (que não são incomuns na qualidade do ajuste, portanto vale a pena mencionar), há três situações relacionadas à taxa de rejeição que se pode considerar:

1) a taxa de rejeição ao simular a partir do nulo (como você parece estar fazendo na sua pergunta)

α

2) a taxa de rejeição ao simular a partir de alguma alternativa

Aqui a taxa de rejeição deve aumentar à medida que n aumenta.

3) a taxa de rejeição para alguma coleta de dados reais

Na prática, o nulo nunca é realmente verdadeiro e os dados reais terão uma mistura de quantidades de não normalidade (conforme medido pela estatística do teste). Se o grau de não normalidade não estiver relacionado ao tamanho da amostra, a taxa de rejeição deve aumentar à medida que n aumenta.

Portanto, de fato, em nenhuma dessas situações devemos ver a taxa de rejeição diminuir com o tamanho da amostra.

Glen_b -Reinstate Monica
fonte