Shapiro-Wilk Test in R: shapiro.test(), Its Sample-Size Limits, and When Not to Use It

Almost everyone who types shapiro.test() into R is using it to make one decision: parametric test or non-parametric test. That is the one decision the Shapiro-Wilk test should not be making for you, and the evidence against it is specific rather than stylistic. This page covers how to run it correctly in R, what R’s output does and does not tell you, and what to do at the two sample sizes where the test is most likely to be actively misleading.

Start here: what the p-value is not evidence for

The Shapiro-Wilk null hypothesis is that the sample came from a normally distributed population. A large p-value therefore means you failed to reject normality, not that normality holds. That distinction is not pedantry — it is the whole reason the test misbehaves at both ends of the sample-size range:

• Small n. The test has little power, so genuinely skewed data routinely pass. Rochon, Gondan and Kieser put it flatly: “for small samples, the Shapiro-Wilk test lacks power to detect deviations from normality.” A clean shapiro.test() on n = 12 is close to uninformative.
• Large n. The test becomes sensitive enough to flag deviations too small to matter. Ghasemi and Zahediasl note that in large samples “significant results would be derived even in the case of a small deviation from normality” — a deviation that typically will not affect the validity of the parametric test you were about to abandon.

So the test is weakest exactly where you need it (small samples, where robustness arguments are thinnest) and hyperactive exactly where you do not (large samples, where the central limit theorem has already done the work). Jeffrey Franc, writing in Prehospital and Disaster Medicine in 2025, calls this the sample-size paradox and argues that normality tests should stop being used as gatekeepers at all.

The two-stage procedure inflates conditional Type I error

The specific harm has been simulated. Rochon and colleagues modelled the standard workflow — Strategy I: “the two-sample t test was conducted if both samples had passed the preliminary Shapiro-Wilk test for normality; otherwise, we applied Mann-Whitney’s U test” — across equally sized samples of n = 10, 20, 30, 40 and 50, at pretest levels of .005, .010, .050 and .100.

For exponentially distributed data at n = 30, the Type I error of the t test conditional on having passed the pretest was 10.8% at α_pre = .005 and 17.0% at α_pre = .100, against a nominal 5%. The pretest does not screen out the bad cases; it selects the samples that happen to look normal, which is a biased subset.

The honest caveat, from the same paper: the unconditional Type I error of the whole two-stage procedure “can be considered robust,” staying near 5%. Their conclusion balances both findings: “From a formal perspective, preliminary testing for normality is incorrect and should therefore be avoided… From a practical perspective, however, preliminary testing does not seem to cause much harm, at least for the cases we have investigated. The worst that can be said is that preliminary testing is unnecessary… If the application of the t test is doubtful, the unconditional use of nonparametric tests seems to be the best choice.”

Read that as: decide your analysis from the design and the measurement scale, before you see the data. If you have reason to doubt a mean comparison, use a rank-based test unconditionally — see what the Mann-Whitney U test actually assumes before treating it as a free pass — rather than letting one p-value flip you between two analyses.

Running it in R

shapiro.test() lives in stats, which is attached by default. There is nothing to install.

set.seed(1)
x <- rnorm(100, mean = 5, sd = 3)
shapiro.test(x)

#   Shapiro-Wilk normality test
#
# data:  x
# W = ..., p-value = ...

It takes exactly one argument, a numeric vector. There is no data argument and no formula interface, so you pass a column: shapiro.test(df$weight), not shapiro.test(weight ~ group, data = df). It returns an object of class "htest" containing statistic (named W), p.value, method and data.name, so you can pull values out programmatically:

res <- shapiro.test(x)
res$statistic   # the W statistic
res$p.value     # the p-value

The R source is short enough to be worth knowing verbatim. Before computing anything it sorts the data and drops missing values with complete.cases(), so NAs are handled for you — but note that the effective n is the number of non-missing values, not nrow().

W and what it measures

W compares the order statistics of your sample against the order statistics expected under normality. It is bounded above by 1, and values close to 1 indicate close agreement with a normal distribution. Because W is a goodness-of-fit measure rather than a magnitude of practical deviation, W = 0.98 with p < 0.001 in a large sample and W = 0.98 with p = 0.4 in a small one describe roughly the same amount of departure from normality — the p-value is the thing that changed, and it changed because n changed. Report W alongside n, always.

The three ways R will refuse

All three are hard stop() calls in the source, not warnings:

• sample size must be between 3 and 5000 — the number of non-missing values must satisfy 3 ≤ n ≤ 5000. The upper bound is not arbitrary: R’s implementation follows Royston’s 1995 Remark AS R94 algorithm, which is specified for n up to 5000.
• all 'x' values are identical — triggered when the range is exactly zero. R also rescales internally when the range is below 1e-10, which is a numerical guard, not a statistical one.
• A non-numeric input fails the is.numeric(x) check. Factors and character columns are the usual culprits; a “numeric” column read from a CSV with a stray comma or unit suffix will have been coerced to character.

The p-value itself is approximated, and R’s documentation warns that it is “adequate for p.value < 0.1”. The underlying algorithm is exact at n = 3 and uses separate approximations for 4 ≤ n ≤ 11 and for n ≥ 12. In practice this means you should not read fine distinctions among large p-values — p = 0.62 versus p = 0.41 is not a meaningful comparison, and neither is evidence of normality.

Hitting the 5000 limit

The sample size must be between 3 and 5000 error is the single most common reason people land on this page. The reflex is to find a test that will run on 50,000 rows. That is the wrong reflex: at n = 50,000, any normality test will reject on deviations so small they are invisible in a histogram, and a t test or linear regression at that sample size is robust to exactly those deviations anyway.

Lumley, Diehr, Emerson and Chen made this argument directly in the Annual Review of Public Health: the belief that the t test and linear regression require normally distributed outcomes is a misconception; their real usefulness is that in large samples they are valid for any distribution. The authors demonstrated this by simulation on extremely non-normal data. If R is refusing to run shapiro.test() because your sample is too big, the sample is big enough that you did not need the test.

If you still want a formal test at n > 5000, the nortest package (version 1.0-4, maintained by Uwe Ligges) provides five omnibus normality tests with no such ceiling:

install.packages("nortest")
library(nortest)

ad.test(x)       # Anderson-Darling
cvm.test(x)      # Cramer-von Mises
lillie.test(x)   # Lilliefors (Kolmogorov-Smirnov)
pearson.test(x)  # Pearson chi-square
sf.test(x)       # Shapiro-Francia

One trap to avoid: do not reach for base R’s ks.test(x, "pnorm", mean(x), sd(x)). Estimating the mean and SD from the same data invalidates the reference distribution the test uses, making the p-value too large. That is what lillie.test() corrects for, and it is covered in detail in the guide to the Kolmogorov-Smirnov test and the estimated-parameter trap.

For models, test the residuals — not the outcome, and not the predictor

Linear regression and ANOVA assume normality of the errors, which you estimate with the residuals. They do not assume that y is normally distributed marginally, and they never assume anything about the distribution of x. Running shapiro.test(df$y) on a variable that is bimodal because it has two groups in it will reject normality every time, while the residuals within group may be perfectly well behaved.

fit <- lm(y ~ x + group, data = df)
shapiro.test(residuals(fit))

Even here the pretest logic does not hold up. Rochon’s Strategy II tested exactly this — the t test was run “if the residuals (x_i−m_X), (y_i−m_Y) from both samples had passed the pretest” — and it distorted conditional error rates as well.

Knief and Forstmeier’s 2021 Monte Carlo study in Behavior Research Methods is the useful counterweight to panic here: they found that p-values from Gaussian models fitted to non-normal data are generally reliable if either the dependent variable or the predictor is normally distributed, that bias appears only when both are heavily skewed, and that judging significance at α = 0.05 remains safe unless the sample size is very small. That is a much better guide to whether you have a problem than a residual p-value.

If you do conclude the residuals are badly non-normal, the response is a modelling decision, not a test switch: a different error family, a transformation with its own costs, robust standard errors, or a bootstrap.

Testing by group

Because shapiro.test() has no formula interface, group-wise testing needs a split. In base R:

tapply(df$y, df$group, shapiro.test)

# or, extracting just the p-values
sapply(split(df$y, df$group), function(v) shapiro.test(v)$p.value)

In the tidyverse, rstatix::shapiro_test() takes a data frame and unquoted variable names and returns a tidy table, which composes with dplyr::group_by():

library(rstatix)
iris %>% shapiro_test(Sepal.Length)
iris %>% shapiro_test(Sepal.Length, Petal.Width)

# multivariate normality
mshapiro_test(iris[, 1:3])

Two cautions. First, a group with fewer than three non-missing values throws the 3-to-5000 error and will abort a whole tapply() call. Second, testing six groups gives you six chances to reject by accident — and the more groups you test, the more likely at least one “fails,” which under the two-stage logic would push you to change the analysis for all of them. That multiplicity is rarely acknowledged and is another reason the pretest workflow is fragile.

Look at the plot, and prefer it

A Q-Q plot tells you the shape of the departure — right skew, heavy tails, discreteness, one outlier — which is the information that actually changes what you do. The p-value tells you only whether n was large enough to detect something.

qqnorm(residuals(fit)); qqline(residuals(fit))
hist(residuals(fit), breaks = 30)

Base R’s qqnorm() gives no sense of how much scatter is normal, which is why beginners over-read wiggles in the tails. car::qqPlot() fixes that by drawing a point-wise confidence envelope, on by default:

library(car)
qqPlot(x)                                    # default method: envelope = TRUE
qqPlot(fit, envelope = .99)                  # lm method

For the default method the reference distribution is "norm" and the comparison line is drawn through the quartiles. For an lm object the method plots studentized residuals against a t distribution by default, with a robust line and a simulated envelope. Points that stay inside the envelope are consistent with normal sampling variation, however wiggly the pattern looks.

Data that should never be Shapiro-Wilk tested

The test asks whether a sample came from a continuous normal distribution. Some variables cannot have come from any continuous distribution, so the answer is known before you run it:

• Likert items and other ordinal scales. A five-point item has five possible values. Testing it for normality answers a question that does not apply.
• Counts, proportions and bounded scores. These are better handled with an appropriate model family than with a normality decision.
• Heavily rounded measurements. Ties are the signature; R only errors when every value is identical, so a badly discretised variable will still return a p-value that looks legitimate.
• Aggregated data with structure. Repeated measures or clustered observations violate the independence the test assumes, before normality is even in question.

The conceptual background — what the normality assumption attaches to and how to check it — is covered in the guide to the normality of distribution assumption, which is worth reading before you decide the test result means anything.

How to report it

If a normality check belongs in your write-up at all, report it as a diagnostic, not as a decision rule:

• Give W, the p-value and n: “Shapiro-Wilk W = 0.97, p = 0.08, n = 64.” W and n without each other are uninterpretable.
• State what you tested — raw values, model residuals, or within-group values. These are different claims.
• Never write “the data were normally distributed (p > .05).” Write that the assumption was assessed and no substantial departure was evident, and say what the Q-Q plot showed.
• State your analysis choice and its justification independently of the test result. “We pre-specified a Mann-Whitney U test because the outcome is an ordinal score” is defensible; “Shapiro-Wilk was significant so we switched” is the two-stage procedure the simulations warn about.

If you are still deciding which test to run at all, work from the design and the measurement scale rather than from a normality p-value. That is the order in which the choice is defensible.

Frequently asked questions

How do I fix “Error in shapiro.test(x) : sample size must be between 3 and 5000”?

You cannot raise the limit — it is a hard stop() in R’s source, matching the range of Royston’s 1995 algorithm. If n > 5000, use nortest::ad.test() or another nortest function, or better, stop testing: at that sample size a t test or linear model is valid for essentially any outcome distribution. If n < 3, there is nothing to test. Do not sample 5000 rows at random to get under the ceiling — the result then depends on which rows you drew.

What p-value means my data are normal?

None. The test can only fail to reject normality; it cannot confirm it. A p-value above 0.05 in a small sample usually reflects low power rather than a normal population. Judge the Q-Q plot for shape and use the p-value, if at all, as a rough flag.

Should I run shapiro.test() on my raw variable or on the residuals?

On the residuals, for any regression or ANOVA. Those models assume normal errors, not a normal marginal outcome and never a normal predictor. For a simple two-group comparison, testing within each group is closer to the assumption than testing the pooled outcome, which will look bimodal whenever the groups differ.

Shapiro-Wilk or Kolmogorov-Smirnov in R?

Shapiro-Wilk, in almost all cases. Ghasemi and Zahediasl are blunt that the K-S test “has low power and it should not be seriously considered for testing normality,” and that Shapiro-Wilk “provides better power than the K-S test even after the Lilliefors correction.” If you use ks.test() for normality with parameters estimated from the same data, the p-value is wrong in the anti-conservative direction; use nortest::lillie.test() instead.

My reviewer asked for a normality test. What do I do?

Run it, report W, p and n, and pair it with a Q-Q plot — then justify your analysis on design grounds regardless of the result. The literature supports this: Franc’s 2025 editorial argues normality tests should not gatekeep methods, and Rochon and colleagues conclude that preliminary testing is at best unnecessary. If the outcome is ordinal or the mean is not the estimand you care about, say so; that is a stronger justification than any p-value.

Does shapiro.test() handle missing values?

Yes. The source removes them with complete.cases() before testing, so no na.rm argument is needed. Be aware that the sample-size limits apply to the count after removal — a column with 5200 rows and 300 NAs will run, and a group with two non-missing values will error.

Is this the same test SPSS runs?

Both implement the Shapiro-Wilk W, but SPSS reaches it through Analyze > Descriptive Statistics > Explore rather than a function call, and reports it beside the Kolmogorov-Smirnov statistic with a Lilliefors correction. The interpretation traps are identical; the menu route and output layout are covered in the guide to testing normality in SPSS.

References

• R Core Team. shapiro.test: Shapiro-Wilk Normality Test, R stats package documentation.
• R source code for shapiro.test: src/library/stats/R/shapiro.test.R (R sources mirror).
• Royston, P. (1995). Remark AS R94: A remark on Algorithm AS 181: The W test for normality. Applied Statistics, 44, 547–551. doi:10.2307/2986146
• Rochon, J., Gondan, M. & Kieser, M. (2012). To test or not to test: preliminary assessment of normality when comparing two independent samples. BMC Medical Research Methodology, 12, 81. PMC3444333
• Lumley, T., Diehr, P., Emerson, S. & Chen, L. (2002). The importance of the normality assumption in large public health data sets. Annual Review of Public Health, 23, 151–169. PMID 11910059
• Ghasemi, A. & Zahediasl, S. (2012). Normality tests for statistical analysis: a guide for non-statisticians. International Journal of Endocrinology and Metabolism. PMC3693611
• Knief, U. & Forstmeier, W. (2021). Violating the normality assumption may be the lesser of two evils. Behavior Research Methods, 53, 2576–2590. doi:10.3758/s13428-021-01587-5
• Franc, J. M. (2025). The misuse of normality tests as gatekeepers for research in prehospital and disaster medicine. Prehospital and Disaster Medicine, 40(5), 241–242. doi:10.1017/S1049023X25101465
• Ligges, U. nortest: Tests for Normality, CRAN, version 1.0-4.
• Fox, J. & Weisberg, S. qqPlot: Quantile-Comparison Plot, car package reference manual.
• Kassambara, A. shapiro_test: Shapiro-Wilk Normality Test, rstatix package reference manual.