What is the best way to assess normality in a large dataset?

Question

I have a regression equation whose residuals are homoscedastic but have a slight negative skew due to several outlier values. The Q-Q plot is pretty close as to what would be expected from a normal distribution, but has a slight deviation in the upper quantile that seems to be related to these values. I am trying to test whether these outlier values cause the residuals of the regression to depart from normality enough that it is inappropriate to use the regression equation to predict the values of new data.

I originally tested for normality using a Kolmogorov-Smirnov test using the estimates mean and standard deviation of the sample as the comparative sample and that said the results did not depart from near normality, but as I read more about statistics I heard several people recommend that the Kolmogorov-Smirnov test should never be used. So I tried a Shapiro-Wilk test instead and the result came back as significant.

However, I suspect the reason the Shapiro-Wilk test produced a significant result was because of the large size (~350) of the dataset, as I've read that at larger sample sizes even small deviations from normality will cause the Shapiro-Wilk test to give a significant result, even if the data approximate normality. This is supported by the fact that even when these outlier values are removed the Shapiro-Wilk test still finds the residuals to significantly depart from normality. This is the reason I originally used a Kolmogorov-Smirnov test, because it is supposedly more robust to large sample sizes.

I realize my sample size is probably large enough that I don't have to be concerned about non-normality, but it would be nice to find a better way to say that the residuals are close enough to normal to work than simply invoking the Central Limit Theorem. Given this, what would be the best way to test if my data are "nearly normal"?

Answer 1

I agree with @Dave on on the graphical examination, but if you really want something quantitative, you could try Kullback-Leibler divergence of your observed distribution vs a normal distribution with the empirical mean and standard deviation plugged in. I suppose you could do some inference using the bootstrap.

Zero is perfect agreement. To interpret non-zero values, you may need to learn about nats (natural units of information)

I have never tried this.

Answer 2

I have usually been happy with results from the Shapiro-Wilk test. In large part, decisions have to be based on what qualifies as 'nearly normal' for the purpose at hand, but I have found that a Shapiro-Wilk test often rejects just when I judge the normal probability plot to show 'noticeable' departure from linearity.

Distributions $\mathsf{Beta}(n,n)$ tend to normal as $n$ increases. Below I look at samples of size $n = 5000$ for $n=10, 20,$ where for my taste $n=10$ shows peculiar tails in a probability plot, while $n=20$ seems tolerably close to normal. (Of course one must anticipate extra variability in the tails of all probability plots.)

Below show Kolmogorov-Smirnov one and two sample tests, where the beta samples are scaled to have $\mu=0,\sigma=1.$

[In R, $n=5000$ is the largest sample size accommodated by K-S tests. I have seen K-S tests deprecated for poor power and also for being 'too fussy' about minor departures from normality. Maybe there circumstances in which each complaint is valid. Overall, I have found K-S tests not to be as reliable as S-W tests--in situations where both apply.]

You can judge for yourself from tests and probability plots how closely my version of 'nearly normal' matches yours.

Noticeably nonnormal (S-W Rejects).

set.seed(2030)
n = 5000
z.10 = (rbeta(n,10,10)-.5)/sqrt(1/84)
z = rnorm(n)      # normal sample for comparison

shapiro.test(z.10)

        Shapiro-Wilk normality test

data:  z.10
W = 0.9988, p-value = 0.0009713

ks.test(z.10, z)

        Two-sample Kolmogorov-Smirnov test

data:  z.10 and z
D = 0.0122, p-value = 0.8508
alternative hypothesis: two-sided

ks.test(z.10, pnorm)

        One-sample Kolmogorov-Smirnov test

data:  z.10
D = 0.014173, p-value = 0.2677
alternative hypothesis: two-sided

More nearly normal (S-W Accepts).

z.20 = (rbeta(n,20,20)-.5)/sqrt(1/184)
shapiro.test(z.20)

        Shapiro-Wilk normality test

data:  z.20
W = 0.99947, p-value = 0.177

ks.test(z.20, z)

        Two-sample Kolmogorov-Smirnov test

data:  z.20 and z
D = 0.0182, p-value = 0.3791
alternative hypothesis: two-sided

ks.test(z.20, pnorm)

        One-sample Kolmogorov-Smirnov test

data:  z.20
D = 0.018738, p-value = 0.05973
alternative hypothesis: two-sided

Probability plots.

par(mfrow=c(1,3))
qqnorm(z.10, main="Scaled BETA(10,10)"); qqline(z.10, col="red", lwd=2)
qqnorm(z.20, main="Scaled BETA(20,20)"); qqline(z.20, col="blue", lwd=2)
qqnorm(z, main="NORM(1,1)"); qqline(z, col="green2", lwd=2)
par(mfrow=c(1,1))