Ist Shapiro-Wilk der beste Normalitätstest? Warum ist es vielleicht besser als andere Tests wie Anderson-Darling?

Ich habe irgendwo in der Literatur gelesen, dass der Shapiro-Wilk-Test als der beste Normalitätstest angesehen wird, weil bei einem gegebenen Signifikanzniveau, , die Wahrscheinlichkeit, die Nullhypothese abzulehnen, wenn sie falsch ist, höher ist als im Fall des anderen Normalitätstests.$\alpha$

Können Sie mir bitte mit mathematischen Argumenten erklären, wie genau dies im Vergleich zu einigen anderen Normalitätstests (etwa dem Anderson-Darling-Test) funktioniert?

Zunächst ein allgemeiner Kommentar: Beachten Sie, dass der Anderson-Darling-Test für vollständig spezifizierte Verteilungen gilt, während der Shapiro-Wilk für Normalen mit jedem Mittelwert und jeder Varianz gilt. Wie in D'Agostino & Stephens [ 1 ] erwähnt, passt sich der Anderson-Darling auf sehr bequeme Weise an den Schätzungsfall an, ist mit dem Lilliefors-Test vergleichbar (aber konvergiert schneller und wird auf eine Weise modifiziert, die einfacher zu handhaben ist als dieser der Fall Kolmogorov-Smirnov). Insbesondere sind bei dem normalen Wert von n = 5 Tabellen des asymptotischen Wertes von A ∗ = A 2 ( 1 + 4$^{[1]}$$n=5$kann verwendet werden (keine Prüfung der Anpassungsgüte für n <5).$A^*=A^2\left(1+\frac{4}{n}-\frac{25}{n^2}\right)$

Ich habe irgendwo in der Literatur gelesen, dass der Shapiro-Wilk-Test als der beste Normalitätstest angesehen wird, weil bei einem gegebenen Signifikanzniveau α die Wahrscheinlichkeit, die Nullhypothese zu verwerfen, wenn sie falsch ist, höher ist als bei der anderen Normalität Tests.

Generell ist dies falsch.

Welche Normalitätstests "besser" sind, hängt davon ab, für welche Klassen von Alternativen Sie sich interessieren. Ein Grund für die Beliebtheit des Shapiro-Wilk ist, dass er in der Regel eine sehr gute Leistung bei einer Vielzahl nützlicher Alternativen aufweist. Es taucht in vielen Machtstudien auf und funktioniert normalerweise sehr gut, aber es ist nicht allgemein am besten.

Es ist ziemlich einfach, Alternativen zu finden, unter denen es weniger mächtig ist.

Zum Beispiel hat es bei Alternativen mit leichtem Schwanz oft weniger Leistung als der Schülerspektrum (Vergleichen Sie sie zum Beispiel mit einem Normalitätstest für einheitliche Daten - bein=30hat ein Test auf der Basis vonueine Leistung von ungefähr 63% im Vergleich zu etwas mehr als 38% für den Shapiro Wilk).$u=\frac{\max(x)−\min(x)}{sd(x)}$$n=30$$u$

Der Anderson-Darling (angepasst an die Parameterschätzung) schneidet bei der doppelten Exponentialfunktion besser ab. Moment-Skewness ist besser gegen einige Skew-Alternativen.

Können Sie mir bitte mit mathematischen Argumenten erklären, wie genau dies im Vergleich zu einigen anderen Normalitätstests (etwa dem Anderson-Darling-Test) funktioniert?

Ich werde es allgemein erklären (wenn Sie genauere Details wünschen, sind die Originalarbeiten und einige der späteren Arbeiten, in denen sie besprochen werden, die beste Wahl):

Betrachten Sie einen einfacheren, aber eng verwandten Test, den Shapiro-Francia. Es ist effektiv eine Funktion der Korrelation zwischen der Ordnungsstatistik und der erwarteten Ordnungsstatistik unter Normalbedingungen (und als solche ein ziemlich direktes Maß dafür, wie gerade die Linie im normalen QQ-Diagramm ist). Wie ich mich erinnere, ist Shapiro-Wilk leistungsfähiger, da es auch die Kovarianzen zwischen den Ordnungsstatistiken berücksichtigt und einen besten linearen Schätzer von aus dem QQ-Plot erzeugt, der dann mit s skaliert wird . Wenn die Verteilung nicht normal ist, liegt das Verhältnis nicht nahe bei 1.$\sigma$$s$

Zum Vergleich: Anderson-Darling basiert wie Kolmogorov-Smirnov und Cramér-von Mises auf der empirischen CDF. Insbesondere basiert es auf gewichteten Abweichungen zwischen dem ECDF und dem theoretischen ECDF (die Varianzgewichtung macht es empfindlicher für Abweichungen im Heck).

Der Test von Shapiro und Chen [ 2 ] (1995) (basierend auf Abständen zwischen Ordnungsstatistiken) zeigt oft etwas mehr Kraft als der Shapiro-Wilk (aber nicht immer); Sie verhalten sich oft sehr ähnlich.$^{[2]}$

-

Verwenden Sie den Shapiro Wilk, weil er häufig leistungsstark und allgemein verfügbar ist und vielen Menschen bekannt ist (ohne die Notwendigkeit, detailliert zu erklären, was er ist, wenn Sie ihn in einem Papier verwenden). Verwenden Sie ihn nur nicht unter der Illusion, dass er ist "der beste Normalitätstest". Es gibt keinen einzigen besten Normalitätstest.

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

[2]: Chen, L. und Shapiro, S. (1995)
"Ein alternativer Test auf Normalität basierend auf normalisierten Abständen."
Journal of Statistical Computation and Simulation 53 , 269-287.

Der Vergleich, den Sie gelesen haben, enthielt eindeutig nicht SnowsPenultimateNormalityTest ( http://cran.r-project.org/web/packages/TeachingDemos/TeachingDemos.pdf ), da er über alle Alternativen hinweg die höchstmögliche Leistung bietet. Es sollte also als "am besten" bewertet werden, wenn nur die Leistung in Betracht gezogen wird.

Ich stimme jedoch dem Kommentar von Nick Cox zu, dass der beste Test eher eine Handlung als ein formaler Test ist, da die Frage nach "Normal genug" viel wichtiger ist als "Genau normal". Wenn Sie einen aussagekräftigen Test wünschen, würde ich vorschlagen, den qq-Plot mit der in diesem Artikel beschriebenen Methode zu kombinieren:

Buja, A., Cook, D. Hofmann, H., Lawrence, M. Lee, E.-K., Swayne, DF und Wickham, H. (2009) Statistical Inference für die explorative Datenanalyse und Modelldiagnose Phil. Dr. Trans. R. Soc. A 2009 367, 4361-4383 doi: 10.1098 / rsta.2009.0120

Eine Implementierung davon ist die vis.testFunktion im TeachingDemos-Paket für R (dasselbe Paket wie SnowsPenultimateNormalityTest).

Ich komme zu spät zur Party, werde aber mit Verweisen auf die veröffentlichten, von Experten überprüften Forschungsergebnisse antworten. Der Grund, warum ich die Frage von OP nicht mit Ja / Nein beantworte, ist, dass sie komplizierter ist, als es scheint. Es gibt keinen Test, der für Stichproben mit oder ohne Ausreißer am leistungsfähigsten wäre. Ausreißer können die Leistung eines Tests erheblich verringern und für einen anderen Test erhöhen. Einige Tests funktionieren besser, wenn die Probe aus einer symmetrischen Verteilung stammt usw.

• Henry C. Thode, Testen auf Normalität , 2002 - Dies ist das umfassendste Buch zu diesem Thema. Wenn ich es auf eine einfache Antwort beschränken musste, ist SW in allen Fällen nicht mächtiger als AD. Hier zwei Auszüge für Ihr Lesevergnügen.

Ab Abschnitt 7.1.5: Auf der Grundlage der Leistung steht die Wahl der Prüfung in direktem Zusammenhang mit den verfügbaren Informationen oder den Annahmen, die bezüglich der Alternative getroffen wurden. Je spezifischer die Alternative ist, desto spezifischer und leistungsfähiger wird der Test in der Regel sein. Dies führt auch zu den zuverlässigsten Empfehlungen.

und

$K_s^2$$A^2$

• Romao, Xavier, Raimundo Delgado und Anibal Costa. "Ein empirischer Leistungsvergleich von univariaten Anpassungstests auf Normalität." Journal of Statistical Computation and Simulation 80.5 (2010): 545-591. Dies ist die neueste veröffentlichte Forschung zu dem Thema, von dem ich weiß.

Die Studie befasst sich mit der Durchführung von 33 Normalitätstests für verschiedene Stichprobengrößen unter Berücksichtigung verschiedener Signifikanzniveaus und für eine Reihe von symmetrischen, asymmetrischen und modifizierten Normalverteilungen. Die aus der Studie resultierenden allgemeinen Empfehlungen für Normalitätstests richten sich nach der Art der Nichtnormalität

Wenn Sie ihre Forschung wirklich auf Ja / Nein reduzieren möchten, lautet die Antwort JA. Shapiro-Wilks-Test scheint in den meisten Fällen etwas leistungsfähiger zu sein als Anderson-Darling. Sie empfehlen den Shapiro Wilk-Test, wenn Sie keine bestimmte alternative Verteilung im Auge haben. Wenn Sie sich jedoch für dieses Thema interessieren, ist das Papier lesenswert. Schauen Sie sich wenigstens die Tische an.

• Edith Seier, Normalitätstests: Leistungsvergleich , in International Encyclopedia of Statistical Science, 2014 - Ein Überblick über veröffentlichte Forschungsergebnisse zu diesem Thema. Auch hier hängt die Antwort von der Stichprobe und Ihrem Wissen über die alternative Verteilung ab, aber die trivialisierte Antwort wäre JA, Shapiro-Wilk ist normalerweise mächtiger, aber nicht immer.

• Henry C. Thode, Normalitätstests , in International Encyclopedia of Statistical Science, 2014 - Beschreibung der gängigen Normalitätstests. Seine Empfehlung:

$A^2$

Hier ging es um univariate Tests. Die Thode (2002) hat auch multivariate Tests, zensierte Daten, normale Gemische, Tests in Gegenwart von Ausreißern und vieles mehr.

Eine ernstere Antwort auf diese Frage und insbesondere das anhaltende Interesse von @ silverfish. Ein Ansatz zur Beantwortung derartiger Fragen besteht darin, einige zu vergleichende Simulationen durchzuführen. Nachfolgend finden Sie einen R-Code, der Daten unter verschiedenen Alternativen simuliert und mehrere der Normalitätstests durchführt und die Leistung vergleicht (sowie ein Konfidenzintervall für die Leistung, da die Leistung durch Simulation geschätzt wird). Ich habe die Stichprobengröße etwas angepasst, da es nicht interessant war, wenn viele der Potenzen nahe 100% oder 5% lagen. Ich fand runde Zahlen, die Potenzen nahe 80% ergaben. Jeder, der interessiert ist, kann diesen Code leicht nehmen und für andere Annahmen, andere Alternativen usw. ändern.

Sie sehen, dass es Alternativen gibt, bei denen einige Tests besser und andere schlechter abschneiden. Die wichtige Frage ist dann, welche Alternativen für Ihre wissenschaftlichen Fragestellungen am realistischsten sind. Dies sollte wirklich durch eine Simulation der Auswirkung der Arten von interessierender Nichtnormalität auf andere durchgeführte Tests ergänzt werden. Einige dieser Arten von Nicht-Normalität wirken sich stark auf andere Tests auf normaler Basis aus, andere wirken sich nicht sehr auf sie aus.

> 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