¿Es Shapiro-Wilk la mejor prueba de normalidad? ¿Por qué podría ser mejor que otras pruebas como Anderson-Darling?

He leído en alguna parte de la literatura que la prueba de Shapiro-Wilk se considera la mejor prueba de normalidad porque, para un nivel de significancia dado, , la probabilidad de rechazar la hipótesis nula si es falsa es mayor que en el caso de la otra normalidad. pruebas$\alpha$

¿Podría explicarme, utilizando argumentos matemáticos si es posible, cómo funciona exactamente en comparación con algunas de las otras pruebas de normalidad (digamos la prueba de Anderson-Darling)?

Primero un comentario general: tenga en cuenta que la prueba de Anderson-Darling es para distribuciones completamente especificadas, mientras que Shapiro-Wilk es para normales con cualquier media y varianza. Sin embargo, como se señaló en D'Agostino & Stephens [ 1 ], Anderson-Darling se adapta de una manera muy conveniente al caso de estimación, similar a (pero converge más rápido y se modifica de una manera más fácil de manejar que la prueba de Lilliefors para el caso Kolmogorov-Smirnov). Específicamente, en la normal, por n = 5 , tablas del valor asintótico de A ∗ = A 2 ( 1 + 4$^{[1]}$$n=5$puede usarse (no pruebe la bondad de ajuste para n <5).$A^*=A^2\left(1+\frac{4}{n}-\frac{25}{n^2}\right)$

He leído en alguna parte de la literatura que la prueba de Shapiro-Wilk se considera la mejor prueba de normalidad porque, para un nivel de significancia dado, α, la probabilidad de rechazar la hipótesis nula si es falsa es mayor que en el caso de la otra normalidad. pruebas

Como una declaración general, esto es falso.

Las pruebas de normalidad que sean "mejores" dependen de las clases de alternativas que le interesen. Una de las razones por las que Shapiro-Wilk es popular es que tiende a tener muy buen poder en una amplia gama de alternativas útiles. Aparece en muchos estudios de poder, y generalmente funciona muy bien, pero no es universalmente el mejor.

Es bastante fácil encontrar alternativas bajo las cuales es menos poderoso.

Por ejemplo, frente a las alternativas de cola ligera, a menudo tiene menos potencia que el rango estudiado (compárelos en una prueba de normalidad en datos uniformes, por ejemplo, enn=30, una prueba basada enutiene un poder de aproximadamente el 63% en comparación con un poco más del 38% para el Shapiro Wilk).$u=\frac{\max(x)−\min(x)}{sd(x)}$$n=30$$u$

El Anderson-Darling (ajustado para la estimación de parámetros) funciona mejor en el doble exponencial. La asimetría de momento funciona mejor frente a algunas alternativas de asimetría.

¿Podría explicarme, utilizando argumentos matemáticos si es posible, cómo funciona exactamente en comparación con algunas de las otras pruebas de normalidad (digamos la prueba de Anderson-Darling)?

Explicaré en términos generales (si desea detalles más específicos, los documentos originales y algunos de los documentos posteriores que los discuten serían su mejor opción):

Considere una prueba más simple pero estrechamente relacionada, la Shapiro-Francia; es efectivamente una función de la correlación entre las estadísticas de orden y las estadísticas de orden esperadas bajo normalidad (y como tal, una medida bastante directa de "qué tan recta es la línea" en el gráfico QQ normal). Como recuerdo, el Shapiro-Wilk es más poderoso porque también tiene en cuenta las covarianzas entre las estadísticas de orden, produciendo un mejor estimador lineal de del gráfico QQ, que luego se escala por s . Cuando la distribución está lejos de ser normal, la relación no es cercana a 1.$\sigma$$s$

En comparación, el Anderson-Darling, como el Kolmogorov-Smirnov y el Cramér-von Mises, se basa en el CDF empírico. Específicamente, se basa en desviaciones ponderadas entre ECDF y ECDF teórico (la ponderación por varianza lo hace más sensible a las desviaciones en la cola).

La prueba de Shapiro y Chen [ 2 ] (1995) (basada en los espacios entre las estadísticas de orden) a menudo exhibe un poco más de poder que el Shapiro-Wilk (pero no siempre); a menudo se desempeñan de manera muy similar.$^{[2]}$

-

Use Shapiro Wilk porque a menudo es poderoso, está ampliamente disponible y muchas personas están familiarizadas con él (eliminando la necesidad de explicar en detalle qué es si lo usa en un documento), simplemente no lo use bajo la ilusión de que es "La mejor prueba de normalidad". No hay una mejor prueba de normalidad.

[1]: D'Agostino, RB y Stephens, MA (1986)
Goodness of Fit Techniques ,
Marcel Dekker, Nueva York.

[2]: Chen, L. y Shapiro, S. (1995)
"Una prueba alternativa para la normalidad basada en espaciamientos normalizados".
Revista de cálculo y simulación estadística 53 , 269-287.

Claramente, la comparación que leyó no incluyó SnowsPenultimateNormalityTest ( http://cran.r-project.org/web/packages/TeachingDemos/TeachingDemos.pdf ) ya que tiene el mayor poder posible en todas las alternativas. Por lo tanto, debe considerarse "Mejor" si el poder es la única consideración (tenga en cuenta que mis opiniones están claramente sesgadas, pero documentadas en el enlace / documentación).

Sin embargo, estoy de acuerdo con el comentario de Nick Cox de que la mejor prueba es una trama en lugar de una prueba formal, ya que la cuestión de "lo suficientemente normal" es mucho más importante que "exactamente normal". Si desea una prueba significativa, sugeriría combinar la gráfica qq con la metodología en este documento:

Buja, A., Cook, D. Hofmann, H., Lawrence, M. Lee, E.-K., Swayne, DF y Wickham, H. (2009) Inferencia estadística para análisis de datos exploratorios y diagnóstico de modelos Phil. Trans. R. Soc. A 2009 367, 4361-4383 doi: 10.1098 / rsta.2009.0120

Una implementación de eso es la vis.testfunción en el paquete TeachingDemos para R (mismo paquete que SnowsPenultimateNormalityTest).

Llego tarde a la fiesta, pero responderé con referencias a la investigación publicada revisada por pares. La razón por la que no respondo Sí / No a la pregunta de OP es que es más complicado de lo que parece. No hay una prueba que sea la más poderosa para muestras provenientes de cualquier distribución con o sin valores atípicos. Los valores atípicos pueden disminuir severamente el poder de una prueba y aumentar para otra. Algunas pruebas funcionan mejor cuando la muestra proviene de distribución simétrica, etc.

• Henry C. Thode, Testing for Normality , 2002 - Este es el libro más completo sobre el tema. Si tuviera que simplificarlo a una respuesta simple, entonces SW no es más poderoso que AD en todos los casos. Aquí hay dos extractos para su placer de lectura.

De la sección 7.1.5: Sobre la base del poder, la elección de la prueba está directamente relacionada con la información disponible o las suposiciones relacionadas con la alternativa. Cuanto más específica sea la alternativa, más específica y poderosa será la prueba; Esto también dará como resultado las recomendaciones más confiables.

y

$K_s^2$$A^2$

• Romao, Xavier, Raimundo Delgado y Anibal Costa. "Una comparación de poder empírico de pruebas univariadas de bondad de ajuste para la normalidad". Journal of Statistical Computation and Simulation 80.5 (2010): 545-591. Esta es la investigación publicada más reciente sobre el tema que conozco.

El estudio aborda el desempeño de 33 pruebas de normalidad, para varios tamaños de muestra, considerando varios niveles de significancia y para una serie de distribuciones normales simétricas, asimétricas y modificadas. Las recomendaciones generales para las pruebas de normalidad resultantes del estudio se definen de acuerdo con la naturaleza de la no normalidad.

Si realmente quiere reducir su investigación a sí / no, entonces la respuesta es SÍ. La prueba de Shapiro-Wilks parece ser un poco más poderosa en la mayoría de los casos que Anderson-Darling. Recomiendan la prueba de Shapiro Wilk cuando no tiene una distribución alternativa particular en mente. Sin embargo, si está interesado en este tema, vale la pena leer el documento. Al menos mira las mesas.

• Edith Seier, Pruebas de normalidad: Comparación de potencia , en International Encyclopedia of Statistical Science, 2014 - Una encuesta de investigación publicada sobre el tema. Una vez más, la respuesta depende de la muestra y de su conocimiento sobre la distribución alternativa, pero la respuesta trivializada sería SÍ, Shapiro-Wilk suele ser más poderoso, pero no siempre.

• Henry C. Thode, Pruebas de normalidad , en International Encyclopedia of Statistical Science, 2014 - Descripción de las pruebas de normalidad populares. Su recomendación:

$A^2$

Ahora, esto se trataba de pruebas univariadas. Thode (2002) también tiene pruebas multivariadas, datos censurados, mezclas normales, pruebas en presencia de valores atípicos y mucho más.

Una respuesta más seria para avanzar en esta pregunta y especialmente el continuo interés de @ silverfish. Un enfoque para responder preguntas como esta es ejecutar algunas simulaciones para comparar. A continuación se muestra un código R que simula datos bajo varias alternativas y realiza varias de las pruebas de normalidad y compara la potencia (y un intervalo de confianza en la potencia ya que la potencia se estima mediante simulación). Ajusté un poco los tamaños de muestra porque no era interesante cuando muchas de las potencias estaban cerca del 100% o 5%, encontré números redondos que daban potencias cerca del 80%. Cualquier persona interesada podría tomar fácilmente este código y modificarlo para diferentes supuestos, diferentes alternativas, etc.

Puede ver que existen alternativas para las cuales algunas de las pruebas funcionan mejor y otras en las que lo hacen peor. La pregunta importante es, entonces, qué alternativas son más realistas para sus preguntas / área científica. Esto realmente debe seguirse con una simulación del efecto de los tipos de no normalidad de interés en otras pruebas que se realizan. Algunos de estos tipos de no normalidad afectan en gran medida a otras pruebas basadas en la normalidad, otras no las afectan demasiado.

> library(nortest)
> 
> simfun1 <- function(fun=function(n) rnorm(n), n=250) {
+   x <- fun(n)
+   c(sw=shapiro.test(x)$p.value, sf=sf.test(x)$p.value, ad=ad.test(x)$p.value,
+     cvm=cvm.test(x)$p.value, lillie=lillie.test(x)$p.value, 
+     pearson=pearson.test(x)$p.value, snow=0)
+ }
> 
> ### Test size using null hypothesis near true
> 
> out1 <- replicate(10000, simfun1())
> apply(out1, 1, function(x) mean(x<=0.05))
     sw      sf      ad     cvm  lillie pearson    snow 
 0.0490  0.0520  0.0521  0.0509  0.0531  0.0538  1.0000 
> apply(out1, 1, function(x) prop.test(sum(x<=0.05),length(x))$conf.int)  #$
             sw         sf         ad        cvm     lillie    pearson      snow
[1,] 0.04489158 0.04776981 0.04786582 0.04671398 0.04882619 0.04949870 0.9995213
[2,] 0.05345887 0.05657820 0.05668211 0.05543493 0.05772093 0.05844785 1.0000000
> 
> ### Test again with mean and sd different
> 
> out2 <- replicate(10000, simfun1(fun=function(n) rnorm(n,100,5)))
> apply(out2, 1, function(x) mean(x<=0.05))
     sw      sf      ad     cvm  lillie pearson    snow 
 0.0482  0.0513  0.0461  0.0477  0.0515  0.0506  1.0000 
> apply(out2, 1, function(x) prop.test(sum(x<=0.05),length(x))$conf.int)  #$
             sw         sf         ad        cvm     lillie    pearson      snow
[1,] 0.04412478 0.04709785 0.04211345 0.04364569 0.04728982 0.04642612 0.9995213
[2,] 0.05262633 0.05585073 0.05043938 0.05210583 0.05605860 0.05512303 1.0000000
> 
> #### now for the power under different forms of non-normality
> 
> ## heavy tails, t(3)
> rt3 <- function(n) rt(n, df=3)
> 
> out3 <- replicate(10000, simfun1(fun=rt3, n=75))
There were 50 or more warnings (use warnings() to see the first 50)
> round(apply(out3, 1, function(x) mean(x<=0.05, na.rm=TRUE)),3)
     sw      sf      ad     cvm  lillie pearson    snow 
  0.788   0.831   0.756   0.726   0.624   0.440   1.000 
> round(apply(out3, 1, function(x){ 
+ prop.test(sum(x<=0.05,na.rm=TRUE),sum(!is.na(x)))$conf.int),3) }  #$
        sw    sf    ad   cvm lillie pearson snow
[1,] 0.780 0.824 0.748 0.717  0.614   0.431    1
[2,] 0.796 0.838 0.765 0.734  0.633   0.450    1
> 
> 
> ## light tails, uniform
> u <- function(n) runif(n)
> 
> out4 <- replicate(10000, simfun1(fun=u, n=65))
> round(apply(out4, 1, function(x) mean(x<=0.05, na.rm=TRUE)),3)
     sw      sf      ad     cvm  lillie pearson    snow 
  0.906   0.712   0.745   0.591   0.362   0.270   1.000 
> round(apply(out4, 1, function(x){ 
+ prop.test(sum(x<=0.05,na.rm=TRUE),sum(!is.na(x)))$conf.int),3) }  #$
        sw    sf    ad   cvm lillie pearson snow
[1,] 0.900 0.703 0.737 0.581  0.353   0.261    1
[2,] 0.911 0.720 0.754 0.600  0.372   0.279    1
> 
> ## double exponential, Laplace
> de <- function(n) sample(c(-1,1), n, replace=TRUE) * rexp(n)
> 
> out5 <- replicate(10000, simfun1(fun=de, n=100))
> round(apply(out5, 1, function(x) mean(x<=0.05, na.rm=TRUE)),3)
     sw      sf      ad     cvm  lillie pearson    snow 
  0.796   0.844   0.824   0.820   0.706   0.477   1.000 
> round(apply(out5, 1, function(x){ 
+ prop.test(sum(x<=0.05,na.rm=TRUE),sum(!is.na(x)))$conf.int),3) }  #$
        sw    sf    ad   cvm lillie pearson snow
[1,] 0.788 0.837 0.817 0.813  0.697   0.467    1
[2,] 0.804 0.851 0.832 0.828  0.715   0.486    1
> 
> ## skewed, gamma(2,2)
> g22 <- function(n) rgamma(n,2,2)
> 
> out6 <- replicate(10000, simfun1(fun=g22, n=50))
Warning message:
In cvm.test(x) :
  p-value is smaller than 7.37e-10, cannot be computed more accurately
> round(apply(out6, 1, function(x) mean(x<=0.05, na.rm=TRUE)),3)
     sw      sf      ad     cvm  lillie pearson    snow 
  0.954   0.930   0.893   0.835   0.695   0.656   1.000 
> round(apply(out6, 1, function(x){ 
+ prop.test(sum(x<=0.05,na.rm=TRUE),sum(!is.na(x)))$conf.int),3) }  #$
        sw    sf    ad   cvm lillie pearson snow
[1,] 0.950 0.925 0.886 0.827  0.686   0.646    1
[2,] 0.958 0.935 0.899 0.842  0.704   0.665    1
> 
> ## skewed, gamma(2,2)
> g99 <- function(n) rgamma(n,9,9)
> 
> out7 <- replicate(10000, simfun1(fun=g99, n=150))
> round(apply(out7, 1, function(x) mean(x<=0.05, na.rm=TRUE)),3)
     sw      sf      ad     cvm  lillie pearson    snow 
  0.844   0.818   0.724   0.651   0.526   0.286   1.000 
> round(apply(out7, 1, function(x){ 
+ prop.test(sum(x<=0.05,na.rm=TRUE),sum(!is.na(x)))$conf.int),3) }  #$
        sw    sf    ad   cvm lillie pearson snow
[1,] 0.836 0.810 0.715 0.642  0.516   0.277    1
[2,] 0.851 0.826 0.732 0.660  0.536   0.294    1
> 
> ## tails normal, middle not
> mid <- function(n) {
+   x <- rnorm(n)
+   x[ x > -0.5 & x < 0.5 ] <- 0
+   x
+ }
> 
> out9 <- replicate(10000, simfun1(fun=mid, n=30))
Warning messages:
1: In cvm.test(x) :
  p-value is smaller than 7.37e-10, cannot be computed more accurately
2: In cvm.test(x) :
  p-value is smaller than 7.37e-10, cannot be computed more accurately
> round(apply(out9, 1, function(x) mean(x<=0.05, na.rm=TRUE)),3)
     sw      sf      ad     cvm  lillie pearson    snow 
  0.374   0.371   0.624   0.739   0.884   0.948   1.000 
> round(apply(out9, 1, function(x){ 
+ prop.test(sum(x<=0.05,na.rm=TRUE),sum(!is.na(x)))$conf.int),3) }  #$
        sw    sf    ad   cvm lillie pearson snow
[1,] 0.365 0.362 0.614 0.730  0.878   0.943    1
[2,] 0.384 0.381 0.633 0.747  0.890   0.952    1
> 
> ## mixture on variance
> mv <- function(n, p=0.1, sd=3) {
+   rnorm(n,0, ifelse(runif(n)<p, sd, 1))
+ }
> 
> out10 <- replicate(10000, simfun1(fun=mv, n=100))
Warning message:
In cvm.test(x) :
  p-value is smaller than 7.37e-10, cannot be computed more accurately
> round(apply(out10, 1, function(x) mean(x<=0.05, na.rm=TRUE)),3)
     sw      sf      ad     cvm  lillie pearson    snow 
  0.800   0.844   0.682   0.609   0.487   0.287   1.000 
> round(apply(out10, 1, function(x){ 
+ prop.test(sum(x<=0.05,na.rm=TRUE),sum(!is.na(x)))$conf.int),3) }  #$
        sw    sf    ad   cvm lillie pearson snow
[1,] 0.792 0.837 0.673 0.599  0.477   0.278    1
[2,] 0.808 0.851 0.691 0.619  0.497   0.296    1
> 
> ## mixture on mean
> mm <- function(n, p=0.3, mu=2) {
+   rnorm(n, ifelse(runif(n)<p, mu, 0), 1)
+ }
> 
> out11 <- replicate(10000, simfun1(fun=mm, n=400))
> round(apply(out11, 1, function(x) mean(x<=0.05, na.rm=TRUE)),3)
     sw      sf      ad     cvm  lillie pearson    snow 
  0.776   0.710   0.808   0.788   0.669   0.354   1.000 
> round(apply(out11, 1, function(x){ 
+ prop.test(sum(x<=0.05,na.rm=TRUE),sum(!is.na(x)))$conf.int),3) }  #$
        sw    sf    ad   cvm lillie pearson snow
[1,] 0.768 0.701 0.801 0.780  0.659   0.344    1
[2,] 0.784 0.719 0.816 0.796  0.678   0.363    1