This guide explains how to check whether a variable follows a normal distribution. It covers the background of the normal curve, then three ways to assess normality (descriptive statistics, graphs and formal tests), and applies each step in R with a worked example.

Part 1. Background and Methods

1.1 The Normal Distribution

Short history. Abraham de Moivre found the curve in 1733 as an approximation to the binomial distribution. In the early 1800s, Laplace and Gauss developed it further. Gauss used it in 1809 to model measurement errors in astronomy, and Laplace linked it to the central limit theorem in 1810. (As a side note, Gauss is often called the “Prince of Mathematicians”.)

Description

# Draw the standard normal curve (mean = 0, SD = 1)
# dnorm() gives the density (height of the curve) for each z value
curve(dnorm(x), from = -4, to = 4, lwd = 2,
      main = "The Standard Normal Distribution",
      xlab = "z", ylab = "Density")

# Add a dashed vertical line at the mean (z = 0)
abline(v = 0, lty = 2)

The probability density function of the normal distribution is

\[f(x) = \frac{1}{\sigma\sqrt{2\pi}} \, e^{-\frac{(x-\mu)^2}{2\sigma^2}}\]

where \(\mu\) is the mean and \(\sigma\) is the standard deviation.

  • In a normal distribution, we see that the observations are densely concentrated around the mean, while the tails are sparse.
  • In other words, the probability of observing values near the mean is higher than in the tails.
  • The curve peaks in the middle and flattens out toward the tails.

Characteristics

  • Symmetric and bell-shaped
  • Mean, median and mode are equal
  • Defined by two parameters, the mean and the standard deviation
  • Asymptotic, the tails approach the x-axis but never touch it
  • The total area under the curve equals 1

Common misconception. The y-axis shows density, not probability. Density can even be greater than 1 when the standard deviation is small. Probability is the area under the curve between two values. For a continuous variable, the probability of one exact value is zero.

1.2 Areas Under the Curve and Confidence Levels

  • The 68-95-99.7 rule (within ±1, ±2 and ±3 standard deviations)
  • Critical z values, 1.96 for 95% and 2.58 for 99%
  • These critical values are used later for the skewness and kurtosis z-ratios

1.3 Interactive Exploration

Open the interactive normal-app to change the mean and standard deviation and watch the intervals.

1.4 Why Normality Matters

Many widely used procedures are parametric and rely on the normality assumption to some degree.

Which models require it. Normality belongs mainly to the general linear model family (t-test, ANOVA, ANCOVA, ordinary regression), where the residuals are assumed to be normal. In t-tests and ANOVA this means normality within each group. Generalized linear models, such as logistic or Poisson regression, assume other distributions instead. Many nonlinear regression models still assume normal errors. Nonparametric tests (Mann-Whitney U, Wilcoxon, Kruskal-Wallis) do not require normality, though they still assume independence and, for some interpretations, similarly shaped distributions.

Why it matters. The p-values and confidence intervals of these procedures come from the t and F distributions, which are correct only when the assumption holds. Violations can inflate or deflate the Type I error rate, reduce power and distort interval coverage.

The role of sample size. Thanks to the central limit theorem, parametric tests are fairly robust in large samples. Heavy skewness, extreme outliers and small or unequal groups can still distort the results.

When normality is not met. Options include transformations, handling outliers, nonparametric alternatives and bootstrapping. Whatever the choice, report it transparently.

1.5 Outliers and Common Causes of Non-normality

  • Outliers. A few extreme values can pull the mean, inflate the standard deviation and create skewness. They may be genuine but rare cases, or signs of a problem, so each one should be investigated before any decision is made.
  • Ceiling and floor effects. When a test is too easy, many scores pile up near the maximum (ceiling); when it is too hard, they pile up near the minimum (floor). Either way, one tail is cut off and the distribution becomes skewed.
  • Mixed subgroups. When a sample combines two different populations, for example two age groups or two levels of expertise, their scores can produce a flattened or bimodal shape even if each group is normal on its own.
  • Data entry errors. Typing mistakes, wrong codes or missing-value codes such as 999 left in the data can create impossible values that look like outliers. Checking minimum and maximum values is a simple first safeguard.

1.6 Methods for Checking Normality

  1. Descriptive statistics
    • Mean, median and mode (the mode is unstable for continuous scores)
    • Skewness and kurtosis divided by their standard errors, compared with ±1.96 (or ±2.58). These coefficients are sensitive to sample size.
  2. Standardized z-scores
    • If almost all z-scores fall between −3 and +3, the distribution is reasonably close to normal
    • Values beyond ±3 are possible outliers
  3. Graphical methods
    • Histogram with a normal curve, boxplot, Q-Q plot
  4. Formal hypothesis tests
    • Shapiro-Wilk
    • Kolmogorov-Smirnov with the Lilliefors correction, which provides a safeguard against Type I error when the mean and SD are estimated from the sample

1.7 Cautions

  • Large samples make tests and z-ratios significant even for trivial departures (high power)
  • Small samples miss real problems (low power, Type II error)
  • Never rely on a single method, use the evidence together

Part 2. Application in R

2.1 Setup and the Data Set

The sat.act data in the psych package come from the SAPA online project. It has 700 people and six variables, gender (1 = male, 2 = female), education (0 to 5), age, ACT composite, SATV and SATQ. Scores are self-reported, and SATQ has some missing values.

install.packages(c("psych", "nortest", "DescTools"))  # run once
library(psych)    # sat.act data, describe(), skew(), kurtosi()
library(nortest)  # lillie.test(): K-S test with Lilliefors correction

data(sat.act)       # load the dataset
head(sat.act)       # first six rows
##       gender education age ACT SATV SATQ
## 29442      2         3  19  24  500  500
## 29457      2         3  23  35  600  500
## 29498      2         3  20  21  480  470
## 29503      1         4  27  26  550  520
## 29504      1         2  33  31  600  550
## 29518      1         5  26  28  640  640
str(sat.act)        # structure of the variables
## 'data.frame':    700 obs. of  6 variables:
##  $ gender   : int  2 2 2 1 1 1 2 1 2 2 ...
##  $ education: int  3 3 3 4 2 5 5 3 4 5 ...
##  $ age      : int  19 23 20 27 33 26 30 19 23 40 ...
##  $ ACT      : int  24 35 21 26 31 28 36 22 22 35 ...
##  $ SATV     : int  500 600 480 550 600 640 610 520 400 730 ...
##  $ SATQ     : int  500 500 470 520 550 640 500 560 600 800 ...
describe(sat.act)   # descriptive statistics
##           vars   n   mean     sd median trimmed    mad min max range  skew kurtosis   se
## gender       1 700   1.65   0.48      2    1.68   0.00   1   2     1 -0.61    -1.62 0.02
## education    2 700   3.16   1.43      3    3.31   1.48   0   5     5 -0.68    -0.07 0.05
## age          3 700  25.59   9.50     22   23.86   5.93  13  65    52  1.64     2.42 0.36
## ACT          4 700  28.55   4.82     29   28.84   4.45   3  36    33 -0.66     0.53 0.18
## SATV         5 700 612.23 112.90    620  619.45 118.61 200 800   600 -0.64     0.33 4.27
## SATQ         6 687 610.22 115.64    620  617.25 118.61 200 800   600 -0.59    -0.02 4.41

The output confirms 700 cases, with 687 valid values for SATQ. Note that gender and education are coded categories, so their skewness and kurtosis values are not meaningful.

2.2 Descriptive Check

x <- sat.act$SATV
mean(x); median(x)
## [1] 612.2343
## [1] 620
# Mode
as.numeric(names(sort(table(x), decreasing = TRUE))[1])
## [1] 600
# Skewness and kurtosis with standard errors
n  <- length(x)
sk <- skew(x, type = 2)      # type 2 matches SPSS and the SE formula below
ku <- kurtosi(x, type = 2)
se_sk <- sqrt(6 * n * (n - 1) / ((n - 2) * (n + 1) * (n + 3)))
se_ku <- 2 * se_sk * sqrt((n^2 - 1) / ((n - 3) * (n + 5)))
c(z_skew = sk / se_sk, z_kurt = ku / se_ku)
##    z_skew    z_kurt 
## -6.998793  1.873844
# Ready-made alternative: DescTools
library(DescTools)
Skew(x, method = 2, conf.level = 0.95, ci.type = "classic")  # skewness with 95% CI
##   skewness     lwr.ci     upr.ci 
## -0.6465796 -0.1810702  0.1810702
Kurt(x, method = 2, conf.level = 0.95, ci.type = "classic")  # excess kurtosis with 95% CI
##   kurtosis     lwr.ci     upr.ci 
##  0.3457394 -0.3616293  0.3616293

Interpreting the output. The mean of SATV (612.23) is lower than the median (620), which points to a negative (left) skew. The skewness coefficient is −0.65, and its z ratio (−7.00) is well beyond −2.58, so the skewness is significant at the 0.01 level. The excess kurtosis is 0.35, with a z ratio of 1.87, which stays within ±1.96, so the kurtosis is not significant. In the DescTools output with ci.type = "classic", the interval is centered at zero (±1.96 × SE) and works as an acceptance region: the skewness (−0.65) falls outside [−0.18, 0.18], while the kurtosis (0.35) falls just inside [−0.36, 0.36]. Both approaches lead to the same conclusion. SATV is moderately left-skewed, consistent with the ceiling at 800, while its tails are close to normal.

2.3 Standardized z-Scores and Outliers

z <- as.numeric(scale(x))
range(z)
## [1] -3.651239  1.663077
sum(abs(z) > 3)
## [1] 4
which(abs(z) > 3)
## [1]  99 220 251 440

Interpreting the output. The z-scores range from −3.65 to 1.66. Four cases (rows 99, 220, 251 and 440) fall below −3, and none exceed +3. This is 0.6% of the sample, slightly more than the 0.3% expected under normality, and all of them are on the low side. That pattern matches the negative skewness found in 2.2: the long tail is on the left, while high scores are cut off by the ceiling at 800. Each flagged case should still be checked with sat.act[c(99, 220, 251, 440), ] to see whether it is a genuine low score or an error.

2.4 Graphical Check

hist(x, freq = FALSE, main = "SAT Verbal", xlab = "SATV")
curve(dnorm(x, mean(sat.act$SATV), sd(sat.act$SATV)), add = TRUE, lwd = 2)

The SATV histogram. The bars show a long left tail stretching down to 200, with only a few cases below 400. Most scores are concentrated between 550 and 700. On the right side, the bars between 750 and 800 sit clearly above the normal curve, and the distribution stops abruptly at 800. This sharp cut-off is the ceiling effect: the test cannot record scores above 800, so high scorers pile up at the maximum. Together, the long left tail and the ceiling explain the negative skewness found in 2.2.

boxplot(x, horizontal = TRUE, main = "SAT Verbal")

The SATV boxplot. The box runs from about 550 (Q1) to 700 (Q3), with the median at 620, slightly closer to the upper edge of the box. The left whisker is much longer, reaching down to about 330, while the right whisker stops at 800, the ceiling. Several points beyond the left whisker (around 200, 240 and 300) are flagged as potential outliers, all on the low side. Note that the boxplot rule (1.5 × IQR) flags a few more cases than the |z| > 3 rule in 2.3, since the two criteria differ. Both, however, point to the same pattern: a long left tail and no unusual high scores.

qqnorm(x); qqline(x, lwd = 2)

The SATV Q-Q plot. In the middle range (theoretical quantiles from about −1 to +1.5), the points follow the reference line closely. At the lower end, they fall well below the line, meaning the lowest scores are more extreme than a normal distribution would predict; this is the long left tail. At the upper end, from about +2 onward, the points flatten into a horizontal row at 800, where many students share the maximum score; this is the ceiling effect. The step-like pattern throughout reflects the fact that SAT scores are reported in steps of 10, so many cases share the same value. Taken together, the three graphs tell the same story as the descriptive statistics: SATV is moderately left-skewed, with a ceiling at 800.

2.5 Formal Hypothesis Tests

shapiro.test(x)
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.96708, p-value = 1.934e-11
lillie.test(x)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  x
## D = 0.10971, p-value < 2.2e-16

Interpreting the output. Both tests take normality as the null hypothesis. If p > 0.05, the null is retained, meaning there is no evidence against normality. If p < 0.05, the distribution differs significantly from normal. The W statistic of Shapiro-Wilk ranges up to 1, and values closer to 1 indicate closer agreement with normality; for the Lilliefors test, a smaller D indicates closer agreement.

For SATV. Both tests reject normality. The Shapiro-Wilk test gives W = 0.967 with p < 0.001, and the Lilliefors test gives D = 0.110 with p < 0.001. Although W is fairly close to 1, the large sample (n = 700) gives the tests enough power to detect even modest departures. In this case, however, the departure is not trivial: the descriptive statistics and all three graphs confirm a moderate left skew and a ceiling at 800. The formal tests therefore agree with the other evidence rather than flagging a negligible deviation.

2.6 By Group (Gender)

by(sat.act$SATV, sat.act$gender, shapiro.test)
## sat.act$gender: 1
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.96258, p-value = 4.76e-06
## 
## ------------------------------------------------------------------------------------------ 
## sat.act$gender: 2
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.9669, p-value = 1.387e-08
describeBy(sat.act$SATV, group = sat.act$gender)
## 
##  Descriptive statistics by group 
## group: 1
##    vars   n   mean     sd median trimmed    mad min max range  skew kurtosis   se
## X1    1 247 615.11 114.16    630  622.07 118.61 200 800   600 -0.63     0.13 7.26
## ------------------------------------------------------------------------------------------ 
## group: 2
##    vars   n   mean     sd median trimmed    mad min max range  skew kurtosis   se
## X1    1 453 610.66 112.31    620  617.91 103.78 200 800   600 -0.65     0.42 5.28
par(mfrow = c(1, 2))
for (g in 1:2) {
  y <- sat.act$SATV[sat.act$gender == g]
  qqnorm(y, main = paste("Gender", g)); qqline(y)
}

par(mfrow = c(1, 1))

2.7 Interpretation

All the evidence points in the same direction. The descriptive statistics show a significant negative skew with near-normal kurtosis, the z-scores flag four low outliers and no high ones, the histogram, boxplot and Q-Q plot all show a long left tail and a ceiling at 800, and both formal tests reject normality. SATV is therefore moderately left-skewed. The ceiling at 800 is a real cause of this skew, so the departure reflects the measurement itself rather than data errors.