24
$\begingroup$

I have read somewhere in the literature that the Shapiro–Wilk test is considered to be the best normality test because for a given significance level, $\alpha$, the probability of rejecting the null hypothesis if it's false is higher than in the case of the other normality tests.

Could you please explain to me, using mathematical arguments if possible, how exactly it works compared to some of the other normality tests (say the Anderson–Darling test)?

$\endgroup$
  • 4
    $\begingroup$ Note that the power depends on the way in which the null hypothesis is false, which for a general-purpose goodness-of-fit test can be any of innumerable ways. Without having checked I'd still bet that each of the common normality tests is most powerful against certain alternatives. $\endgroup$ – Scortchi Mar 20 '14 at 11:34
  • 5
    $\begingroup$ Not the answer you seek, perhaps, but I'd say that the best normality test is a normal probability plot, i.e. a quantile-quantile plot of observed values versus normal quantiles. The Shapiro-Wilk test is indeed often commended, but it can't tell you exactly how your data differ from a normal. Often unimportant differences are flagged by the test, because they do qualify as significant for large sample sizes, and the opposite problem can also bite you. $\endgroup$ – Nick Cox Mar 20 '14 at 11:51
  • $\begingroup$ Maybe you read it in Normality Tests for Statistical Analysis: A Guide for Non-Statisticians $\endgroup$ – lithic Jun 17 '16 at 16:29
17
+100
$\begingroup$

First a general comment: Note that the Anderson-Darling test is for completely specified distributions, while the Shapiro-Wilk is for normals with any mean and variance. However, as noted in D'Agostino & Stephens$^{[1]}$ the Anderson-Darling adapts in a very convenient way to the estimation case, akin to (but converges faster and is modified in a way that's simpler to deal with than the Lilliefors test for the Kolmogorov-Smirnov case). Specifically, at the normal, by $n=5$, tables of the asymptotic value of $A^*=A^2\left(1+\frac{4}{n}-\frac{25}{n^2}\right)$ may be used (don't be testing goodness of fit for n<5).

I have read somewhere in the literature that the Shapiro–Wilk test is considered to be the best normality test because for a given significance level, α, the probability of rejecting the null hypothesis if it's false is higher than in the case of the other normality tests.

As a general statement this is false.

Which normality tests are "better" depends on which classes of alternatives you're interested in. One reason the Shapiro-Wilk is popular is that it tends to have very good power under a broad range of useful alternatives. It comes up in many studies of power, and usually performs very well, but it's not universally best.

It's quite easy to find alternatives under which it's less powerful.

For example, against light tailed alternatives it often has less power than the studentized range $u=\frac{\max(x)−\min(x)}{sd(x)}$ (compare them on a test of normality on uniform data, for example - at $n=30$, a test based on $u$ has power of about 63% compared to a bit over 38% for the Shapiro Wilk).

The Anderson-Darling (adjusted for parameter estimation) does better at the double exponential. Moment-skewness does better against some skew alternatives.

Could you please explain to me, using mathematical arguments if possible, how exactly it works compared to some of the other normality tests (say the Anderson–Darling test)?

I will explain in general terms (if you want more specific details the original papers and some of the later papers that discuss them would be your best bet):

Consider a simpler but closely related test, the Shapiro-Francia; it's effectively a function of the correlation between the order statistics and the expected order statistics under normality (and as such, a pretty direct measure of "how straight the line is" in the normal Q-Q plot). As I recall, the Shapiro-Wilk is more powerful because it also takes into account the covariances between the order statistics, producing a best linear estimator of $\sigma$ from the Q-Q plot, which is then scaled by $s$. When the distribution is far from normal, the ratio isn't close to 1.

By comparison the Anderson-Darling, like the Kolmogorov-Smirnov and the Cramér-von Mises, is based on the empirical CDF. Specifically, it's based on weighted deviations between ECDF and theoretical ECDF (the weighting-for-variance makes it more sensitive to deviations in the tail).

The test by Shapiro and Chen$^{[2]}$ (1995) (based on spacings between order statistics) often exhibits slightly more power than the Shapiro-Wilk (but not always); they often perform very similarly.

--

Use the Shapiro Wilk because it's often powerful, widely available and many people are familiar with it (removing the need to explain in detail what it is if you use it in a paper) -- just don't use it under the illusion that it's "the best normality test". There isn't one best normality test.

[1]: D’Agostino, R. B. and Stephens, M. A. (1986)
Goodness of Fit Techniques,
Marcel Dekker, New York.

[2]: Chen, L. and Shapiro, S. (1995)
"An Alternative test for normality based on normalized spacings."
Journal of Statistical Computation and Simulation 53, 269-287.

$\endgroup$
  • $\begingroup$ My classmate told me:"If sample size >50,you should use Kolmogorov-Smirnov." Is that correct? $\endgroup$ – kittygirl Apr 4 at 3:18
  • $\begingroup$ No. To my recollection the original 1965 paper by Shapiro and Wilk only gave the required constants ($a_i$) used in the linear estimate of $\sigma$ for $n$ up to $50$ but that was over half a century ago. Things have moved on a little since then. Even without that, the Shapiro Francia or the Anderson-Darling (also adjusted for parameter estimation) are usually better choices; these often have considerably lower power against typically interesting alternatives. (& if you're estimating mean and sd from the sample, you're not strictly doing a Kolmogorov-Smirnov, but rather a Lilliefors test) $\endgroup$ – Glen_b Apr 4 at 4:20
  • $\begingroup$ In short, there was a brief period of a few years post 1967 (the initial publication of Lilliefors' work) where it might have been a justifiable piece of advice, but not for a long time since $\endgroup$ – Glen_b Apr 4 at 4:24
  • $\begingroup$ When sample size>5000, run shapiro.test in R will get error sample size must be between 3 and 5000.Then, what else test should be used? $\endgroup$ – kittygirl Apr 4 at 9:25
  • $\begingroup$ 1. At large n you'll almost always reject any simple distributional model (even when it's quite a suitable approximation); it may be more advisable to do something else (why are you testing normality?) 2. It's not really a matter of "should" about it; there's no single goodness of test that's always better than any other. It just happens that the Shapiro Wilk is reasonably good. However, a suitable alternative at large n is the Shapiro-Francia test. If you can find an implementation of the Chen-Shapiro test at large n (assuming there's a good reason to test at all), consider that instead. $\endgroup$ – Glen_b Apr 4 at 12:01
12
$\begingroup$

Clearly the comparison that you read did not include SnowsPenultimateNormalityTest (http://cran.r-project.org/web/packages/TeachingDemos/TeachingDemos.pdf) since it has the highest possible power across all alternatives. So it should be considered "Best" if power is the only consideration (Note that my opinions are clearly biased, but documented in the link/documentation).

However, I agree with Nick Cox's comment that the best test is a plot rather than a formal test since the question of "Normal enough" is much more important than "Exactly normal". If you want a meaningful test then I would suggest combining the qq plot with the methodology in this paper:

Buja, A., Cook, D. Hofmann, H., Lawrence, M. Lee, E.-K., Swayne, D.F and Wickham, H. (2009) Statistical Inference for exploratory data analysis and model diagnostics Phil. Trans. R. Soc. A 2009 367, 4361-4383 doi: 10.1098/rsta.2009.0120

One implementation of that is the vis.test function in the TeachingDemos package for R (same package as SnowsPenultimateNormalityTest).

$\endgroup$
  • 1
    $\begingroup$ I agree with the sentiment, but gazing at QQ plots isn't much a solution if one wants to do much of anything with more than 10-20 variables. $\endgroup$ – Andrew M Dec 20 '14 at 0:19
  • $\begingroup$ We print QQ plots in addition to normality tests. They're not exclusive, but complementary tools. $\endgroup$ – Aksakal Dec 20 '14 at 2:12
  • 2
    $\begingroup$ @Andrew M So, is looking through the results of 100 or 1000 or so Shapiro-Wilk or other tests what you prefer? It's not difficult to automate looking at Q-Q plots, say 25 at a time. Often a glance is enough to see (literally) a real problem. $\endgroup$ – Nick Cox Dec 20 '14 at 13:51
10
$\begingroup$

I'm late to the party, but will answer with references to the published peer-reviewed research. The reason why I don't answer Yes/No to OP's question is that it is more complicated than it may seem. There isn't one test which would be the most powerful for samples coming from any distribution with or without outliers. Outliers may severely diminish power of one test and increase for another. Some test work better when the sample comes from symmetrical distribution etc.

  • Henry C. Thode, Testing for Normality, 2002 - This is a the most comprehensive book on the subject. If I had to dumb it down to a simple answer, then SW is not more powerful than AD in all cases. Here are two excerpt for your reading pleasure.

From section 7.1.5: On the basis of power, the choice of test is directly related to the information available or the assumptions made concerning the alternative. The more specific the alternative, the more specific and more powerful the test will usually be; this will also result in the most reliable recommendations.

and

A joint skewness and kurtosis test such as $K_s^2$ provides high power against a wide range of alternatives, as does the Anderson-Darling $A^2$. Wilk-Shapiro W showed relatively high power among skewed and short-tailed symmetric alternatives when compared to other tests, and respectable power for long-tailed symmetric alternatives.

  • Romao, Xavier, Raimundo Delgado, and Anibal Costa. "An empirical power comparison of univariate goodness-of-fit tests for normality." Journal of Statistical Computation and Simulation 80.5 (2010): 545-591. This is the most recent published research on the subject I know of.

The study addresses the performance of 33 normality tests, for various sample sizes, considering several significance levels and for a number of symmetric, asymmetric and modified normal distributions. General recommendations for normality testing resulting from the study are defined according to the nature of the non-normality

If you really want to boil down their research to yes/no, then the answer is YES. Shapiro-Wilks test seems to be a little bit more powerful in most cases than Anderson-Darling. They recommend Shapiro Wilk test when you don't have a particular alternative distribution in mind. However, if you're interested in this subject, the paper is worth reading. At least look at the tables.

  • Edith Seier, Normality Tests: Power Comparison, in International Encyclopedia of Statistical Science, 2014 - A survey of published research on the subject. Again, the answer depends on the sample and your knowledge about the alternative distribution, but trivialized answer would be YES, Shapiro-Wilk is usually more powerful, but not always.

  • Henry C. Thode, Normality Tests, in International Encyclopedia of Statistical Science, 2014 - Description of popular normality tests. His recommendation:

As indicated previously, the number of normality tests is large, too large for even the majority of them to be mentioned here. Overall the best tests appear to be the moment tests, Shapiro–Wilk W, Anderson–Darling $A^2$ (see Anderson-Darling Tests of Goodness-of-Fit), and the Jarque–Bera test. Specifics on these and many other normality tests and their characteristics can be found in Thode (2002) and on general goodness of t issues, including normality tests, in D’Agostino and Stephens (1986).

Now, this was all about univariate tests. The Thode (2002) also has multivariate test, censored data, normal mixtures, testing in the presence of outliers, and much more.

$\endgroup$
9
$\begingroup$

A more serious answer to further this question and especially @silverfish's continued interest. One approach to answering questions like this is to run some simulations to compare. Below is some R code that simulates data under various alternatives and does several of the normality tests and compares the power (and a confidence interval on the power since power is estimated via simulation). I tweaked the sample sizes somewhat because it was not interesting when many of the powers were close to 100% or 5%, I found round numbers that gave powers near 80%. Anyone interested could easily take this code and modify it for different assumptions, different alternatives, etc.

You can see that there are alternatives for which some of the tests do better and others where they do worse. The important question is then which alternatives are most realistic for your scientific questions/area. This really should be followed up with a simulation of the effect of the types of non-normality of interest on other tests being done. Some of these types of non-normality greatly affect other normal based tests, others don't affect them much.

> 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
$\endgroup$
  • $\begingroup$ Thanks for adding this answer which is a nice complement to your earlier one. No disparagement was intended to your other answer, which is one of my favourites on SE! $\endgroup$ – Silverfish Dec 19 '14 at 22:23
  • 1
    $\begingroup$ I took the liberty of editing your code to add code formatting & fit everything within the window, Greg. I think it will be easier to read this way, but if you don't like it, roll it back w/ my apologies. $\endgroup$ – gung Dec 20 '14 at 0:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.