When checking for normality, would you prefer numerical methods or graphical methods?

Question

I would like to know what is your preferred way to decide when the deviation from the normal distribution is relevant. when saying "relevant", I not saying necessarily "statistically significant".

For example, let's take the following residuals from an ANOVA mdoel:

res1 <-  c(-0.7890, -0.7090, -0.2390, 0.6610, -0.2390, 0.2610, -0.4390, -0.1390,  0.1610, -0.1390,  0.4700, -0.2300,  0.0200, -0.1300, -0.5300,  0.1700, 
-0.2300, -0.3300,  0.1700, -0.5300,  0.2205, -0.1795, -0.1295, -0.1995, -0.2395,  0.2305, -0.0695, -0.1295, -0.2995,  0.0905,  0.6600,  0.6600, 
 0.5600, -0.3400,  0.5600, -0.4400, -0.0400, -0.9400, -1.3400, -0.4400, -0.3390, -0.3390, -0.0390,  0.6610, -0.2390,  0.6610,  0.1610,  0.1610, 
0.6610,  0.2610,  0.4700, -0.1300,  0.0700, -0.0800, -0.0300,  0.4700, -0.0300, -0.0300,  0.4700, -0.0300,  0.3705, -0.1295, -0.1295, -0.1295,
-0.1295,  0.3705, -0.0295,  0.2705, -0.1295,  0.3705,  0.7600,  0.6600,  0.6600, -0.1400,  0.5600,  0.1600,  0.1600, -0.5400, -0.8400, -0.3400)

library(nortest)
ad.test(res1)

library("ggpubr")
ggqqplot(res1, ylab = "res1",
         ggtheme = theme_minimal())

In this example, the AD statistic is significant at 0.05, however, in the qqplot, most of the point fall inside the confidence interval. Would you say the numerical test and the qq graphic contradict each other. Or I'm misinterpreting the qqplot (observations form several horizontal lines, does this mean anything?).

Answer 1

The process of looking for non-normality, besides assuming you have a fairly large sample size, has a problem with distortion of statistical inference at the end (optimistically narrow confidence intervals, etc.). The entire process was never well founded in statistical theory, and semiparametrics and Bayesian methods are superior in almost every way. Start with John Krushke's BEST method and look at semiparametric ordinal regression, which is transformation-invariant with respect to $Y$. The latter is contained in a case study in my RMS book and course notes.

Tests of normality are essentially misleading because (1) they don't have a power of 1.0 so you'll get false negatives and (2) they tempt one to not pre-specify the analysis, which leads to problems mentioned above. Think about the Bayesian t-test for example: it can always allow for unequal variance and always allow for non-normality. Kruschke does the latter by using a t-distribution rather than a normal distribution for the raw data.

Answer 2

Comment: The procedures I would use and the results I would require would depend on the source of the data and the application at hand.

Not very demanding: Presumably, your data are residuals from some sort of linear model, and you have good reasons for expecting the residuals to be nearly normal. In this situation, I would want to look at the following results from R:

res1 <-  c(-0.7890, -0.7090, -0.2390,  0.6610, -0.2390,  0.2610, -0.4390,
           -0.1390,  0.1610, -0.1390,  0.4700, -0.2300,  0.0200, -0.1300, 
           -0.5300,  0.1700, -0.2300, -0.3300,  0.1700, -0.5300,  0.2205,
           -0.1795, -0.1295, -0.1995, -0.2395,  0.2305, -0.0695, -0.1295, 
           -0.2995,  0.0905,  0.6600,  0.6600,  0.5600, -0.3400,  0.5600,
           -0.4400, -0.0400, -0.9400, -1.3400, -0.4400, -0.3390, -0.3390, 
           -0.0390,  0.6610, -0.2390,  0.6610,  0.1610,  0.1610,  0.6610,
            0.2610,  0.4700, -0.1300,  0.0700, -0.0800, -0.0300,  0.4700,
           -0.0300, -0.0300,  0.4700, -0.0300,  0.3705, -0.1295, -0.1295, 
           -0.1295, -0.1295,  0.3705, -0.0295,  0.2705, -0.1295,  0.3705,  
            0.7600,  0.6600,  0.6600, -0.1400,  0.5600,  0.1600,  0.1600, 
           -0.5400, -0.8400, -0.3400)

par(mfrow=c(1,2))
 hist(res1, prob=T, col="skyblue2")
  curve(dnorm(x, mean(res1), sd(res1)), add=T, col="red")
 qqnorm(res1)
  abline(a=mean(res1), b=sd(res1), col="red")
par(mfrow=c(1,1))

At left, the normal curve that matches the mean and SD of the residuals does not show a very good fit to their histogram. This is not unusual for a sample as small as 80.

At right, except for bad behavior in the upper tail, where it departs markedly from a straight line, the normal probability plot of the residuals does not seem catastrophically bad.

Generally, one cannot expect residuals to look 'precisely' normal, and I have gotten what I judged to be useful results from linear models for which the residuals were no closer to normal than these. So using these graphical methods to assess normality, I would say the residuals seem close enough to normal for most practical purposes.

More demanding: However, if these observations were simulated from a program that claims to produce normally distributed observations, I would want to see results from several well-regarded formal tests of normality.

You have already seen that an Anderson-Darling test rejects the null hypothesis of normality. Below we see that a Shapiro-Wilk test also rejects. Furthermore, by the standards I would expect from a program to generate normal data, the graphical procedures shown above give disappointing results.

shapiro.test(res1)

        Shapiro-Wilk normality test

data:  res1
W = 0.96538, p-value = 0.02895

In this situation I would want to assess several larger samples. For example, looking at samples of sizes 100, 1000, and 5000, generated using the R function rnorm, I got Shapiro-Wilk p-values 0.17, 0.46, and 0.57, respectively. All consistent with normal.

From what I could infer from your current sample, I would not be willing to generate normal samples using a program that gives such results. I know that there are other programs that will do better and I would want to use one of them instead.