Testes de hipóteses como Anderson-Darling ou Shapiro-Wilk verificam a normalidade de uma distribuição. No entanto, se o tamanho da amostra for muito grande, o teste é extremamente “preciso”, mas praticamente inútil porque o intervalo de confiança é muito pequeno. Eles sempre rejeitarão a hipótese nula, mesmo se a distribuição for razoavelmente normal o suficiente.

Como devemos testar a normalidade quando o tamanho da amostra é muito grande, além de visualizar a densidade estimada?

A motivação é que queremos automatizar a verificação da normalidade de um conjunto de dados grande em uma plataforma de software, onde tudo precisa ser automatizado, não visualizado e inspecionado manualmente por humanos.

Uma descrição do conjunto de dados e os rótulos de nome de variáveis podem ser encontrados aqui. O código abaixo cria o conjunto de dados e então exibe as primeiras 6 observações.

set.seed(201)
N = 10000 # tamanho da amostra
dados = data.frame(dados=c(rnorm(N,0.5,0.078),rbeta(N,20,20)),Distr=c(rep("Normal",N),rep("Beta",N)))
library(janitor) # contém a função tabyl() para tabulação
library(dplyr) # operações em conjuntos de dados, por exemplo %>%
options(dplyr.summarise.inform = FALSE) # para suprimir mensagens de aviso 
library(arsenal) # função tableby()
library(kableExtra) # exibe a formatação da tabela
kbl(dados[c(1:6,(2*N-5):(2*N)), 1:2], 
    caption = "Tabela 1: Mostrando as primeiras 6 observações.") %>%
  kable_styling(bootstrap_options = "striped", full_width = TRUE, position = "left")
Tabela 1: Mostrando as primeiras 6 observações.
dados Distr
1 0.5223127 Normal
2 0.4731254 Normal
3 0.5254096 Normal
4 0.3676667 Normal
5 0.3997324 Normal
6 0.4941475 Normal
19995 0.4719642 Beta
19996 0.5772815 Beta
19997 0.5708642 Beta
19998 0.5942047 Beta
19999 0.5634953 Beta
20000 0.5565097 Beta

1 Estatísticas descritivas

A sintaxe abaixo usa as funções do pacote dplyr para calcular o mínimo, máximo, média, desvio padrão e distância interquantilica de cada amostra. Os resultados são apresentados na Tabela 2.

table2 <- dados %>% 
    group_by(Distr) %>% 
    summarise(Minimum = min(dados),
              Maximum = max(dados),
              Mean = mean(dados), 
              SD = sd(dados),
              IQR = diff(quantile(dados, c(1, 3)/4)))
names(table2)[1] = c("Distribuição")

kbl(table2,digits = 2, 
    caption = "Table 2: Estatísticas descritivas amostrais segundo a distribuição geradora.") %>%
  kable_styling(bootstrap_options = "striped", full_width = TRUE, position = "left")
Table 2: Estatísticas descritivas amostrais segundo a distribuição geradora.
Distribuição Minimum Maximum Mean SD IQR
Beta 0.24 0.76 0.5 0.08 0.11
Normal 0.20 0.77 0.5 0.08 0.11


Apresentação gráfica das amostras.

par(mfrow=c(1,2),mar=c(3,3,1,1))
x = seq(0,1,by=0.01)
plot(x,dnorm(x,0.5,0.078),type="l",lty=1,col="red")
densidade = density(dados[dados$Distr=="Normal",1])
lines(densidade$x,densidade$y,col="violet")
grid()
legend("right",legend=c("Normal","Estimativa"),col = c("red","violet"),lty=1)
rug(dados[dados$Distr=="Normal",1])
plot(x,dbeta(x,20,20),type="l",lty=1,col="red")
densidade = density(dados[dados$Distr=="Beta",1])
lines(densidade$x,densidade$y,col="violet")
grid()
legend("right",legend=c("Beta","Estimativa"),col = c("red","violet"),lty=1)
rug(dados[dados$Distr=="Beta",1])


2 Testes de normalidade

Vejamos o que acontece se aplicarmos testes de normalidade às amostras completas.

dados.normal = dados[dados$Distr=="Normal",1]
dados.normal = (dados.normal-mean(dados.normal))/sd(dados.normal)
ks.test(dados.normal,"pnorm")
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dados.normal
## D = 0.0063948, p-value = 0.8081
## alternative hypothesis: two-sided
dados.beta = dados[dados$Distr=="Beta",1]
dados.beta = (dados.beta-mean(dados.beta))/sd(dados.beta)
ks.test(dados.beta,"pnorm")
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dados.beta
## D = 0.008949, p-value = 0.3998
## alternative hypothesis: two-sided


Lembremos que a hipótese nula é \(H_0 \, : \, F(x)=F_0(x)\), para todo \(x\) contra a alternativa geral \(H_1 \, : F(x)\neq F_0(x)\), para algum \(x\). A distribuição \(F_0(x)\) é complemtamente especificada. Acontece que a função ks.test calcula o \(p\)-valor exato se \(n\leq 100\), caso contrário utiliza a distribuição assintótica. Nada é informado acerca de limites no tamanho de amostra.

A resposta obtida é contraditória, aceitamos ambas amostas, tanto a gerada da distribuição Normal(0.5,0.078) é considerada normal quanto a amostra gerada pela distribuição Beta(20,20).

nortest::lillie.test(dados[dados$Distr=="Normal",1])
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  dados[dados$Distr == "Normal", 1]
## D = 0.0063948, p-value = 0.4133
nortest::lillie.test(dados[dados$Distr=="Beta",1])
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  dados[dados$Distr == "Beta", 1]
## D = 0.008949, p-value = 0.06881


Neste caso a resposta também é contraditória. Aceita-se a normalidade dos dados N(0.5,0.078) e não rejeitamos a normalidade dos dados Beta(20,20). Informa-se muito acerca da forma de obtensão do \(p\)-valor e nada acerca do tamanho da amostra permitido.

Lembremos que, se \(X\) tiver distribuição Beta\((\alpha,\beta)\) então, \(\mbox{E}(X)=\alpha/(\alpha+\beta)\), o que significa que nesta amostra a esperança é 0.5 e a variância é \(\mbox{Var}(X)=\alpha\beta/\big((\alpha+\beta)^2(\alpha+\beta-1)\big)\), logo 0.006410256 e desvio padrão 0.08006407.

x = dados[dados$Distr=="Normal",1]
goftest::cvm.test(x,"pnorm",mean=mean(x),sd=sd(x),estimated = TRUE)
## 
##  Cramer-von Mises test of goodness-of-fit
##  Braun's adjustment using 100 groups
##  Null hypothesis: Normal distribution
##  with parameters mean = 0.500642926430135, sd = 0.0779263329113196
##  Parameters assumed to have been estimated from data
## 
## data:  x
## omega2max = 0.52116, p-value = 0.9716
x = dados[dados$Distr=="Beta",1]
goftest::cvm.test(x,"pnorm",mean=mean(x),sd=sd(x),estimated = TRUE)
## 
##  Cramer-von Mises test of goodness-of-fit
##  Braun's adjustment using 100 groups
##  Null hypothesis: Normal distribution
##  with parameters mean = 0.498814202242074, sd = 0.0773226882183705
##  Parameters assumed to have been estimated from data
## 
## data:  x
## omega2max = 1.0029, p-value = 0.2092

Este resultado também é desconsertante! continuamos aceitamos ambas amostras como normais.

DescTools::JarqueBeraTest(dados[dados$Distr=="Normal",1])
## 
##  Robust Jarque Bera Test
## 
## data:  dados[dados$Distr == "Normal", 1]
## X-squared = 1.2778, df = 2, p-value = 0.5279
DescTools::JarqueBeraTest(dados[dados$Distr=="Beta",1])
## 
##  Robust Jarque Bera Test
## 
## data:  dados[dados$Distr == "Beta", 1]
## X-squared = 8.9491, df = 2, p-value = 0.0114


Este teste está de acordo com nossas amostras.

Dispoibilizamos aqui um artigo que faz um comparativo do poder de decisão de 40 testes de normalidade segundo o tamanho da amostra até 1024 e considerando a simetria da amostra.

3 Amostras grandes

Uma proposta de teste de normalidade em amostras grandes seria utilizar um dos testes de normalidade, desenvolvidos para amostras medianas, e aplicá-los a diversas amostras pequenas extraídas da amostra original, considerada grande.

3.1 Teste Shapiro-Wilks

O teste de Shapiro–Wilk testa a hipótese nula de que uma amostra \(X_1,\cdots,X_n\) veio de uma população normalmente distribuída. Foi publicado em 1965 por Samuel Sanford Shapiro e Martin Wilk. A estatística do teste é \[ W=\dfrac{\displaystyle \Big(\sum_{i=1}^n a_i x_{(i)}\Big)^2}{\displaystyle \sum_{i=1}^n (x_i-\overline{x})^2}, \] onde \(x_{(i)}\) são as estatísticas amostrais de ordem.

Os coeficientes \(a_i\) são dados como \[ (a_1,\cdots,a_n)=\dfrac{m^\top V^{-1}}{C}, \] sendo \(C\) a norma do vetor \[ C=|| V^{-1}m|| = (m^\top V^{-1}V^{-1}m)^{1/2} \] e o vetor \(m=(m_1,\cdots,m_n)^\top\) são os valores esperados das estatísticas de ordem de variáveis aleatórias independentes e distribuídas de forma idêntica, amostradas da distribuição normal padrão; finalmente, \(V\) é a matriz de covariância dessas estatísticas de ordem normal.

Neste caso, o tamanho da amostra é limitado entre 3 e 5000.

Verifiquemos a normalidade da amostra original, de tamanho 10000. Para isso extraímos um milhão de amostras de tamanho 1000 e calculamos a taxa de rejeição. Consideramos que se a taxa de rejeição for menor de 0.05 então a amostra a original deve ser Normal, caso contrário, rejeitamos a hipótese nula e a distribuição dos dados não é Normal.

start_time <- Sys.time()
pv.normal <- replicate( 10^6, shapiro.test(sample(dados[dados$Distr=="Normal",1],1000))$p.value )
end_time <- Sys.time()
end_time - start_time # minutos
## Time difference of 7.965397 mins
mean(pv.normal <= 0.05)
## [1] 0.030723


Com estes resultados aceitamos a hipótese de normalidade da amostra normal original.

start_time <- Sys.time()
pv.beta <- replicate( 10^6, shapiro.test(sample(dados[dados$Distr=="Beta",1],1000))$p.value )
end_time <- Sys.time()
end_time - start_time # minutos
## Time difference of 7.935529 mins
mean(pv.beta <= 0.05)
## [1] 0.058009


Com estes resultados rejeitamos a hipótese de normalidade da amostra beta original.