How to calculate critical values ($W_{\alpha}$) for the Shapiro-Wilk test?

Question

When performing the Shapiro-Wilk test, we can obtain the critical $W_{\alpha}$ statistic from tables, given $\alpha$ is 0.05 or 0.01.

If $\alpha$ is nonstandard, say 0.07 or 0.02 for example, how we can calculate the value of $W_{\alpha}$?

Answer 1

I think you're asking how to calculate a critical value for the Shapiro-Wilk test at significance levels other than the usual tabulated ones (I mention this because it's possible you're really interested in how to compute a p-value, which is a closely related issue, though now Chris and I have edited your question, my doubt looks a bit odd.)

If your desired significance level is between tabulated ones, such as wanting 2% and having 5% and 1%, you could use interpolation, which will at least be approximately suitable. It looks to me like $\Phi(1-p)$ is close to linear* in $\log(n(1-W))$ for $n\ge 12$, and even down at $n=5$ it's locally close enough that linear interpolation should work fine if done on those scales.

More generally, computer software is the most obvious choice. Some software offers direct calculation of p-values for the Shapiro-Wilk that may avoid the need to use critical values at all.

Finally, simulation is an option; one can simulate the statistic and hence obtain simulations from the distribution under the null; this allows one to compute estimates of quantiles of the distribution.

* Edit: looking at the p-value code in R, that shift in nearness-to-linearity between n=11 and n=12 is because R is using a different approximation to compute p-values for $n$ below 12; that shouldn't affect the suitability of linear interpolation at say 12, but it does suggest to me that it appearing to be very close to linear is more that the approximation is close to linear for $n\ge 12$; the actual transformed distribution is probably somewhat less linear down that low. [Now I look, the help even gives the information that a different approximation is used below 12.]

R also offers the following references, in case you want to write your own code:

[1] Patrick Royston (1982) An extension of Shapiro and Wilk's W test for normality to large samples. Applied Statistics, 31, 115-124.

[2] Patrick Royston (1982) Algorithm AS 181: The W test for Normality. Applied Statistics, 31, 176-180.

[3] Patrick Royston (1995) Remark AS R94: A remark on Algorithm AS 181: The W test for normality. Applied Statistics, 44, 547-551.

Answer 2

I think this depends quite a lot on what you intend to do, and to what detail you would like to define your answer. The $P$-value of the Shapiro-Wilk test is not a trivial thing to compute unfortunately. I'll give a simple answer to identify specific levels of $W_{\alpha}$ with a small degree of error, then we'll trace how you can actually calculate the $P$-value yourself.

Quick and dirty method

Assuming that you are working in R, we can simply generate a number of data sets and test them until we find the actual value of $W$ that gives a particular $\alpha$ value. Note that $W_{\alpha}$ is dependant upon the sample size that you are using, so you will have to perform this method with your own sample size $n$.

We're going to basically generate a normally distributed random sequence of numbers of size $n$ with the function rnorm(n), then calculate the $W$ statistic from R until we find one that we like.

find.W <- function(alpha = 0.05, error = 0.000001, n = 100){
  not.done <- TRUE
  while(not.done){
    a <- shapiro.test(rnorm(n))
    if(a$p.value < alpha+error && a$p.value > alpha-error){
      not.done = FALSE
      W <- a$statistic
    }
  }
  return(W)
}

Doing this once we find that for sample size 100, we find that

w <- find.W(alpha = 0.07, error = 0.001, n = 100); w
0.9765172 

This is of course stochastic. If we repeat this 1000 times to get a better idea of the real number:

w <- vector()
for(i in 1:1000){
  w[i] <- find.W(alpha = 0.07, error = 0.001, n = 100)
}
mean(w); sd(w)

[1] 0.9764422
[1] 4.544819e-05

We find that for $\alpha = 0.07$, $W_{\alpha} \approx 0.9764$. If we look into it a little bit further, we find that this is only an approximation of the true $W$ statistic. To find the real value $W$ value, let's trace how this one is calculated.

Theoretical value of $W$

If we look into the source code for R, we can pretty easily find how the $W$ statistic is calculated. Looking at the source of the shapiro.test() function, we find that it calls a C file named S_Wilks.c. Looking into it, we find the source code here [1]. Inside the code, a paper gets referenced with the theoretical aspect. After a little digging, it appears that the approximation method comes from this [2,3] paper. You can read this yourself if you would like to know the actually theoretical aspects to manually calculated $P$-values. The methods used in the original Shapiro & Wilk 1965 paper are from an unpublished manuscript that has since been lost to time (at least, I can't find it).

[1] Source code which R uses: http://www.matrixscience.com/msparser/help/group___shapiro_wilk_source_code.html

[2] Approximation of $W$ paper: Some Techniques for Assessing Multivarate Normality Based on the Shapiro- Wilk W, J. P. Royston, Journal of the Royal Statistical Society. Series C (Applied Statistics), Vol. 32, No. 2 (1983), pp. 121-133

[3] Approximation of $W$ technique: Statistical Algorithms: Algorithm AS, J. P. Royston, Journal of the Royal Statistical Society. Series C (Applied Statistics), Vol. 32, No. 2 (1983), pp. 176-180

Archive on JSTOR found here for that particular issue.